Rewrote dip::GetOptimalDFTSize and dip::OptimalFourierTransformSize.

changelogs/diplib_next.md

@@ -45,6 +45,15 @@ title: "Changes DIPlib 3.x.x"
 - `dip::ImageRead()` throws a `dip::RunTimeError` instead of a `dip::ParameterError` if the file doesn't exist
   or is of an unrecognized type.
 
+- `dip::GetOptimalDFTSize()` has a new argument `maxFactor` that determines what the largest factor for the output
+  integer can be. Valid values are 5, 7 and 11. A new function `dip::MaxFactor()` will give the appropriate value
+  depending on whether a real or complex-valued transform is being computed, and what FFT implementation is being
+  used. When using *PocketFFT*, this function can now also return much larger values than before.
+
+- `dip::OptimalFourierTransformSize()`, which depends on `dip::GetOptimalDFTSize()` (see the bullet point above
+  this one) has a new argument `purpose` that can be either `"real"` or `"complex"`, and uses the new
+  `dip::MaxFactor()` to fill in the `maxFactor` argument.
+
 - Minimum required version of *CMake* is now 3.12.
 
 ### Bug fixes

include/diplib/dft.h

@@ -1,5 +1,5 @@
 /*
- * (c)2017-2022, Cris Luengo.
+ * (c)2017-2024, Cris Luengo.
  *
  * Encapsulates code Taken from OpenCV 3.1: (c)2000, Intel Corporation.
  * (see src/transform/opencv_dxt.cpp)
@@ -33,7 +33,10 @@
 namespace dip {
 
 
-/// \addtogroup transform
+/// \group transform_lowlevel Low-level transform support
+/// \ingroup transform
+/// \brief Low-level functionality for computing the Discrete Fourier Transform.
+/// \addtogroup
 
 
 namespace Option {
@@ -124,6 +127,7 @@ class DFT {
       /// buffers with this length. If `inverse` is `true`, an inverse transform will be computed.
       ///
       /// `options` determines some properties for the algorithm that will compute the DFT.
+      ///
       ///  - \ref Option::DFTOption::InPlace means the input and output pointers passed to \ref Apply must be the same.
       ///  - \ref Option::DFTOption::TrashInput means that the algorithm is free to overwrite the input array.
       ///    Ignored when working in place.
@@ -154,7 +158,7 @@ class DFT {
       /// regions of memory. When using FFTW, the `inplace` parameter to the constructor or \ref Initialize
       /// must be `true` if the two pointers here are the same, or `false` if they are different.
       ///
-      /// The input array is not marked `const`. If \ref Option::DFTOption::TrashInput` is given when planning,
+      /// The input array is not marked `const`. If \ref Option::DFTOption::TrashInput is given when planning,
       /// the input array can be overwritten with intermediate data, but otherwise will be left intact.
       ///
       /// `scale` is a real scalar that the output values are multiplied by. It is typically set to `1/size` for
@@ -280,6 +284,7 @@ class RDFT {
       /// The complex buffer has a size of `size/2+1`.
       ///
       /// `options` determines some properties for the algorithm that will compute the DFT.
+      ///
       ///  - \ref Option::DFTOption::InPlace means the input and output pointers passed to \ref Apply must be the same.
       ///    Do note that the complex array has one or two floats more than the real array, the buffer must be large
       ///    enough.
@@ -345,15 +350,21 @@ class RDFT {
 
 
 /// \brief Returns a size equal or larger to `size0` that is efficient for the DFT implementation.
+/// The value returned is a product of small primes.
 ///
 /// Set `larger` to false to return a size equal or smaller instead.
 ///
-/// Returns 0 if `size0` is too large for the DFT implementation.
+/// `maxFactor` should be 5, 7 or 11, and gives the largest integer that should be considered a "small prime".
+/// When using FFTW, set `maxFactor` to 7, as by default it works most efficiently for sizes that can be factored
+/// into primes smaller or equal to 7. When using PocketFFT, set it to 5 for complex-to-real or real-to-complex
+/// transforms, and to 11 for complex-to-complex transforms. See \ref dip::MaxFactor.
+///
+/// If there is no suitable size that the implementation can use (see \ref maximumDFTSize), this function
+/// will return 0.
 ///
 /// Prefer to use \ref dip::OptimalFourierTransformSize in your applications, it will throw an error if
-/// the transform size is too large.
-DIP_EXPORT dip::uint GetOptimalDFTSize( dip::uint size0, bool larger = true );
-
+/// the transform size is too large, and determines `maxFactor` for you.
+DIP_EXPORT dip::uint GetOptimalDFTSize( dip::uint size0, bool larger = true, dip::uint maxFactor = 5 );
 
 /// \brief The largest size supported by \ref DFT and \ref FourierTransform. Is equal to 2^31^-1 when using FFTW,
 /// or 2^64^-1 when using PocketFFT.
@@ -362,6 +373,9 @@ DIP_EXPORT extern dip::uint const maximumDFTSize;
 /// \brief Is `true` if `dip::DFT` and `dip::RDFT` use the FFTW library, or false if they use PocketFFT.
 DIP_EXPORT extern bool const usingFFTW;
 
+/// \brief The `maxFactor` parameter for \ref GetOptimalDFTSize. `complex` determines whether the transform to
+/// be computed is complex-to-complex or not.
+DIP_EXPORT dip::uint MaxFactor( bool complex = true );
 
 /// \endgroup
 

include/diplib/library/stringparams.h

@@ -173,6 +173,9 @@ constexpr char const* REAL = "real";
 constexpr char const* SYMMETRIC = "symmetric";
 constexpr char const* CORNER = "corner";
 //constexpr char const* FAST = "fast";
+constexpr char const* LARGER = "larger";
+constexpr char const* SMALLER = "smaller";
+constexpr char const* COMPLEX = "complex";
 
 // Distance transforms
 //constexpr char const* FAST = "fast";

include/diplib/linear.h

@@ -150,7 +150,7 @@ DIP_NODISCARD inline Image SeparableConvolution(
 /// then a periodic boundary condition is imposed. This is the natural boundary condition for the method,
 /// the image will be Fourier transformed as-is. For other boundary conditions, the image will be padded before
 /// the transform is applied. The padding will extend the image by at least half the size of `filter` in all
-/// dimensions, and the padding will make the image size a multiple of 2, 3 and/or 5 so that it is cheaper
+/// dimensions, and the padding will make the image size a multiple of small integers so that it is cheaper
 /// to compute the Fourier transform (see \ref dip::OptimalFourierTransformSize).
 /// The output image will be cropped to the size of the input.
 ///

include/diplib/transform.h

@@ -1,5 +1,5 @@
 /*
- * (c)2017-2021, Cris Luengo.
+ * (c)2017-2024, Cris Luengo.
  * Based on original DIPlib code: (c)1995-2014, Delft University of Technology.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
@@ -66,7 +66,7 @@ namespace dip {
 ///   transform to be computed.
 /// - "real": assumes that the (complex) input is conjugate symmetric, and returns a real-valued
 ///   result. Can only be used together with "inverse".
-/// - "fast": pads the input to a "nice" size, multiple of 2, 3 and 5, which can be processed faster.
+/// - "fast": pads the input to a "nice" size (multiple of small integers), which can be processed faster.
 /// - "corner": sets the origin to the top-left corner of the image (both in the spatial and the
 ///   frequency domain). This yields a standard DFT (Discrete Fourier Transform).
 /// - "symmetric": the normalization is made symmetric, where both forward and inverse transforms
@@ -131,18 +131,26 @@ DIP_NODISCARD inline Image InverseFourierTransform(
    return out;
 }
 
-/// \brief Returns the next larger (or smaller) multiple of {2, 3, 5}, an image of this size is more
+/// \brief Returns the next larger (or smaller) multiple of small integers. An image of this size is more
 /// efficient for FFT computations.
 ///
-/// The largest value that can be returned is 2125764000
-/// (smaller than 2^31^-1, the largest possible value of an `int` on most platforms).
+/// Pad an image with zeros to the next larger size or crop the image to the next smaller size to
+/// improve FFT performance.
+///
+/// The largest value that can be returned depends on the FFT implementation used. When using *FFTW*,
+/// the largest FFT that can be computed is 2^31^-1 (the maximum value for an `int`). Otherwise any
+/// `dip::uint` can be returned. If the result is larger than the maximum FFT size, an exception
+/// will be thrown.
 ///
 /// By default, `which` is `"larger"`, in which case it returns the next larger value. Set it
 /// to `"smaller"` to obtain the next smaller value instead.
 ///
-/// Pad an image with zeros to the next larger size or crop the image to the next smaller size to
-/// improve FFT performance.
-DIP_EXPORT dip::uint OptimalFourierTransformSize( dip::uint size, dip::String const& which = "larger" );
+/// `purpose` is either `"real"` (the default) or `"complex"`, and should be set according to the
+/// type of transform that will be computed. This determines what are considered "small primes".
+/// For *PocketFFT* this changes depending on whether the transform is complex-to-complex
+/// (complex-valued transform), or is real-to-complex or complex-to-real (real-valued transform).
+/// See \ref dip::GetOptimalDFTSize.
+DIP_EXPORT dip::uint OptimalFourierTransformSize( dip::uint size, dip::String const& which = S::LARGER, dip::String const& purpose = S::REAL );
 
 
 /// \brief Computes the Riesz transform of a scalar image.

pydip/src/analysis.cpp

@@ -169,7 +169,7 @@ void init_analysis( py::module& m ) {
           "in"_a, "options"_a = dip::StringSet{}, "process"_a = dip::BooleanArray{} );
    m.def( "InverseFourierTransform", py::overload_cast< dip::Image const&, dip::Image&, dip::StringSet, dip::BooleanArray >( &dip::InverseFourierTransform ),
           "in"_a, py::kw_only(), "out"_a, "options"_a = dip::StringSet{}, "process"_a = dip::BooleanArray{} );
-   m.def( "OptimalFourierTransformSize", &dip::OptimalFourierTransformSize, "size"_a, "which"_a = "larger" );
+   m.def( "OptimalFourierTransformSize", &dip::OptimalFourierTransformSize, "size"_a, "which"_a = dip::S::LARGER, "purpose"_a = dip::S::REAL );
    m.def( "RieszTransform", py::overload_cast< dip::Image const&, dip::String const&, dip::String const&, dip::BooleanArray >( &dip::RieszTransform ),
           "in"_a, "inRepresentation"_a = dip::S::SPATIAL, "outRepresentation"_a = dip::S::SPATIAL, "process"_a = dip::BooleanArray{} );
    m.def( "RieszTransform", py::overload_cast< dip::Image const&, dip::Image&, dip::String const&, dip::String const&, dip::BooleanArray >( &dip::RieszTransform ),

src/DIPlib_sources.cmake

@@ -284,6 +284,7 @@ support/matrix.cpp
 support/numeric.cpp
 support/thin_plate_spline.cpp
 transform/dft.cpp
+transform/dft_optimal_size.cpp
 transform/fftw3api.h
 transform/fourier.cpp
 transform/haar.cpp

src/deconvolution/common_deconv_utility.h

@@ -56,12 +56,13 @@ inline void FourierTransformImageAndKernel(
    DIP_THROW_IF( pad && isOtf, E::ILLEGAL_FLAG_COMBINATION );
    if( pad || ( powersOfTwo > 0 )) {
       dip::uint multiple = static_cast< dip::uint >( std::pow( 2, powersOfTwo ));
+      dip::String purpose = in.DataType().IsComplex() ? S::COMPLEX : S::REAL;
       dip::UnsignedArray sizes = in.Sizes();
       for( dip::uint ii = 0; ii < nDims; ++ii ) {
          if( pad ) {
             sizes[ ii ] += 2 * psf.Size( ii );
          }
-         sizes[ ii ] = OptimalFourierTransformSize( div_ceil( sizes[ ii ], multiple )) * multiple;
+         sizes[ ii ] = OptimalFourierTransformSize( div_ceil( sizes[ ii ], multiple ), S::LARGER, purpose ) * multiple;
       }
       Image tmp = ExtendImageToSize( in, sizes, S::CENTER );
       DIP_STACK_TRACE_THIS( FourierTransform( tmp, G ));

src/deconvolution/richardson_lucy.cpp

@@ -71,7 +71,7 @@ void RichardsonLucy(
    if( pad ) {
       dip::UnsignedArray sizes = in.Sizes();
       for( dip::uint ii = 0; ii < nDims; ++ii ) {
-         sizes[ ii ] += OptimalFourierTransformSize( sizes[ ii ] + psf.Size( ii ) - 1 );
+         sizes[ ii ] += OptimalFourierTransformSize( sizes[ ii ] + psf.Size( ii ) - 1, S::LARGER, S::REAL );
       }
       g = ExtendImageToSize( in, sizes, S::CENTER );
    } else {

src/linear/convolution.cpp

@@ -440,17 +440,18 @@ void ConvolveFT(
 
    // Prepare input image
    bool inPadding = ( inSpatial && filterSpatial && outSpatial && !boundaryCondition.empty() );
-   UnsignedArray inSizes = in.Sizes();
+   UnsignedArray const& inSizes = in.Sizes();
    bool real = true;
    Image inFT;
    bool reuseInFT = false;
    if( inSpatial ) {
       real &= in.DataType().IsReal();
+      dip::String purpose = real ? S::COMPLEX : S::REAL;
       if( inPadding ) {
          // Pad the input image with at least the size of `filter`, but make it larger so it's a nice size
          UnsignedArray sizes = inSizes;
          for( dip::uint ii = 0; ii < sizes.size(); ++ii ) {
-            sizes[ ii ] = OptimalFourierTransformSize( sizes[ ii ] + filter.Size( ii ) - 1 );
+            sizes[ ii ] = OptimalFourierTransformSize( sizes[ ii ] + filter.Size( ii ) - 1, S::LARGER, purpose );
          }
          ExtendImageToSize( in, inFT, sizes, S::CENTER, boundaryCondition );
          DIP_STACK_TRACE_THIS( FourierTransform( inFT, inFT ));