ankafilters/src/ufo-anka-backproject-task.c

/*
 * Copyright (C) 2011-2014 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/>.
 */
#include <stdio.h>
#include <math.h>
#include <glib.h>
#include <glib/gprintf.h>

#ifdef __APPLE__
#include <OpenCL/cl.h>
#else
#include <CL/cl.h>
#endif

#include "ufo-anka-backproject-task.h"

/* Copy only neccessary projection region */
/* TODO: make this a parameter? */
#define COPY_PROJECTION_REGION 1
#define EXTRACT_INT(region, index) g_value_get_int (g_value_array_get_nth ((region), (index)))
#define EXTRACT_FLOAT(region, index) g_value_get_float (g_value_array_get_nth ((region), (index)))
#define REGION_SIZE(region) ((EXTRACT_INT ((region), 2)) == 0) ? 0 : \
                            ((EXTRACT_INT ((region), 1) - EXTRACT_INT ((region), 0) - 1) /\
                            EXTRACT_INT ((region), 2) + 1)
#define PAD_TO_DIVIDE(dividend, divisor) ((dividend) + (divisor) - (dividend) % (divisor))


/**
 * SECTION:ufo-anka-backproject-task
 * @Short_description: Backproject projection by projection
 * @Title: anka_backproject
 *
 */

typedef enum {
    PARAM_Z,
    PARAM_CENTER,
    PARAM_LAMINO
} Param;

struct _UfoAnkaBackprojectTaskPrivate {
    /* private */
    gboolean generated;
    guint count;
    /* sine and cosine table size based on BURST */
    gsize table_size;

    /* OpenCL */
    cl_context context;
    cl_kernel vector_kernel;
    cl_kernel scalar_kernel;
    cl_sampler sampler;
    /* Buffered images for invoking backprojection on BURST projections at once.
     * We potentially don't need to copy the last image and can use the one from
     * framework directly but it seems to have no performance effects. */
    cl_mem images[BURST];

    /* properties */
    GValueArray *x_region;
    GValueArray *y_region;
    GValueArray *region;
    GValueArray *center;
    GValueArray *projection_offset;
    float sines[BURST], cosines[BURST];
    guint num_projections;
    gfloat overall_angle;
    gfloat tomo_angle;
    gfloat lamino_angle;
    gfloat z;
    Param parameter;
};

static void ufo_task_interface_init (UfoTaskIface *iface);

G_DEFINE_TYPE_WITH_CODE (UfoAnkaBackprojectTask, ufo_anka_backproject_task, UFO_TYPE_TASK_NODE,
                         G_IMPLEMENT_INTERFACE (UFO_TYPE_TASK,
                                                ufo_task_interface_init))

#define UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE(obj) (G_TYPE_INSTANCE_GET_PRIVATE((obj), UFO_TYPE_ANKA_BACKPROJECT_TASK, UfoAnkaBackprojectTaskPrivate))

enum {
    PROP_0,
    PROP_X_REGION,
    PROP_Y_REGION,
    PROP_Z,
    PROP_REGION,
    PROP_PROJECTION_OFFSET,
    PROP_CENTER,
    PROP_NUM_PROJECTIONS,
    PROP_OVERALL_ANGLE,
    PROP_TOMO_ANGLE,
    PROP_LAMINO_ANGLE,
    PROP_PARAMETER,
    N_PROPERTIES
};

static GParamSpec *properties[N_PROPERTIES] = { NULL, };

static inline void
swap (gint *first, gint *second) {
    gint tmp;

    tmp = *first;
    *first = *second;
    *second = tmp;
}

/**
 * Determine the left and right column to read from a projection at a given
 * tomographic angle.
 */
static void
determine_x_extrema (gfloat extrema[2], GValueArray *x_extrema, GValueArray *y_extrema,
                     gfloat tomo_angle, gfloat x_center)
{
    gfloat sin_tomo, cos_tomo;
    gint x_min, x_max, y_min, y_max;

    sin_tomo = sin (tomo_angle);
    cos_tomo = cos (tomo_angle);
    x_min = EXTRACT_INT (x_extrema, 0);
    /* The interval is right-opened when OpenCL indices for both x and y are generated, */
    /* so the last index doesn't count */
    x_max = EXTRACT_INT (x_extrema, 1) - 1;
    y_min = EXTRACT_INT (y_extrema, 0);
    y_max = EXTRACT_INT (y_extrema, 1) - 1;

    if (sin_tomo < 0) {
        swap (&y_min, &y_max);
    }
    if (cos_tomo < 0) {
        swap (&x_min, &x_max);
    }

    extrema[0] = cos_tomo * x_min + sin_tomo * y_min + x_center;
    /* +1 because extrema[1] will be accessed by interpolation
     * but the region in copying is right-open */
    extrema[1] = cos_tomo * x_max + sin_tomo * y_max + x_center + 1;
}

/**
 * Determine the top and bottom row to read from a projection at given
 * tomographic and laminographic angles.
 */
static void
determine_y_extrema (gfloat extrema[2], GValueArray *x_extrema, GValueArray *y_extrema,
                     gfloat z_extrema[2], gfloat tomo_angle, gfloat lamino_angle,
                     gfloat y_center)
{
    gfloat sin_tomo, cos_tomo, sin_lamino, cos_lamino;
    gint x_min, x_max, y_min, y_max;

    sin_tomo = sin (tomo_angle);
    cos_tomo = cos (tomo_angle);
    sin_lamino = sin (lamino_angle);
    cos_lamino = cos (lamino_angle);
    x_min = EXTRACT_INT (x_extrema, 0);
    x_max = EXTRACT_INT (x_extrema, 1) - 1;
    y_min = EXTRACT_INT (y_extrema, 0);
    y_max = EXTRACT_INT (y_extrema, 1) - 1;

    if (sin_tomo < 0) {
        swap (&x_min, &x_max);
    }
    if (cos_tomo > 0) {
        swap (&y_min, &y_max);
    }

    extrema[0] = sin_tomo * x_min - cos_tomo * y_min;
    extrema[1] = sin_tomo * x_max - cos_tomo * y_max;
    extrema[0] = extrema[0] * cos_lamino + z_extrema[0] * sin_lamino + y_center;
    extrema[1] = extrema[1] * cos_lamino + z_extrema[1] * sin_lamino + y_center + 1;
}

/**
 * clip:
 * @result: resulting clipped extrema
 * @extrema: (min, max)
 * @maximum: projection width or height
 *
 * Clip extrema to an allowed interval [0, projection width/height)
 */
static void
clip (gint result[2], gfloat extrema[2], gint maximum)
{
    result[0] = (gint) floorf (extrema[0]);
    result[1] = (gint) ceilf (extrema[1]);

    if (result[0] < 0) {
        result[0] = 0;
    }
    if (result[0] > maximum) {
        result[0] = maximum;
    }
    if (result[1] < 0) {
        result[1] = 0;
    }
    if (result[1] > maximum) {
        result[1] = maximum;
    }

    if (result[0] == result[1]) {
        if (result[1] < maximum) {
            result[1]++;
        } else if (result[0] > 0) {
            result[0]--;
        } else {
            g_warning ("Cannot extend");
        }
    } else if (result[0] > result[1]) {
        g_warning ("Invalid extrema: minimum larger than maximum");
    }
}

/**
 * Determine the left and right column to read from a projection at a given
 * tomographic angle. The result is bound to [0, projection width)
 */
static void
determine_x_region (gint result[2], GValueArray *x_extrema, GValueArray *y_extrema, gfloat tomo_angle,
                    gfloat x_center, gint width)
{
    gfloat extrema[2];

    determine_x_extrema (extrema, x_extrema, y_extrema, tomo_angle, x_center);

    clip (result, extrema, width);
}

/**
 * Determine the top and bottom column to read from a projection at given
 * tomographic and laminographic angles. The result is bound to
 * [0, projection height).
 */
static void
determine_y_region (gint result[2], GValueArray *x_extrema, GValueArray *y_extrema, gfloat z_extrema[2],
                    gfloat tomo_angle, gfloat lamino_angle, gfloat y_center, gint height)
{
    gfloat extrema[2];

    determine_y_extrema (extrema, x_extrema, y_extrema, z_extrema, tomo_angle,
                         lamino_angle, y_center);
    clip (result, extrema, height);
}

static void
set_region (GValueArray *src, GValueArray **dst)
{
    if (EXTRACT_INT (src, 0) > EXTRACT_INT (src, 1)) {
        g_log ("Ufo", G_LOG_LEVEL_CRITICAL,
               "Error <%s:%i>: Invalid region [\"from\", \"to\", \"step\"]: [%d, %d, %d], "\
               "\"from\" has to be less than or equal to \"to\"",
               __FILE__, __LINE__,
               EXTRACT_INT (src, 0), EXTRACT_INT (src, 1), EXTRACT_INT (src, 2));
    }
    else {
        g_value_array_free (*dst);
        *dst = g_value_array_copy (src);
    }
}

static void
copy_to_image (UfoBuffer *input,
               cl_mem output_image,
               cl_command_queue cmd_queue,
               size_t origin[3],
               size_t region[3],
               gint in_width)
{
    const UfoBufferLocation location = ufo_buffer_get_location (input);
    cl_mem input_data;
    gfloat *input_data_host;
    size_t src_offset;

    switch (location) {
        case UFO_BUFFER_LOCATION_HOST:
            input_data_host = ufo_buffer_get_host_array (input, NULL);
            src_offset = origin[1] * in_width + origin[0];
            UFO_RESOURCES_CHECK_CLERR (clEnqueueWriteImage (cmd_queue,
                                                            output_image,
                                                            CL_TRUE,
                                                            origin,
                                                            region,
                                                            0,
                                                            0,
                                                            input_data_host + src_offset,
                                                            0,
                                                            NULL,
                                                            NULL));
            break;
        case UFO_BUFFER_LOCATION_DEVICE:
            input_data = ufo_buffer_get_device_array (input, cmd_queue);
            src_offset = (origin[1] * in_width + origin[0]) * sizeof (cl_float);
            UFO_RESOURCES_CHECK_CLERR (clEnqueueCopyBufferToImage (cmd_queue,
                                                                   input_data,
                                                                   output_image,
                                                                   src_offset,
                                                                   origin,
                                                                   region,
                                                                   0,
                                                                   NULL,
                                                                   NULL));
            break;
        case UFO_BUFFER_LOCATION_DEVICE_IMAGE:
            input_data = ufo_buffer_get_device_image (input, cmd_queue);
            UFO_RESOURCES_CHECK_CLERR (clEnqueueCopyImage (cmd_queue,
                                                           input_data,
                                                           output_image,
                                                           origin,
                                                           origin,
                                                           region,
                                                           0,
                                                           NULL,
                                                           NULL));
            break;
        default:
            g_warning ("Invalid input buffer location");
            break;
    }
}


UfoNode *
ufo_anka_backproject_task_new (void)
{
    return UFO_NODE (g_object_new (UFO_TYPE_ANKA_BACKPROJECT_TASK, NULL));
}

static void
ufo_anka_backproject_task_setup (UfoTask *task,
                                 UfoResources *resources,
                                 GError **error)
{
    UfoAnkaBackprojectTaskPrivate *priv;
    cl_int cl_error;
    gint i;
    gchar *vector_kernel_name, *kernel_filename;

    vector_kernel_name = g_strdup_printf ("backproject_burst_%d", BURST);
    if (!vector_kernel_name) {
        g_warning ("Error making burst kernel name");
    }

    priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE (task);
    priv->context = ufo_resources_get_context (resources);
    switch (priv->parameter) {
        case PARAM_Z:
            kernel_filename = g_strdup ("z_kernel.cl");
            break;
        case PARAM_CENTER:
            kernel_filename = g_strdup ("center_kernel.cl");
            break;
        case PARAM_LAMINO:
            kernel_filename = g_strdup ("lamino_kernel.cl");
            break;
        default:
            g_warning ("Unkown varying parameter");
            break;
    }
    priv->vector_kernel = ufo_resources_get_kernel (resources, kernel_filename,
                                                    vector_kernel_name, error);
    priv->scalar_kernel = ufo_resources_get_kernel (resources, kernel_filename,
                                                    "backproject_burst_1", error);
    priv->sampler = clCreateSampler (priv->context,
                                     (cl_bool) FALSE,
                                     CL_ADDRESS_CLAMP,
                                     CL_FILTER_LINEAR,
                                     &cl_error);

    UFO_RESOURCES_CHECK_CLERR (clRetainContext (priv->context));
    UFO_RESOURCES_CHECK_CLERR (cl_error);
    if (priv->vector_kernel) {
        UFO_RESOURCES_CHECK_CLERR (clRetainKernel (priv->vector_kernel));
    }
    if (priv->scalar_kernel) {
        UFO_RESOURCES_CHECK_CLERR (clRetainKernel (priv->scalar_kernel));
    }

    for (i = 0; i < BURST; i++) {
        priv->images[i] = NULL;
    }

    switch (BURST) {
        case 1: priv->table_size = sizeof (cl_float); break;
        case 2: priv->table_size = sizeof (cl_float2); break;
        case 4: priv->table_size = sizeof (cl_float4); break;
        case 8: priv->table_size = sizeof (cl_float8); break;
        case 16: priv->table_size = sizeof (cl_float16); break;
        default: g_warning ("Unsupported vector size"); break;
    }

    g_free (vector_kernel_name);
    g_free (kernel_filename);
}

static void
ufo_anka_backproject_task_get_requisition (UfoTask *task,
                                           UfoBuffer **inputs,
                                           UfoRequisition *requisition)
{
    UfoAnkaBackprojectTaskPrivate *priv;
    gfloat start, stop, step;

    priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE (task);
    start = EXTRACT_FLOAT (priv->region, 0);
    stop = EXTRACT_FLOAT (priv->region, 1);
    step = EXTRACT_FLOAT (priv->region, 2);

    if (!priv->num_projections) {
        g_warning ("Number of projections has not been set");
    }
    if (step == 0.0f) {
        g_warning ("Step in region is 0");
    }

    requisition->n_dims = 3;
    requisition->dims[0] = REGION_SIZE (priv->x_region);
    requisition->dims[1] = REGION_SIZE (priv->y_region);
    requisition->dims[2] = (gint) ceil ((stop - start) / step);
}

static guint
ufo_anka_backproject_task_get_num_inputs (UfoTask *task)
{
    return 1;
}

static guint
ufo_anka_backproject_task_get_num_dimensions (UfoTask *task,
                                              guint input)
{
    g_return_val_if_fail (input == 0, 0);

    return 3;
}

static gboolean
ufo_anka_backproject_task_equal_real (UfoNode *n1,
                                      UfoNode *n2)
{
    g_return_val_if_fail (UFO_IS_ANKA_BACKPROJECT_TASK (n1) && UFO_IS_ANKA_BACKPROJECT_TASK (n2), FALSE);

    return UFO_ANKA_BACKPROJECT_TASK (n1)->priv->vector_kernel == UFO_ANKA_BACKPROJECT_TASK (n2)->priv->vector_kernel;
}

static UfoTaskMode
ufo_anka_backproject_task_get_mode (UfoTask *task)
{
    return UFO_TASK_MODE_REDUCTOR | UFO_TASK_MODE_GPU;
}

static gboolean
ufo_anka_backproject_task_process (UfoTask *task,
                                   UfoBuffer **inputs,
                                   UfoBuffer *output,
                                   UfoRequisition *requisition)
{
    UfoAnkaBackprojectTaskPrivate *priv;
    UfoRequisition in_req;
    UfoGpuNode *node;
    UfoProfiler *profiler;
    gfloat tomo_angle, *sines, *cosines;
    gint i, index;
    gint cumulate;
    gsize table_size;
    gboolean scalar;
    /* regions stripped off the "to" value */
    gfloat x_region[2], y_region[2], z_region[2], x_center[2], z_ends[2], lamino_angles[2],
           y_center, sin_lamino, cos_lamino, norm_factor;
    gint x_copy_region[2], y_copy_region[2];
    cl_kernel kernel;
    cl_command_queue cmd_queue;
    cl_mem out_mem;
    cl_int cl_error;
    /* image creation and copying */
    cl_image_format image_fmt;
    size_t origin[3];
    size_t region[3];
    /* keep the warp size satisfied but make sure the local grid is localized
     * around a point in 3D for efficient caching */
    const gint real_size[4] = {requisition->dims[0], requisition->dims[1], requisition->dims[2], 0};
    const gsize local_work_size[] = {16, 8, 8};
    gsize global_work_size[3];

    global_work_size[0] = requisition->dims[0] % local_work_size[0] ?
                          PAD_TO_DIVIDE (requisition->dims[0], local_work_size[0]) :
                          requisition->dims[0];
    global_work_size[1] = requisition->dims[1] % local_work_size[1] ?
                          PAD_TO_DIVIDE (requisition->dims[1], local_work_size[1]) :
                          requisition->dims[1];
    global_work_size[2] = requisition->dims[2] % local_work_size[2] ?
                          PAD_TO_DIVIDE (requisition->dims[2], local_work_size[2]) :
                          requisition->dims[2];

    priv = UFO_ANKA_BACKPROJECT_TASK (task)->priv;
    node = UFO_GPU_NODE (ufo_task_node_get_proc_node (UFO_TASK_NODE (task)));
    cmd_queue = ufo_gpu_node_get_cmd_queue (node);
    out_mem = ufo_buffer_get_device_array (output, cmd_queue);
    ufo_buffer_get_requisition (inputs[0], &in_req);

    index = priv->count % BURST;
    tomo_angle = priv->tomo_angle > -G_MAXFLOAT ? priv->tomo_angle :
                 priv->overall_angle * priv->count / priv->num_projections;
    norm_factor = priv->overall_angle / priv->num_projections;
    priv->sines[index] = sin (tomo_angle);
    priv->cosines[index] = cos (tomo_angle);
    x_region[0] = (gfloat) EXTRACT_INT (priv->x_region, 0);
    x_region[1] = (gfloat) EXTRACT_INT (priv->x_region, 2);
    y_region[0] = (gfloat) EXTRACT_INT (priv->y_region, 0);
    y_region[1] = (gfloat) EXTRACT_INT (priv->y_region, 2);
    if (priv->parameter == PARAM_Z) {
        z_ends[0] = z_region[0] = EXTRACT_FLOAT (priv->region, 0);
        z_region[1] = EXTRACT_FLOAT (priv->region, 2);
        z_ends[1] = EXTRACT_FLOAT (priv->region, 1);
    } else {
        z_ends[0] = priv->z;
        z_ends[1] = priv->z + 1.0f;
    }
    if (priv->parameter == PARAM_CENTER) {
        x_center[0] = EXTRACT_FLOAT (priv->region, 0) - EXTRACT_INT (priv->projection_offset, 0);
        x_center[1] = EXTRACT_FLOAT (priv->region, 2);
    } else {
        x_center[0] = x_center[1] = EXTRACT_FLOAT (priv->center, 0) - EXTRACT_INT (priv->projection_offset, 0);
    }
    y_center = EXTRACT_FLOAT (priv->center, 1) - EXTRACT_INT (priv->projection_offset, 1);
    if (priv->parameter == PARAM_LAMINO) {
        lamino_angles[0] = EXTRACT_FLOAT (priv->region, 0);
        lamino_angles[1] = EXTRACT_FLOAT (priv->region, 2);
    } else {
        lamino_angles[0] = lamino_angles[1] = priv->lamino_angle;
    }
    sin_lamino = sinf (priv->lamino_angle);
    cos_lamino = cosf (priv->lamino_angle);
    scalar = priv->count >= priv->num_projections / BURST * BURST ? 1 : 0;

    /* If COPY_PROJECTION_REGION is True we copy only the part necessary  */
    /* for a given tomographic and laminographic angle */
    /* TODO: Extend the region determination to be able to handle PARAM_LAMINO */
    if (COPY_PROJECTION_REGION && priv->parameter != PARAM_LAMINO) {
        determine_x_region (x_copy_region, priv->x_region, priv->y_region, tomo_angle,
                            EXTRACT_FLOAT (priv->center, 0), in_req.dims[0]);
        determine_y_region (y_copy_region, priv->x_region, priv->y_region, z_ends,
                            tomo_angle, priv->lamino_angle, EXTRACT_FLOAT (priv->center, 1),
                            in_req.dims[1]);
        origin[0] = x_copy_region[0];
        origin[1] = y_copy_region[0];
        origin[2] = 0;
        region[0] = x_copy_region[1] - x_copy_region[0];
        region[1] = y_copy_region[1] - y_copy_region[0];
    } else {
        origin[0] = origin[1] = origin[2] = 0;
        region[0] = in_req.dims[0];
        region[1] = in_req.dims[1];
    }
    region[2] = 1;

    if (priv->images[index] == NULL) {
        /* TODO: dangerous, don't rely on the ufo-buffer */
        image_fmt.image_channel_order = CL_R;
        image_fmt.image_channel_data_type = CL_FLOAT;
        /* TODO: what with the "other" API? */
        priv->images[index] = clCreateImage2D (priv->context,
                                               CL_MEM_READ_ONLY,
                                               &image_fmt,
                                               in_req.dims[0],
                                               in_req.dims[1],
                                               0,
                                               NULL,
                                               &cl_error);
        UFO_RESOURCES_CHECK_CLERR (cl_error);
    }

    copy_to_image (inputs[0], priv->images[index], cmd_queue, origin, region, in_req.dims[0]);

    if (scalar) {
        kernel = priv->scalar_kernel;
        cumulate = priv->count;
        table_size = sizeof (cl_float);
        sines = &priv->sines[index];
        cosines = &priv->cosines[index];
        i = 1;
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, 0, sizeof (cl_mem), &priv->images[index]));
    } else {
        kernel = priv->vector_kernel;
        cumulate = priv->count + 1 == BURST ? 0 : 1;
        table_size = priv->table_size;
        sines = priv->sines;
        cosines = priv->cosines;
        i = BURST;
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, index, sizeof (cl_mem), &priv->images[index]));
    }

    if (scalar || index == BURST - 1) {
        /* Execute the kernel after BURST images have arrived, i.e. we use more
         * projections at one invocation, so the number of read/writes to the
         * result is reduced by a factor of BURST. If there are not enough
         * projecttions left, execute the scalar kernel */
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_mem), &out_mem));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_sampler), &priv->sampler));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_int3), real_size));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float2), x_center));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float), (cl_float *) &y_center));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float2), x_region));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float2), y_region));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float2), z_region));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float2), lamino_angles));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float), &sin_lamino));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float), &cos_lamino));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, table_size, sines));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, table_size, cosines));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_float), &norm_factor));
        UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i, sizeof (cl_int), (cl_int *) &cumulate));

        profiler = ufo_task_node_get_profiler (UFO_TASK_NODE (task));
        ufo_profiler_call (profiler, cmd_queue, kernel, 3, global_work_size, local_work_size);
    }

    priv->count++;

    return TRUE;
}

static gboolean
ufo_anka_backproject_task_generate (UfoTask *task,
                                    UfoBuffer *output,
                                    UfoRequisition *requisition)
{
    UfoAnkaBackprojectTaskPrivate *priv;

    priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE (task);

    if (priv->generated) {
        return FALSE;
    }

    priv->generated = TRUE;

    return TRUE;
}

static void
ufo_anka_backproject_task_finalize (GObject *object)
{
    UfoAnkaBackprojectTaskPrivate *priv;
    gint i;

    priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE (object);
    g_value_array_free (priv->x_region);
    g_value_array_free (priv->y_region);
    g_value_array_free (priv->region);
    g_value_array_free (priv->projection_offset);
    g_value_array_free (priv->center);

    if (priv->vector_kernel) {
        UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->vector_kernel));
        priv->vector_kernel = NULL;
    }
    if (priv->scalar_kernel) {
        UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->scalar_kernel));
        priv->scalar_kernel = NULL;
    }
    if (priv->context) {
        UFO_RESOURCES_CHECK_CLERR (clReleaseContext (priv->context));
        priv->context = NULL;
    }
    if (priv->sampler) {
        UFO_RESOURCES_CHECK_CLERR (clReleaseSampler (priv->sampler));
        priv->sampler = NULL;
    }

    for (i = 0; i < BURST; i++) {
        if (priv->images[i] != NULL) {
            UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (priv->images[i]));
            priv->images[i] = NULL;
        }
    }

    G_OBJECT_CLASS (ufo_anka_backproject_task_parent_class)->finalize (object);
}

static void
ufo_task_interface_init (UfoTaskIface *iface)
{
    iface->setup = ufo_anka_backproject_task_setup;
    iface->get_requisition = ufo_anka_backproject_task_get_requisition;
    iface->get_num_inputs = ufo_anka_backproject_task_get_num_inputs;
    iface->get_num_dimensions = ufo_anka_backproject_task_get_num_dimensions;
    iface->get_mode = ufo_anka_backproject_task_get_mode;
    iface->process = ufo_anka_backproject_task_process;
    iface->generate = ufo_anka_backproject_task_generate;
}

static void
ufo_anka_backproject_task_set_property (GObject *object,
                                        guint property_id,
                                        const GValue *value,
                                        GParamSpec *pspec)
{
    UfoAnkaBackprojectTaskPrivate *priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE (object);
    GValueArray *array;


    switch (property_id) {
        case PROP_X_REGION:
            array = (GValueArray *) g_value_get_boxed (value);
            set_region (array, &priv->x_region);
            break;
        case PROP_Y_REGION:
            array = (GValueArray *) g_value_get_boxed (value);
            set_region (array, &priv->y_region);
            break;
        case PROP_Z:
            priv->z = g_value_get_float (value);
            break;
        case PROP_REGION:
            array = (GValueArray *) g_value_get_boxed (value);
            g_value_array_free (priv->region);
            priv->region = g_value_array_copy (array);
            break;
        case PROP_PROJECTION_OFFSET:
            array = (GValueArray *) g_value_get_boxed (value);
            g_value_array_free (priv->projection_offset);
            priv->projection_offset = g_value_array_copy (array);
            break;
        case PROP_CENTER:
            array = (GValueArray *) g_value_get_boxed (value);
            g_value_array_free (priv->center);
            priv->center = g_value_array_copy (array);
            break;
        case PROP_NUM_PROJECTIONS:
            priv->num_projections = g_value_get_uint (value);
            break;
        case PROP_OVERALL_ANGLE:
            priv->overall_angle = g_value_get_float (value);
            break;
        case PROP_TOMO_ANGLE:
            priv->tomo_angle = g_value_get_float (value);
            break;
        case PROP_LAMINO_ANGLE:
            priv->lamino_angle = g_value_get_float (value);
            break;
        case PROP_PARAMETER:
            if (!g_strcmp0 (g_value_get_string (value), "z")) {
                priv->parameter = PARAM_Z;
            } else if (!g_strcmp0 (g_value_get_string (value), "x-center")) {
                priv->parameter = PARAM_CENTER;
            } else if (!g_strcmp0 (g_value_get_string (value), "lamino-angle")) {
                priv->parameter = PARAM_LAMINO;
            }
            break;
        default:
            G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
            break;
    }
}

static void
ufo_anka_backproject_task_get_property (GObject *object,
                                        guint property_id,
                                        GValue *value,
                                        GParamSpec *pspec)
{
    UfoAnkaBackprojectTaskPrivate *priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE (object);

    switch (property_id) {
        case PROP_X_REGION:
            g_value_set_boxed (value, priv->x_region);
            break;
        case PROP_Y_REGION:
            g_value_set_boxed (value, priv->y_region);
            break;
        case PROP_Z:
            g_value_set_float (value, priv->z);
            break;
        case PROP_REGION:
            g_value_set_boxed (value, priv->region);
            break;
        case PROP_PROJECTION_OFFSET:
            g_value_set_boxed (value, priv->projection_offset);
            break;
        case PROP_CENTER:
            g_value_set_boxed (value, priv->center);
            break;
        case PROP_NUM_PROJECTIONS:
            g_value_set_uint (value, priv->num_projections);
            break;
        case PROP_OVERALL_ANGLE:
            g_value_set_float (value, priv->overall_angle);
            break;
        case PROP_TOMO_ANGLE:
            g_value_set_float (value, priv->tomo_angle);
            break;
        case PROP_LAMINO_ANGLE:
            g_value_set_float (value, priv->lamino_angle);
            break;
        case PROP_PARAMETER:
            switch (priv->parameter) {
                case PARAM_Z:
                    g_value_set_string (value, "z");
                    break;
                case PARAM_CENTER:
                    g_value_set_string (value, "x-center");
                    break;
                case PARAM_LAMINO:
                    g_value_set_string (value, "lamino-angle");
                    break;
            }
            break;
        default:
            G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
            break;
    }
}

static void
ufo_anka_backproject_task_class_init (UfoAnkaBackprojectTaskClass *klass)
{
    GObjectClass *oclass;
    UfoNodeClass *node_class;
    
    oclass = G_OBJECT_CLASS (klass);
    node_class = UFO_NODE_CLASS (klass);

    oclass->finalize = ufo_anka_backproject_task_finalize;
    oclass->set_property = ufo_anka_backproject_task_set_property;
    oclass->get_property = ufo_anka_backproject_task_get_property;


    GParamSpec *region_vals = g_param_spec_int ("region_values",
                                                "Region values",
                                                "Elements in regions",
                                                G_MININT,
                                                G_MAXINT,
                                                (gint) 0,
                                                G_PARAM_READWRITE);

    GParamSpec *float_region_vals = g_param_spec_float ("float_region_values",
                                                        "Float Region values",
                                                        "Elements in float regions",
                                                        -G_MAXFLOAT,
                                                        G_MAXFLOAT,
                                                        0.0f,
                                                        G_PARAM_READWRITE);

    properties[PROP_X_REGION] =
        g_param_spec_value_array ("x-region",
                                  "X region for reconstruction as (from, to, step)",
                                  "X region for reconstruction as (from, to, step)",
                                  region_vals,
                                  G_PARAM_READWRITE);

    properties[PROP_Y_REGION] =
        g_param_spec_value_array ("y-region",
                                  "Y region for reconstruction as (from, to, step)",
                                  "Y region for reconstruction as (from, to, step)",
                                  region_vals,
                                  G_PARAM_READWRITE);

    properties[PROP_Z] =
        g_param_spec_float ("z",
                            "Z coordinate of the reconstructed slice",
                            "Z coordinate of the reconstructed slice",
                            -G_MAXFLOAT,
                            G_MAXFLOAT,
                            0.0f,
                            G_PARAM_READWRITE);

    properties[PROP_REGION] =
        g_param_spec_value_array ("region",
                                  "Region for the parameter along z-axis as (from, to, step)",
                                  "Region for the parameter along z-axis as (from, to, step)",
                                  float_region_vals,
                                  G_PARAM_READWRITE);

    properties[PROP_PROJECTION_OFFSET] =
        g_param_spec_value_array ("projection-offset",
                                  "Offset to projection data as (x, y)",
                                  "Offset to projection data as (x, y) for the case input data \
                                  is cropped to the necessary range of interest",
                                  region_vals,
                                  G_PARAM_READWRITE);

    properties[PROP_CENTER] =
        g_param_spec_value_array ("center",
                                  "Center of the volume with respect to projections (x, y)",
                                  "Center of the volume with respect to projections (x, y), (rotation axes)",
                                  float_region_vals,
                                  G_PARAM_READWRITE);

    properties[PROP_OVERALL_ANGLE] =
        g_param_spec_float ("overall-angle",
                            "Angle covered by all projections",
                            "Angle covered by all projections (can be negative for negative steps "
                            "in case only num-projections is specified",
                            -G_MAXFLOAT,
                            G_MAXFLOAT,
                            G_PI,
                            G_PARAM_READWRITE);

    properties[PROP_NUM_PROJECTIONS] =
        g_param_spec_uint ("num-projections",
                          "Number of projections",
                          "Number of projections",
                          0,
                          16384,
                          0,
                          G_PARAM_READWRITE);

    properties[PROP_TOMO_ANGLE] =
        g_param_spec_float ("tomo-angle",
                            "Tomographic rotation angle in radians",
                            "Tomographic rotation angle in radians (used for acquiring projections)",
                            -G_MAXFLOAT,
                            G_MAXFLOAT,
                            0.0f,
                            G_PARAM_READWRITE);

    properties[PROP_LAMINO_ANGLE] =
        g_param_spec_float ("lamino-angle",
                            "Absolute laminogrpahic angle in radians",
                            "Absolute laminogrpahic angle in radians determining the sample tilt",
                            0.0f,
                            (float) G_PI / 2,
                            0.0f,
                            G_PARAM_READWRITE);

    properties[PROP_PARAMETER] =
        g_param_spec_string ("parameter",
                             "Which paramter will be varied along the z-axis",
                             "Which paramter will be varied along the z-axis, from \"z\", \"x-center\", \"lamino-angle\"",
                             "z",
                             G_PARAM_READWRITE);

    for (guint i = PROP_0 + 1; i < N_PROPERTIES; i++)
        g_object_class_install_property (oclass, i, properties[i]);

    node_class->equal = ufo_anka_backproject_task_equal_real;

    g_type_class_add_private (klass, sizeof(UfoAnkaBackprojectTaskPrivate));
}

static void
ufo_anka_backproject_task_init(UfoAnkaBackprojectTask *self)
{
    UfoAnkaBackprojectTaskPrivate *priv;
    self->priv = priv = UFO_ANKA_BACKPROJECT_TASK_GET_PRIVATE(self);
    guint i;
    GValue int_zero = G_VALUE_INIT;
    GValue float_zero = G_VALUE_INIT;

    g_value_init (&int_zero, G_TYPE_INT);
    g_value_init (&float_zero, G_TYPE_FLOAT);
    g_value_set_int (&int_zero, 0);
    g_value_set_float (&float_zero, 0.0f);
    self->priv->x_region = g_value_array_new (3);
    self->priv->y_region = g_value_array_new (3);
    self->priv->region = g_value_array_new (3);
    self->priv->z = 0.0f;
    self->priv->projection_offset = g_value_array_new (2);
    self->priv->center = g_value_array_new (2);

    for (i = 0; i < 3; i++) {
        g_value_array_insert (self->priv->x_region, i, &int_zero);
        g_value_array_insert (self->priv->y_region, i, &int_zero);
        g_value_array_insert (self->priv->region, i, &float_zero);
        if (i < 2) {
            g_value_array_insert (self->priv->projection_offset, i, &int_zero);
            g_value_array_insert (self->priv->center, i, &float_zero);
        }
    }

    self->priv->num_projections = 0;
    self->priv->overall_angle = G_PI;
    self->priv->tomo_angle = -G_MAXFLOAT;
    self->priv->lamino_angle = 0.0f;
    self->priv->parameter = PARAM_Z;
    self->priv->count = 0;
    self->priv->generated = FALSE;
}