Examples with CustomTricubicFunction uk.ac.sussex.gdsc.core.math.interpolation.CustomTricubicFunction used on opensource projects

Example 1 with CustomTricubicFunction

use of uk.ac.sussex.gdsc.core.math.interpolation.CustomTricubicFunction 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);
}

Example 2 with CustomTricubicFunction

use of uk.ac.sussex.gdsc.core.math.interpolation.CustomTricubicFunction in project GDSC-SMLM by aherbert.

the class CubicSplineDataTest method canExternaliseFunction.

private static void canExternaliseFunction(RandomSeed seed, boolean singlePrecision) throws IOException {
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    final int x = 6;
    final int y = 5;
    final int z = 4;
    final int size = x * y;
    final CustomTricubicFunction[][] splines = new CustomTricubicFunction[z][x * y];
    final double[] a = new double[64];
    for (int zz = 0; zz < z; zz++) {
        for (int i = 0; i < size; i++) {
            for (int j = 0; j < 64; j++) {
                a[j] = r.nextDouble();
            }
            splines[zz][i] = CustomTricubicFunctionUtils.create(a);
            if (singlePrecision) {
                splines[zz][i] = splines[zz][i].toSinglePrecision();
            }
        }
    }
    final CubicSplineData f1 = new CubicSplineData(x, y, splines);
    final ByteArrayOutputStream b = new ByteArrayOutputStream();
    f1.write(b);
    final byte[] bytes = b.toByteArray();
    final CubicSplineData f2 = CubicSplineData.read(new ByteArrayInputStream(bytes));
    final double[] exp = new double[64];
    final double[] obs = new double[64];
    for (int zz = 0; zz < z; zz++) {
        for (int i = 0; i < size; i++) {
            f1.splines[zz][i].getCoefficients(exp);
            f2.splines[zz][i].getCoefficients(obs);
            Assertions.assertArrayEquals(exp, obs);
        }
    }
}

Example 3 with CustomTricubicFunction

use of uk.ac.sussex.gdsc.core.math.interpolation.CustomTricubicFunction in project GDSC-SMLM by aherbert.

the class CubicSplineManager method createCubicSpline.

/**
 * Creates the cubic spline.
 *
 * @param imagePsf the image PSF details
 * @param image the image
 * @param singlePrecision Set to true to use single precision (float values) to store the cubic
 *        spline coefficients
 * @return the cubic spline PSF
 */
public static CubicSplinePsf createCubicSpline(ImagePSFOrBuilder imagePsf, ImageStack image, final boolean singlePrecision) {
    final int maxx = image.getWidth();
    final int maxy = image.getHeight();
    final int maxz = image.getSize();
    final float[][] psf = new float[maxz][];
    for (int z = 0; z < maxz; z++) {
        psf[z] = ImageJImageConverter.getData(image.getPixels(z + 1), null);
    }
    // We reduce by a factor of 3
    final int maxi = (maxx - 1) / 3;
    final int maxj = (maxy - 1) / 3;
    final int maxk = (maxz - 1) / 3;
    final int size = maxi * maxj;
    final CustomTricubicFunction[][] splines = new CustomTricubicFunction[maxk][size];
    final int threadCount = Prefs.getThreads();
    final Ticker ticker = ImageJUtils.createTicker((long) maxi * maxj * maxk, threadCount);
    final ExecutorService threadPool = Executors.newFixedThreadPool(threadCount);
    final LocalList<Future<?>> futures = new LocalList<>(maxk);
    // spline node along each dimension, i.e. dimension length = n*3 + 1 with n the number of nodes.
    for (int k = 0; k < maxk; k++) {
        final int kk = k;
        futures.add(threadPool.submit(() -> {
            final CubicSplineCalculator calc = new CubicSplineCalculator();
            final double[] value = new double[64];
            final int zz = 3 * kk;
            for (int j = 0, index = 0; j < maxj; j++) {
                // 4x4 block origin in the XY data
                int index0 = 3 * j * maxx;
                for (int i = 0; i < maxi; i++, index++) {
                    ticker.tick();
                    int count = 0;
                    for (int z = 0; z < 4; z++) {
                        final float[] data = psf[zz + z];
                        for (int y = 0; y < 4; y++) {
                            for (int x = 0, ii = index0 + y * maxx; x < 4; x++) {
                                value[count++] = data[ii++];
                            }
                        }
                    }
                    splines[kk][index] = CustomTricubicFunctionUtils.create(calc.compute(value));
                    if (singlePrecision) {
                        splines[kk][index] = splines[kk][index].toSinglePrecision();
                    }
                    index0 += 3;
                }
            }
        }));
    }
    ticker.stop();
    threadPool.shutdown();
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    // Normalise
    double maxSum = 0;
    for (int k = 0; k < maxk; k++) {
        double sum = 0;
        for (int i = 0; i < size; i++) {
            sum += splines[k][i].value000();
        }
        if (maxSum < sum) {
            maxSum = sum;
        }
    }
    if (maxSum == 0) {
        throw new IllegalStateException("The cubic spline has no maximum signal");
    }
    final double scale = 1.0 / maxSum;
    for (int k = 0; k < maxk; k++) {
        for (int i = 0; i < size; i++) {
            splines[k][i] = splines[k][i].scale(scale);
        }
    }
    // Create on an integer scale
    final CubicSplineData f = new CubicSplineData(maxi, maxj, splines);
    // Create a new info with the PSF details
    final ImagePSF.Builder b = ImagePSF.newBuilder();
    b.setImageCount(imagePsf.getImageCount());
    // Reducing the image has the effect of enlarging the pixel size
    b.setPixelSize(imagePsf.getPixelSize() * 3.0);
    b.setPixelDepth(imagePsf.getPixelDepth() * 3.0);
    // The centre has to be moved as we reduced the image size by 3.
    // In the ImagePSF the XY centre puts 0.5 at the centre of the pixel.
    // The spline puts 0,0 at the centre of each pixel for convenience.
    double cx = maxi / 2.0;
    if (imagePsf.getXCentre() != 0) {
        cx = (imagePsf.getXCentre() - 0.5) / 3;
    }
    double cy = maxj / 2.0;
    if (imagePsf.getYCentre() != 0) {
        cy = (imagePsf.getYCentre() - 0.5) / 3;
    }
    double cz = maxk / 2.0;
    if (imagePsf.getZCentre() != 0) {
        cz = imagePsf.getZCentre() / 3;
    } else if (imagePsf.getCentreImage() != 0) {
        cz = (imagePsf.getCentreImage() - 1) / 3.0;
    }
    b.setXCentre(cx);
    b.setYCentre(cy);
    b.setZCentre(cz);
    return new CubicSplinePsf(b.build(), f);
}

Example 4 with CustomTricubicFunction

use of uk.ac.sussex.gdsc.core.math.interpolation.CustomTricubicFunction 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();
}

Example 5 with CustomTricubicFunction

use of uk.ac.sussex.gdsc.core.math.interpolation.CustomTricubicFunction in project GDSC-SMLM by aherbert.

the class CubicSplineData method read.

/**
 * Read a tricubic spline from the input stream.
 *
 * <p>Note: For best performance a buffered input stream should be used.
 *
 * @param inputStream the input stream
 * @param progress the progress
 * @return the tricubic spline data
 * @throws IOException Signals that an I/O exception has occurred.
 */
public static CubicSplineData read(InputStream inputStream, TrackProgress progress) throws IOException {
    // Read dimensions
    final DataInput in = new DataInputStream(inputStream);
    final int maxx = in.readInt();
    final int maxy = in.readInt();
    final int maxz = in.readInt();
    final Ticker ticker = Ticker.createStarted(progress, (long) maxx * maxy * maxz, false);
    // Read precision
    final boolean singlePrecision = in.readBoolean();
    final SplineReader reader = (singlePrecision) ? new FloatSplineReader() : new DoubleSplineReader();
    final int size = maxx * maxy;
    final CustomTricubicFunction[][] splines = new CustomTricubicFunction[maxz][maxx * maxy];
    for (int z = 0; z < maxz; z++) {
        for (int i = 0; i < size; i++) {
            splines[z][i] = reader.read(in);
            ticker.tick();
        }
    }
    ticker.stop();
    // Skip validation
    return new CubicSplineData(maxx, maxy, splines, false);
}