/*
* 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 .
*/
#include "config.h"
#include
#ifdef __APPLE__
#include
#else
#include
#endif
#include "ufo-retrieve-phase-task.h"
#define IS_POW_OF_2(x) !(x & (x - 1))
typedef enum {
METHOD_TIE = 0,
METHOD_CTF,
METHOD_CTF_MULTI,
METHOD_QP,
METHOD_QP2,
N_METHODS
} Method;
static GEnumValue method_values[] = {
{ METHOD_TIE, "METHOD_TIE", "tie" },
{ METHOD_CTF, "METHOD_CTF", "ctf" },
{ METHOD_CTF_MULTI, "METHOD_CTF_MULTI", "ctf_multidistance" },
{ METHOD_QP, "METHOD_QP", "qp" },
{ METHOD_QP2, "METHOD_QP2", "qp2" },
{ 0, NULL, NULL}
};
struct _UfoRetrievePhaseTaskPrivate {
Method method;
gfloat energy;
GValueArray *distance;
gfloat distance_x, distance_y;
gfloat pixel_size;
gfloat regularization_rate;
gfloat binary_filter;
gfloat frequency_cutoff;
gboolean output_filter;
gfloat prefac[2];
cl_kernel *kernels;
cl_kernel mult_by_value_kernel, ctf_multi_apply_dist_kernel;
cl_context context;
UfoBuffer *filter_buffer;
};
static void ufo_task_interface_init (UfoTaskIface *iface);
G_DEFINE_TYPE_WITH_CODE (UfoRetrievePhaseTask, ufo_retrieve_phase_task, UFO_TYPE_TASK_NODE,
G_IMPLEMENT_INTERFACE (UFO_TYPE_TASK,
ufo_task_interface_init))
#define UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE(obj) (G_TYPE_INSTANCE_GET_PRIVATE((obj), UFO_TYPE_RETRIEVE_PHASE_TASK, UfoRetrievePhaseTaskPrivate))
enum {
PROP_0,
PROP_METHOD,
PROP_ENERGY,
PROP_DISTANCE,
PROP_DISTANCE_X,
PROP_DISTANCE_Y,
PROP_PIXEL_SIZE,
PROP_REGULARIZATION_RATE,
PROP_BINARY_FILTER_THRESHOLDING,
PROP_FREQUENCY_CUTOFF,
PROP_OUTPUT_FILTER,
N_PROPERTIES
};
static GParamSpec *properties[N_PROPERTIES] = { NULL, };
UfoNode *
ufo_retrieve_phase_task_new (void)
{
return UFO_NODE (g_object_new (UFO_TYPE_RETRIEVE_PHASE_TASK, NULL));
}
static void
ufo_retrieve_phase_task_setup (UfoTask *task,
UfoResources *resources,
GError **error)
{
UfoRetrievePhaseTaskPrivate *priv;
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);
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->n_values > 0) {
priv->prefac[0] = priv->prefac[1] = tmp * g_value_get_double (g_value_array_get_nth (priv->distance, 0));
} 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 specified");
return;
}
if (priv->distance->n_values > 1 && priv->method != METHOD_CTF_MULTI) {
g_set_error (error, UFO_TASK_ERROR, UFO_TASK_ERROR_SETUP,
"When multiple distances are speicified method must be set to \"ctf_multidistance\"");
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);
priv->kernels[METHOD_CTF_MULTI] = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "ctf_multidistance_square", NULL, error);
priv->kernels[METHOD_QP] = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "qp_method", NULL, error);
priv->kernels[METHOD_QP2] = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "qp2_method", NULL, error);
priv->mult_by_value_kernel = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "mult_by_value", NULL, error);
priv->ctf_multi_apply_dist_kernel = ufo_resources_get_kernel (resources, "phase-retrieval.cl", "ctf_multidistance_apply_distance", NULL, error);
UFO_RESOURCES_CHECK_SET_AND_RETURN (clRetainContext(priv->context), error);
if (priv->filter_buffer == NULL) {
UfoRequisition requisition;
requisition.n_dims = 2;
requisition.dims[0] = 1;
requisition.dims[1] = 1;
priv->filter_buffer = ufo_buffer_new(&requisition, priv->context);
}
for (int i = 0; i < N_METHODS; i++) {
if (priv->kernels[i] != NULL) {
UFO_RESOURCES_CHECK_SET_AND_RETURN (clRetainKernel(priv->kernels[i]), error);
}
}
if (priv->mult_by_value_kernel != NULL) {
UFO_RESOURCES_CHECK_SET_AND_RETURN (clRetainKernel (priv->mult_by_value_kernel), error);
}
if (priv->ctf_multi_apply_dist_kernel != NULL) {
UFO_RESOURCES_CHECK_SET_AND_RETURN (clRetainKernel (priv->ctf_multi_apply_dist_kernel), error);
}
}
static void
ufo_retrieve_phase_task_get_requisition (UfoTask *task,
UfoBuffer **inputs,
UfoRequisition *requisition,
GError **error)
{
UfoRetrievePhaseTaskPrivate *priv;
priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE (task);
ufo_buffer_get_requisition (inputs[0], requisition);
if (priv->output_filter) {
requisition->dims[0] >>= 1;
}
if (!IS_POW_OF_2 (requisition->dims[0]) || !IS_POW_OF_2 (requisition->dims[1])) {
g_set_error_literal (error, UFO_TASK_ERROR, UFO_TASK_ERROR_GET_REQUISITION,
"Please, perform zeropadding of your dataset along both directions (width, height) up to length of power of 2 (e.g. 256, 512, 1024, 2048, etc.)");
}
}
static guint
ufo_filter_task_get_num_inputs (UfoTask *task)
{
UfoRetrievePhaseTaskPrivate *priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE (task);
return MAX(1, priv->distance->n_values);
}
static guint
ufo_filter_task_get_num_dimensions (UfoTask *task, guint input)
{
return 2;
}
static UfoTaskMode
ufo_filter_task_get_mode (UfoTask *task)
{
return UFO_TASK_MODE_PROCESSOR | UFO_TASK_MODE_GPU;
}
static gboolean
ufo_retrieve_phase_task_process (UfoTask *task,
UfoBuffer **inputs,
UfoBuffer *output,
UfoRequisition *requisition)
{
UfoRetrievePhaseTaskPrivate *priv;
UfoGpuNode *node;
UfoProfiler *profiler;
gsize global_work_size[3];
cl_int cl_err;
guint i;
gfloat *distances, lambda, distance;
gfloat fill_pattern = 0.0f;
cl_mem current_in_mem, in_sum_mem, in_mem, out_mem, filter_mem, distances_mem;
cl_kernel method_kernel;
cl_command_queue cmd_queue;
priv = UFO_RETRIEVE_PHASE_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);
profiler = ufo_task_node_get_profiler (UFO_TASK_NODE (task));
global_work_size[0] = requisition->dims[0];
global_work_size[1] = requisition->dims[1];
if (!priv->output_filter) {
/* Filter is real as opposed to the complex input, so the width is only half of the interleaved input */
global_work_size[0] >>= 1;
}
if (ufo_buffer_cmp_dimensions (priv->filter_buffer, requisition) != 0) {
ufo_buffer_resize (priv->filter_buffer, requisition);
filter_mem = ufo_buffer_get_device_array (priv->filter_buffer, cmd_queue);
method_kernel = priv->kernels[(gint)priv->method];
if (priv->method == METHOD_CTF_MULTI) {
lambda = 6.62606896e-34 * 299792458 / (priv->energy * 1.60217733e-16);
distances = g_malloc0 (priv->distance->n_values * sizeof (gfloat));
for (i = 0; i < priv->distance->n_values; i++) {
distances[i] = g_value_get_double (g_value_array_get_nth (priv->distance, i));
}
distances_mem = clCreateBuffer (priv->context,
CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
priv->distance->n_values * sizeof(float),
distances,
&cl_err);
UFO_RESOURCES_CHECK_CLERR (cl_err);
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 0, sizeof (cl_mem), &distances_mem));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 1, sizeof (guint), &priv->distance->n_values));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 2, sizeof (gfloat), &lambda));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 3, sizeof (gfloat), &priv->pixel_size));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 4, sizeof (gfloat), &priv->regularization_rate));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 5, sizeof (cl_mem), &filter_mem));
g_free (distances);
} else {
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));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (method_kernel, 4, sizeof (cl_mem), &filter_mem));
}
ufo_profiler_call (profiler, cmd_queue, method_kernel, requisition->n_dims, global_work_size, NULL);
if (priv->method == METHOD_CTF_MULTI) {
UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (distances_mem));
}
}
else {
filter_mem = ufo_buffer_get_device_array (priv->filter_buffer, cmd_queue);
}
if (priv->output_filter) {
ufo_buffer_copy (priv->filter_buffer, output);
}
else {
if (priv->method == METHOD_CTF_MULTI) {
/* Sum input frequencies first, then proceed wth complex division */
in_sum_mem = clCreateBuffer (priv->context,
CL_MEM_READ_WRITE,
requisition->dims[0] * requisition->dims[1] * sizeof(float),
NULL,
&cl_err);
UFO_RESOURCES_CHECK_CLERR (cl_err);
UFO_RESOURCES_CHECK_CLERR (clEnqueueFillBuffer (cmd_queue, in_sum_mem, &fill_pattern, sizeof (cl_float),
0, requisition->dims[0] * requisition->dims[1] * sizeof(float),
0, NULL, NULL));
for (i = 0; i < priv->distance->n_values; i++) {
distance = g_value_get_double (g_value_array_get_nth (priv->distance, i));
current_in_mem = ufo_buffer_get_device_array (inputs[i], cmd_queue);
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 0, sizeof (cl_mem), ¤t_in_mem));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 1, sizeof (gfloat), &distance));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 2, sizeof (guint), &priv->distance->n_values));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 3, sizeof (gfloat), &lambda));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 4, sizeof (gfloat), &priv->pixel_size));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 5, sizeof (gfloat), &priv->regularization_rate));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->ctf_multi_apply_dist_kernel, 6, sizeof (cl_mem), &in_sum_mem));
ufo_profiler_call (profiler, cmd_queue, priv->ctf_multi_apply_dist_kernel, requisition->n_dims, global_work_size, NULL);
}
in_mem = in_sum_mem;
} else {
in_mem = ufo_buffer_get_device_array (inputs[0], cmd_queue);
}
out_mem = ufo_buffer_get_device_array (output, cmd_queue);
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->mult_by_value_kernel, 0, sizeof (cl_mem), &in_mem));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->mult_by_value_kernel, 1, sizeof (cl_mem), &filter_mem));
UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->mult_by_value_kernel, 2, sizeof (cl_mem), &out_mem));
ufo_profiler_call (profiler, cmd_queue, priv->mult_by_value_kernel, requisition->n_dims, requisition->dims, NULL);
if (priv->method == METHOD_CTF_MULTI) {
UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (in_sum_mem));
}
}
return TRUE;
}
static void
ufo_retrieve_phase_task_get_property (GObject *object,
guint property_id,
GValue *value,
GParamSpec *pspec)
{
UfoRetrievePhaseTaskPrivate *priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE (object);
switch (property_id) {
case PROP_METHOD:
g_value_set_enum (value, priv->method);
break;
case PROP_ENERGY:
g_value_set_float (value, priv->energy);
break;
case PROP_DISTANCE:
g_value_set_boxed (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;
case PROP_REGULARIZATION_RATE:
g_value_set_float (value, priv->regularization_rate);
break;
case PROP_BINARY_FILTER_THRESHOLDING:
g_value_set_float (value, priv->binary_filter);
break;
case PROP_FREQUENCY_CUTOFF:
g_value_set_float (value, priv->frequency_cutoff);
break;
case PROP_OUTPUT_FILTER:
g_value_set_boolean (value, priv->output_filter);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
break;
}
}
static void
ufo_retrieve_phase_task_set_property (GObject *object,
guint property_id,
const GValue *value,
GParamSpec *pspec)
{
UfoRetrievePhaseTaskPrivate *priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE (object);
GValueArray *array;
switch (property_id) {
case PROP_METHOD:
priv->method = g_value_get_enum (value);
break;
case PROP_ENERGY:
priv->energy = g_value_get_float (value);
break;
case PROP_DISTANCE:
array = (GValueArray *) g_value_get_boxed (value);
if (array) {
g_value_array_free (priv->distance);
priv->distance = g_value_array_copy (array);
}
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;
case PROP_REGULARIZATION_RATE:
priv->regularization_rate = g_value_get_float (value);
break;
case PROP_BINARY_FILTER_THRESHOLDING:
priv->binary_filter = g_value_get_float (value);
break;
case PROP_FREQUENCY_CUTOFF:
priv->frequency_cutoff = g_value_get_float (value);
break;
case PROP_OUTPUT_FILTER:
priv->output_filter = g_value_get_boolean (value);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
break;
}
}
static void
ufo_retrieve_phase_task_finalize (GObject *object)
{
UfoRetrievePhaseTaskPrivate *priv;
priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE (object);
if (priv->kernels) {
for (int i = 0; i < N_METHODS; i++) {
if (priv->kernels[i]) {
UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->kernels[i]));
priv->kernels[i] = NULL;
}
}
}
g_free(priv->kernels);
if (priv->mult_by_value_kernel) {
UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->mult_by_value_kernel));
priv->mult_by_value_kernel = NULL;
}
if (priv->ctf_multi_apply_dist_kernel) {
UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->ctf_multi_apply_dist_kernel));
priv->ctf_multi_apply_dist_kernel = NULL;
}
if (priv->context) {
UFO_RESOURCES_CHECK_CLERR (clReleaseContext (priv->context));
priv->context = NULL;
}
if (priv->filter_buffer) {
g_object_unref(priv->filter_buffer);
}
g_value_array_free (priv->distance);
G_OBJECT_CLASS (ufo_retrieve_phase_task_parent_class)->finalize (object);
}
static void
ufo_task_interface_init (UfoTaskIface *iface)
{
iface->setup = ufo_retrieve_phase_task_setup;
iface->get_requisition = ufo_retrieve_phase_task_get_requisition;
iface->get_num_inputs = ufo_filter_task_get_num_inputs;
iface->get_num_dimensions = ufo_filter_task_get_num_dimensions;
iface->get_mode = ufo_filter_task_get_mode;
iface->process = ufo_retrieve_phase_task_process;
}
static void
ufo_retrieve_phase_task_class_init (UfoRetrievePhaseTaskClass *klass)
{
GObjectClass *gobject_class = G_OBJECT_CLASS (klass);
gobject_class->set_property = ufo_retrieve_phase_task_set_property;
gobject_class->get_property = ufo_retrieve_phase_task_get_property;
gobject_class->finalize = ufo_retrieve_phase_task_finalize;
GParamSpec *double_region_vals = g_param_spec_double ("double-region-values",
"Double Region values",
"Elements in double regions",
-INFINITY,
INFINITY,
0.0,
G_PARAM_READWRITE);
properties[PROP_METHOD] =
g_param_spec_enum ("method",
"Method name",
"Method name",
g_enum_register_static ("ufo_retrieve_phase_method", method_values),
METHOD_TIE,
G_PARAM_READWRITE);
properties[PROP_ENERGY] =
g_param_spec_float ("energy",
"Energy value",
"Energy value.",
0, G_MAXFLOAT, 20,
G_PARAM_READWRITE);
properties[PROP_DISTANCE] =
g_param_spec_value_array ("distance",
"Distance value",
"Distance value.",
double_region_vals,
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] =
g_param_spec_float ("pixel-size",
"Pixel size",
"Pixel size.",
0, G_MAXFLOAT, 0.75e-6,
G_PARAM_READWRITE);
properties[PROP_REGULARIZATION_RATE] =
g_param_spec_float ("regularization-rate",
"Regularization rate value",
"Regularization rate value.",
0, G_MAXFLOAT, 2.5,
G_PARAM_READWRITE);
properties[PROP_BINARY_FILTER_THRESHOLDING] =
g_param_spec_float ("thresholding-rate",
"Binary thresholding rate value",
"Binary thresholding rate value.",
0, G_MAXFLOAT, 0.1,
G_PARAM_READWRITE);
properties[PROP_FREQUENCY_CUTOFF] =
g_param_spec_float ("frequency-cutoff",
"Cut-off frequency in radians",
"Cut-off frequency in radians",
0, G_MAXFLOAT, G_MAXFLOAT,
G_PARAM_READWRITE);
properties[PROP_OUTPUT_FILTER] =
g_param_spec_boolean ("output-filter",
"Output the frequency filter, not the result of the filtering",
"Output the frequency filter, not the result of the filtering",
FALSE,
G_PARAM_READWRITE);
for (guint i = PROP_0 + 1; i < N_PROPERTIES; i++)
g_object_class_install_property (gobject_class, i, properties[i]);
g_type_class_add_private (gobject_class, sizeof(UfoRetrievePhaseTaskPrivate));
}
static void
ufo_retrieve_phase_task_init(UfoRetrievePhaseTask *self)
{
UfoRetrievePhaseTaskPrivate *priv;
self->priv = priv = UFO_RETRIEVE_PHASE_TASK_GET_PRIVATE(self);
priv->method = METHOD_TIE;
priv->energy = 20.0f;
priv->distance = g_value_array_new (1);
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;
priv->frequency_cutoff = G_MAXFLOAT;
priv->kernels = (cl_kernel *) g_malloc0(N_METHODS * sizeof(cl_kernel));
priv->filter_buffer = NULL;
priv->output_filter = FALSE;
}