Remove FFTConvolution class

[ctrueden]

The FFTConvolution class was copied from ImgLib2 wholesale. This class was taken from imglib2-algorithms-gpl, which is licensed under the GPL. Hence, this class is incompatible with ImageJ OPS's BSD-2 licensing.

To address this issue in the short to medium term, we can split the imagej.ops.convolve package into its own component, ij-ops-gpl, which depends on imglib2-algorithms-gpl. But ideally, we don't want any part of OPS to be licensed under the GPL; better would be to find a way to do convolution without that code.

See also imglib/imglib2#61

[ctrueden]

Pushed from 0.1.0 to 0.2.0, because this isn't really a blocker for the initial release announcement. The announcement just needs to warn that the development is still very active and subject to restructuring within the next weeks.

[bnorthan]

Hi Curtis

My understanding, after talking to you about it at the hackathon, is that there are 2 parts to this.

1. FFTConvolution is GPL. I think we would want a slightly different class hierarchy anyway. Convolution could be built from the fft op. High level convolution is not to hard to write (ffts and a complex multiply).
2. Mines FFT is CPL. This can be solved by using a different 1d FFT routine in the imglib2 FFTMethods class. ND-FFT is implemented by applying Mines FFT line by line. There are only a few calls to Mines that would have to be replaced.

I looked at JTransforms and the only thing I couldn't figure out is how to get the efficient sizes. FFTs work fastest at certain sizes (power of 2 is best, but other sizes work also) and you need to know these sizes to set up the padding correctly. I could not figure out how to do this with the JTransforms API.

@StephanPreibisch do you have any plans for this??

[dscho]

I looked at JTransforms and the only thing I couldn't figure out is how to get the efficient sizes.

AFAICT JTransforms only handles powers of two really efficiently, everything else is handled by the exact same code:

https://github.com/wendykierp/JTransforms/blob/ede249c9824262bd9b5f571e56f7a0fa596f6f20/src/main/java/org/jtransforms/fft/DoubleFFT_2D.java#L117

There is another issue we might have to hash out with @wendykierp at some stage: due to the static way nthreads is determined the code is not cluster-safe (where you would need to specify on a per-Context level how many threads are permissible, as opposed to setting it globally).

[wendykierp]

There are three FFT algorithms implemented in JTransforms: split-radix (only for powers of two), mixed-radix (2-, 3-, 4-, 5-, n-) and Bluestein (used when the reminder from the integer division of input size by 2, 3, 4 and 5 is greater or equal 211):

https://github.com/wendykierp/JTransforms/blob/ede249c9824262bd9b5f571e56f7a0fa596f6f20/src/main/java/org/jtransforms/fft/DoubleFFT_1D.java#L117

[bnorthan]

Thanks @dscho and @wendykierp. So then can we use this approach http://ltfat.sourceforge.net/notes/ltfatnote017.pdf to find the efficient sizes in between powers of 2??

[dscho]

@bnorthan interesting, but I think a table lookup with 2^20 entries is overkill.

Given a number N and a set of prime numbers p_i you should be able to set up a linear programming algorithm that minimizes \sum_i (a_i \cdot\log p_i) - \log N under the constraint that it has to be non-negative.

[bnorthan]

I was thinking of a smaller lookup table... but that equation looks like it would work as well.

[dscho]

@bnorthan see also 5.7 of http://www.math.leidenuniv.nl/~psh/ANTproc/09andrew.pdf

[dscho]

@bnorthan I ran out of time, but I did implement something you might find useful:

package org.scijava.math;

import java.util.Arrays;

/**
 * A class to determine the next <i>k</i>-smooth number (i.e. a number divisible
 * only by prime numbers up to <i>k</i>), given a lower bound <i>n</i>. Based on
 * A. Granville, <i>Smooth numbers: computational number theory and beyond</i>:
 * <blockquote> 5.7. <b>Finding smooth numbers computationally.</b> The obvious
 * way to find y-smooth numbers in (x; x + z) with z &le; x is to initialize an
 * array a[i] := 0 for 1 &le; i &le; z (where a[i] corresponds to x + i). For
 * each successively larger prime power p<sup>j</sup> &le; x + z with p &le; y,
 * determine the smallest i such that p<sup>j</sup> divides x + i and then add
 * log(p) to a[i], a[i + p<sup>j</sup>], a[i + 2p<sup>j</sup>] and so on, up
 * until the end of the array. When we’ve finished, if any a[i] &ge; log(x),
 * then x + i is y-smooth. </blockquote>
 * 
 * @author Johannes Schindelin
 */
public class NextSmoothNumber {
    public static int nextSmooth(final int y, final int x, final int z) {
        final double[] a = new double[z];
        handlePrime(2, x, a);
        handlePrime(3, x, a);
        handlePrime(5, x, a);
        handlePrime(7, x, a);
        double log = Math.log(x);
        for (int i = 0; i < a.length; i++) {
            if (a[i] >= log) return x + i;
        }
        System.err.println(Arrays.toString(a));
        System.err.println(a[32805 - x]);
        return -1;
    }

    private static void handlePrime(final int p, final int x, final double[] a) {
        double log = Math.log(p);
        for (int power = p; power <= x + a.length; power *= p) {
            int j = x % power;
            if (j > 0) j = power - j;
            while (j < a.length) {
                a[j] += log;
                j += power;
            }
        }
    }

    private static int log2(int x) {
        for (int j = 0, k = 1; ; j++, k *= 2) {
            if (k >= x) return j;
        }
    }

    public static void main(String... args) {
        int x = args.length == 0 ? 32769 : Integer.parseInt(args[0]);
        System.err.println(log2(x));
        System.out.println(nextSmooth(7, x, 4 * log2(x)));
    }
}

What I wanted to do is to run a loop from 32000 to 33000 and verify that the produced numbers are 7-smooth...

[bnorthan]

Nice. I'll take a look at it and see if I can finish off the test.