stripe-filter: add vertical sigma parameter

Which allows filtering of rings which are not perfectly straight or
which are broader. vertical-sigma parameter roughly specifies the
inclination from "perfectly vertical" stripes.

3 files changed, 56 insertions, 22 deletions

diff --git a/docs/filters.rst b/docs/filters.rst
index 2906401..5227f3b 100644
--- a/docs/filters.rst
+++ b/docs/filters.rst
@@ -930,17 +930,27 @@ Stripe filtering
 
     Filter vertical stripes. The input and output are in 2D frequency domain.
     The filter multiplies horizontal frequencies (for frequency ky=0) with a
-    Gaussian profile centered at 0 frequency.
+    Gaussian profile centered at 0 frequency if ``vertical-sigma`` is 0.
+    Otherwise it applies also a vertical Gaussian profile (1 - Gaussian), which
+    enables filtering of not perfectly vertical stripes, which is useful for
+    broader stripes and stripes which are not perfectly straight.
 
     Example usage::
 
         $ ufo-launch read path=sino.tif ! fft dimensions=2 ! filter-stripes sigma=1 ! ifft dimensions=2 ! write filename=sino-filtered.tif
 
-    .. gobj:prop:: sigma:float
+    .. gobj:prop:: horizontal-sigma:float
+
+        Horizontal filter strength, which is the sigma of the Gaussian. Small
+        values, e.g. 1e-7 cause only the zero frequency to remain in the
+        signal, i.e. stronger filtering. Values around 1 are a good starting
+        point.
+
+    .. gobj:prop:: vertical-sigma:float
 
-        Filter strength, which is the sigma of the gaussian. Small values, e.g.
-        1e-7 cause only the zero frequency to remain in the signal, i.e.
-        stronger filtering.
+        Vertical filter strength, which is the sigma of the Gaussian. The larger
+        the value, the more non-vertical frequencies are removed. Value around 4
+        is a good starting point.
 
 
 1D stripe filtering
diff --git a/src/kernels/filter.cl b/src/kernels/filter.cl
index ddc91f6..186b2e6 100644
--- a/src/kernels/filter.cl
+++ b/src/kernels/filter.cl
@@ -31,7 +31,8 @@ filter (global float *input,
 kernel void
 stripe_filter (global float *input,
                global float *output,
-               const float sigma)
+               const float horizontal_sigma,
+               const float vertical_sigma)
 {
     const int idx = get_global_id(0);
     const int idy = get_global_id(1);
@@ -41,11 +42,17 @@ stripe_filter (global float *input,
     /* Swap frequencies for interleaved complex input */
     /* No need to negate the second half of frequencies for Gaussian */
     const float x = idx >= width >> 1 ? (width >> 1) - (idx >> 1) : idx >> 1;
-    const float weight = exp (- x * x / (2 * sigma * sigma));
+    const float x_weight = exp (- x * x / (2 * horizontal_sigma * horizontal_sigma));
 
-    if (idy == 0) {
-        output[index] = input[index] * weight;
+    if (vertical_sigma == 0.0f) {
+        if (idy == 0) {
+            output[index] = input[index] * x_weight;
+        } else {
+            output[index] = input[index];
+        }
     } else {
-        output[index] = input[index];
+        const float y = idy >= height >> 1 ? height - idy : idy;
+        const float y_weight = exp (- y * y / (2 * vertical_sigma * vertical_sigma));
+        output[index] = input[index] * (1 - y_weight * (1 - x_weight));
     }
 }
diff --git a/src/ufo-filter-stripes-task.c b/src/ufo-filter-stripes-task.c
index bd56b7e..bb87498 100644
--- a/src/ufo-filter-stripes-task.c
+++ b/src/ufo-filter-stripes-task.c
@@ -29,7 +29,8 @@
 
 
 struct _UfoFilterStripesTaskPrivate {
-    gfloat sigma;
+    gfloat horizontal_sigma;
+    gfloat vertical_sigma;
     cl_kernel kernel;
 };
 
@@ -43,7 +44,8 @@ G_DEFINE_TYPE_WITH_CODE (UfoFilterStripesTask, ufo_filter_stripes_task, UFO_TYPE
 
 enum {
     PROP_0,
-    PROP_SIGMA,
+    PROP_HORIZONTAL_SIGMA,
+    PROP_VERTICAL_SIGMA,
     N_PROPERTIES
 };
 
@@ -76,7 +78,8 @@ ufo_filter_stripes_task_process (UfoTask *task,
 
     UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 0, sizeof (cl_mem), &in_mem));
     UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 1, sizeof (cl_mem), &out_mem));
-    UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 2, sizeof (cl_float), &priv->sigma));
+    UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 2, sizeof (cl_float), &priv->horizontal_sigma));
+    UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 3, sizeof (cl_float), &priv->vertical_sigma));
 
     profiler = ufo_task_node_get_profiler (UFO_TASK_NODE (task));
     ufo_profiler_call (profiler, cmd_queue, priv->kernel, 2, requisition->dims, NULL);
@@ -163,8 +166,11 @@ ufo_filter_stripes_task_set_property (GObject *object,
     UfoFilterStripesTaskPrivate *priv = UFO_FILTER_STRIPES_TASK_GET_PRIVATE (object);
 
     switch (property_id) {
-        case PROP_SIGMA:
-            priv->sigma = g_value_get_float (value);
+        case PROP_HORIZONTAL_SIGMA:
+            priv->horizontal_sigma = g_value_get_float (value);
+            break;
+        case PROP_VERTICAL_SIGMA:
+            priv->vertical_sigma = g_value_get_float (value);
             break;
         default:
             G_OBJECT_WARN_INVALID_PROPERTY_ID(object, property_id, pspec);
@@ -181,8 +187,11 @@ ufo_filter_stripes_task_get_property (GObject *object,
     UfoFilterStripesTaskPrivate *priv = UFO_FILTER_STRIPES_TASK_GET_PRIVATE (object);
 
     switch (property_id) {
-        case PROP_SIGMA:
-            g_value_set_float (value, priv->sigma);
+        case PROP_HORIZONTAL_SIGMA:
+            g_value_set_float (value, priv->horizontal_sigma);
+            break;
+        case PROP_VERTICAL_SIGMA:
+            g_value_set_float (value, priv->vertical_sigma);
             break;
         default:
             G_OBJECT_WARN_INVALID_PROPERTY_ID(object, property_id, pspec);
@@ -201,13 +210,20 @@ ufo_filter_stripes_task_class_init (UfoFilterStripesTaskClass *klass)
     oclass->set_property = ufo_filter_stripes_task_set_property;
     oclass->get_property = ufo_filter_stripes_task_get_property;
 
-    properties[PROP_SIGMA] =
-        g_param_spec_float ("sigma",
-            "Sigma of the gaussian window",
-            "Sigma of the gaussian window",
+    properties[PROP_HORIZONTAL_SIGMA] =
+        g_param_spec_float ("horizontal-sigma",
+            "Sigma of the gaussian window in the horizontal direction",
+            "Sigma of the gaussian window in the horizontal direction",
             0.0, G_MAXFLOAT, 1e-7,
             G_PARAM_READWRITE);
 
+    properties[PROP_VERTICAL_SIGMA] =
+        g_param_spec_float ("vertical-sigma",
+            "Sigma of the gaussian window in the vertical direction",
+            "Sigma of the gaussian window in the vertical direction",
+            0.0, G_MAXFLOAT, 0.0f,
+            G_PARAM_READWRITE);
+
     for (guint i = PROP_0 + 1; i < N_PROPERTIES; i++)
         g_object_class_install_property (oclass, i, properties[i]);
 
@@ -220,5 +236,6 @@ ufo_filter_stripes_task_init (UfoFilterStripesTask *self)
     UfoFilterStripesTaskPrivate *priv;
     self->priv = priv = UFO_FILTER_STRIPES_TASK_GET_PRIVATE (self);
     priv->kernel = NULL;
-    priv->sigma = 1e-7;
+    priv->horizontal_sigma = 1e-7;
+    priv->vertical_sigma = 0.0f;
 }