Add Cooley-Tukey discrete FFT implementation utility

to later analyze the quality of sampling algorithms for path tracing
by their resulting image frequency spectrum.

Ideally, the resulting image should not contain
lower-frequency components, because the human eye
is very sensitive to those.

For example, blue noise contains mostly high-frequency components.

src/org/lwjgl/demo/util/FFT.java

@@ -0,0 +1,134 @@
+/*
+ * Copyright LWJGL. All rights reserved.
+ * License terms: https://www.lwjgl.org/license
+ */
+package org.lwjgl.demo.util;
+
+import javax.imageio.ImageIO;
+import java.awt.image.BufferedImage;
+import java.awt.image.WritableRaster;
+
+/**
+ * Implementation of Cooley-Tukey radix-2 Decimation-In-Time (DIT) discrete Fourier transform (DFT) algorithm
+ * to analyze the frequency spectrum of an image, such as those produced by sampling in path tracing.
+ * <p>
+ * This can be used to assess the quality of a sampling algorithm by analyzing the frequency spectrum of the resulting
+ * image. Ideally, it should not contain any low-frequency components, which the human eye is very sensitive to.
+ *
+ * @author Kai Burjack
+ */
+public class FFT {
+
+    private static double[][][] transpose(double[][][] m) {
+        int N = m.length;
+        double[][][] res = new double[N][N][];
+        for (int i = 0; i < N; i++)
+            for (int j = 0; j < N; j++)
+                res[i][j] = m[j][i];
+        return res;
+    }
+
+    private static double[][][] dfft2D(double[][][] input) {
+        int N = input.length;
+        for (int i = 0; i < N; i++)
+            input[i] = dfft1D(input[i]);
+        double[][][] t = transpose(input);
+        for (int i = 0; i < N; i++)
+            t[i] = dfft1D(t[i]);
+        return transpose(t);
+    }
+
+    private static double[][] dfft1D(double[][] input) {
+        int N = input.length;
+        if (N == 1) return input;
+        double[][] e = new double[N>>>1][], o = new double[N>>>1][];
+        for (int i = 0; i < N>>>1; i++) {
+            e[i] = input[i<<1];
+            o[i] = input[(i<<1) + 1];
+        }
+        double[][] q = dfft1D(e), v = dfft1D(o), r = new double[N][];
+        for (int i = 0; i < N>>>1; i++) {
+            double k = (i<<1) * Math.PI / N;
+            double[] c = new double[]{org.joml.Math.cos(k), org.joml.Math.sin(k)};
+            r[i] = new double[]{q[i][0] + c[0] * v[i][0] - c[1] * v[i][1], q[i][1] + c[0] * v[i][1] + c[1] * v[i][0]};
+            r[i + (N>>>1)] = new double[]{q[i][0] - (c[0] * v[i][0] - c[1] * v[i][1]), q[i][1] - (c[0] * v[i][1] + c[1] * v[i][0])};
+        }
+        return r;
+    }
+
+    private static double[][][] toComplexArray(double[][] input) {
+        int N = input.length;
+        double[][][] output = new double[N][N][];
+        for (int i = 0; i < N; i++)
+            for (int j = 0; j < N; j++)
+                output[i][j] = new double[]{input[i][j], 0};
+        return output;
+    }
+
+    private static double[][][] dfftShifted(double[][][] input) {
+        int N = input.length;
+        double[][][] res = new double[N][N][];
+        int hN = N>>>1;
+        for (int i = 0; i < hN; i++)
+            for (int j = 0; j < hN; j++) {
+                res[i + hN][j + hN] = input[i][j];
+                res[i][j] = input[i + hN][j + hN];
+                res[i + hN][j] = input[i][j + hN];
+                res[i][j + hN] = input[i + hN][j];
+            }
+        return res;
+    }
+
+    /**
+     * Compute the frequency spectrum of the given input image, which is the magnitude of the 2D DFT of the input image.
+     *
+     * @param input
+     *         the input image
+     * @return the frequency spectrum of the input image
+     */
+    public static byte[][] frequencySpectrum(double[][] input) {
+        double[][][] shifted = dfftShifted(dfft2D(toComplexArray(input)));
+        int N = shifted.length;
+        double max = 0;
+        double[][] mag = new double[N][N];
+        for (int i = 0; i < N; i++)
+            for (int j = 0; j < N; j++) {
+                mag[i][j] = Math.hypot(shifted[i][j][0], shifted[i][j][1]);
+                if (mag[i][j] > max)
+                    max = mag[i][j];
+            }
+        byte[][] res = new byte[N][N];
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                int value = (int) (255 * (Math.log1p(mag[i][j]) / Math.log1p(max)));
+                res[i][j] = (byte) (value & 0xFF);
+            }
+        }
+        return res;
+    }
+
+    // Test
+
+    private static double[][] loadFromResource(String resource) throws Exception {
+        BufferedImage image = ImageIO.read(FFT.class.getClassLoader().getResourceAsStream(resource));
+        int N = image.getWidth();
+        double[][] output = new double[N][N];
+        for (int i = 0; i < N; i++)
+            for (int j = 0; j < N; j++)
+                output[i][j] = (image.getRGB(i, j) & 0xFF) / 255.0;
+        return output;
+    }
+
+    public static void main(String[] args) throws Exception {
+        double[][] input = loadFromResource("org/lwjgl/demo/opengl/raytracing/tutorial4_2/blueNoise.png");
+        byte[][] spectrum = frequencySpectrum(input);
+        int N = spectrum.length;
+        BufferedImage res = new BufferedImage(N, N, BufferedImage.TYPE_BYTE_GRAY);
+        WritableRaster raster = res.getRaster();
+        for (int i = 0; i < N; i++)
+            for (int j = 0; j < N; j++)
+                raster.setSample(i, j, 0, spectrum[i][j]);
+        javax.imageio.ImageIO.write(res, "PNG", new java.io.File("spectrum.png"));
+    }
+
+}