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"));
+ }
+
+}