Separate morphology filter into X and Y processing phases.

Gives a massive performance boost for large radii.

(bzr r9508.1.36)

2 files changed, 70 insertions, 171 deletions

diff --git a/src/display/cairo-templates.h b/src/display/cairo-templates.h
index 88f3bf41a..900ca2d54 100644
--- a/src/display/cairo-templates.h
+++ b/src/display/cairo-templates.h
@@ -358,7 +358,7 @@ struct SurfaceSynth {
         cairo_surface_flush(surface);
     }
 
-    guint32 pixelAt(int x, int y) {
+    guint32 pixelAt(int x, int y) const {
         if (_alpha) {
             unsigned char *px = _px + y*_stride + x;
             return *px << 24;
@@ -367,7 +367,7 @@ struct SurfaceSynth {
             return *reinterpret_cast<guint32*>(px);
         }
     }
-    guint32 alphaAt(int x, int y) {
+    guint32 alphaAt(int x, int y) const {
         if (_alpha) {
             unsigned char *px = _px + y*_stride + x;
             return *px;
@@ -379,7 +379,7 @@ struct SurfaceSynth {
     }
 
     // retrieve a pixel value with bilinear interpolation
-    guint32 pixelAt(double x, double y) {
+    guint32 pixelAt(double x, double y) const {
         if (_alpha) {
             return alphaAt(x, y) << 24;
         }
@@ -416,7 +416,7 @@ struct SurfaceSynth {
     }
 
     // retrieve an alpha value with bilinear interpolation
-    guint32 alphaAt(double x, double y) {
+    guint32 alphaAt(double x, double y) const {
         double xf = floor(x), yf = floor(y);
         int xi = xf, yi = yf;
         guint32 xif = round((x - xf) * 255), yif = round((y - yf) * 255);
@@ -441,7 +441,7 @@ struct SurfaceSynth {
     }
 
     // compute surface normal at given coordinates using 3x3 Sobel gradient filter
-    NR::Fvector surfaceNormalAt(int x, int y, double scale) {
+    NR::Fvector surfaceNormalAt(int x, int y, double scale) const {
         // Below there are some multiplies by zero. They will be optimized out.
         // Do not remove them, because they improve readability.
         // NOTE: fetching using alphaAt is slightly lazy.
diff --git a/src/display/nr-filter-morphology.cpp b/src/display/nr-filter-morphology.cpp
index 51433da14..5be2c8b29 100644
--- a/src/display/nr-filter-morphology.cpp
+++ b/src/display/nr-filter-morphology.cpp
@@ -31,100 +31,71 @@ FilterPrimitive * FilterMorphology::create() {
 FilterMorphology::~FilterMorphology()
 {}
 
-struct MorphologyErode : public SurfaceSynth {
-    MorphologyErode(cairo_surface_t *in, double xradius, double yradius)
-        : SurfaceSynth(in)
-        , _xr(round(xradius))
-        , _yr(round(yradius))
-    {}
-    guint32 operator()(int x, int y) {
-        int startx = std::max(x - _xr, 0), endx = std::min(x + _xr + 1, _w);
-        int starty = std::max(y - _yr, 0), endy = std::min(y + _yr + 1, _h);
+enum MorphologyOp {
+    ERODE,
+    DILATE
+};
 
-        guint32 ao = 255;
-        guint32 ro = 255;
-        guint32 go = 255;
-        guint32 bo = 255;
+namespace {
 
-        if (_alpha) {
-            ao = 0xff000000;
-            for (int i = starty; i < endy; ++i) {
-                for (int j = startx; j < endx; ++j) {
-                    guint32 px = pixelAt(j, i);
-                    ao = std::min(ao, px & 0xff000000);
-                }
-            }
-            return ao;
-        } else {
-            for (int i = starty; i < endy; ++i) {
-                for (int j = startx; j < endx; ++j) {
-                    guint32 px = pixelAt(j, i);
-                    EXTRACT_ARGB32(px, a,r,g,b);
-                    if (a) {
-                        r = unpremul_alpha(r, a);
-                        g = unpremul_alpha(g, a);
-                        b = unpremul_alpha(b, a);
-                        ao = std::min(ao, a);
-                        ro = std::min(ro, r);
-                        go = std::min(go, g);
-                        bo = std::min(bo, b);
-                    } else {
-                        // zero pixel is guaranteed to be the minimum
-                        ao = 0; ro = 0; go = 0; bo = 0;
-                        goto end_loop;
-                    }
-                }
-            }
-            end_loop:
+template <MorphologyOp oper> guint32 extreme(guint32 a, guint32 b);
+template <> guint32 extreme<ERODE>(guint32 a, guint32 b) { return std::min(a, b); }
+template <> guint32 extreme<DILATE>(guint32 a, guint32 b) { return std::max(a, b); }
 
-            ro = premul_alpha(ro, ao);
-            go = premul_alpha(go, ao);
-            bo = premul_alpha(bo, ao);
-            ASSEMBLE_ARGB32(pxout, ao,ro,go,bo)
-            return pxout;
-        }
-    }
-private:
-    int _xr, _yr;
-};
-
-struct MorphologyDilate : public SurfaceSynth {
-    MorphologyDilate(cairo_surface_t *in, double xradius, double yradius)
+/* 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. We can do this, because the morphology operation
+ * is separable just like Gaussian blur. */
+template <MorphologyOp OP, Geom::Dim2 axis>
+struct Morphology : public SurfaceSynth {
+    Morphology(cairo_surface_t *in, double xradius)
         : SurfaceSynth(in)
-        , _xr(round(xradius))
-        , _yr(round(yradius))
+        , _radius(round(xradius))
     {}
     guint32 operator()(int x, int y) {
-        int startx = std::max(x - _xr, 0), endx = std::min(x + _xr + 1, _w);
-        int starty = std::max(y - _yr, 0), endy = std::min(y + _yr + 1, _h);
+        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);
+        }
 
-        guint32 ao = 0;
-        guint32 ro = 0;
-        guint32 go = 0;
-        guint32 bo = 0;
+        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) {
-            for (int i = starty; i < endy; ++i) {
-                for (int j = startx; j < endx; ++j) {
-                    guint32 px = pixelAt(j, i);
-                    ao = std::max(ao, px & 0xff000000);
-                }
+            ao = (OP == DILATE ? 0 : 0xff000000);
+            for (int i = start; i < end; ++i) {
+                guint32 px = (axis == Geom::X ? pixelAt(i, y) : pixelAt(x, i));
+                ao = extreme<OP>(ao, px & 0xff000000);
             }
             return ao;
         } else {
-            for (int i = starty; i < endy; ++i) {
-                for (int j = startx; j < endx; ++j) {
-                    guint32 px = pixelAt(j, i);
-                    EXTRACT_ARGB32(px, a,r,g,b)
-                    if (a == 0) continue; // this cannot affect the maximum
-
+            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);
+                if (a) {
                     r = unpremul_alpha(r, a);
                     g = unpremul_alpha(g, a);
                     b = unpremul_alpha(b, a);
-                    ao = std::max(ao, a);
-                    ro = std::max(ro, r);
-                    go = std::max(go, g);
-                    bo = std::max(bo, b);
+
+                    ao = extreme<OP>(ao, a);
+                    ro = extreme<OP>(ro, r);
+                    go = extreme<OP>(go, g);
+                    bo = extreme<OP>(bo, b);
+                } else {
+                    if (OP == DILATE) {
+                        continue; // zero pixel will not affect the maximum
+                    } else {
+                        // zero pixel is guaranteed to be the minimum
+                        ao = 0; ro = 0; go = 0; bo = 0;
+                        break;
+                    }
                 }
             }
 
@@ -136,112 +107,40 @@ struct MorphologyDilate : public SurfaceSynth {
         }
     }
 private:
-    int _xr, _yr;
+    int _radius;
 };
 
+} // end anonymous namespace
+
 void FilterMorphology::render_cairo(FilterSlot &slot)
 {
     cairo_surface_t *input = slot.getcairo(_input);
-    cairo_surface_t *out = ink_cairo_surface_create_identical(input);
 
     Geom::Matrix p2pb = slot.get_units().get_matrix_primitiveunits2pb();
     double xr = xradius * p2pb.expansionX();
     double yr = yradius * p2pb.expansionY();
 
+    cairo_surface_t *interm = ink_cairo_surface_create_identical(input);
+
     if (Operator == MORPHOLOGY_OPERATOR_DILATE) {
-        ink_cairo_surface_synthesize(out, MorphologyDilate(input, xr, yr));
+        ink_cairo_surface_synthesize(interm, Morphology<DILATE, Geom::X>(input, xr));
     } else {
-        ink_cairo_surface_synthesize(out, MorphologyErode(input, xr, yr));
-    }
-
-    slot.set(_output, out);
-    cairo_surface_destroy(out);
-}
-
-/*
-int FilterMorphology::render(FilterSlot &slot, FilterUnits const &units) {
-    NRPixBlock *in = slot.get(_input);
-    if (!in) {
-        g_warning("Missing source image for feMorphology (in=%d)", _input);
-        return 1;
-    }
-
-    NRPixBlock *out = new NRPixBlock;
-
-    // this primitive is defined for premultiplied RGBA values,
-    // thus convert them to that format
-    bool free_in_on_exit = false;
-    if (in->mode != NR_PIXBLOCK_MODE_R8G8B8A8P) {
-        NRPixBlock *original_in = in;
-        in = new NRPixBlock;
-        nr_pixblock_setup_fast(in, NR_PIXBLOCK_MODE_R8G8B8A8P,
-                               original_in->area.x0, original_in->area.y0,
-                               original_in->area.x1, original_in->area.y1,
-                               true);
-        nr_blit_pixblock_pixblock(in, original_in);
-        free_in_on_exit = true;
+        ink_cairo_surface_synthesize(interm, Morphology<ERODE, Geom::X>(input, xr));
     }
 
-    Geom::Matrix p2pb = units.get_matrix_primitiveunits2pb();
-    int const xradius = (int)round(this->xradius * p2pb.expansionX());
-    int const yradius = (int)round(this->yradius * p2pb.expansionY());
-
-    int x0=in->area.x0;
-    int y0=in->area.y0;
-    int x1=in->area.x1;
-    int y1=in->area.y1;
-    int w=x1-x0, h=y1-y0;
-    int x, y, i, j;
-    int rmax,gmax,bmax,amax;
-    int rmin,gmin,bmin,amin;
-
-    nr_pixblock_setup_fast(out, in->mode, x0, y0, x1, y1, true);
-
-    unsigned char *in_data = NR_PIXBLOCK_PX(in);
-    unsigned char *out_data = NR_PIXBLOCK_PX(out);
-
-    for(x = 0 ; x < w ; x++){
-        for(y = 0 ; y < h ; y++){
-            rmin = gmin = bmin = amin = 255;
-            rmax = gmax = bmax = amax = 0;
-            for(i = x - xradius ; i < x + xradius ; i++){
-                if (i < 0 || i >= w) continue;
-                for(j = y - yradius ; j < y + yradius ; j++){
-                    if (j < 0 || j >= h) continue;
-                    if(in_data[4*(i + w*j)]>rmax) rmax = in_data[4*(i + w*j)];
-                    if(in_data[4*(i + w*j)+1]>gmax) gmax = in_data[4*(i + w*j)+1];
-                    if(in_data[4*(i + w*j)+2]>bmax) bmax = in_data[4*(i + w*j)+2];
-                    if(in_data[4*(i + w*j)+3]>amax) amax = in_data[4*(i + w*j)+3];
+    cairo_surface_t *out = ink_cairo_surface_create_identical(input);
 
-                    if(in_data[4*(i + w*j)]<rmin) rmin = in_data[4*(i + w*j)];
-                    if(in_data[4*(i + w*j)+1]<gmin) gmin = in_data[4*(i + w*j)+1];
-                    if(in_data[4*(i + w*j)+2]<bmin) bmin = in_data[4*(i + w*j)+2];
-                    if(in_data[4*(i + w*j)+3]<amin) amin = in_data[4*(i + w*j)+3];
-                }
-            }
-            if (Operator==MORPHOLOGY_OPERATOR_DILATE){
-                out_data[4*(x + w*y)]=rmax;
-                out_data[4*(x + w*y)+1]=gmax;
-                out_data[4*(x + w*y)+2]=bmax;
-                out_data[4*(x + w*y)+3]=amax;
-            } else {
-                out_data[4*(x + w*y)]=rmin;
-                out_data[4*(x + w*y)+1]=gmin;
-                out_data[4*(x + w*y)+2]=bmin;
-                out_data[4*(x + w*y)+3]=amin;
-            }
-        }
+    if (Operator == MORPHOLOGY_OPERATOR_DILATE) {
+        ink_cairo_surface_synthesize(out, Morphology<DILATE, Geom::Y>(interm, yr));
+    } else {
+        ink_cairo_surface_synthesize(out, Morphology<ERODE, Geom::Y>(interm, yr));
     }
 
-    if (free_in_on_exit) {
-        nr_pixblock_release(in);
-        delete in;
-    }
+    cairo_surface_destroy(interm);
 
-    out->empty = FALSE;
     slot.set(_output, out);
-    return 0;
-}*/
+    cairo_surface_destroy(out);
+}
 
 void FilterMorphology::area_enlarge(NRRectL &area, Geom::Matrix const &trans)
 {