image: generate Gaussian kernels.

src/image/filter/blur.cpp

@@ -34,98 +34,6 @@
 #define TEMPLATE_FILE "filter/blur.c"
 #include "pipi-template.h"
 
-/* Any standard deviation below this value will be rounded up, in order
- * to avoid ridiculously low values. exp(-1/(2*0.2*0.2)) is < 10^-5 so
- * there is little chance that any value below 0.2 will be useful. */
-#define BLUR_EPSILON 0.2
-
-pipi_image_t *pipi_gaussian_blur(pipi_image_t *src, float radius)
-{
-    return pipi_gaussian_blur_ext(src, radius, radius, 0.0, 0.0, 0.0);
-}
-
-pipi_image_t *pipi_gaussian_blur_ext(pipi_image_t *src, float rx, float ry,
-                                     float angle, float dx, float dy)
-{
-    pipi_image_t *ret;
-    double *kernel;
-    double Kx, Ky, t = 0.0, sint, cost, bbx, bby;
-    int i, j, krx, kry, m, n;
-
-    if(rx < BLUR_EPSILON) rx = BLUR_EPSILON;
-    if(ry < BLUR_EPSILON) ry = BLUR_EPSILON;
-
-    sint = sin(angle * M_PI / 180.);
-    cost = cos(angle * M_PI / 180.);
-
-    /* Compute the final ellipse's bounding box */
-    bbx = sqrt(rx * rx * cost * cost + ry * ry * sint * sint);
-    bby = sqrt(ry * ry * cost * cost + rx * rx * sint * sint);
-
-    /* FIXME: the kernel becomes far too big with large values of dx, because
-     * we grow both left and right. Fix the growing direction. */
-    krx = (int)(3. * bbx + .99999 + ceil(abs(dx)));
-    m = 2 * krx + 1;
-    Kx = -1. / (2. * rx * rx);
-
-    kry = (int)(3. * bby + .99999 + ceil(abs(dy)));
-    n = 2 * kry + 1;
-    Ky = -1. / (2. * ry * ry);
-
-    kernel = malloc(m * n * sizeof(double));
-
-    for(j = -kry; j <= kry; j++)
-    {
-        for(i = -krx; i <= krx; i++)
-        {
-            /* FIXME: this level of interpolation sucks. We should
-             * interpolate on the full NxN grid for better quality. */
-            static double const samples[] =
-            {
-                 .0,  .0, 1,
-                -.40, -.40, 0.8,
-                -.30,  .0,  0.9,
-                -.40,  .40, 0.8,
-                 .0,   .30, 0.9,
-                 .40,  .40, 0.8,
-                 .30,  .0,  0.9,
-                 .40, -.40, 0.8,
-                 .0,  -.30, 0.9,
-            };
-            double u, v, ex, ey, d = 0.;
-            unsigned int k;
-
-            for(k = 0; k < sizeof(samples) / sizeof(*samples) / 3; k++)
-            {
-                u = ((double)i + samples[k * 3]) * cost
-                     - ((double)j + samples[k * 3 + 1]) * sint + dx;
-                v = ((double)i + samples[k * 3]) * sint
-                     + ((double)j + samples[k * 3 + 1]) * cost + dy;
-                ex = Kx * u * u;
-                ey = Ky * v * v;
-                d += samples[k * 3 + 2] * exp(ex + ey);
-
-                /* Do not interpolate if this is a standard gaussian. */
-                if(!dx && !dy && !angle)
-                    break;
-            }
-
-            kernel[(j + kry) * m + i + krx] = d;
-            t += d;
-        }
-    }
-
-    for(j = 0; j < n; j++)
-        for(i = 0; i < m; i++)
-            kernel[j * m + i] /= t;
-
-    ret = pipi_convolution(src, m, n, kernel);
-
-    free(kernel);
-
-    return ret;
-}
-
 pipi_image_t *pipi_box_blur(pipi_image_t *src, int size)
 {
     return pipi_box_blur_ext(src, size, size);

src/image/stock.cpp

@@ -277,5 +277,84 @@ Array2D<float> Image::EdiffKernel(EdiffAlgorithm algorithm)
     return ret;
 }
 
+/* Any standard deviation below this value will be rounded up, in order
+ * to avoid ridiculously low values. exp(-1/(2*0.2*0.2)) is < 10^-5 so
+ * there is little chance that any value below 0.2 will be useful. */
+#define BLUR_EPSILON 0.2f
+
+Array2D<float> Image::GaussianKernel(vec2 radius, float angle, vec2 delta)
+{
+    Array2D<float> kernel;
+
+    if (radius.x < BLUR_EPSILON)
+        radius.x = BLUR_EPSILON;
+    if (radius.y < BLUR_EPSILON)
+        radius.y = BLUR_EPSILON;
+
+    float const sint = lol::sin(angle);
+    float const cost = lol::cos(angle);
+
+    /* Compute the final ellipse's bounding box */
+    float const bbx = lol::sqrt(sq(radius.x * cost) + sq(radius.y * sint));
+    float const bby = lol::sqrt(sq(radius.y * cost) + sq(radius.x * sint));
+
+    /* FIXME: the kernel becomes far too big with large values of dx, because
+     * we grow both left and right. Fix the growing direction. */
+    int const krx = (int)(3.f * bbx + .99999f + lol::ceil(lol::abs(delta.x)));
+    int const kry = (int)(3.f * bby + .99999f + lol::ceil(lol::abs(delta.y)));
+    ivec2 size(2 * krx + 1, 2 * kry + 1);
+    float const Kx = -1.f / (2.f * radius.x * radius.x);
+    float const Ky = -1.f / (2.f * radius.y * radius.y);
+
+    kernel.SetSize(size);
+
+    float t = 0.f;
+
+    for (int j = -kry; j <= kry; j++)
+    {
+        for (int i = -krx; i <= krx; i++)
+        {
+            /* FIXME: this level of interpolation sucks. We should
+             * interpolate on the full NxN grid for better quality. */
+            static vec3 const samples[] =
+            {
+                vec3( 0.0f,  0.0f, 1.0f),
+                vec3(-0.4f, -0.4f, 0.8f),
+                vec3(-0.3f,  0.0f, 0.9f),
+                vec3(-0.4f,  0.4f, 0.8f),
+                vec3( 0.0f,  0.3f, 0.9f),
+                vec3( 0.4f,  0.4f, 0.8f),
+                vec3( 0.3f,  0.0f, 0.9f),
+                vec3( 0.4f, -0.4f, 0.8f),
+                vec3( 0.0f, -0.3f, 0.9f),
+            };
+
+            float d = 0.f;
+
+            for (auto p : samples)
+            {
+                float u = (i + p.x) * cost - (j + p.y) * sint + delta.x;
+                float v = (i + p.x) * sint + (j + p.y) * cost + delta.y;
+                float ex = Kx * u * u;
+                float ey = Ky * v * v;
+                d += p.z * lol::exp(ex + ey);
+
+                /* Do not interpolate if this is a standard gaussian. */
+                if (!delta.x && !delta.y && !angle)
+                    break;
+            }
+
+            kernel[i + krx][j + kry] = d;
+            t += d;
+        }
+    }
+
+    for (int j = 0; j < size.y; j++)
+        for (int i = 0; i < size.x; i++)
+            kernel[i][j] *= (1.f / t);
+
+    return kernel;
+}
+
 } /* namespace lol */
 

src/lol/image/image.h

@@ -95,6 +95,7 @@ public:
     static Array2D<float> HalftoneKernel(ivec2 size);
     static Array2D<float> EdiffKernel(EdiffAlgorithm algorithm);
     static Array2D<float> NormalizeKernel(Array2D<float> const &kernel);
+    static Array2D<float> GaussianKernel(vec2 radius, float angle, vec2 delta);
 
     /* Rendering */
     bool Stock(char const *desc);