dip::Reconstruction() is now significantly faster.

changelogs/diplib_3.1.0.md

@@ -42,12 +42,11 @@ title: "Changes DIPlib 3.1.0"
   using it (including `dip::MeasurementTool::Measure()`), `dip::GetImageChainCodes()`, `dip::GetObjectLabels()`,
   `dip::Relabel()`, `dip::ChordLength()` and `dip::PairCorrelation()`.
 
-- Improved speed of `dip::Measurement::AddObjectIDs()` and
-  `dip::Measurement::operator+()`.
-
-- `dip::Label` is slightly more efficient for 3D and higher-dimensional images; added tests.
-
-- Sped up `dip::MomentAccumulator`.
+- The speed of the following functions has been improved:
+    - `dip::Measurement::AddObjectIDs()` and `dip::Measurement::operator+()`.
+    - `dip::Label` (slightly more efficient for 3D and higher-dimensional images).
+    - `dip::MomentAccumulator`.
+    - `dip::MorphologicalReconstruction` and functions that depend on it.
 
 - `dip::OptimalFourierTransformSize()` has a new option to return a smaller or equal size, rather than
   a larger or equal size, so we can crop an image for efficient FFT instead of padding.
@@ -99,7 +98,7 @@ title: "Changes DIPlib 3.1.0"
 - There was a strange rounding error when creating disk-shaped filter kernels and structuring elements,
   for some even integer sizes, which caused these kernels to be not symmetric over 90 degree rotations.
 
-- Fixed bug in `dip::MeanAbs()` and `dip::SumAbs()` for complex inputs.
+- `dip::MeanAbs()` and `dip::SumAbs()` could produce wrong results for complex inputs.
 
 
 

src/morphology/areaopening.cpp

@@ -2,7 +2,7 @@
  * DIPlib 3.0
  * This file contains the area opening and related functions.
  *
- * (c)2008, 2017, Cris Luengo.
+ * (c)2008, 2017-2021 Cris Luengo.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.

src/morphology/reconstruction.cpp

@@ -2,7 +2,7 @@
  * DIPlib 3.0
  * This file contains the morphological reconstruction and related functions.
  *
- * (c)2017-2018, Cris Luengo.
+ * (c)2017-2021, Cris Luengo.
  * Based on original DIPlib code: (c)1995-2014, Delft University of Technology.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
@@ -27,7 +27,7 @@
 #include "diplib/neighborlist.h"
 #include "diplib/iterators.h"
 #include "diplib/overload.h"
-#include "watershed_support.h"
+#include "diplib/border.h"
 
 namespace dip {
 
@@ -68,55 +68,49 @@ void MorphologicalReconstructionInternal(
    TPI minval = *static_cast< TPI const* >( c_minval.Origin() );
 
    // Walk over the entire image & put all the pixels larger than minval on the heap.
-   // Also set done=DIP_FALSE for all pixels and reduce the value of c_out where c_out > c_in.
-   JointImageIterator< TPI, TPI, bin > it( { c_in, c_out, c_done } );
-   if( dilation ) {
-      do {
-         if( it.Out() > it.In() ) {
-            it.Out() = it.In();
-            it.template Sample< 2 >() = true;
-         } else {
-            it.template Sample< 2 >() = false;
-         }
-         if( it.Out() > minval ) {
-            Q.push( Qitem< TPI >{ it.Out(), it.template Offset< 2 >() } ); // offset pushed is that of `done`, so we can test it fast.
-            //std::cout << " - Pushed " << it.Out();
-         }
-      } while( ++it );
-   } else {
-      do {
-         if( it.Out() < it.In() ) {
-            it.Out() = it.In();
-            it.template Sample< 2 >() = true;
-         } else {
-            it.template Sample< 2 >() = false;
-         }
-         if( it.Out() < minval ) {
-            Q.push( Qitem< TPI >{ it.Out(), it.template Offset< 2 >() } ); // offset pushed is that of `done`, so we can test it fast.
-            //std::cout << " - Pushed " << it.Out();
-         }
-      } while( ++it );
-   }
+   // TODO: Many of these heap items are not necessary, we will enqueue the same pixel multiple times until it
+   //       is actually processed. We cannot update items in the queue, we can only drop them when popping and
+   //       seeing that the pixel was already processed. If there was a cheap way to determine here if a pixel
+   //       needs to be enqueued at this point or not, it might help speed up things a bit. We'd need to enqueue
+   //       only pixels that:
+   //         - have a neighbor that has a lower value in `out`, and doesn't have the lower bit of `done` set, and
+   //         - don't have a neighbor with a higher value that will propagate into it.
+   //       Not an easy test!
+   JointImageIterator< bin, TPI  > it( { c_done, c_out } );
+   do {
+      if( it.Out() != minval ) {
+         Q.push( Qitem< TPI >{ it.Out(), it.template Offset< 0 >() } ); // offset pushed is that of `done`, so we can test it fast.
+         //std::cout << " - Pushed " << it.Out();
+      }
+   } while( ++it );
 
    // Start processing pixels
    TPI* in = static_cast< TPI* >( c_in.Origin() );
    TPI* out = static_cast< TPI* >( c_out.Origin() );
-   bin* done = static_cast< bin* >( c_done.Origin() );
+   uint8* done = static_cast< uint8* >( c_done.Origin() ); // It's binary, but we use other bit planes too.
+   // To wit: the bottom bit is set if he pixel value will no longer change (set in caller, updated here)
+   //         the 2nd bit is set if we already propagated out from this pixel
+   //         the 3rd bit is set if this is a border pixel (set in caller)
    auto coordinatesComputer = c_done.OffsetToCoordinatesComputer();
    BooleanArray skipar( nNeigh );
    while( !Q.empty() ) {
       dip::sint offsetDone = Q.top().offset;
       Q.pop();
+      //std::cout << " - Popped: " << done[ offsetDone ];
+      if( done[ offsetDone ] & 2u ) {
+         // 2nd bit is set if we already propagated out from this pixel
+         continue;
+      }
       UnsignedArray coords = coordinatesComputer( offsetDone );
       dip::sint offsetIn = c_in.Offset( coords );
       dip::sint offsetOut = c_out.Offset( coords );
       //std::cout << " - Popped " << offsetDone << " (" << out[ offsetOut ] << ")";
       // Iterate over all neighbors
       auto lit = neighborList.begin();
       for( dip::uint jj = 0; jj < nNeigh; ++jj, ++lit ) {
-         if( lit.IsInImage( coords, imsz )) {
+         if( !( done[ offsetDone ] & 4u ) || lit.IsInImage( coords, imsz )) { // test IsInImage only for border pixels
             // Propagate this pixel's value to its unfinished neighbours
-            if( !done[ offsetDone + neighborOffsetsDone[ jj ]] ) {
+            if( !( done[ offsetDone + neighborOffsetsDone[ jj ]] & 1u )) { // if the bottom bit is set, we don't need to propagate
                TPI newval = in[ offsetIn + neighborOffsetsIn[ jj ]];
                if( dilation ) {
                   newval = std::min( newval, out[ offsetOut ] );
@@ -139,7 +133,7 @@ void MorphologicalReconstructionInternal(
          }
       }
       // Mark this pixels done
-      done[ offsetDone ] = true;
+      done[ offsetDone ] = 3; // bottom two bits both set
    }
 }
 
@@ -176,10 +170,12 @@ void MorphologicalReconstruction (
       out.Strip(); // We can work in-place if c_marker and out are the same image, but c_in must be separate from out.
    }
    DIP_STACK_TRACE_THIS( Convert( marker, out, in.DataType() ));
+   DIP_STACK_TRACE_THIS( dilation ? Infimum( in, out, out ) : Supremum( in, out, out ));
    Image minval = dilation ? Minimum( out ) : Maximum( out ); // same data type as `out`
 
    // Intermediate image
-   Image done( out.Sizes(), 1, DT_BIN );
+   Image done = in == out; // Sets the bottom bit of the dip::bin byte
+   detail::ProcessBorders< bin >( done, []( bin* ptr, dip::sint ){ *reinterpret_cast< uint8* >( ptr ) += 4; } ); // Set the 3rd bit to indicate a border pixel
 
    // Create array with offsets to neighbours
    NeighborList neighborList( { Metric::TypeCode::CONNECTED, connectivity }, nDims );