Morphological filters linear in the number of output pixels.

(bzr r11079)

1 files changed, 131 insertions, 66 deletions

diff --git a/src/display/nr-filter-morphology.cpp b/src/display/nr-filter-morphology.cpp
index b6aea1b06..18f99cdcd 100644
--- a/src/display/nr-filter-morphology.cpp
+++ b/src/display/nr-filter-morphology.cpp
@@ -11,6 +11,8 @@
 
 #include <cmath>
 #include <algorithm>
+#include <deque>
+#include <functional>
 #include "display/cairo-templates.h"
 #include "display/cairo-utils.h"
 #include "display/nr-filter-morphology.h"
@@ -40,72 +42,118 @@ namespace {
 
 /* This performs one "half" of the morphology operation by calculating 
  * the componentwise extreme in the specified axis with the given radius.
- * Performing the operation one axis at a time gives us a MASSIVE performance boost
- * at large morphology radii. Extreme of row extremes is equal to the extreme
- * of components, so this doesn't change the result. */
-template <MorphologyOp OP, Geom::Dim2 axis>
-struct Morphology : public SurfaceSynth {
-    Morphology(cairo_surface_t *in, double xradius)
-        : SurfaceSynth(in)
-        , _radius(round(xradius))
-    {}
-    guint32 operator()(int x, int y) {
-        int start, end;
-        if (axis == Geom::X) {
-            start = std::max(0, x - _radius);
-            end = std::min(x + _radius + 1, _w);
-        } else {
-            start = std::max(0, y - _radius);
-            end = std::min(y + _radius + 1, _h);
+ * Extreme of row extremes is equal to the extreme of components, so this
+ * doesn't change the result.
+ * The algorithm is due to: Petr Dokládal, Eva Dokládalová (2011), "Computationally efficient, one-pass algorithm for morphological filters"
+ * TODO: Currently only the 1D algorithm is implemented, but it should not be too difficult (and at the very least more memory efficient) to implement the full 2D algorithm.
+ *       One problem with the 2D algorithm is that it is harder to parallelize.
+ */
+template <typename Comparison, Geom::Dim2 axis, int BPP>
+void morphologicalFilter1D(cairo_surface_t * const input, cairo_surface_t * const out, double radius) {
+    Comparison comp;
+
+    int w = cairo_image_surface_get_width(out);
+    int h = cairo_image_surface_get_height(out);
+    if (axis == Geom::Y) std::swap(w,h);
+
+    int stridein = cairo_image_surface_get_stride(input);
+    int strideout = cairo_image_surface_get_stride(out);
+
+    unsigned char *in_data = cairo_image_surface_get_data(input);
+    unsigned char *out_data = cairo_image_surface_get_data(out);
+
+    int ri = round(radius); // TODO: Support fractional radii?
+    int wi = 2*ri+1;
+
+    #if HAVE_OPENMP
+    int limit = w * h;
+    Inkscape::Preferences *prefs = Inkscape::Preferences::get();
+    int numOfThreads = prefs->getIntLimited("/options/threading/numthreads", omp_get_num_procs(), 1, 256);
+    #pragma omp parallel for if(limit > OPENMP_THRESHOLD) num_threads(numOfThreads)
+    #endif // HAVE_OPENMP
+    for (int i = 0; i < h; ++i) {
+        // TODO: Store position and value in one 32 bit integer? 24 bits should be enough for a position, it would be quite strange to have an image with a width/height of more than 16 million(!).
+        std::deque< std::pair<int,unsigned char> > vals[BPP]; // In my tests it was actually slightly faster to allocate it here than allocate it once for all threads and retrieving the correct set based on the thread id.
+
+        // Initialize with transparent black
+        for(int p = 0; p < BPP; ++p) {
+            vals[p].push_back(std::pair<int,unsigned char>(-1,0)); // TODO: Only do this when performing an erosion?
         }
 
-        guint32 ao = (OP == DILATE ? 0 : 255);
-        guint32 ro = (OP == DILATE ? 0 : 255);
-        guint32 go = (OP == DILATE ? 0 : 255);
-        guint32 bo = (OP == DILATE ? 0 : 255);
-
-        if (_alpha) {
-            ao = (OP == DILATE ? 0 : 0xff000000);
-            for (int i = start; i < end; ++i) {
-                guint32 px = (axis == Geom::X ? pixelAt(i, y) : pixelAt(x, i));
-                if (OP == DILATE) {
-                    if (px > ao) ao = px;
-                } else {
-                    if (px < ao) ao = px;
+        // Process "row"
+        unsigned char *in_p = in_data + i * (axis == Geom::X ? stridein : BPP);
+        unsigned char *out_p = out_data + i * (axis == Geom::X ? strideout : BPP);
+        /* This is the "short but slow" version, which might be easier to follow, the longer but faster version follows.
+        for (int j = 0; j < w+ri; ++j) {
+            for(int p = 0; p < BPP; ++p) { // Iterate over channels
+                // Push new value onto FIFO, erasing any previous values that are "useless" (see paper) or out-of-range
+                if (!vals[p].empty() && vals[p].front().first+wi <= j) vals[p].pop_front(); // out-of-range
+                if (j < w) {
+                    while(!vals[p].empty() && !comp(vals[p].back().second, *in_p)) vals[p].pop_back(); // useless
+                    vals[p].push_back(std::make_pair(j, *in_p));
+                    ++in_p;
+                } else if (j == w) { // Transparent black beyond the image. TODO: Only do this when performing an erosion?
+                    while(!vals[p].empty() && !comp(vals[p].back().second, 0)) vals[p].pop_back();
+                    vals[p].push_back(std::make_pair(j, 0));
                 }
-            }
-            return ao;
-        } else {
-            for (int i = start; i < end; ++i) {
-                guint32 px = (axis == Geom::X ? pixelAt(i, y) : pixelAt(x, i));
-                EXTRACT_ARGB32(px, a,r,g,b);
-
-                // this will be compiled to conditional moves;
-                // the operator comparison will be evaluated at compile time.
-                // therefore there will be no branching in this loop
-                if (OP == DILATE) {
-                    if (a > ao) ao = a;
-                    if (r > ro) ro = r;
-                    if (g > go) go = g;
-                    if (b > bo) bo = b;
-                } else {
-                    if (a < ao) ao = a;
-                    if (r < ro) ro = r;
-                    if (g < go) go = g;
-                    if (b < bo) bo = b;
+                // Set output
+                if (j >= ri) {
+                    *out_p = vals[p].front().second;
+                    ++out_p;
                 }
-
-                // TODO: verify whether this check gives any speedup.
-                if (OP == ERODE && a == 0) break;
             }
-
-            ASSEMBLE_ARGB32(pxout, ao,ro,go,bo)
-            return pxout;
+            if (axis == Geom::Y && j < w  ) in_p += stridein - BPP;
+            if (axis == Geom::Y && j >= ri) out_p += strideout - BPP;
+        }*/
+        for (int j = 0; j < std::min(ri,w); ++j) {
+            for(int p = 0; p < BPP; ++p) { // Iterate over channels
+                // Push new value onto FIFO, erasing any previous values that are "useless" (see paper) or out-of-range
+                if (!vals[p].empty() && vals[p].front().first <= j) vals[p].pop_front(); // out-of-range
+                while(!vals[p].empty() && !comp(vals[p].back().second, *in_p)) vals[p].pop_back(); // useless
+                vals[p].push_back(std::make_pair(j+wi, *in_p));
+                ++in_p;
+            }
+            if (axis == Geom::Y) in_p += stridein - BPP;
+        }
+        // We have now done all preparatory work.
+        // If w<=ri, then the following loop does nothing (which is as it should).
+        for (int j = ri; j < w; ++j) {
+            for(int p = 0; p < BPP; ++p) { // Iterate over channels
+                // Push new value onto FIFO, erasing any previous values that are "useless" (see paper) or out-of-range
+                if (!vals[p].empty() && vals[p].front().first <= j) vals[p].pop_front(); // out-of-range
+                while(!vals[p].empty() && !comp(vals[p].back().second, *in_p)) vals[p].pop_back(); // useless
+                vals[p].push_back(std::make_pair(j+wi, *in_p));
+                ++in_p;
+                // Set output
+                *out_p = vals[p].front().second;
+                ++out_p;
+            }
+            if (axis == Geom::Y) {
+                in_p += stridein - BPP;
+                out_p += strideout - BPP;
+            }
+        }
+        // We have now done all work which involves both input and output.
+        // The following loop makes sure that the border is handled correctly.
+        for(int p = 0; p < BPP; ++p) { // Iterate over channels
+            while(!vals[p].empty() && !comp(vals[p].back().second, 0)) vals[p].pop_back();
+            vals[p].push_back(std::make_pair(w+wi, 0));
+        }
+        // Now we just have to finish the output.
+        for (int j = std::max(w,ri); j < w+ri; ++j) {
+            for(int p = 0; p < BPP; ++p) { // Iterate over channels
+                // Remove out-of-range values
+                if (!vals[p].empty() && vals[p].front().first <= j) vals[p].pop_front(); // out-of-range
+                // Set output
+                *out_p = vals[p].front().second;
+                ++out_p;
+            }
+            if (axis == Geom::Y) out_p += strideout - BPP;
         }
     }
-private:
-    int _radius;
-};
+
+    cairo_surface_mark_dirty(out);
+}
 
 } // end anonymous namespace
 
@@ -122,23 +170,40 @@ void FilterMorphology::render_cairo(FilterSlot &slot)
     }
 
     Geom::Affine p2pb = slot.get_units().get_matrix_primitiveunits2pb();
-    double xr = xradius * p2pb.expansionX();
-    double yr = yradius * p2pb.expansionY();
+    double xr = fabs(xradius * p2pb.expansionX());
+    double yr = fabs(yradius * p2pb.expansionY());
+    int bpp = cairo_image_surface_get_format(input) == CAIRO_FORMAT_A8 ? 1 : 4;
 
     cairo_surface_t *interm = ink_cairo_surface_create_identical(input);
 
     if (Operator == MORPHOLOGY_OPERATOR_DILATE) {
-        ink_cairo_surface_synthesize(interm, Morphology<DILATE, Geom::X>(input, xr));
+        if (bpp == 1) {
+            morphologicalFilter1D< std::greater<unsigned char>, Geom::X, 1 >(input, interm, xr);
+        } else {
+            morphologicalFilter1D< std::greater<unsigned char>, Geom::X, 4 >(input, interm, xr);
+        }
     } else {
-        ink_cairo_surface_synthesize(interm, Morphology<ERODE, Geom::X>(input, xr));
+        if (bpp == 1) {
+            morphologicalFilter1D< std::less<unsigned char>, Geom::X, 1 >(input, interm, xr);
+        } else {
+            morphologicalFilter1D< std::less<unsigned char>, Geom::X, 4 >(input, interm, xr);
+        }
     }
 
-    cairo_surface_t *out = ink_cairo_surface_create_identical(input);
+    cairo_surface_t *out = ink_cairo_surface_create_identical(interm);
 
     if (Operator == MORPHOLOGY_OPERATOR_DILATE) {
-        ink_cairo_surface_synthesize(out, Morphology<DILATE, Geom::Y>(interm, yr));
+        if (bpp == 1) {
+            morphologicalFilter1D< std::greater<unsigned char>, Geom::Y, 1 >(interm, out, yr);
+        } else {
+            morphologicalFilter1D< std::greater<unsigned char>, Geom::Y, 4 >(interm, out, yr);
+        }
     } else {
-        ink_cairo_surface_synthesize(out, Morphology<ERODE, Geom::Y>(interm, yr));
+        if (bpp == 1) {
+            morphologicalFilter1D< std::less<unsigned char>, Geom::Y, 1 >(interm, out, yr);
+        } else {
+            morphologicalFilter1D< std::less<unsigned char>, Geom::Y, 4 >(interm, out, yr);
+        }
     }
 
     cairo_surface_destroy(interm);