Merge pull request #181 from ufo-kit/add-horizontal-interpolator

Add HorizontalInterpolate task

6 files changed, 294 insertions, 0 deletions

diff --git a/docs/filters.rst b/docs/filters.rst
index d962b31..bd7044d 100644
--- a/docs/filters.rst
+++ b/docs/filters.rst
@@ -543,6 +543,23 @@ Non-local-means denoising
         Sigma influencing the Gaussian weighting.
 
 
+Horizontal interpolation
+------------------------
+
+.. gobj:class:: horizontal-interpolate
+
+    Interpolate masked values in rows of an image. For all pixels equal to one
+    in the mask, find the closest pixel where mask is zero to the left and right
+    and linearly interpolate the value in the current pixel based on the found
+    left and right values. If the mask goes to the left or right border of the
+    image and on the other side there are at least two non-masked pixels
+    :math:`x_1` and :math:`x_2`, compute the value in the current pixel
+    :math:`x` by (in case the mask goes to the right border, left is analogous)
+    :math:`f(x) = f(x_2) + (x - x_2) * (f(x_2) - f(x_1))`. In case there is only
+    one valid pixel on one of the borders and all the others are masked, use
+    that pixel's value in all the remaining ones.
+
+
 Stream transformations
 ======================
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ad55660..d4bc70c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -39,6 +39,7 @@ set(ufofilter_SRCS
     ufo-get-dup-circ-task.c
     ufo-gradient-task.c
     ufo-general-backproject-task.c
+    ufo-horizontal-interpolate-task.c
     ufo-ifft-task.c
     ufo-interpolate-task.c
     ufo-interpolate-stream-task.c
diff --git a/src/kernels/interpolator.cl b/src/kernels/interpolator.cl
index fa3d765..6600af9 100644
--- a/src/kernels/interpolator.cl
+++ b/src/kernels/interpolator.cl
@@ -49,3 +49,60 @@ interpolate_horizontally (global float *a,
     output[idy * width + idx + offset] = (1.0f - alpha) * a[idy * left_width + idx + offset] +
                                          alpha * weight * b[idy * right_width + idx];
 }
+
+kernel void
+interpolate_mask_horizontally (global float *input,
+                               global float *mask,
+                               global float *output)
+{
+   const int idx = get_global_id (0);
+   const int idy = get_global_id (1);
+   const int width = get_global_size (0);
+   const int offset = idy * width;
+   int left = idx, right = idx;
+   float span, diff;
+
+   if (mask[offset + idx]) {
+        while (left > -1 && mask[offset + left] > 0) {
+            left--;
+        }
+        while (right < width && mask[offset + right] > 0) {
+            right++;
+        }
+
+        if (left < 0) {
+            /* Mask spans to the left border, use only the right value */
+            if (right < width - 1) {
+                /* There are two valid pixels on the right, use gradient */
+                diff = input[offset + right] - input[offset + right + 1];
+                output[offset + idx] = input[offset + right] + diff * (right - idx);
+            } else if (right == width - 1) {
+                /* There is only one valid pixel on the right, which is the only
+                 * valid pixel in this row, so use it's value everywhere */
+                output[offset + idx] = input[offset + right];
+            } else {
+                /* All pixels in this row are invalid, just copy the input */
+                output[offset + idx] = input[offset + idx];
+            }
+        } else if (right >= width) {
+            /* Mask spans to the right border, use only the left value */
+            if (left > 0) {
+                /* There are two valid pixels on the left, use gradient */
+                diff = input[offset + left] - input[offset + left - 1];
+                output[offset + idx] = input[offset + left] + diff * (idx - left);
+            } else if (left == 0) {
+                /* There is only one valid pixel on the left, which is the only
+                 * valid pixel in this row, so use it's value everywhere */
+                output[offset + idx] = input[offset + left];
+            }
+        } else {
+            /* There are valid pixels on both sides, use standard linear
+             * interpolation */
+            span = (float) (right - left);
+            output[offset + idx] = (right - idx) / span * input[offset + left] +
+                                   (idx - left) / span * input[offset + right];
+        }
+   } else {
+       output[offset + idx] = input[offset + idx];
+   }
+}
diff --git a/src/meson.build b/src/meson.build
index e589ea1..eb981b2 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -33,6 +33,7 @@ plugins = [
     'forwardproject',
     'get-dup-circ',
     'gradient',
+    'horizontal-interpolate',
     'interpolate',
     'interpolate-stream',
     'loop',
diff --git a/src/ufo-horizontal-interpolate-task.c b/src/ufo-horizontal-interpolate-task.c
new file mode 100644
index 0000000..d97d951
--- /dev/null
+++ b/src/ufo-horizontal-interpolate-task.c
@@ -0,0 +1,165 @@
+/*
+ * Copyright (C) 2011-2015 Karlsruhe Institute of Technology
+ *
+ * This file is part of Ufo.
+ *
+ * This library is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation, either
+ * version 3 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifdef __APPLE__
+#include <OpenCL/cl.h>
+#else
+#include <CL/cl.h>
+#endif
+
+#include "ufo-horizontal-interpolate-task.h"
+
+
+struct _UfoHorizontalInterpolateTaskPrivate {
+    cl_kernel kernel;
+};
+
+static void ufo_task_interface_init (UfoTaskIface *iface);
+
+G_DEFINE_TYPE_WITH_CODE (UfoHorizontalInterpolateTask, ufo_horizontal_interpolate_task, UFO_TYPE_TASK_NODE,
+                         G_IMPLEMENT_INTERFACE (UFO_TYPE_TASK,
+                                                ufo_task_interface_init))
+
+#define UFO_HORIZONTAL_INTERPOLATE_TASK_GET_PRIVATE(obj) (G_TYPE_INSTANCE_GET_PRIVATE((obj), UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK, UfoHorizontalInterpolateTaskPrivate))
+
+UfoNode *
+ufo_horizontal_interpolate_task_new (void)
+{
+    return UFO_NODE (g_object_new (UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK, NULL));
+}
+
+static void
+ufo_horizontal_interpolate_task_setup (UfoTask *task,
+                       UfoResources *resources,
+                       GError **error)
+{
+    UfoHorizontalInterpolateTaskPrivate *priv = UFO_HORIZONTAL_INTERPOLATE_TASK_GET_PRIVATE (task);
+
+    priv->kernel = ufo_resources_get_kernel (resources, "interpolator.cl", "interpolate_mask_horizontally", NULL, error);
+
+    if (priv->kernel != NULL) {
+        UFO_RESOURCES_CHECK_SET_AND_RETURN (clRetainKernel (priv->kernel), error);
+    }
+}
+
+static void
+ufo_horizontal_interpolate_task_get_requisition (UfoTask *task,
+                                 UfoBuffer **inputs,
+                                 UfoRequisition *requisition,
+                                 GError **error)
+{
+    ufo_buffer_get_requisition (inputs[0], requisition);
+}
+
+static guint
+ufo_horizontal_interpolate_task_get_num_inputs (UfoTask *task)
+{
+    return 2;
+}
+
+static guint
+ufo_horizontal_interpolate_task_get_num_dimensions (UfoTask *task,
+                                                    guint input)
+{
+    return 2;
+}
+
+static UfoTaskMode
+ufo_horizontal_interpolate_task_get_mode (UfoTask *task)
+{
+    return UFO_TASK_MODE_PROCESSOR | UFO_TASK_MODE_GPU;
+}
+
+static gboolean
+ufo_horizontal_interpolate_task_equal_real (UfoNode *n1,
+                                            UfoNode *n2)
+{
+    g_return_val_if_fail (UFO_IS_HORIZONTAL_INTERPOLATE_TASK (n1) && UFO_IS_HORIZONTAL_INTERPOLATE_TASK (n2), FALSE);
+    return UFO_HORIZONTAL_INTERPOLATE_TASK (n1)->priv->kernel == UFO_HORIZONTAL_INTERPOLATE_TASK (n2)->priv->kernel;
+}
+
+static gboolean
+ufo_horizontal_interpolate_task_process (UfoTask *task,
+                                         UfoBuffer **inputs,
+                                         UfoBuffer *output,
+                                         UfoRequisition *requisition)
+{
+    UfoHorizontalInterpolateTaskPrivate *priv;
+    UfoGpuNode *node;
+    UfoProfiler *profiler;
+    cl_command_queue cmd_queue;
+    cl_mem mask_in_mem, in_mem, out_mem;
+
+    priv = UFO_HORIZONTAL_INTERPOLATE_TASK_GET_PRIVATE (task);
+    node = UFO_GPU_NODE (ufo_task_node_get_proc_node (UFO_TASK_NODE (task)));
+    cmd_queue = ufo_gpu_node_get_cmd_queue (node);
+    in_mem = ufo_buffer_get_device_array (inputs[0], cmd_queue);
+    mask_in_mem = ufo_buffer_get_device_array (inputs[1], cmd_queue);
+    out_mem = ufo_buffer_get_device_array (output, cmd_queue);
+    profiler = ufo_task_node_get_profiler (UFO_TASK_NODE (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), &mask_in_mem));
+    UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 2, sizeof (cl_mem), &out_mem));
+    ufo_profiler_call (profiler, cmd_queue, priv->kernel, 2, requisition->dims, NULL);
+
+    return TRUE;
+}
+
+static void
+ufo_horizontal_interpolate_task_finalize (GObject *object)
+{
+    UfoHorizontalInterpolateTaskPrivate *priv = UFO_HORIZONTAL_INTERPOLATE_TASK_GET_PRIVATE (object);
+
+    if (priv->kernel) {
+        UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->kernel));
+        priv->kernel = NULL;
+    }
+
+    G_OBJECT_CLASS (ufo_horizontal_interpolate_task_parent_class)->finalize (object);
+}
+
+static void
+ufo_task_interface_init (UfoTaskIface *iface)
+{
+    iface->setup = ufo_horizontal_interpolate_task_setup;
+    iface->get_num_inputs = ufo_horizontal_interpolate_task_get_num_inputs;
+    iface->get_num_dimensions = ufo_horizontal_interpolate_task_get_num_dimensions;
+    iface->get_mode = ufo_horizontal_interpolate_task_get_mode;
+    iface->get_requisition = ufo_horizontal_interpolate_task_get_requisition;
+    iface->process = ufo_horizontal_interpolate_task_process;
+}
+
+static void
+ufo_horizontal_interpolate_task_class_init (UfoHorizontalInterpolateTaskClass *klass)
+{
+    GObjectClass *oclass = G_OBJECT_CLASS (klass);
+    UfoNodeClass *node_class = UFO_NODE_CLASS (klass);
+
+    oclass->finalize = ufo_horizontal_interpolate_task_finalize;
+    node_class->equal = ufo_horizontal_interpolate_task_equal_real;
+
+    g_type_class_add_private (oclass, sizeof(UfoHorizontalInterpolateTaskPrivate));
+}
+
+static void
+ufo_horizontal_interpolate_task_init(UfoHorizontalInterpolateTask *self)
+{
+    self->priv = UFO_HORIZONTAL_INTERPOLATE_TASK_GET_PRIVATE(self);
+}
diff --git a/src/ufo-horizontal-interpolate-task.h b/src/ufo-horizontal-interpolate-task.h
new file mode 100644
index 0000000..50354dc
--- /dev/null
+++ b/src/ufo-horizontal-interpolate-task.h
@@ -0,0 +1,53 @@
+/*
+ * Copyright (C) 2011-2013 Karlsruhe Institute of Technology
+ *
+ * This file is part of Ufo.
+ *
+ * This library is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation, either
+ * version 3 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef __UFO_HORIZONTAL_INTERPOLATE_TASK_H
+#define __UFO_HORIZONTAL_INTERPOLATE_TASK_H
+
+#include <ufo/ufo.h>
+
+G_BEGIN_DECLS
+
+#define UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK             (ufo_horizontal_interpolate_task_get_type())
+#define UFO_HORIZONTAL_INTERPOLATE_TASK(obj)             (G_TYPE_CHECK_INSTANCE_CAST((obj), UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK, UfoHorizontalInterpolateTask))
+#define UFO_IS_HORIZONTAL_INTERPOLATE_TASK(obj)          (G_TYPE_CHECK_INSTANCE_TYPE((obj), UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK))
+#define UFO_HORIZONTAL_INTERPOLATE_TASK_CLASS(klass)     (G_TYPE_CHECK_CLASS_CAST((klass), UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK, UfoHorizontalInterpolateTaskClass))
+#define UFO_IS_HORIZONTAL_INTERPOLATE_TASK_CLASS(klass)  (G_TYPE_CHECK_CLASS_TYPE((klass), UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK))
+#define UFO_HORIZONTAL_INTERPOLATE_TASK_GET_CLASS(obj)   (G_TYPE_INSTANCE_GET_CLASS((obj), UFO_TYPE_HORIZONTAL_INTERPOLATE_TASK, UfoHorizontalInterpolateTaskClass))
+
+typedef struct _UfoHorizontalInterpolateTask           UfoHorizontalInterpolateTask;
+typedef struct _UfoHorizontalInterpolateTaskClass      UfoHorizontalInterpolateTaskClass;
+typedef struct _UfoHorizontalInterpolateTaskPrivate    UfoHorizontalInterpolateTaskPrivate;
+
+struct _UfoHorizontalInterpolateTask {
+    UfoTaskNode parent_instance;
+
+    UfoHorizontalInterpolateTaskPrivate *priv;
+};
+
+struct _UfoHorizontalInterpolateTaskClass {
+    UfoTaskNodeClass parent_class;
+};
+
+UfoNode  *ufo_horizontal_interpolate_task_new       (void);
+GType     ufo_horizontal_interpolate_task_get_type  (void);
+
+G_END_DECLS
+
+#endif