FFT

CHANGELOG.md

@@ -11,13 +11,29 @@ Added / Changed / Deprecated / Fixed / Removed / Security
 
 > Corresponds to changes in the `develop` branch since the last release
 
+### Added
+
+#### org.ojalgo.data
+
+- There is a new package `org.ojalgo.data.transform` that (so far) contain implementations of the `ZTransform` and the `DiscreteFourierTransform`. Most notably ojAlgo now contains an implementation of the FFT algorithm.
+- `ImageData` has been extended with features to work with Fourier transforms of images, as well as a Gaussian blur function.
+
+#### org.ojalgo.function
+
+- The `PolynomialFunction` interface now extends `org.ojalgo.algebra.Ring` which means it is now possible to `add` and `multiply` polynomials. There's also been some refactoring of the internals of the polynomial implementations.
+- Two new `PrimitiveFunction.Unary` implementations `PeriodicFunction` and `FourierSeries`. It is now possible to take any (periodic) function and, via FFT, get a Fourier series approximation of it (1 method call).
+
+#### org.ojalgo.scalar
+
+- Some minor additions to `ComplexNumber`. It is now possible to do complex number multiplication without creating new instances (calculate the real and imaginary parts separately). Added utilities to get the unit roots.
+
 ## [53.1.1] – 2023-10-16
 
 ### Added
 
 #### org.ojalgo.data
 
-- New `ImageData` class that wraps a `java.awt.image.BufferedImage` and implements `MatrixStore`. Further it adds a few utility methods to simplify working image data - convert to grey scale, re-sample (change size), separate the colour channels...
+- New `ImageData` class that wraps a `java.awt.image.BufferedImage` and implements `MatrixStore`. Further it adds a few utility methods to simplify working with image data - convert to grey scale, re-sample (change size), separate the colour channels...
 
 #### org.ojalgo.matrix
 

src/main/java/org/ojalgo/data/DataProcessors.java

@@ -52,27 +52,27 @@ public class DataProcessors {
     /**
      * Variables centered so that their average will be 0.0
      */
-    public static final Transformation2D<Double> CENTER = DataProcessors.newTransformation2D(ss -> SUBTRACT.by(ss.getMean()));
+    public static final Transformation2D<Double> CENTER = DataProcessors.newColumnsTransformer(ss -> SUBTRACT.by(ss.getMean()));
 
     /**
      * Variables will be centered around 0.0 AND scaled to be [-1.0,1.0]. The minimum value will be
      * transformed to -1.0 and the maximum to +1.0.
      */
     public static final Transformation2D<Double> CENTER_AND_SCALE = DataProcessors
-            .newTransformation2D(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by((ss.getMaximum() - ss.getMinimum()) / TWO)));
+            .newColumnsTransformer(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by((ss.getMaximum() - ss.getMinimum()) / TWO)));
 
     /**
      * Variables scaled to be within [-1.0,1.0] (divide by largest magnitude regardless of sign). If all
      * values are positive the range will within [0.0,1.0]. If all are negative the range will be within
      * [-1.0,0.0]
      */
-    public static final Transformation2D<Double> SCALE = DataProcessors.newTransformation2D(ss -> DIVIDE.by(ss.getLargest()));
+    public static final Transformation2D<Double> SCALE = DataProcessors.newColumnsTransformer(ss -> DIVIDE.by(ss.getLargest()));
 
     /**
      * Will normalise each variable - replace each value with its standard score.
      */
     public static final Transformation2D<Double> STANDARD_SCORE = DataProcessors
-            .newTransformation2D(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by(ss.getStandardDeviation())));
+            .newColumnsTransformer(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by(ss.getStandardDeviation())));
 
     /**
      * Calculate the correlation matrix from a set of variables' samples. Each {@link Access1D} instance
@@ -237,7 +237,20 @@ public static <M extends PhysicalStore<Double>> M covariances(final Factory2D<M>
         return retVal;
     }
 
-    public static Transformation2D<Double> newTransformation2D(final Function<SampleSet, UnaryFunction<Double>> definition) {
+    /**
+     * Creates a {@link Transformation2D} that will apply a {@link UnaryFunction} to each column. That unary
+     * function will be created by the provided {@link Function} using a {@link SampleSet} (of that column) as
+     * input.
+     * <p>
+     * The constants {@link #CENTER}, {@link #SCALE}, {@link #CENTER_AND_SCALE} and {@link #STANDARD_SCORE}
+     * are predefined {@link Transformation2D} instances created by calling this method.
+     * 
+     * @param definition A {@link Function} that will create a {@link UnaryFunction} from a {@link SampleSet}
+     *        to be applied to each column
+     * @return A {@link Transformation2D} that will apply a {@link UnaryFunction} to each column
+     */
+    public static Transformation2D<Double> newColumnsTransformer(final Function<SampleSet, UnaryFunction<Double>> definition) {
+
         return new Transformation2D<>() {
 
             public <T extends Mutate2D.ModifiableReceiver<Double>> void transform(final T transformable) {

src/main/java/org/ojalgo/data/batch/BatchNode.java

@@ -51,6 +51,9 @@
  * A batch processing data node for when there's no way to fit the data in memory.
  * <p>
  * Data is stored in sharded files, and data is written/consumed and processed concurrently.
+ * <p>
+ * The data is processed in batches. Each batch is processed in a single thread. The number of threads is
+ * controlled by {@link #parallelism(IntSupplier)}.
  */
 public final class BatchNode<T> {
 
@@ -182,18 +185,22 @@ private static final class TwoStepWrapper<T> implements TwoStepMapper<T, Boolean
             myActualConsumer = consumerFactory.get();
         }
 
+        @Override
         public void consume(final T item) {
             myActualConsumer.accept(item);
         }
 
+        @Override
         public Boolean getResults() {
             return Boolean.TRUE;
         }
 
+        @Override
         public void merge(final Boolean aggregate) {
             // No need to (not possible to) merge, just continue
         }
 
+        @Override
         public void reset() {
             // No need to (not possible to) reset, just continue
         }

src/main/java/org/ojalgo/data/image/ImageData.java

@@ -22,20 +22,30 @@
 package org.ojalgo.data.image;
 
 import java.awt.Graphics;
+import java.awt.Graphics2D;
+import java.awt.RenderingHints;
 import java.awt.image.BufferedImage;
+import java.awt.image.ConvolveOp;
+import java.awt.image.Kernel;
 import java.io.File;
 import java.io.IOException;
 
 import javax.imageio.ImageIO;
 
+import org.ojalgo.data.transform.DiscreteFourierTransform;
+import org.ojalgo.function.aggregator.Aggregator;
 import org.ojalgo.function.constant.PrimitiveMath;
 import org.ojalgo.function.special.MissingMath;
 import org.ojalgo.matrix.store.MatrixStore;
+import org.ojalgo.matrix.store.PhysicalStore;
 import org.ojalgo.matrix.store.PhysicalStore.Factory;
 import org.ojalgo.matrix.store.Primitive32Store;
+import org.ojalgo.matrix.store.Primitive64Store;
+import org.ojalgo.scalar.ComplexNumber;
 import org.ojalgo.structure.Access2D;
 import org.ojalgo.structure.Mutate2D;
 import org.ojalgo.structure.Structure2D;
+import org.ojalgo.structure.Transformation2D;
 import org.ojalgo.type.NumberDefinition;
 
 /**
@@ -53,7 +63,22 @@
  * it has 256 shades of grey. Each of the sliced colour channels have 256 shades of their respective colour.
  * </ol>
  */
-public class ImageData implements MatrixStore<Double>, Mutate2D.Fillable<Double> {
+public class ImageData implements MatrixStore<Double>, Mutate2D.Receiver<Double> {
+
+    @FunctionalInterface
+    public static interface FrequencyDomainUpdater {
+
+        /**
+         * Used as a callback from {@link #newTransformation(FrequencyDomainUpdater)}.
+         * 
+         * @param distance The distance from the centre. The centre is the zero frequency. The distance is
+         *        scaled so that the radius of the largest circle that can fit in the image rectangle is 100.
+         * @param value The current value
+         * @return The new value
+         */
+        ComplexNumber invoke(double distance, ComplexNumber value);
+
+    }
 
     static final class SingleChannel extends ImageData {
 
@@ -102,6 +127,20 @@ public static ImageData copy(final Access2D<?> values) {
         return retVal;
     }
 
+    /**
+     * Creates a new image, transforming the input (back) from the frequency domain to the spatial domain
+     * using the inverse discrete Fourier transform. The transformation also reverts the shift of the
+     * frequency representation – the input should be in a format as returned by {@link #toFrequencyDomain()}.
+     * The new image will be of the same dimensions as the input matrix.
+     *
+     * @see #toFrequencyDomain()
+     */
+    public static ImageData fromFrequencyDomain(final MatrixStore<ComplexNumber> transformed) {
+        ImageData retVal = ImageData.newGreyScale(transformed);
+        retVal.fillMatching(DiscreteFourierTransform.inverse2D(DiscreteFourierTransform.shift(transformed)));
+        return retVal;
+    }
+
     public static ImageData newColour(final int nbRows, final int nbCols) {
         return new ImageData(new BufferedImage(nbCols, nbRows, BufferedImage.TYPE_INT_ARGB));
     }
@@ -118,6 +157,55 @@ public static ImageData newGreyScale(final Structure2D shape) {
         return ImageData.newGreyScale(shape.getRowDim(), shape.getColDim());
     }
 
+    /**
+     * Creates a new transformation that can be used to transform a matrix of complex numbers in the frequency
+     * domain. The transformation will invoke the updater for each element in the input matrix, passing the
+     * distance from the centre of the matrix as the first argument. The distance is scaled so that the radius
+     * of the largest circle that can fit in the image rectangle is 100.
+     */
+    public static Transformation2D<ComplexNumber> newTransformation(final FrequencyDomainUpdater updater) {
+
+        return new Transformation2D<>() {
+
+            public <T extends Mutate2D.ModifiableReceiver<ComplexNumber>> void transform(final T transformable) {
+
+                int nbRows = transformable.getRowDim();
+                int nbCols = transformable.getColDim();
+
+                double midRow = (nbRows - 1) / 2D;
+                double midCol = (nbCols - 1) / 2D;
+
+                double scale = Math.min(midRow, midCol) / 100D;
+
+                for (int j = 0; j < nbCols; j++) {
+                    for (int i = 0; i < nbRows; i++) {
+                        transformable.set(i, j, updater.invoke(MissingMath.hypot(i - midRow, j - midCol) / scale, transformable.get(i, j)));
+                    }
+                }
+            }
+        };
+    }
+
+    /**
+     * Converts a matrix of complex numbers to an image of its power spectrum (log10 of the squared norms). In
+     * addition there is scaling applied to adapt the values towards the range [0,255].
+     * <p>
+     * This is meant to be used with the output from {@link #toFrequencyDomain()} to get a visual
+     * representation of the frequency content of an image.
+     */
+    public static ImageData ofPowerSpectrum(final Access2D<ComplexNumber> transformed) {
+
+        PhysicalStore<Double> magnitudes = Primitive64Store.getComplexModulus(transformed);
+
+        magnitudes.modifyAll(PrimitiveMath.POWER.parameter(2).andThen(PrimitiveMath.LOG10));
+
+        double average = magnitudes.aggregateAll(Aggregator.AVERAGE).doubleValue();
+
+        magnitudes.modifyAll(PrimitiveMath.MULTIPLY.by(127.5D / average));
+
+        return ImageData.copy(magnitudes);
+    }
+
     public static ImageData read(final File file) {
         try {
             return new ImageData(ImageIO.read(file));
@@ -126,10 +214,37 @@ public static ImageData read(final File file) {
         }
     }
 
+    public static ImageData read(final File directory, final String fileName) {
+        return ImageData.read(new File(directory, fileName));
+    }
+
     public static ImageData wrap(final BufferedImage image) {
         return new ImageData(image);
     }
 
+    private static Kernel getGaussianKernel(final double sigma, final int radius) {
+        int size = radius * 2 + 1;
+        float[] data = new float[size];
+
+        double twoSigmaSquared = 2.0 * sigma * sigma;
+        double sigmaRoot = Math.sqrt(twoSigmaSquared * Math.PI);
+        double total = 0.0;
+
+        for (int i = -radius; i <= radius; i++) {
+            double distance = i * i;
+            int index = i + radius;
+            data[index] = (float) (Math.exp(-distance / twoSigmaSquared) / sigmaRoot);
+            total += data[index];
+        }
+
+        // Normalize the kernel
+        for (int i = 0; i < data.length; i++) {
+            data[i] /= total;
+        }
+
+        return new Kernel(size, 1, data);
+    }
+
     static int toRanged(final int value) {
         return Math.max(0, Math.min(value, 255));
     }
@@ -141,12 +256,37 @@ static int toRanged(final int value) {
         myImage = image;
     }
 
+    /**
+     * Creates a new image (of the same type) blurring the input using a Gaussian blur kernel.
+     * 
+     * @param sigma The standard deviation of the Gaussian blur kernel
+     */
+    public ImageData applyGaussianBlur(final double sigma) {
+
+        int radius = (int) (3.0 * sigma);
+        if (radius % 2 == 0) {
+            radius++; // Ensure radius is odd
+        }
+
+        Kernel kernel = ImageData.getGaussianKernel(sigma, radius);
+        ConvolveOp op = new ConvolveOp(kernel, ConvolveOp.EDGE_NO_OP, null);
+
+        BufferedImage blurredImage = new BufferedImage(myImage.getWidth(), myImage.getHeight(), myImage.getType());
+        Graphics2D g2d = blurredImage.createGraphics();
+        g2d.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
+        g2d.drawImage(myImage, op, 0, 0);
+        g2d.dispose();
+
+        return new ImageData(blurredImage);
+    }
+
     @Override
     public byte byteValue(final int row, final int col) {
         return (byte) this.intValue(row, col);
     }
 
     /**
+     * Creates a new image, converting the input to the specified image type in the process.
      * {@link BufferedImage#TYPE_BYTE_GRAY}, {@link BufferedImage#TYPE_INT_ARGB} or any other valid type...
      */
     public ImageData convertTo(final int imageType) {
@@ -207,6 +347,13 @@ public int getColDim() {
         return myImage.getWidth();
     }
 
+    /**
+     * @return The underlying {@link BufferedImage}
+     */
+    public BufferedImage getImage() {
+        return myImage;
+    }
+
     @Override
     public int getRowDim() {
         return myImage.getHeight();
@@ -337,6 +484,20 @@ public ImageData sliceRedChannel() {
         return new SingleChannel(myImage, MASK_RED, SHIFT_RED);
     }
 
+    /**
+     * Transforms the spatial representation of the image to its frequency representation using the discrete
+     * Fourier transform.
+     * <p>
+     * In addition the frequency representation is shifted so that the zero frequency is in the centre of the
+     * image.
+     * <p>
+     * When you're done processing the frequency representation, you can transform it back to a spatial
+     * representation using {@link #fromFrequencyDomain(MatrixStore)}.
+     */
+    public PhysicalStore<ComplexNumber> toFrequencyDomain() {
+        return DiscreteFourierTransform.shift(DiscreteFourierTransform.transform2D(this)).copy();
+    }
+
     /**
      * The file format is derived from the file name ending (png, jpg...)
      */
@@ -350,6 +511,10 @@ public void writeTo(final File file) {
         }
     }
 
+    public void writeTo(final File directory, final String fileName) {
+        this.writeTo(new File(directory, fileName));
+    }
+
     private double doubleValue(final double row, final double col) {
         return this.doubleValue(MissingMath.roundToInt(row), MissingMath.roundToInt(col));
     }