Color matrix filter. Fix arithmetic operator in feComposite

(bzr r9508.1.23)

5 files changed, 278 insertions, 28 deletions

diff --git a/src/display/cairo-templates.h b/src/display/cairo-templates.h
index 79e461cf8..c64ad78c1 100644
--- a/src/display/cairo-templates.h
+++ b/src/display/cairo-templates.h
@@ -146,7 +146,7 @@ void ink_cairo_surface_blend(cairo_surface_t *in1, cairo_surface_t *in2, cairo_s
                 #pragma omp parallel for num_threads(num_threads)
                 #endif
                 for (int i = 0; i < h; ++i) {
-                    guint8  *in1_p = reinterpret_cast<guint8*>(in1_data) + i * stride1;
+                    guint8 *in1_p = reinterpret_cast<guint8*>(in1_data) + i * stride1;
                     guint8 *in2_p = reinterpret_cast<guint8*>(in2_data) + i * stride2;
                     guint8 *out_p = reinterpret_cast<guint8*>(out_data) + i * strideout;
                     for (int j = 0; j < w; ++j) {
@@ -164,22 +164,29 @@ void ink_cairo_surface_blend(cairo_surface_t *in1, cairo_surface_t *in2, cairo_s
     cairo_surface_mark_dirty(out);
 }
 
-#if 0
 template <typename Filter>
-ink_cairo_surface_filter(cairo_surface_t *in, cairo_surface_t *out, Filter filter)
+void ink_cairo_surface_filter(cairo_surface_t *in, cairo_surface_t *out, Filter filter)
 {
     cairo_surface_flush(in);
 
     // ASSUMPTIONS
     // 1. Cairo ARGB32 surface strides are always divisible by 4
     // 2. We can only receive CAIRO_FORMAT_ARGB32 or CAIRO_FORMAT_A8 surfaces
-    // 3. Surfaces have the same dimensions and pixel formats
+    // 3. Surfaces have the same dimensions
+    // 4. Output surface is A8 if input is A8
 
     int w = cairo_image_surface_get_width(in);
     int h = cairo_image_surface_get_height(in);
     int stridein   = cairo_image_surface_get_stride(in);
     int strideout = cairo_image_surface_get_stride(out);
-    int bpp = cairo_image_surface_get_format(in) == CAIRO_FORMAT_A8 ? 1 : 4;
+    int bppin = cairo_image_surface_get_format(in) == CAIRO_FORMAT_A8 ? 1 : 4;
+    int bppout = cairo_image_surface_get_format(out) == CAIRO_FORMAT_A8 ? 1 : 4;
+    int limit = w * h;
+
+    // Check whether we can loop over pixels without taking stride into account.
+    bool fast_path = true;
+    fast_path &= (stridein == w * bppin);
+    fast_path &= (strideout == w * bppout);
 
     guint32 *const in_data = (guint32*) cairo_image_surface_get_data(in);
     guint32 *const out_data = (guint32*) cairo_image_surface_get_data(out);
@@ -189,17 +196,86 @@ ink_cairo_surface_filter(cairo_surface_t *in, cairo_surface_t *out, Filter filte
     int num_threads = prefs->getIntLimited("/options/threading/numthreads", omp_get_num_procs(), 1, 256);
     #endif
 
-    if (bpp == 4) {
-        
+    if (bppin == 4) {
+        if (bppout == 4) {
+            // bppin == 4, bppout == 4
+            if (fast_path) {
+                #if HAVE_OPENMP
+                #pragma omp parallel for num_threads(num_threads)
+                #endif
+                for (int i = 0; i < limit; ++i) {
+                    *(out_data + i) = filter(*(in_data + i));
+                }
+            } else {
+                #if HAVE_OPENMP
+                #pragma omp parallel for num_threads(num_threads)
+                #endif
+                for (int i = 0; i < h; ++i) {
+                    guint32 *in_p = in_data + i * stridein/4;
+                    guint32 *out_p = out_data + i * strideout/4;
+                    for (int j = 0; j < w; ++j) {
+                        *out_p = filter(*in_p);
+                        ++in_p; ++out_p;
+                    }
+                }
+            }
+        } else {
+            // bppin == 4, bppout == 1
+            // we use this path with COLORMATRIX_LUMINANCETOALPHA
+            #if HAVE_OPENMP
+            #pragma omp parallel for num_threads(num_threads)
+            #endif
+            for (int i = 0; i < h; ++i) {
+                guint32 *in_p = in_data + i * stridein/4;
+                guint8 *out_p = reinterpret_cast<guint8*>(out_data) + i * strideout;
+                for (int j = 0; j < w; ++j) {
+                    guint32 out_px = filter(*in_p);
+                    *out_p = out_px >> 24;
+                    ++in_p; ++out_p;
+                }
+            }
+        }
     } else {
-        #if HAVE_OPENMP
-        #pragma omp parallel for num_threads(num_threads)
-        #endif
+        // bppin == 1, bppout == 1
+        // Note: there is no path for bppin == 1, bppout == 4 because it is useless
+        if (fast_path) {
+            #if HAVE_OPENMP
+            #pragma omp parallel for num_threads(num_threads)
+            #endif
+            for (int i = 0; i < limit; ++i) {
+                guint8 *in_p = reinterpret_cast<guint8*>(in_data) + i;
+                guint8 *out_p = reinterpret_cast<guint8*>(out_data) + i;
+                guint32 in_px = *in_p; in_px <<= 24;
+                guint32 out_px = filter(in_px);
+                *out_p = out_px >> 24;
+            }
+        } else {
+            #if HAVE_OPENMP
+            #pragma omp parallel for num_threads(num_threads)
+            #endif
+            for (int i = 0; i < h; ++i) {
+                guint8 *in_p = reinterpret_cast<guint8*>(in_data) + i * stridein;
+                guint8 *out_p = reinterpret_cast<guint8*>(out_data) + i * strideout;
+                for (int j = 0; j < w; ++j) {
+                    guint32 in_px = *in_p; in_px <<= 24;
+                    guint32 out_px = filter(in_px);
+                    *out_p = out_px >> 24;
+                    ++in_p; ++out_p;
+                }
+            }
+        }
     }
 }
-#endif
 
-// helper macros for pixel extraction
+// Some helpers for pixel manipulation
+
+G_GNUC_CONST inline gint32
+pxclamp(gint32 v, gint32 low, gint32 high) {
+    if (v < low) return low;
+    if (v > high) return high;
+    return v;
+}
+
 #define EXTRACT_ARGB32(px,a,r,g,b) \
     guint32 a, r, g, b; \
     a = (px & 0xff000000) >> 24; \
diff --git a/src/display/nr-filter-colormatrix.cpp b/src/display/nr-filter-colormatrix.cpp
index 0b24649a9..d0926ebcc 100644
--- a/src/display/nr-filter-colormatrix.cpp
+++ b/src/display/nr-filter-colormatrix.cpp
@@ -10,11 +10,13 @@
  * Released under GNU GPL, read the file 'COPYING' for more information
  */
 
+#include <math.h>
+#include <algorithm>
+#include "display/cairo-templates.h"
+#include "display/cairo-utils.h"
 #include "display/nr-filter-colormatrix.h"
 #include "display/nr-filter-units.h"
 #include "display/nr-filter-utils.h"
-#include "libnr/nr-blit.h"
-#include <math.h>
 
 namespace Inkscape {
 namespace Filters {
@@ -30,6 +32,7 @@ FilterPrimitive * FilterColorMatrix::create() {
 FilterColorMatrix::~FilterColorMatrix()
 {}
 
+#if 0
 int FilterColorMatrix::render(FilterSlot &slot, FilterUnits const &/*units*/) {
     NRPixBlock *in = slot.get(_input);
     if (!in) {
@@ -194,6 +197,169 @@ int FilterColorMatrix::render(FilterSlot &slot, FilterUnits const &/*units*/) {
     slot.set(_output, out);
     return 0;
 }
+#endif
+
+struct ColorMatrixMatrix {
+    ColorMatrixMatrix(std::vector<double> const &values) {
+        unsigned limit = std::min(20ul, values.size());
+        for (unsigned i = 0; i < limit; ++i) {
+            if (i % 5 == 4) {
+                _v[i] = round(values[i]*255*255);
+            } else {
+                _v[i] = round(values[i]*255);
+            }
+        }
+        for (unsigned i = limit; i < 20; ++i) {
+            _v[i] = 0;
+        }
+    }
+
+    static inline guint32 premul_alpha(guint32 color, guint32 alpha)
+    {
+        guint32 temp = alpha * color + 128;
+        return (temp + (temp >> 8)) >> 8;
+    }
+
+    guint32 operator()(guint32 in) {
+        EXTRACT_ARGB32(in, a, r, g, b)
+        // we need to un-premultiply alpha values for this type of matrix
+        // TODO: unpremul can be ignored if there is an identity mapping on the alpha channel
+        if (a != 0) {
+            r = (r * 255 + a/2) / a;
+            b = (b * 255 + a/2) / a;
+            g = (g * 255 + a/2) / a;
+        }
+
+        gint32 ro = r*_v[0]  + g*_v[1]  + b*_v[2]  + a*_v[3]  + _v[4];
+        gint32 go = r*_v[5]  + g*_v[6]  + b*_v[7]  + a*_v[8]  + _v[9];
+        gint32 bo = r*_v[10] + g*_v[11] + b*_v[12] + a*_v[13] + _v[14];
+        gint32 ao = r*_v[15] + g*_v[16] + b*_v[17] + a*_v[18] + _v[19];
+        ro = (pxclamp(ro, 0, 255*255) + 127) / 255;
+        go = (pxclamp(go, 0, 255*255) + 127) / 255;
+        bo = (pxclamp(bo, 0, 255*255) + 127) / 255;
+        ao = (pxclamp(ao, 0, 255*255) + 127) / 255;
+
+        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:
+    gint32 _v[20];
+};
+
+struct ColorMatrixSaturate {
+    ColorMatrixSaturate(double v_in) {
+        // clamp parameter instead of clamping color values
+        double v = CLAMP(v_in, 0.0, 1.0);
+        _v[0] = 0.213+0.787*v; _v[1] = 0.715-0.715*v; _v[2] = 0.072-0.072*v;
+        _v[3] = 0.213-0.213*v; _v[4] = 0.715+0.285*v; _v[5] = 0.072-0.072*v;
+        _v[6] = 0.213-0.213*v; _v[7] = 0.715-0.715*v; _v[8] = 0.072+0.928*v;
+    }
+
+    guint32 operator()(guint32 in) {
+        EXTRACT_ARGB32(in, a, r, g, b)
+
+        // Note: this cannot be done in fixed point, because the loss of precision
+        //       causes overflow for some values of v
+        guint32 ro = r*_v[0] + g*_v[1] + b*_v[2] + 0.5;
+        guint32 go = r*_v[3] + g*_v[4] + b*_v[5] + 0.5;
+        guint32 bo = r*_v[6] + g*_v[7] + b*_v[8] + 0.5;
+
+        ASSEMBLE_ARGB32(pxout, a, ro, go, bo)
+        return pxout;
+    }
+private:
+    double _v[9];
+};
+
+struct ColorMatrixHueRotate {
+    ColorMatrixHueRotate(double v) {
+        double sinhue, coshue;
+        sincos(v * M_PI/180.0, &sinhue, &coshue);
+
+        _v[0] = round((0.213 +0.787*coshue -0.213*sinhue)*255);
+        _v[1] = round((0.715 -0.715*coshue -0.715*sinhue)*255);
+        _v[2] = round((0.072 -0.072*coshue +0.928*sinhue)*255);
+
+        _v[3] = round((0.213 -0.213*coshue +0.143*sinhue)*255);
+        _v[4] = round((0.715 +0.285*coshue +0.140*sinhue)*255);
+        _v[5] = round((0.072 -0.072*coshue -0.283*sinhue)*255);
+
+        _v[6] = round((0.213 -0.213*coshue -0.787*sinhue)*255);
+        _v[7] = round((0.715 -0.715*coshue +0.715*sinhue)*255);
+        _v[8] = round((0.072 +0.928*coshue +0.072*sinhue)*255);
+    }
+    guint32 operator()(guint32 in) {
+        EXTRACT_ARGB32(in, a, r, g, b)
+        gint32 maxpx = a*255;
+        gint32 ro = r*_v[0] + g*_v[1] + b*_v[2];
+        gint32 go = r*_v[3] + g*_v[4] + b*_v[5];
+        gint32 bo = r*_v[6] + g*_v[7] + b*_v[8];
+        ro = (pxclamp(ro, 0, maxpx) + 127) / 255;
+        go = (pxclamp(go, 0, maxpx) + 127) / 255;
+        bo = (pxclamp(bo, 0, maxpx) + 127) / 255;
+
+        ASSEMBLE_ARGB32(pxout, a, ro, go, bo)
+        return pxout;
+    }
+private:
+    gint32 _v[9];
+};
+
+struct ColorMatrixLuminanceToAlpha {
+    guint32 operator()(guint32 in) {
+        // original computation in double: r*0.2125 + g*0.7154 + b*0.0721
+        EXTRACT_ARGB32(in, a, r, g, b)
+        // unpremultiply color values
+        if (a != 0) {
+            r = (r * 255 + a/2) / a;
+            b = (b * 255 + a/2) / a;
+            g = (g * 255 + a/2) / a;
+        }
+        guint32 ao = r*54 + g*182 + b*18;
+        return ((ao + 127) / 255) << 24;
+    }
+};
+
+void FilterColorMatrix::render_cairo(FilterSlot &slot)
+{
+    cairo_surface_t *input = slot.getcairo(_input);
+    cairo_surface_t *out = NULL;
+    if (type == COLORMATRIX_LUMINANCETOALPHA) {
+        out = ink_cairo_surface_create_same_size(input, CAIRO_CONTENT_ALPHA);
+    } else {
+        out = ink_cairo_surface_create_identical(input);
+    }
+
+    switch (type) {
+    case COLORMATRIX_MATRIX:
+        ink_cairo_surface_filter(input, out, ColorMatrixMatrix(values));
+        break;
+    case COLORMATRIX_SATURATE:
+        ink_cairo_surface_filter(input, out, ColorMatrixSaturate(value));
+        break;
+    case COLORMATRIX_HUEROTATE:
+        ink_cairo_surface_filter(input, out, ColorMatrixHueRotate(value));
+        break;
+    case COLORMATRIX_LUMINANCETOALPHA:
+        ink_cairo_surface_filter(input, out, ColorMatrixLuminanceToAlpha());
+        break;
+    case COLORMATRIX_ENDTYPE:
+    default:
+        break;
+    }
+
+    slot.set(_output, out);
+    cairo_surface_destroy(out);
+}
+
+bool FilterColorMatrix::can_handle_affine(Geom::Matrix const &)
+{
+    return true;
+}
 
 void FilterColorMatrix::area_enlarge(NRRectL &/*area*/, Geom::Matrix const &/*trans*/)
 {
@@ -207,7 +373,7 @@ void FilterColorMatrix::set_value(gdouble v){
         value = v;
 }
 
-void FilterColorMatrix::set_values(std::vector<gdouble> &v){
+void FilterColorMatrix::set_values(std::vector<gdouble> const &v){
         values = v;
 }
 
diff --git a/src/display/nr-filter-colormatrix.h b/src/display/nr-filter-colormatrix.h
index 47b454c53..e3beb943d 100644
--- a/src/display/nr-filter-colormatrix.h
+++ b/src/display/nr-filter-colormatrix.h
@@ -34,11 +34,13 @@ public:
     static FilterPrimitive *create();
     virtual ~FilterColorMatrix();
 
-    virtual int render(FilterSlot &slot, FilterUnits const &units);
+    virtual void render_cairo(FilterSlot &slot);
+    virtual bool can_handle_affine(Geom::Matrix const &);
+
     virtual void area_enlarge(NRRectL &area, Geom::Matrix const &trans);
     virtual void set_type(FilterColorMatrixType type);
     virtual void set_value(gdouble value);
-    virtual void set_values(std::vector<gdouble> &values);
+    virtual void set_values(std::vector<gdouble> const &values);
 private:
     std::vector<gdouble> values;
     gdouble value;
diff --git a/src/display/nr-filter-composite.cpp b/src/display/nr-filter-composite.cpp
index 5c00898d4..de9905161 100644
--- a/src/display/nr-filter-composite.cpp
+++ b/src/display/nr-filter-composite.cpp
@@ -33,8 +33,8 @@ FilterPrimitive * FilterComposite::create() {
 FilterComposite::~FilterComposite()
 {}
 
-struct BlendArithmetic {
-    BlendArithmetic(double k1, double k2, double k3, double k4)
+struct ComposeArithmetic {
+    ComposeArithmetic(double k1, double k2, double k3, double k4)
         : _k1(round(k1 * 255))
         , _k2(round(k2 * 255*255))
         , _k3(round(k3 * 255*255))
@@ -44,16 +44,21 @@ struct BlendArithmetic {
         EXTRACT_ARGB32(in1, aa, ra, ga, ba)
         EXTRACT_ARGB32(in2, ab, rb, gb, bb)
 
-        guint32 ao = _k1*aa*ab + _k2*aa + _k3*ab + _k4; ao = (ao + 255*255) / (255*255);
-        guint32 ro = _k1*ra*rb + _k2*ra + _k3*rb + _k4; ro = (ro + 255*255) / (255*255);
-        guint32 go = _k1*ga*gb + _k2*ga + _k3*gb + _k4; go = (go + 255*255) / (255*255);
-        guint32 bo = _k1*ba*bb + _k2*ba + _k3*bb + _k4; bo = (bo + 255*255) / (255*255);
+        gint32 ao = _k1*aa*ab + _k2*aa + _k3*ab + _k4;
+        gint32 ro = _k1*ra*rb + _k2*ra + _k3*rb + _k4;
+        gint32 go = _k1*ga*gb + _k2*ga + _k3*gb + _k4;
+        gint32 bo = _k1*ba*bb + _k2*ba + _k3*bb + _k4;
+
+        ao = (pxclamp(ao, 0, 255*255*255) + (255*255/2)) / (255*255);
+        ro = (pxclamp(ro, 0, 255*255*255) + (255*255/2)) / (255*255);
+        go = (pxclamp(go, 0, 255*255*255) + (255*255/2)) / (255*255);
+        bo = (pxclamp(bo, 0, 255*255*255) + (255*255/2)) / (255*255);
 
         ASSEMBLE_ARGB32(pxout, ao, ro, go, bo)
         return pxout;
     }
 private:
-    guint32 _k1, _k2, _k3, _k4;
+    gint32 _k1, _k2, _k3, _k4;
 };
 
 void FilterComposite::render_cairo(FilterSlot &slot)
@@ -64,7 +69,7 @@ void FilterComposite::render_cairo(FilterSlot &slot)
     cairo_surface_t *out = ink_cairo_surface_create_output(input1, input2);
 
     if (op == COMPOSITE_ARITHMETIC) {
-        ink_cairo_surface_blend(input1, input2, out, BlendArithmetic(k1, k2, k3, k4));
+        ink_cairo_surface_blend(input1, input2, out, ComposeArithmetic(k1, k2, k3, k4));
     } else {
         ink_cairo_surface_blit(input2, out);
         cairo_t *ct = cairo_create(out);
diff --git a/src/display/nr-filter.cpp b/src/display/nr-filter.cpp
index a5b5801b1..667a3cc14 100644
--- a/src/display/nr-filter.cpp
+++ b/src/display/nr-filter.cpp
@@ -186,13 +186,14 @@ int Filter::render(NRArenaItem const *item, cairo_t *bgct, NRRectL const *bgarea
         units.set_automatic_resolution(true);
     }
 
-    /*units.set_paraller(false);
+    units.set_paraller(false);
+    Geom::Matrix pbtrans = units.get_matrix_display2pb();
     for (int i = 0 ; i < _primitive_count ; i++) {
-        if (_primitive[i]->get_input_traits() & TRAIT_PARALLER) {
+        if (!_primitive[i]->can_handle_affine(pbtrans)) {
             units.set_paraller(true);
             break;
         }
-    }*/
+    }
     units.set_paraller(true);
 
     FilterSlot slot(const_cast<NRArenaItem*>(item), bgct, bgarea, cairo_get_target(graphic), area, units);