use of uk.ac.sussex.gdsc.core.math.interpolation.CubicSplinePosition in project GDSC-SMLM by aherbert.
the class CubicSplineCalculatorTest method canComputeCoefficientsForGaussianFunction.
@Test
void canComputeCoefficientsForGaussianFunction() {
final int x = 4;
final int y = 4;
final int z = 4;
final double xscale = 1;
final double yscale = 0.5;
final double zscale = 2.0;
final double[] xval = SimpleArrayUtils.newArray(x, 0, xscale);
final double[] yval = SimpleArrayUtils.newArray(y, 0, yscale);
final double[] zval = SimpleArrayUtils.newArray(z, 0, zscale);
final double[][][] fval = createData(x, y, z, null);
final CustomTricubicInterpolatingFunction f1 = new CustomTricubicInterpolator().interpolate(xval, yval, zval, fval);
final double[] exp = new double[64];
f1.getSplineNode(1, 1, 1).getCoefficients(exp);
final CustomTricubicFunction f = CustomTricubicFunctionUtils.create(exp);
final CubicSplinePosition[] s = new CubicSplinePosition[4];
for (int i = 0; i < 4; i++) {
s[i] = new CubicSplinePosition((double) i / 3);
}
final double[][][] value = new double[4][4][4];
final double[] b = new double[64];
int count = 0;
for (int k = 0; k < 4; k++) {
for (int j = 0; j < 4; j++) {
for (int i = 0; i < 4; i++) {
value[i][j][k] = f.value(s[i], s[j], s[k]);
b[count++] = value[i][j][k];
}
}
}
final CubicSplineCalculator calc = new CubicSplineCalculator();
double[] obs = calc.compute(value);
Assertions.assertArrayEquals(exp, obs, 1e-6);
obs = calc.compute(b);
Assertions.assertArrayEquals(exp, obs, 1e-6);
}
use of uk.ac.sussex.gdsc.core.math.interpolation.CubicSplinePosition in project GDSC-SMLM by aherbert.
the class CubicSplineData method sample.
/**
* Sample the function.
*
* <p>n samples will be taken per node in each dimension. A final sample is taken at the end of
* the sample range thus the final range for each axis will be the current axis range.
*
* <p>The procedure setValue(int,int,int,double) method will be executed in ZYX order.
*
* @param nx the number of samples per spline node in the x dimension
* @param ny the number of samples per spline node in the y dimension
* @param nz the number of samples per spline node in the z dimension
* @param procedure the procedure
* @param progress the progress
* @throws IllegalArgumentException If the number of sample is not positive
*/
public void sample(int nx, int ny, int nz, TrivalueProcedure procedure, TrackProgress progress) {
if (nx < 1 || ny < 1 || nz < 1) {
throw new IllegalArgumentException("Samples must be positive");
}
// We can interpolate all nodes n-times plus a final point at the last node
final int maxx = (getMaxX() - 1) * nx;
final int maxy = (getMaxY() - 1) * ny;
final int maxz = (getMaxZ() - 1) * nz;
if (!procedure.setDimensions(maxx + 1, maxy + 1, maxz + 1)) {
return;
}
final Ticker ticker = Ticker.createStarted(progress, (long) (maxx + 1) * (maxy + 1) * (maxz + 1), false);
// Pre-compute interpolation tables
final CubicSplinePosition[] sx = createCubicSplinePosition(nx);
final CubicSplinePosition[] sy = createCubicSplinePosition(ny);
final CubicSplinePosition[] sz = createCubicSplinePosition(nz);
// Write axis values
// Cache the table and the spline position to use for each interpolation point
final int[] xt = new int[maxx + 1];
final int[] xp = new int[maxx + 1];
for (int x = 0; x <= maxx; x++) {
int xposition = x / nx;
int xtable = x % nx;
if (x == maxx) {
// Final interpolation point
xposition--;
xtable = nx;
}
xt[x] = xtable;
xp[x] = xposition;
procedure.setX(x, xposition + (double) xtable / nx);
}
final int[] yt = new int[maxy + 1];
final int[] yp = new int[maxy + 1];
for (int y = 0; y <= maxy; y++) {
int yposition = y / ny;
int ytable = y % ny;
if (y == maxy) {
// Final interpolation point
yposition--;
ytable = ny;
}
yt[y] = ytable;
yp[y] = yposition;
procedure.setY(y, yposition + (double) ytable / ny);
}
final int[] zt = new int[maxz + 1];
final int[] zp = new int[maxz + 1];
for (int z = 0; z <= maxz; z++) {
int zposition = z / nz;
int ztable = z % nz;
if (z == maxz) {
// Final interpolation point
zposition--;
ztable = nz;
}
zt[z] = ztable;
zp[z] = zposition;
procedure.setZ(z, zposition + (double) ztable / nz);
}
// Write interpolated values
for (int z = 0; z <= maxz; z++) {
final CustomTricubicFunction[] xySplines = splines[zp[z]];
for (int y = 0; y <= maxy; y++) {
final int index = yp[y] * getMaxX();
for (int x = 0; x <= maxx; x++) {
procedure.setValue(x, y, z, xySplines[index + xp[x]].value(sx[xt[x]], sy[yt[y]], sz[zt[z]]));
ticker.tick();
}
}
}
ticker.stop();
}
use of uk.ac.sussex.gdsc.core.math.interpolation.CubicSplinePosition in project GDSC-SMLM by aherbert.
the class CubicSplineData method createCubicSplinePosition.
private static CubicSplinePosition[] createCubicSplinePosition(int n) {
// Use an extra one to have the final x=1 interpolation point.
final int n1 = n + 1;
final double step = 1.0 / n;
final CubicSplinePosition[] s = new CubicSplinePosition[n1];
for (int x = 0; x < n; x++) {
s[x] = new CubicSplinePosition(x * step);
}
// Final interpolation point must be exactly 1
s[n] = new CubicSplinePosition(1);
return s;
}
use of uk.ac.sussex.gdsc.core.math.interpolation.CubicSplinePosition in project GDSC-SMLM by aherbert.
the class CubicSplineCalculatorTest method canComputeCoefficientsForDistanceFunction.
@Test
void canComputeCoefficientsForDistanceFunction() {
final double[] exp = new double[64];
int count = 0;
for (int k = 0; k < 4; k++) {
for (int j = 0; j < 4; j++) {
for (int i = 0; i < 4; i++) {
exp[count++] = Math.sqrt(i * i + j * j + k * k);
}
}
}
final CustomTricubicFunction f = CustomTricubicFunctionUtils.create(exp);
final CubicSplinePosition[] s = new CubicSplinePosition[4];
for (int i = 0; i < 4; i++) {
s[i] = new CubicSplinePosition((double) i / 3);
}
final double[][][] value = new double[4][4][4];
final double[] b = new double[64];
count = 0;
for (int k = 0; k < 4; k++) {
for (int j = 0; j < 4; j++) {
for (int i = 0; i < 4; i++) {
value[i][j][k] = f.value(s[i], s[j], s[k]);
b[count++] = value[i][j][k];
}
}
}
final CubicSplineCalculator calc = new CubicSplineCalculator();
double[] obs = calc.compute(value);
Assertions.assertArrayEquals(exp, obs, 1e-6);
obs = calc.compute(b);
Assertions.assertArrayEquals(exp, obs, 1e-6);
}
Aggregations