Examples with FhtFilter uk.ac.sussex.gdsc.smlm.ij.filters.FhtFilter used on opensource projects

Example 1 with FhtFilter

use of uk.ac.sussex.gdsc.smlm.ij.filters.FhtFilter in project GDSC-SMLM by aherbert.

the class JTransformsTest method canComputeUsingFft.

private static void canComputeUsingFft(boolean convolution) {
    final int size = 16;
    final int ex = 5;
    final int ey = 7;
    final int ox = 1;
    final int oy = 2;
    final FloatProcessor fp1 = createProcessor(size, ex, ey, 4, 4, null);
    final FloatProcessor fp2 = createProcessor(size, size / 2 + ox, size / 2 + oy, 4, 4, null);
    final float[] input1 = (float[]) fp1.getPixels();
    final float[] input2 = (float[]) fp2.getPixels();
    final FhtFilter ff = new FhtFilter(input2.clone(), size, size);
    ff.setOperation((convolution) ? Operation.CONVOLUTION : Operation.CORRELATION);
    final float[] e = input1.clone();
    ff.filter(e, size, size);
    // Do the same with JTransforms
    final float[] data1 = new float[input1.length * 2];
    final FloatFFT_2D fft = new FloatFFT_2D(size, size);
    System.arraycopy(input1, 0, data1, 0, input1.length);
    final float[] data2 = new float[data1.length];
    System.arraycopy(input2, 0, data2, 0, input2.length);
    fft.realForwardFull(data1);
    fft.realForwardFull(data2);
    // https://en.wikipedia.org/wiki/Complex_number#Multiplication_and_division
    for (int i = 0; i < data2.length; i += 2) {
        final float a = data1[i];
        final float b = data1[i + 1];
        final float c = data2[i];
        // Get the conjugate for correlation
        final float d = (convolution) ? data2[i + 1] : -data2[i + 1];
        data1[i] = a * c - b * d;
        data1[i + 1] = b * c + a * d;
    }
    fft.complexInverse(data1, true);
    final float[] o = new float[e.length];
    for (int i = 0, j = 0; i < o.length; i++, j += 2) {
        o[i] = data1[j];
    }
    Fht.swapQuadrants(new FloatProcessor(size, size, o));
    Assertions.assertArrayEquals(e, o, 1e-3f);
}

Example 2 with FhtFilter

use of uk.ac.sussex.gdsc.smlm.ij.filters.FhtFilter in project GDSC-SMLM by aherbert.

the class ImageKernelFilter method run.

@Override
public void run(ImageProcessor ip) {
    final float[] data = (float[]) ip.getPixels();
    final int w = ip.getWidth();
    final int h = ip.getHeight();
    if (settings.method == Settings.METHOD_SPATIAL) {
        kf.convolve(data, w, h, settings.border);
    } else {
        // Use a clone for thread safety
        final FhtFilter f = (ticker.getTotal() > 1) ? ff.copy() : ff;
        f.filter(data, w, h, settings.border);
    }
    if (ticker.getTotal() == 1) {
        ip.resetMinAndMax();
    }
    ticker.tick();
}

Example 3 with FhtFilter

use of uk.ac.sussex.gdsc.smlm.ij.filters.FhtFilter in project GDSC-SMLM by aherbert.

the class ImageKernelFilter method setNPasses.

@Override
public void setNPasses(int passes) {
    // Create the kernel from the image
    boolean build = kernelImp.getID() != lastId || settings.method != lastMethod || settings.filter != lastFilter;
    build = build || (settings.method == Settings.METHOD_SPATIAL && kf == null);
    build = build || (settings.method == Settings.METHOD_FHT && ff == null);
    if (build) {
        final Operation operation = Operation.forOrdinal(settings.filter);
        FloatProcessor fp = kernelImp.getProcessor().toFloat(0, null);
        if (settings.method == Settings.METHOD_SPATIAL) {
            if (kf == null || kernelImp.getID() != lastId || settings.zero != lastZero) {
                fp = pad(fp);
                final int kw = fp.getWidth();
                final int kh = fp.getHeight();
                final float[] kernel = (float[]) fp.getPixels();
                kf = (settings.zero) ? new ZeroKernelFilter(kernel, kw, kh) : new KernelFilter(kernel, kw, kh);
            }
            switch(operation) {
                case CONVOLUTION:
                    kf.setConvolution(true);
                    break;
                case CORRELATION:
                    kf.setConvolution(false);
                    break;
                case DECONVOLUTION:
                default:
                    // Spatial filtering does not support anything other than convolution or correlation.
                    ImageJUtils.log("Unsupported operation (%s), default to correlation", operation.getName());
                    kf.setConvolution(false);
                    break;
            }
        } else {
            if (ff == null || kernelImp.getID() != lastId) {
                final int kw = fp.getWidth();
                final int kh = fp.getHeight();
                final float[] kernel = (float[]) fp.getPixels();
                ff = new FhtFilter(kernel, kw, kh);
                ff.initialiseKernel(dataImp.getWidth(), dataImp.getHeight());
            }
            ff.setOperation(operation);
        }
        lastId = kernelImp.getID();
        lastMethod = settings.method;
        lastFilter = settings.filter;
        lastZero = settings.zero;
    }
    ticker = ImageJUtils.createTicker(passes, passes);
}

Example 4 with FhtFilter

use of uk.ac.sussex.gdsc.smlm.ij.filters.FhtFilter in project GDSC-SMLM by aherbert.

the class FhtFilterTest method canWindow.

@Test
void canWindow() {
    final int size = 16;
    final float[] in = SimpleArrayUtils.newFloatArray(size * size, 1);
    final FhtFilter f = new FhtFilter(new float[1], 1, 1);
    for (int i = 1; i < 5; i++) {
        final double[] wx = ImageWindow.tukeyEdge(size, i);
        final float[] e = ImageWindow.applyWindowSeparable(in, size, size, wx, wx);
        final float[] o = in.clone();
        f.applyBorder(o, size, size, i);
        Assertions.assertArrayEquals(e, o);
    }
}

Example 5 with FhtFilter

use of uk.ac.sussex.gdsc.smlm.ij.filters.FhtFilter in project GDSC-SMLM by aherbert.

the class FhtFilterTest method canFilter.

private static void canFilter(RandomSeed seed, Operation operation) {
    final int size = 16;
    final int ex = 5;
    final int ey = 7;
    final int ox = 1;
    final int oy = 2;
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    final FloatProcessor fp1 = createProcessor(size, ex, ey, 4, 4, r);
    // This is offset from the centre
    final FloatProcessor fp2 = createProcessor(size, size / 2 + ox, size / 2 + oy, 4, 4, r);
    final float[] input1 = ((float[]) fp1.getPixels()).clone();
    final float[] input2 = ((float[]) fp2.getPixels()).clone();
    final FHT fht1 = new FHT(fp1);
    fht1.transform();
    final FHT fht2 = new FHT(fp2);
    fht2.transform();
    FHT fhtE;
    switch(operation) {
        case CONVOLUTION:
            fhtE = fht1.multiply(fht2);
            break;
        case CORRELATION:
            fhtE = fht1.conjugateMultiply(fht2);
            break;
        case DECONVOLUTION:
            fhtE = fht1.divide(fht2);
            break;
        default:
            throw new RuntimeException();
    }
    fhtE.inverseTransform();
    fhtE.swapQuadrants();
    final float[] e = (float[]) fhtE.getPixels();
    if (operation == Operation.CORRELATION) {
        // Test the max correlation position
        final int max = SimpleArrayUtils.findMaxIndex(e);
        final int x = max % 16;
        final int y = max / 16;
        Assertions.assertEquals(ex, x + ox);
        Assertions.assertEquals(ey, y + oy);
    }
    // Test verses a spatial domain filter in the middle of the image
    if (operation != Operation.DECONVOLUTION) {
        double sum = 0;
        float[] i2 = input2;
        if (operation == Operation.CONVOLUTION) {
            i2 = i2.clone();
            KernelFilter.rotate180(i2);
        }
        for (int i = 0; i < input1.length; i++) {
            sum += input1[i] * i2[i];
        }
        // double exp = e[size / 2 * size + size / 2];
        // logger.fine(() -> String.format("Sum = %f vs [%d] %f", sum, size / 2 * size + size / 2,
        // exp);
        Assertions.assertEquals(sum, sum, 1e-3);
    }
    // Test the FHT filter
    final FhtFilter ff = new FhtFilter(input2, size, size);
    ff.setOperation(operation);
    ff.filter(input1, size, size);
    // There may be differences due to the use of the JTransforms library
    final double error = (operation == Operation.DECONVOLUTION) ? 5e-2 : 1e-4;
    final FloatFloatBiPredicate predicate = TestHelper.floatsAreClose(error, 0);
    // This tests everything and can fail easily depending on the random generator
    // due to edge artifacts.
    // TestAssertions.assertArrayTest(e, input1, TestHelper.almostEqualFloats(error, 0));
    // This tests the centre to ignore edge differences
    final int min = size / 4;
    final int max = size - min;
    int repeats = 0;
    for (int y = min; y < max; y++) {
        for (int x = min; x < max; x++) {
            repeats++;
        }
    }
    // Use a fail counter for a 'soft' test that detects major problems
    final int failureLimit = TestCounter.computeFailureLimit(repeats, 0.1);
    final TestCounter failCounter = new TestCounter(failureLimit);
    final IndexSupplier msg = new IndexSupplier(2);
    for (int y = min; y < max; y++) {
        msg.set(1, y);
        for (int x = min; x < max; x++) {
            final int xx = x;
            final int i = y * size + x;
            failCounter.run(() -> {
                TestAssertions.assertTest(e[i], input1[i], predicate, msg.set(0, xx));
            });
        }
    }
}