Create half Gaussian kernel (#206)

Expose internally MakeHalfGaussian and update CreateGauss to optionally only generate half of a Gaussian kernel

Author: VisualFox

changelogs/diplib_next.md

@@ -38,6 +38,9 @@ date: 2020-00-00
   The low-level functions `dip::SymmetricEigenDecomposition2()` and `dip::SymmetricEigenDecomposition3()` add
   the same argument, but as an enum rather than a string.
 
+- `dip::CreateGauss()` has a new argument, `extent`, which defaults to `"full"`. When set to `"half"`, the output
+  is the first half (in 1D) or quadrant (in 2D) of the Gaussian kernel.
+
 ### Changed functionality
 
 - All functions that compute a percentile (`dip::Percentile()`, `dip::PercentilePosition()`,

dipimage/derivative.m

@@ -36,7 +36,7 @@
 %
 % DIPlib:
 %  This function calls the DIPlib function <a href="https://diplib.org/diplib-docs/linear.html#dip-Derivative-Image-CL-Image-L-UnsignedArray--FloatArray--String-CL-StringArray-CL-dfloat-">dip::Derivative</a> (which calls <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussFIR-Image-CL-Image-L-FloatArray--UnsignedArray--StringArray-CL-dfloat-">dip::GaussFIR</a>,
-%  <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussIIR-Image-CL-Image-L-FloatArray--UnsignedArray--StringArray-CL-UnsignedArray--String-CL-dfloat-">dip::GaussIIR</a>, <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussFT-Image-CL-Image-L-FloatArray--UnsignedArray--dfloat--String-CL-String-CL-StringArray-CL">dip::GaussFT</a> and <a href="https://diplib.org/diplib-docs/linear.html#dip-FiniteDifference-Image-CL-Image-L-UnsignedArray--String-CL-StringArray-CL-BooleanArray-">dip::FiniteDifference</a>) and <a href="https://diplib.org/diplib-docs/generation_test.html#dip-CreateGauss-Image-L-FloatArray-CL-UnsignedArray--dfloat--UnsignedArray-">dip::CreateGauss</a>.
+%  <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussIIR-Image-CL-Image-L-FloatArray--UnsignedArray--StringArray-CL-UnsignedArray--String-CL-dfloat-">dip::GaussIIR</a>, <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussFT-Image-CL-Image-L-FloatArray--UnsignedArray--dfloat--String-CL-String-CL-StringArray-CL">dip::GaussFT</a> and <a href="https://diplib.org/diplib-docs/linear.html#dip-FiniteDifference-Image-CL-Image-L-UnsignedArray--String-CL-StringArray-CL-BooleanArray-">dip::FiniteDifference</a>) and <a href="https://diplib.org/diplib-docs/generation_test.html#dip-CreateGauss-Image-L-FloatArray-CL-UnsignedArray--dfloat--UnsignedArray--String-CL">dip::CreateGauss</a>.
 
 % (c)2017-2018, Cris Luengo.
 % Based on original DIPlib code: (c)1995-2014, Delft University of Technology.

dipimage/gaussf.m

@@ -32,7 +32,7 @@
 %
 % DIPlib:
 %  This function calls the DIPlib function <a href="https://diplib.org/diplib-docs/linear.html#dip-Derivative-Image-CL-Image-L-UnsignedArray--FloatArray--String-CL-StringArray-CL-dfloat-">dip::Derivative</a> (which calls <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussFIR-Image-CL-Image-L-FloatArray--UnsignedArray--StringArray-CL-dfloat-">dip::GaussFIR</a>,
-%  <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussIIR-Image-CL-Image-L-FloatArray--UnsignedArray--StringArray-CL-UnsignedArray--String-CL-dfloat-">dip::GaussIIR</a> and <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussFT-Image-CL-Image-L-FloatArray--UnsignedArray--dfloat--String-CL-String-CL-StringArray-CL">dip::GaussFT</a>) and <a href="https://diplib.org/diplib-docs/generation_test.html#dip-CreateGauss-Image-L-FloatArray-CL-UnsignedArray--dfloat--UnsignedArray-">dip::CreateGauss</a>.
+%  <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussIIR-Image-CL-Image-L-FloatArray--UnsignedArray--StringArray-CL-UnsignedArray--String-CL-dfloat-">dip::GaussIIR</a> and <a href="https://diplib.org/diplib-docs/linear.html#dip-GaussFT-Image-CL-Image-L-FloatArray--UnsignedArray--dfloat--String-CL-String-CL-StringArray-CL">dip::GaussFT</a>) and <a href="https://diplib.org/diplib-docs/generation_test.html#dip-CreateGauss-Image-L-FloatArray-CL-UnsignedArray--dfloat--UnsignedArray--String-CL">dip::CreateGauss</a>.
 
 % (c)2017-2018, Cris Luengo.
 % Based on original DIPlib code: (c)1995-2014, Delft University of Technology.

include/diplib/generation.h

@@ -628,6 +628,19 @@ DIP_NODISCARD inline Image CreateDelta( UnsignedArray const& sizes, String const
    return out;
 }
 
+// Create a half Gaussian, used by MakeGaussian in linear/gauss.cpp (where it is defined)
+// Not documented because not a public function. We suggest you use `CreateGauss()` instead.
+// Length will be given by truncation and sigma, but limited to meaningful values.
+//
+// !!! warning
+//     The second half of the gaussian will need to be scaled by -1.0 for odd derivative order (ex 1, 3, ...)
+// Defined in src/linear/gauss.cpp
+DIP_EXPORT std::vector< dfloat > MakeHalfGaussian(
+      dfloat sigma,
+      dip::uint derivativeOrder = 0,
+      dfloat truncation = 3.0,
+      DataType dt = DT_DFLOAT  // if not DT_DFLOAT, assumed to be DT_SFLOAT
+);
 
 // Create 1D Gaussian, used in linear/gauss.cpp (where it is defined) and in nonlinear/bilateral.cpp
 // Not documented because not a public function. We suggest you use `CreateGauss()` instead.
@@ -654,6 +667,10 @@ DIP_EXPORT std::vector< dfloat > MakeGaussian(
 /// By setting `exponents` to a positive value for each dimension, the created kernel will be multiplied by
 /// the coordinates to the power of `exponents`.
 ///
+/// `extent` defaults to `"full"`. Set it to `"half"` to generate only the first half (along each dimension)
+/// of the kernel.
+/// The second half of the gaussian will need to be scaled by -1.0 for odd derivative order (ex 1, 3, ...)
+///
 /// !!! warning
 ///     Convolving an image with the result of this function is much less efficient than calling \ref Gauss.
 // Defined in src/linear/gauss.cpp
@@ -662,16 +679,18 @@ DIP_EXPORT void CreateGauss(
       FloatArray const& sigmas,
       UnsignedArray derivativeOrder = { 0 },
       dfloat truncation = 3.0,
-      UnsignedArray exponents = { 0 }
+      UnsignedArray exponents = { 0 },
+      String const& extent = "full"
 );
 DIP_NODISCARD inline Image CreateGauss(
       FloatArray const& sigmas,
       UnsignedArray derivativeOrder = { 0 },
       dfloat truncation = 3.0,
-      UnsignedArray exponents = { 0 }
+      UnsignedArray exponents = { 0 },
+      String const& extent = "full"
 ) {
    Image out;
-   CreateGauss( out, sigmas, std::move( derivativeOrder ), truncation, std::move( exponents ));
+   CreateGauss( out, sigmas, std::move( derivativeOrder ), truncation, std::move( exponents ), extent );
    return out;
 }
 

pydip/src/documentation_strings.h

@@ -147,6 +147,7 @@ constexpr char const* dip·viewer·Manager·createWindow·WindowPtr· = "Create
 constexpr char const* dip·viewer·Manager·activeWindows = "Returns the number of managed windows.";
 constexpr char const* dip·viewer·Manager·destroyWindows = "Destroys all windows.";
 constexpr char const* dip·viewer·Manager·processEvents = "Processes event queue.";
+constexpr char const* dip·viewer·Manager·screenSize·C = "Returns the size of the screen in pixels.";
 constexpr char const* dip·viewer·Manager·swapBuffers·Window·P = "Swap display buffers.";
 constexpr char const* dip·viewer·Manager·setWindowTitle·Window·P·char·CP = "Sets a Window's title.";
 constexpr char const* dip·viewer·Manager·refreshWindow·Window·P = "Refresh a Window's contents.";
@@ -1970,7 +1971,7 @@ constexpr char const* dip·GaussianLineClip·Image·CL·Image·L·Image·Pixel·
 constexpr char const* dip·FillDelta·Image·L·String·CL = "Fills an image with a delta function.";
 constexpr char const* dip·CreateDelta·Image·L·UnsignedArray·CL·String·CL = "Creates a delta function image.";
 constexpr char const* dip·CreateDelta·UnsignedArray·CL·String·CL = "Overload for the function above, which takes image sizes instead of an image.";
-constexpr char const* dip·CreateGauss·Image·L·FloatArray·CL·UnsignedArray··dfloat··UnsignedArray· = "Creates a Gaussian kernel.";
+constexpr char const* dip·CreateGauss·Image·L·FloatArray·CL·UnsignedArray··dfloat··UnsignedArray··String·CL = "Creates a Gaussian kernel.";
 constexpr char const* dip·CreateGabor·Image·L·FloatArray·CL·FloatArray·CL·dfloat· = "Creates a Gabor kernel.";
 constexpr char const* dip·FTEllipsoid·Image·L·FloatArray··dfloat· = "Generates the Fourier transform of an ellipsoid.";
 constexpr char const* dip·FTEllipsoid·UnsignedArray·CL·FloatArray··dfloat· = "Overload for the function above, which takes image sizes instead of an image.";

pydip/src/documentation_urls.py

@@ -148,6 +148,7 @@
     ('dip.viewer.Manager.activeWindows', 'dip-viewer-Manager.html#dip-viewer-Manager-activeWindows'),
     ('dip.viewer.Manager.destroyWindows', 'dip-viewer-Manager.html#dip-viewer-Manager-destroyWindows'),
     ('dip.viewer.Manager.processEvents', 'dip-viewer-Manager.html#dip-viewer-Manager-processEvents'),
+    ('dip.viewer.Manager.screenSize', 'dip-viewer-Manager.html#dip-viewer-Manager-screenSize-C'),
     ('dip.viewer.Manager.swapBuffers', 'dip-viewer-Manager.html#dip-viewer-Manager-swapBuffers-Window-P'),
     ('dip.viewer.Manager.setWindowTitle', 'dip-viewer-Manager.html#dip-viewer-Manager-setWindowTitle-Window-P-char-CP'),
     ('dip.viewer.Manager.refreshWindow', 'dip-viewer-Manager.html#dip-viewer-Manager-refreshWindow-Window-P'),
@@ -2012,7 +2013,7 @@
     ('dip.FillDelta', 'generation_test.html#dip-FillDelta-Image-L-String-CL'),
     ('dip.CreateDelta', 'generation_test.html#dip-CreateDelta-Image-L-UnsignedArray-CL-String-CL'),
     ('dip.CreateDelta', 'generation_test.html#dip-CreateDelta-UnsignedArray-CL-String-CL'),
-    ('dip.CreateGauss', 'generation_test.html#dip-CreateGauss-Image-L-FloatArray-CL-UnsignedArray--dfloat--UnsignedArray-'),
+    ('dip.CreateGauss', 'generation_test.html#dip-CreateGauss-Image-L-FloatArray-CL-UnsignedArray--dfloat--UnsignedArray--String-CL'),
     ('dip.CreateGabor', 'generation_test.html#dip-CreateGabor-Image-L-FloatArray-CL-FloatArray-CL-dfloat-'),
     ('dip.FTEllipsoid', 'generation_test.html#dip-FTEllipsoid-Image-L-FloatArray--dfloat-'),
     ('dip.FTEllipsoid', 'generation_test.html#dip-FTEllipsoid-UnsignedArray-CL-FloatArray--dfloat-'),

pydip/src/generation.cpp

@@ -105,10 +105,10 @@ void init_generation( py::module& m ) {
           "sizes"_a, "origin"_a = "", doc_strings::dip·CreateDelta·UnsignedArray·CL·String·CL );
    m.def( "CreateDelta", py::overload_cast< dip::Image&, dip::UnsignedArray const&, dip::String const& >( &dip::CreateDelta ),
           py::kw_only(), "out"_a, "sizes"_a, "origin"_a = "", doc_strings::dip·CreateDelta·Image·L·UnsignedArray·CL·String·CL );
-   m.def( "CreateGauss", py::overload_cast< dip::FloatArray const&, dip::UnsignedArray, dip::dfloat, dip::UnsignedArray >( &dip::CreateGauss ),
-          "sigmas"_a, "order"_a = dip::UnsignedArray{ 0 }, "truncation"_a = 3.0, "exponents"_a = dip::UnsignedArray{ 0 }, doc_strings::dip·CreateGauss·Image·L·FloatArray·CL·UnsignedArray··dfloat··UnsignedArray· );
-   m.def( "CreateGauss", py::overload_cast< dip::Image&, dip::FloatArray const&, dip::UnsignedArray, dip::dfloat, dip::UnsignedArray >( &dip::CreateGauss ),
-          py::kw_only(), "out"_a, "sigmas"_a, "order"_a = dip::UnsignedArray{ 0 }, "truncation"_a = 3.0, "exponents"_a = dip::UnsignedArray{ 0 }, doc_strings::dip·CreateGauss·Image·L·FloatArray·CL·UnsignedArray··dfloat··UnsignedArray· );
+   m.def( "CreateGauss", py::overload_cast< dip::FloatArray const&, dip::UnsignedArray, dip::dfloat, dip::UnsignedArray, dip::String const& >( &dip::CreateGauss ),
+          "sigmas"_a, "order"_a = dip::UnsignedArray{ 0 }, "truncation"_a = 3.0, "exponents"_a = dip::UnsignedArray{ 0 }, "extent"_a = "full", doc_strings::dip·CreateGauss·Image·L·FloatArray·CL·UnsignedArray··dfloat··UnsignedArray··String·CL );
+   m.def( "CreateGauss", py::overload_cast< dip::Image&, dip::FloatArray const&, dip::UnsignedArray, dip::dfloat, dip::UnsignedArray, dip::String const& >( &dip::CreateGauss ),
+          py::kw_only(), "out"_a, "sigmas"_a, "order"_a = dip::UnsignedArray{ 0 }, "truncation"_a = 3.0, "exponents"_a = dip::UnsignedArray{ 0 }, "extent"_a = "full", doc_strings::dip·CreateGauss·Image·L·FloatArray·CL·UnsignedArray··dfloat··UnsignedArray··String·CL );
    m.def( "CreateGabor", py::overload_cast< dip::FloatArray const&, dip::FloatArray const&, dip::dfloat >( &dip::CreateGabor ),
           "sigmas"_a, "frequencies"_a, "truncation"_a = 3.0, doc_strings::dip·CreateGabor·Image·L·FloatArray·CL·FloatArray·CL·dfloat· );
    m.def( "CreateGabor", py::overload_cast< dip::Image&, dip::FloatArray const&, dip::FloatArray const&, dip::dfloat >( &dip::CreateGabor ),

src/linear/gauss.cpp

@@ -50,8 +50,7 @@ inline dip::uint HalfGaussianSize(
    return clamp_cast< dip::uint >( std::ceil( truncation * sigma ));
 }
 
-// Creates a half Gaussian kernel, with the x=0 at the right end (last element) of the output array.
-std::vector< dfloat > MakeHalfGaussian(
+std::vector< dfloat > MakeHalfGaussianInternal(
       dfloat sigma,
       dip::uint derivativeOrder,
       dfloat truncation,
@@ -148,6 +147,23 @@ std::vector< dfloat > MakeHalfGaussian(
 
 } // namespace
 
+// Creates a half Gaussian kernel, with the x=0 at the right end (last element) of the output array.
+std::vector< dfloat > MakeHalfGaussian(
+      dfloat sigma,
+      dip::uint derivativeOrder,
+      dfloat truncation,
+      DataType dt
+) {
+   // Handle sigma == 0.0
+   if( sigma == 0.0 ) {
+      return { 1.0 };
+   }
+   // Create half Gaussian
+   std::vector< dfloat > gaussian;
+   DIP_STACK_TRACE_THIS( gaussian = MakeHalfGaussianInternal( sigma, derivativeOrder, truncation, dt ));
+   return gaussian;
+}
+
 // Create 1D full Gaussian
 std::vector< dfloat > MakeGaussian(
       dfloat sigma,
@@ -161,7 +177,7 @@ std::vector< dfloat > MakeGaussian(
    }
    // Create half Gaussian
    std::vector< dfloat > gaussian;
-   DIP_STACK_TRACE_THIS( gaussian = MakeHalfGaussian( sigma, derivativeOrder, truncation, dt ));
+   DIP_STACK_TRACE_THIS( gaussian = MakeHalfGaussianInternal( sigma, derivativeOrder, truncation, dt ));
    dip::uint halfFilterSize = gaussian.size() - 1;
    // Complete the Gaussian
    gaussian.resize( halfFilterSize * 2 + 1 );
@@ -177,22 +193,27 @@ void CreateGauss(
       FloatArray const& sigmas,
       UnsignedArray orders,
       dfloat truncation,
-      UnsignedArray exponents
+      UnsignedArray exponents,
+      String const& extent
 ) {
    // Verify dimensionality
    dip::uint nDims = sigmas.size();
    DIP_STACK_TRACE_THIS( ArrayUseParameter( orders, nDims, dip::uint( 0 )));
    DIP_STACK_TRACE_THIS( ArrayUseParameter( exponents, nDims, dip::uint( 0 )));
 
+   bool full{};
+   DIP_STACK_TRACE_THIS( full = BooleanFromString( extent, "full", "half" ));
+
    // Create 1D gaussian for each dimension
    std::vector< std::vector< dfloat >> gaussians( nDims );
    UnsignedArray outSizes( nDims );
    UnsignedArray centers( nDims );
+
    for( dip::uint ii = 0; ii < nDims; ++ii ) {
-      DIP_STACK_TRACE_THIS( gaussians[ ii ] = MakeGaussian( sigmas[ ii ], orders[ ii ], truncation, DT_DFLOAT ));
+      DIP_STACK_TRACE_THIS( gaussians[ ii ] = full ? MakeGaussian( sigmas[ ii ], orders[ ii ], truncation, DT_DFLOAT ) : MakeHalfGaussian( sigmas[ ii ], orders[ ii ], truncation, DT_DFLOAT ));
       dip::uint gaussianLength = gaussians[ ii ].size();
       outSizes[ ii ] = gaussianLength;
-      centers[ ii ] = ( gaussianLength - 1 ) / 2;
+      centers[ ii ] = full ? ( gaussianLength - 1 ) / 2 : gaussianLength - 1;
    }
 
    // Create output image