1 files changed, 264 insertions, 0 deletions
diff --git a/src/display/nr-filter-morphology.cpp b/src/display/nr-filter-morphology.cpp
new file mode 100644
index 0000000..700d69f
--- /dev/null
+++ b/src/display/nr-filter-morphology.cpp
@@ -0,0 +1,264 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * feMorphology filter primitive renderer
+ *
+ * Authors:
+ *   Felipe Corrêa da Silva Sanches <juca@members.fsf.org>
+ *
+ * Copyright (C) 2007 authors
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"  // only include where actually required!
+#endif
+
+#include <cmath>
+#include <algorithm>
+#include <deque>
+#include <functional>
+#include "display/cairo-templates.h"
+#include "display/cairo-utils.h"
+#include "display/nr-filter-morphology.h"
+#include "display/nr-filter-slot.h"
+#include "display/nr-filter-units.h"
+
+namespace Inkscape {
+namespace Filters {
+
+FilterMorphology::FilterMorphology()
+= default;
+
+FilterPrimitive * FilterMorphology::create() {
+    return new FilterMorphology();
+}
+
+FilterMorphology::~FilterMorphology()
+= default;
+
+enum MorphologyOp {
+    ERODE,
+    DILATE
+};
+
+namespace {
+
+/* This performs one "half" of the morphology operation by calculating 
+ * the componentwise extreme in the specified axis with the given radius.
+ * 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);
+    (void) numOfThreads; // suppress unused variable warning
+    #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?
+        }
+
+        // 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));
+                }
+                // Set output
+                if (j >= ri) {
+                    *out_p = vals[p].front().second;
+                    ++out_p;
+                }
+            }
+            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;
+        }
+    }
+
+    cairo_surface_mark_dirty(out);
+}
+
+} // end anonymous namespace
+
+void FilterMorphology::render_cairo(FilterSlot &slot)
+{
+    cairo_surface_t *input = slot.getcairo(_input);
+
+    if (xradius == 0.0 || yradius == 0.0) {
+        // output is transparent black
+        cairo_surface_t *out = ink_cairo_surface_create_identical(input);
+        copy_cairo_surface_ci(input, out);
+        slot.set(_output, out);
+        cairo_surface_destroy(out);
+        return;
+    }
+
+    int device_scale = slot.get_device_scale();
+    Geom::Affine p2pb = slot.get_units().get_matrix_primitiveunits2pb();
+    double xr = fabs(xradius * p2pb.expansionX()) * device_scale;
+    double yr = fabs(yradius * p2pb.expansionY()) * device_scale;
+    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) {
+        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 {
+        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(interm);
+
+    // color_interpolation_filters for out same as input. See spec (DisplacementMap).
+    copy_cairo_surface_ci(input, out);
+
+    if (Operator == MORPHOLOGY_OPERATOR_DILATE) {
+        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 {
+        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);
+
+    slot.set(_output, out);
+    cairo_surface_destroy(out);
+}
+
+void FilterMorphology::area_enlarge(Geom::IntRect &area, Geom::Affine const &trans)
+{
+    int enlarge_x = ceil(xradius * trans.expansionX());
+    int enlarge_y = ceil(yradius * trans.expansionY());
+
+    area.expandBy(enlarge_x, enlarge_y);
+}
+
+double FilterMorphology::complexity(Geom::Affine const &trans)
+{
+    int enlarge_x = ceil(xradius * trans.expansionX());
+    int enlarge_y = ceil(yradius * trans.expansionY());
+    return enlarge_x * enlarge_y;
+}
+
+void FilterMorphology::set_operator(FilterMorphologyOperator &o){
+    Operator = o;
+}
+
+void FilterMorphology::set_xradius(double x){
+    xradius = x;
+}
+
+void FilterMorphology::set_yradius(double y){
+    yradius = y;
+}
+
+} /* namespace Filters */
+} /* namespace Inkscape */
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :