Merge pull request #205 from ufo-kit/directional-phaseretrieval

phase retrieval: enable separate x and y distance

3 files changed, 64 insertions, 12 deletions

diff --git a/docs/filters.rst b/docs/filters.rst
index ff901bb..2906401 100644
--- a/docs/filters.rst
+++ b/docs/filters.rst
@@ -1235,6 +1235,10 @@ Phase retrieval
 
     Computes and applies a fourier filter to correct phase-shifted data.
     Expects frequencies as an input and produces frequencies as an output.
+    Propagation distance can be specified for both x and y directions together
+    by the :gobj:prop:`distance` parameter or separately by
+    :gobj:prop:`distance-x` and :gobj:prop:`distance-y`, which is useful e.g.
+    when pixel size is not symmetrical.
 
     .. gobj:prop:: method:enum
 
@@ -1248,6 +1252,14 @@ Phase retrieval
 
         Distance in meters.
 
+    .. gobj:prop:: distance-x:float
+
+        Distance in x-direction in meters.
+
+    .. gobj:prop:: distance-y:float
+
+        Distance in y-direction in meters.
+
     .. gobj:prop:: pixel-size:float
 
         Pixel size in meters.
diff --git a/src/kernels/phase-retrieval.cl b/src/kernels/phase-retrieval.cl
index 6f43efb..fa6db33 100644
--- a/src/kernels/phase-retrieval.cl
+++ b/src/kernels/phase-retrieval.cl
@@ -29,14 +29,14 @@
     float n_idy = (idy >= height >> 1) ? idy - height : idy; \
     n_idx = n_idx / width; \
     n_idy = n_idy / height; \
-    float sin_arg = prefac * (n_idy * n_idy + n_idx * n_idx); \
+    float sin_arg = prefac.x * (n_idx * n_idx) + prefac.y * (n_idy * n_idy); \
 
 #define COMMON_SETUP    \
     COMMON_SETUP_TIE;   \
     float sin_value = sin(sin_arg);
 
 kernel void
-tie_method(float prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
+tie_method(float2 prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
 {
     COMMON_SETUP_TIE;
     if (sin_arg >= frequency_cutoff)
@@ -46,7 +46,7 @@ tie_method(float prefac, float regularize_rate, float binary_filter_rate, float
 }
 
 kernel void
-ctf_method(float prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
+ctf_method(float2 prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
 {
     COMMON_SETUP;
 
@@ -57,7 +57,7 @@ ctf_method(float prefac, float regularize_rate, float binary_filter_rate, float
 }
 
 kernel void
-qp_method(float prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
+qp_method(float2 prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
 {
     COMMON_SETUP;
 
@@ -70,7 +70,7 @@ qp_method(float prefac, float regularize_rate, float binary_filter_rate, float f
 }
 
 kernel void
-qp2_method(float prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
+qp2_method(float2 prefac, float regularize_rate, float binary_filter_rate, float frequency_cutoff, global float *output)
 {
     COMMON_SETUP;
     float cacl_filter_value = 0.5f / (sin_value + pow(10, -regularize_rate));
diff --git a/src/ufo-retrieve-phase-task.c b/src/ufo-retrieve-phase-task.c
index 2e98b3d..a9a2035 100644
--- a/src/ufo-retrieve-phase-task.c
+++ b/src/ufo-retrieve-phase-task.c
@@ -47,14 +47,14 @@ static GEnumValue method_values[] = {
 struct _UfoRetrievePhaseTaskPrivate {
     Method method;
     gfloat energy;
-    gfloat distance;
+    gfloat distance, distance_x, distance_y;
     gfloat pixel_size;
     gfloat regularization_rate;
     gfloat binary_filter;
     gfloat frequency_cutoff;
     gboolean output_filter;
 
-    gfloat prefac;
+    gfloat prefac[2];
     cl_kernel *kernels;
     cl_kernel mult_by_value_kernel;
     cl_context context;
@@ -75,6 +75,8 @@ enum {
     PROP_METHOD,
     PROP_ENERGY,
     PROP_DISTANCE,
+    PROP_DISTANCE_X,
+    PROP_DISTANCE_Y,
     PROP_PIXEL_SIZE,
     PROP_REGULARIZATION_RATE,
     PROP_BINARY_FILTER_THRESHOLDING,
@@ -97,13 +99,23 @@ ufo_retrieve_phase_task_setup (UfoTask *task,
                                GError **error)
 {
     UfoRetrievePhaseTaskPrivate *priv;
-    gfloat lambda;
+    gfloat lambda, tmp;
 
     priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE (task);
     priv->context = ufo_resources_get_context (resources);
 
     lambda = 6.62606896e-34 * 299792458 / (priv->energy * 1.60217733e-16);
-    priv->prefac = G_PI * lambda * priv->distance / (priv->pixel_size * priv->pixel_size);
+    tmp = G_PI * lambda / (priv->pixel_size * priv->pixel_size);
+    if (priv->distance_x != 0.0f && priv->distance_y != 0.0f) {
+        priv->prefac[0] = tmp * priv->distance_x;
+        priv->prefac[1] = tmp * priv->distance_y;
+    } else if (priv->distance != 0.0f) {
+        priv->prefac[0] = priv->prefac[1] = tmp * priv->distance;
+    } else {
+        g_set_error (error, UFO_TASK_ERROR, UFO_TASK_ERROR_SETUP,
+                     "Either both, distance_x and distance_y must be non-zero, or distance must be non-zero");
+        return;
+    }
 
     priv->kernels[METHOD_TIE] = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "tie_method", NULL, error);
     priv->kernels[METHOD_CTF] = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "ctf_method", NULL, error);
@@ -200,7 +212,7 @@ ufo_retrieve_phase_task_process (UfoTask *task,
 
         method_kernel = priv->kernels[(gint)priv->method];
 
-        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 0, sizeof (gfloat), &priv->prefac));
+        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 0, sizeof (cl_float2), &priv->prefac));
         UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 1, sizeof (gfloat), &priv->regularization_rate));
         UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 2, sizeof (gfloat), &priv->binary_filter));
         UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 3, sizeof (gfloat), &priv->frequency_cutoff));
@@ -250,6 +262,12 @@ ufo_retrieve_phase_task_get_property (GObject *object,
         case PROP_DISTANCE:
             g_value_set_float (value, priv->distance);
             break;
+        case PROP_DISTANCE_X:
+            g_value_set_float (value, priv->distance_x);
+            break;
+        case PROP_DISTANCE_Y:
+            g_value_set_float (value, priv->distance_y);
+            break;
         case PROP_PIXEL_SIZE:
             g_value_set_float (value, priv->pixel_size);
             break;
@@ -289,6 +307,12 @@ ufo_retrieve_phase_task_set_property (GObject *object,
         case PROP_DISTANCE:
             priv->distance = g_value_get_float (value);
             break;
+        case PROP_DISTANCE_X:
+            priv->distance_x = g_value_get_float (value);
+            break;
+        case PROP_DISTANCE_Y:
+            priv->distance_y = g_value_get_float (value);
+            break;
         case PROP_PIXEL_SIZE:
             priv->pixel_size = g_value_get_float (value);
             break;
@@ -385,7 +409,21 @@ ufo_retrieve_phase_task_class_init (UfoRetrievePhaseTaskClass *klass)
         g_param_spec_float ("distance",
             "Distance value",
             "Distance value.",
-            0, G_MAXFLOAT, 0.945,
+            0, G_MAXFLOAT, 0.0f,
+            G_PARAM_READWRITE);
+
+    properties[PROP_DISTANCE_X] =
+        g_param_spec_float ("distance-x",
+            "Distance value for x-direction",
+            "Distance value for x-direction",
+            0, G_MAXFLOAT, 0.0f,
+            G_PARAM_READWRITE);
+
+    properties[PROP_DISTANCE_Y] =
+        g_param_spec_float ("distance-y",
+            "Distance value for y-direction",
+            "Distance value for y-direction",
+            0, G_MAXFLOAT, 0.0f,
             G_PARAM_READWRITE);
 
     properties[PROP_PIXEL_SIZE] =
@@ -436,7 +474,9 @@ ufo_retrieve_phase_task_init(UfoRetrievePhaseTask *self)
     self->priv = priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE(self);
     priv->method = METHOD_TIE;
     priv->energy = 20.0f;
-    priv->distance = 0.945f;
+    priv->distance = 0.0f;
+    priv->distance_x = 0.0f;
+    priv->distance_y = 0.0f;
     priv->pixel_size = 0.75e-6f;
     priv->regularization_rate = 2.5f;
     priv->binary_filter = 0.1f;