Merge pull request #179 from ufo-kit/filter-stripes-sigma

Filter stripes sigma

3 files changed, 86 insertions, 4 deletions

diff --git a/docs/filters.rst b/docs/filters.rst
index 96a40ea..d962b31 100644
--- a/docs/filters.rst
+++ b/docs/filters.rst
@@ -801,6 +801,26 @@ Frequency filtering
         Theta parameter of Faris-Byer filter.
 
 
+Stripe filtering
+----------------
+
+.. gobj:class:: filter-stripes
+
+    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.
+
+    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
+
+        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.
+
+
 1D stripe filtering
 -------------------
 
diff --git a/src/kernels/filter.cl b/src/kernels/filter.cl
index 445a735..ddc91f6 100644
--- a/src/kernels/filter.cl
+++ b/src/kernels/filter.cl
@@ -30,16 +30,22 @@ filter (global float *input,
 
 kernel void
 stripe_filter (global float *input,
-               global float *output)
+               global float *output,
+               const float sigma)
 {
     const int idx = get_global_id(0);
     const int idy = get_global_id(1);
     const int width = get_global_size(0);
     const int height = get_global_size(1);
     const int index = idy * get_global_size(0) + idx;
+    /* 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));
 
-    if (((idy < 1) || (idy >= height - 1)) && (idx > 0) && (idx < width - 1))
-        output[index] = 0.0f;
-    else
+    if (idy == 0) {
+        output[index] = input[index] * weight;
+    } else {
         output[index] = input[index];
+    }
 }
diff --git a/src/ufo-filter-stripes-task.c b/src/ufo-filter-stripes-task.c
index 3b6aebd..7b07c5d 100644
--- a/src/ufo-filter-stripes-task.c
+++ b/src/ufo-filter-stripes-task.c
@@ -28,6 +28,7 @@
 
 
 struct _UfoFilterStripesTaskPrivate {
+    gfloat sigma;
     cl_kernel kernel;
 };
 
@@ -41,9 +42,12 @@ G_DEFINE_TYPE_WITH_CODE (UfoFilterStripesTask, ufo_filter_stripes_task, UFO_TYPE
 
 enum {
     PROP_0,
+    PROP_SIGMA,
     N_PROPERTIES
 };
 
+static GParamSpec *properties[N_PROPERTIES] = { NULL, };
+
 UfoNode *
 ufo_filter_stripes_task_new (void)
 {
@@ -71,6 +75,7 @@ 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));
 
     profiler = ufo_task_node_get_profiler (UFO_TASK_NODE (task));
     ufo_profiler_call (profiler, cmd_queue, priv->kernel, 2, requisition->dims, NULL);
@@ -149,12 +154,62 @@ ufo_task_interface_init (UfoTaskIface *iface)
 }
 
 static void
+ufo_filter_stripes_task_set_property (GObject *object,
+                                      guint property_id,
+                                      const GValue *value,
+                                      GParamSpec *pspec)
+{
+    UfoFilterStripesTaskPrivate *priv = UFO_FILTER_STRIPES_TASK_GET_PRIVATE (object);
+
+    switch (property_id) {
+        case PROP_SIGMA:
+            priv->sigma = g_value_get_float (value);
+            break;
+        default:
+            G_OBJECT_WARN_INVALID_PROPERTY_ID(object, property_id, pspec);
+            break;
+    }
+}
+
+static void
+ufo_filter_stripes_task_get_property (GObject *object,
+                                      guint property_id,
+                                      GValue *value,
+                                      GParamSpec *pspec)
+{
+    UfoFilterStripesTaskPrivate *priv = UFO_FILTER_STRIPES_TASK_GET_PRIVATE (object);
+
+    switch (property_id) {
+        case PROP_SIGMA:
+            g_value_set_float (value, priv->sigma);
+            break;
+        default:
+            G_OBJECT_WARN_INVALID_PROPERTY_ID(object, property_id, pspec);
+            break;
+    }
+}
+
+static void
 ufo_filter_stripes_task_class_init (UfoFilterStripesTaskClass *klass)
 {
     GObjectClass *oclass;
 
     oclass = G_OBJECT_CLASS (klass);
+
     oclass->finalize = ufo_filter_stripes_task_finalize;
+    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",
+            0.0, G_MAXFLOAT, 1e-7,
+            G_PARAM_READWRITE);
+
+    for (guint i = PROP_0 + 1; i < N_PROPERTIES; i++)
+        g_object_class_install_property (oclass, i, properties[i]);
+
     g_type_class_add_private(klass, sizeof(UfoFilterStripesTaskPrivate));
 }
 
@@ -164,4 +219,5 @@ 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;
 }