summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
authorTomas Farago <sensej007@email.cz>2018-05-30 09:25:58 +0200
committerTomas Farago <sensej007@email.cz>2018-06-29 09:36:11 +0200
commit1da4a2349032596cea6a93dee6f8500aab852c6f (patch)
tree6979c2a8d740137aede4ad69e59657a4b988a7bf /src
parentde2007ebaa0a4146b12bfb7675a32b999e8e36b1 (diff)
downloadufo-filters-1da4a2349032596cea6a93dee6f8500aab852c6f.tar.gz
ufo-filters-1da4a2349032596cea6a93dee6f8500aab852c6f.tar.bz2
ufo-filters-1da4a2349032596cea6a93dee6f8500aab852c6f.tar.xz
ufo-filters-1da4a2349032596cea6a93dee6f8500aab852c6f.zip
Add cone beam reconstructor
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt11
-rw-r--r--src/common/ufo-conebeam.c38
-rw-r--r--src/common/ufo-conebeam.h55
-rw-r--r--src/common/ufo-ctgeometry.c119
-rw-r--r--src/common/ufo-ctgeometry.h57
-rw-r--r--src/common/ufo-scarray.c148
-rw-r--r--src/common/ufo-scarray.h49
-rw-r--r--src/kernels/CMakeLists.txt8
-rw-r--r--src/kernels/conebeam.cl36
-rw-r--r--src/kernels/meson.build7
-rw-r--r--src/kernels/templates/general_bp_body.in33
-rw-r--r--src/kernels/templates/general_bp_definitions.in3
-rw-r--r--src/kernels/templates/general_bp_header_scalar.in22
-rw-r--r--src/kernels/templates/general_bp_header_vector.in22
-rw-r--r--src/meson.build27
-rw-r--r--src/ufo-cone-beam-projection-weight-task.c320
-rw-r--r--src/ufo-cone-beam-projection-weight-task.h53
-rw-r--r--src/ufo-general-backproject-task.c2310
-rw-r--r--src/ufo-general-backproject-task.h53
19 files changed, 3370 insertions, 1 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index dbcd45b..21f323e 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,6 +15,7 @@ set(ufofilter_SRCS
ufo-cut-sinogram-task.c
ufo-center-of-rotation-task.c
ufo-concatenate-result-task.c
+ ufo-cone-beam-projection-weight-task.c
ufo-contrast-task.c
ufo-correlate-stacks-task.c
ufo-dfi-sinc-task.c
@@ -37,6 +38,7 @@ set(ufofilter_SRCS
ufo-forwardproject-task.c
ufo-get-dup-circ-task.c
ufo-gradient-task.c
+ ufo-general-backproject-task.c
ufo-ifft-task.c
ufo-interpolate-task.c
ufo-interpolate-stream-task.c
@@ -114,6 +116,15 @@ set(retrieve_phase_aux_SRCS
set(lamino_backproject_aux_SRCS
lamino-roi.c)
+set(cone_beam_projection_weight_aux_SRCS
+ common/ufo-scarray.c)
+
+set(general_backproject_aux_SRCS
+ common/ufo-math.c
+ common/ufo-conebeam.c
+ common/ufo-scarray.c
+ common/ufo-ctgeometry.c)
+
file(GLOB ufofilter_KERNELS "kernels/*.cl")
#}}}
#{{{ Variables
diff --git a/src/common/ufo-conebeam.c b/src/common/ufo-conebeam.c
new file mode 100644
index 0000000..5f531ae
--- /dev/null
+++ b/src/common/ufo-conebeam.c
@@ -0,0 +1,38 @@
+/*
+ * Copyright (C) 2015-2017 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 "ufo-conebeam.h"
+
+static UfoUniRecoNodeProps node_props[] = {
+ {"GENERIC", 8, 0},
+ {"GEFORCE_GTX_TITAN", 24, 32},
+ {"GEFORCE_GTX_1080_TI", 24, 32},
+};
+
+GHashTable *
+ufo_get_node_props_table (void)
+{
+ guint i;
+ GHashTable *table = g_hash_table_new (g_str_hash, g_str_equal);
+
+ for (i = 0; i < sizeof (node_props) / sizeof (UfoUniRecoNodeProps); i++) {
+ g_hash_table_insert (table, node_props[i].name, &node_props[i]);
+ }
+
+ return table;
+}
diff --git a/src/common/ufo-conebeam.h b/src/common/ufo-conebeam.h
new file mode 100644
index 0000000..58db453
--- /dev/null
+++ b/src/common/ufo-conebeam.h
@@ -0,0 +1,55 @@
+/*
+ * Copyright (C) 2015-2017 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_CONEBEAM_H
+#define UFO_CONEBEAM_H
+
+#include <glib.h>
+#include <glib-object.h>
+
+typedef enum {
+ UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_X,
+ UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Y,
+ UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Z,
+ UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_X,
+ UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Y,
+ UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Z,
+ UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_X,
+ UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Y,
+ UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Z,
+ UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_X,
+ UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Y,
+ UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Z,
+ UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_X,
+ UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Y,
+ UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Z,
+ UFO_UNI_RECO_PARAMETER_CENTER_POSITION_X,
+ UFO_UNI_RECO_PARAMETER_CENTER_POSITION_Z,
+ UFO_UNI_RECO_PARAMETER_Z
+} UfoUniRecoParameter;
+
+typedef struct {
+ gchar *name;
+ guint burst;
+ guint max_regcount;
+} UfoUniRecoNodeProps;
+
+GHashTable *ufo_get_node_props_table (void);
+
+#endif
diff --git a/src/common/ufo-ctgeometry.c b/src/common/ufo-ctgeometry.c
new file mode 100644
index 0000000..589ed82
--- /dev/null
+++ b/src/common/ufo-ctgeometry.c
@@ -0,0 +1,119 @@
+/*
+ * Copyright (C) 2017 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 <math.h>
+#include "ufo-math.h"
+#include "ufo-conebeam.h"
+#include "ufo-ctgeometry.h"
+#include "ufo-scarray.h"
+
+UfoScpoint *
+ufo_scpoint_new (UfoScarray *x, UfoScarray *y, UfoScarray *z)
+{
+ UfoScpoint *point = g_new0 (UfoScpoint, 1);
+
+ point->x = ufo_scarray_copy (x);
+ point->y = ufo_scarray_copy (y);
+ point->z = ufo_scarray_copy (z);
+
+ return point;
+}
+
+UfoScpoint *
+ufo_scpoint_copy (const UfoScpoint *point)
+{
+ return ufo_scpoint_new (point->x, point->y, point->z);
+}
+
+void
+ufo_scpoint_free (UfoScpoint *point)
+{
+ ufo_scarray_free (point->x);
+ ufo_scarray_free (point->y);
+ ufo_scarray_free (point->z);
+ g_free (point);
+}
+
+gboolean
+ufo_scpoint_are_almost_zero (UfoScpoint *point)
+{
+ return ufo_scarray_is_almost_zero (point->x) &&
+ ufo_scarray_is_almost_zero (point->y) &&
+ ufo_scarray_is_almost_zero (point->z);
+}
+
+UfoScvector *
+ufo_scvector_new (UfoScpoint *position, UfoScpoint *angle)
+{
+ UfoScvector *vector = g_new0 (UfoScvector, 1);
+ vector->position = ufo_scpoint_copy (position);
+ vector->angle = ufo_scpoint_copy (angle);
+
+ return vector;
+}
+
+void
+ufo_scvector_free (UfoScvector *vector)
+{
+ ufo_scpoint_free (vector->position);
+ ufo_scpoint_free (vector->angle);
+ g_free (vector);
+}
+
+/**
+ * Create a new computed tomography geometry with parallel beam geometry.
+ */
+UfoCTGeometry *
+ufo_ctgeometry_new (void)
+{
+ UfoCTGeometry *geometry = g_new0 (UfoCTGeometry, 1);
+ UfoScpoint *position, *angle;
+ UfoScarray *one, *one_inf;
+ GValue double_value = G_VALUE_INIT;
+
+ g_value_init (&double_value, G_TYPE_DOUBLE);
+ g_value_set_double (&double_value, -INFINITY);
+ one = ufo_scarray_new (1, G_TYPE_DOUBLE, NULL);
+ one_inf = ufo_scarray_new (1, G_TYPE_DOUBLE, &double_value);
+ position = ufo_scpoint_new (one, one, one);
+ angle = ufo_scpoint_new (one, one, one);
+
+ geometry->source_position = ufo_scpoint_new (one, one_inf, one);
+ geometry->volume_angle = ufo_scpoint_new (one, one, one);
+ geometry->axis = ufo_scvector_new (position, angle);
+ geometry->detector = ufo_scvector_new (position, angle);
+
+ ufo_scarray_free (one);
+ ufo_scarray_free (one_inf);
+ ufo_scpoint_free (position);
+ ufo_scpoint_free (angle);
+ g_value_unset (&double_value);
+
+ return geometry;
+}
+
+void
+ufo_ctgeometry_free (UfoCTGeometry *geometry)
+{
+ ufo_scpoint_free (geometry->source_position);
+ ufo_scpoint_free (geometry->volume_angle);
+ ufo_scvector_free (geometry->axis);
+ ufo_scvector_free (geometry->detector);
+ g_free (geometry);
+}
diff --git a/src/common/ufo-ctgeometry.h b/src/common/ufo-ctgeometry.h
new file mode 100644
index 0000000..60bdbfe
--- /dev/null
+++ b/src/common/ufo-ctgeometry.h
@@ -0,0 +1,57 @@
+/*
+ * Copyright (C) 2017 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_CTGEOMETRY_H
+#define UFO_CTGEOMETRY_H
+
+#include <ufo/ufo.h>
+#include "ufo-scarray.h"
+#include "ufo-conebeam.h"
+
+typedef struct {
+ UfoScarray *x;
+ UfoScarray *y;
+ UfoScarray *z;
+} UfoScpoint;
+
+typedef struct {
+ UfoScpoint *position;
+ UfoScpoint *angle;
+} UfoScvector;
+
+typedef struct {
+ UfoScpoint *source_position;
+ UfoScpoint *volume_angle;
+ UfoScvector *axis;
+ UfoScvector *detector;
+} UfoCTGeometry;
+
+UfoScpoint *ufo_scpoint_new (UfoScarray *x,
+ UfoScarray *y,
+ UfoScarray *z);
+UfoScpoint *ufo_scpoint_copy (const UfoScpoint *point);
+void ufo_scpoint_free (UfoScpoint *point);
+gboolean ufo_scpoint_are_almost_zero (UfoScpoint *point);
+UfoScvector *ufo_scvector_new (UfoScpoint *position,
+ UfoScpoint *angle);
+void ufo_scvector_free (UfoScvector *vector);
+UfoCTGeometry *ufo_ctgeometry_new (void);
+void ufo_ctgeometry_free (UfoCTGeometry *geometry);
+
+#endif
diff --git a/src/common/ufo-scarray.c b/src/common/ufo-scarray.c
new file mode 100644
index 0000000..db338e4
--- /dev/null
+++ b/src/common/ufo-scarray.c
@@ -0,0 +1,148 @@
+/*
+ * Copyright (C) 2017 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 "ufo-math.h"
+#include "ufo-scarray.h"
+
+struct _UfoScarray {
+ GValueArray *array;
+};
+
+UfoScarray *
+ufo_scarray_new (guint num_elements, GType type, GValue *init_value)
+{
+ UfoScarray *scarray;
+ GValue zero = G_VALUE_INIT;
+ GValue value = G_VALUE_INIT;
+
+ if (init_value == NULL) {
+ g_value_init (&zero, G_TYPE_INT);
+ g_value_init (&value, type);
+ g_value_set_int (&zero, 0);
+ g_value_transform (&zero, &value);
+ } else {
+ value = *init_value;
+ }
+
+ scarray = g_new0 (UfoScarray, 1);
+ scarray->array = g_value_array_new (num_elements);
+
+ for (guint i = 0; i < num_elements; i++)
+ g_value_array_append (scarray->array, &value);
+
+ if (init_value == NULL) {
+ g_value_unset (&zero);
+ g_value_unset (&value);
+ }
+
+ return scarray;
+}
+
+UfoScarray *
+ufo_scarray_copy (const UfoScarray *array)
+{
+ UfoScarray *cpy = g_new0 (UfoScarray, 1);
+ cpy->array = g_value_array_copy (array->array);
+
+ return cpy;
+}
+
+void
+ufo_scarray_free (UfoScarray *scarray)
+{
+ g_value_array_free (scarray->array);
+ g_free (scarray);
+}
+
+void
+ufo_scarray_set_value (UfoScarray *scarray,
+ GValue *value)
+{
+ g_value_set_boxed (value, scarray->array);
+}
+
+void
+ufo_scarray_get_value (UfoScarray *scarray,
+ const GValue *value)
+{
+ g_value_array_free (scarray->array);
+ scarray->array = g_value_array_copy (g_value_get_boxed (value));
+}
+
+void
+ufo_scarray_insert (UfoScarray *scarray,
+ guint index,
+ const GValue *value)
+{
+ g_value_array_insert (scarray->array, index, value);
+}
+
+gint
+ufo_scarray_get_int (UfoScarray *scarray,
+ guint index)
+{
+ if (scarray->array->n_values == 1)
+ return g_value_get_float (g_value_array_get_nth (scarray->array, 0));
+
+ g_assert (scarray->array->n_values > index);
+ return g_value_get_int (g_value_array_get_nth (scarray->array, index));
+}
+
+gfloat
+ufo_scarray_get_float (UfoScarray *scarray,
+ guint index)
+{
+ if (scarray->array->n_values == 1)
+ return g_value_get_float (g_value_array_get_nth (scarray->array, 0));
+
+ g_assert (scarray->array->n_values > index);
+ return g_value_get_float (g_value_array_get_nth (scarray->array, index));
+}
+
+gdouble
+ufo_scarray_get_double (UfoScarray *scarray,
+ guint index)
+{
+ if (scarray->array->n_values == 1)
+ return g_value_get_double (g_value_array_get_nth (scarray->array, 0));
+
+ g_assert (scarray->array->n_values > index);
+ return g_value_get_double (g_value_array_get_nth (scarray->array, index));
+}
+
+gboolean
+ufo_scarray_has_n_values (UfoScarray *scarray,
+ guint num)
+{
+ return scarray->array->n_values == num;
+}
+
+gboolean
+ufo_scarray_is_almost_zero (UfoScarray *scarray)
+{
+ guint i;
+
+ for (i = 0; i < scarray->array->n_values; i++) {
+ if (!UFO_MATH_ARE_ALMOST_EQUAL (ufo_scarray_get_double (scarray, i), 0)) {
+ return FALSE;
+ }
+ }
+
+ return TRUE;
+}
diff --git a/src/common/ufo-scarray.h b/src/common/ufo-scarray.h
new file mode 100644
index 0000000..292a136
--- /dev/null
+++ b/src/common/ufo-scarray.h
@@ -0,0 +1,49 @@
+/*
+ * Copyright (C) 2017 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_SCARRAY_H
+#define UFO_SCARRAY_H
+
+#include <ufo/ufo.h>
+
+typedef struct _UfoScarray UfoScarray;
+
+UfoScarray *ufo_scarray_new (guint num_elements,
+ GType type,
+ GValue *init_value);
+UfoScarray *ufo_scarray_copy (const UfoScarray *scarray);
+void ufo_scarray_free (UfoScarray *scarray);
+void ufo_scarray_set_value (UfoScarray *scarray,
+ GValue *value);
+void ufo_scarray_get_value (UfoScarray *scarray,
+ const GValue *value);
+void ufo_scarray_insert (UfoScarray *scarray,
+ guint index,
+ const GValue *value);
+gint ufo_scarray_get_int (UfoScarray *scarray,
+ guint index);
+gfloat ufo_scarray_get_float (UfoScarray *scarray,
+ guint index);
+gdouble ufo_scarray_get_double (UfoScarray *scarray,
+ guint index);
+gboolean ufo_scarray_has_n_values (UfoScarray *scarray,
+ guint num);
+gboolean ufo_scarray_is_almost_zero (UfoScarray *scarray);
+
+#endif
diff --git a/src/kernels/CMakeLists.txt b/src/kernels/CMakeLists.txt
index ec433ef..dcb37ca 100644
--- a/src/kernels/CMakeLists.txt
+++ b/src/kernels/CMakeLists.txt
@@ -7,6 +7,10 @@ set(Z_TEMPLATE ${CMAKE_CURRENT_SOURCE_DIR}/templates/z_template.in)
set(LAMINO_TEMPLATE ${CMAKE_CURRENT_SOURCE_DIR}/templates/lamino_template.in)
set(CENTER_TEMPLATE ${CMAKE_CURRENT_SOURCE_DIR}/templates/center_template.in)
set(ROLL_TEMPLATE ${CMAKE_CURRENT_SOURCE_DIR}/templates/roll_template.in)
+set(GENERAL_BP_DEFINITIONS templates/general_bp_definitions.in)
+set(GENERAL_BP_HEADER_SCALAR templates/general_bp_header_scalar.in)
+set(GENERAL_BP_HEADER_VECTOR templates/general_bp_header_vector.in)
+set(GENERAL_BP_BODY templates/general_bp_body.in)
set(GENERATOR ${CMAKE_CURRENT_SOURCE_DIR}/tools/make_burst_kernels.py)
set(Z_KERNEL z_kernel.cl)
set(LAMINO_KERNEL lamino_kernel.cl)
@@ -32,6 +36,10 @@ install(
${CMAKE_CURRENT_BINARY_DIR}/${CENTER_KERNEL}
${CMAKE_CURRENT_BINARY_DIR}/${LAMINO_KERNEL}
${CMAKE_CURRENT_BINARY_DIR}/${ROLL_KERNEL}
+ ${CMAKE_CURRENT_SOURCE_DIR}/${GENERAL_BP_DEFINITIONS}
+ ${CMAKE_CURRENT_SOURCE_DIR}/${GENERAL_BP_HEADER_SCALAR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/${GENERAL_BP_HEADER_VECTOR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/${GENERAL_BP_BODY}
DESTINATION ${UFO_KERNELDIR}
)
diff --git a/src/kernels/conebeam.cl b/src/kernels/conebeam.cl
new file mode 100644
index 0000000..c508628
--- /dev/null
+++ b/src/kernels/conebeam.cl
@@ -0,0 +1,36 @@
+/*
+ * Copyright (C) 2011-2017 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/>.
+ */
+
+kernel void weight_projection (global float *projection,
+ global float *result,
+ const float2 center,
+ const float source_distance,
+ const float overall_distance,
+ const float magnification_recip,
+ const float cos_theta)
+{
+ int idx = get_global_id (0);
+ int idy = get_global_id (1);
+ int index = idy * get_global_size (0) + idx;
+ float x = (((float) idx) - center.x + 0.5f) * magnification_recip;
+ float y = (((float) idy) - center.y + 0.5f) * magnification_recip;
+
+ /* Based on eq. 28 from "Direct cone beam SPECT reconstruction with camera tilt" */
+ result[index] = projection[index] * source_distance * cos_theta / sqrt (overall_distance * overall_distance + x * x + y * y);
+}
diff --git a/src/kernels/meson.build b/src/kernels/meson.build
index 542b5bb..864a94c 100644
--- a/src/kernels/meson.build
+++ b/src/kernels/meson.build
@@ -5,6 +5,7 @@ kernel_files = [
'bin.cl',
'clip.cl',
'complex.cl',
+ 'conebeam.cl',
'correlate.cl',
'cut.cl',
'cut-sinogram.cl',
@@ -38,7 +39,11 @@ kernel_files = [
'split.cl',
'swap-quadrants.cl',
'transpose.cl',
- 'zeropad.cl'
+ 'zeropad.cl',
+ 'templates/general_bp_definitions.in',
+ 'templates/general_bp_body.in',
+ 'templates/general_bp_header_scalar.in',
+ 'templates/general_bp_header_vector.in',
]
install_data(kernel_files,
diff --git a/src/kernels/templates/general_bp_body.in b/src/kernels/templates/general_bp_body.in
new file mode 100644
index 0000000..5f9409b
--- /dev/null
+++ b/src/kernels/templates/general_bp_body.in
@@ -0,0 +1,33 @@
+{
+ int idx = get_global_id (0);
+ int idy = get_global_id (1);
+ int idz = get_global_id (2);
+ cfloat3 voxel_0, voxel, detector_normal;
+ cfloat project_tmp, coeff, detector_offset, tmp_transformation;
+ rtype result = 0.0;
+ %tmpl%
+
+ if (idx < real_size.x && idy < real_size.y && idz < real_size.z) {
+ voxel_0.x = mad((cfloat) idx, x_region.y, x_region.x);
+ voxel_0.y = mad((cfloat) idy, y_region.y, y_region.x);
+ voxel_0.z = slice_z_position;
+
+ // assign z-coordinate to the chosen parameter i.e. param = region[idz];
+ %tmpl%
+
+ // Start rotation angle independent transformations and temporary assignments
+ %tmpl%
+ // End rotation angle independent transformations and temporary assignments
+
+ // Start rotation angle dependent transformation, pixel fetch and slice weighing
+ %tmpl%
+ // End rotation angle dependent transformation, pixel fetch and slice weighing
+
+ if (iteration) {
+ volume[idz * real_size.x * real_size.y + idy * real_size.x + idx] += %tmpl%;
+ } else {
+ volume[idz * real_size.x * real_size.y + idy * real_size.x + idx] = %tmpl%;
+ }
+ }
+}
+
diff --git a/src/kernels/templates/general_bp_definitions.in b/src/kernels/templates/general_bp_definitions.in
new file mode 100644
index 0000000..1006974
--- /dev/null
+++ b/src/kernels/templates/general_bp_definitions.in
@@ -0,0 +1,3 @@
+#define rotate_x(trig, point) ((cfloat3)(((point).x), mad ((trig).y, (point).y, - ((trig).x * (point).z)), mad ((trig).x, (point).y, ((trig).y * (point).z))))
+#define rotate_y(trig, point) ((cfloat3)(mad ((trig).y, (point).x, ((trig).x * (point).z)), ((point).y), mad (-(trig).x, (point).x, ((trig).y * (point).z))))
+#define rotate_z(trig, point) ((cfloat3)(mad ((trig).y, (point).x, - ((trig).x * (point).y)), mad ((trig).x, (point).x, ((trig).y * (point).y)), ((point).z)))
diff --git a/src/kernels/templates/general_bp_header_scalar.in b/src/kernels/templates/general_bp_header_scalar.in
new file mode 100644
index 0000000..15a5ecb
--- /dev/null
+++ b/src/kernels/templates/general_bp_header_scalar.in
@@ -0,0 +1,22 @@
+kernel void backproject (const sampler_t sampler,
+ const int3 real_size,
+ const cfloat2 x_region,
+ const cfloat2 y_region,
+ const cfloat slice_z_position,
+ cfloat2 axis_angle_x,
+ cfloat2 axis_angle_y,
+ cfloat2 volume_angle_x,
+ cfloat2 volume_angle_y,
+ cfloat2 volume_angle_z,
+ cfloat2 detector_angle_x,
+ cfloat2 detector_angle_y,
+ cfloat2 detector_angle_z,
+ cfloat3 center_position,
+ cfloat3 source_position,
+ cfloat3 detector_position,
+ const cfloat norm_factor,
+ const cfloat2 gray_limit,
+%tmpl%%tmpl%
+ const int iteration,
+ global stype *volume,
+ constant cfloat2 *region)
diff --git a/src/kernels/templates/general_bp_header_vector.in b/src/kernels/templates/general_bp_header_vector.in
new file mode 100644
index 0000000..7804a25
--- /dev/null
+++ b/src/kernels/templates/general_bp_header_vector.in
@@ -0,0 +1,22 @@
+kernel void backproject (const sampler_t sampler,
+ const int3 real_size,
+ const cfloat2 x_region,
+ const cfloat2 y_region,
+ const cfloat slice_z_position,
+ global cfloat2 *axis_angle_x,
+ global cfloat2 *axis_angle_y,
+ global cfloat2 *volume_angle_x,
+ global cfloat2 *volume_angle_y,
+ global cfloat2 *volume_angle_z,
+ global cfloat2 *detector_angle_x,
+ global cfloat2 *detector_angle_y,
+ global cfloat2 *detector_angle_z,
+ global cfloat3 *center_position,
+ global cfloat3 *source_position,
+ global cfloat3 *detector_position,
+ const cfloat norm_factor,
+ const cfloat2 gray_limit,
+%tmpl%%tmpl%
+ const int iteration,
+ global stype *volume,
+ constant cfloat2 *region)
diff --git a/src/meson.build b/src/meson.build
index 08388f0..235d82c 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -142,6 +142,33 @@ foreach plugin: plugins
)
endforeach
+# generalized backproject and conebeam
+
+shared_module('generalbackproject',
+ sources: [
+ 'ufo-general-backproject-task.c',
+ 'common/ufo-conebeam.c',
+ 'common/ufo-ctgeometry.c',
+ 'common/ufo-math.c',
+ 'common/ufo-scarray.c',
+ ],
+ dependencies: deps,
+ name_prefix: 'libufofilter',
+ install: true,
+ install_dir: plugin_install_dir,
+)
+
+shared_module('conebeamprojectionweight',
+ sources: [
+ 'ufo-cone-beam-projection-weight-task.c',
+ 'common/ufo-conebeam.c'
+ ],
+ dependencies: deps,
+ name_prefix: 'libufofilter',
+ install: true,
+ install_dir: plugin_install_dir,
+)
+
# fft plugins
fft_deps = deps
diff --git a/src/ufo-cone-beam-projection-weight-task.c b/src/ufo-cone-beam-projection-weight-task.c
new file mode 100644
index 0000000..87b7271
--- /dev/null
+++ b/src/ufo-cone-beam-projection-weight-task.c
@@ -0,0 +1,320 @@
+/*
+ * Copyright (C) 2011-2017 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 <math.h>
+#include <glib.h>
+#include <glib-object.h>
+#ifdef __APPLE__
+#include <OpenCL/cl.h>
+#else
+#include <CL/cl.h>
+#endif
+
+#include "ufo-cone-beam-projection-weight-task.h"
+#include "common/ufo-scarray.h"
+
+
+struct _UfoConeBeamProjectionWeightTaskPrivate {
+ /* Properties */
+ UfoScarray *center_position_x, *center_position_z, *source_distance, *detector_distance, *axis_angle_x;
+ /* Private */
+ guint count;
+ /* OpenCL */
+ cl_context context;
+ cl_kernel kernel;
+};
+
+static void ufo_task_interface_init (UfoTaskIface *iface);
+
+G_DEFINE_TYPE_WITH_CODE (UfoConeBeamProjectionWeightTask, ufo_cone_beam_projection_weight_task, UFO_TYPE_TASK_NODE,
+ G_IMPLEMENT_INTERFACE (UFO_TYPE_TASK, ufo_task_interface_init))
+
+#define UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_PRIVATE(obj) (G_TYPE_INSTANCE_GET_PRIVATE((obj), UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK, UfoConeBeamProjectionWeightTaskPrivate))
+
+enum {
+ PROP_0,
+ PROP_CENTER_POSITION_X,
+ PROP_CENTER_POSITION_Z,
+ PROP_SOURCE_DISTANCE,
+ PROP_DETECTOR_DISTANCE,
+ PROP_AXIS_ANGLE_X,
+ N_PROPERTIES
+};
+
+static GParamSpec *properties[N_PROPERTIES] = { NULL, };
+
+UfoNode *
+ufo_cone_beam_projection_weight_task_new (void)
+{
+ return UFO_NODE (g_object_new (UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK, NULL));
+}
+
+static void
+ufo_cone_beam_projection_weight_task_setup (UfoTask *task,
+ UfoResources *resources,
+ GError **error)
+{
+ UfoConeBeamProjectionWeightTaskPrivate *priv;
+
+ priv = UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_PRIVATE (task);
+ priv->context = ufo_resources_get_context (resources);
+ priv->kernel = ufo_resources_get_kernel (resources, "conebeam.cl", "weight_projection", NULL, error);
+
+ UFO_RESOURCES_CHECK_CLERR (clRetainContext (priv->context));
+ if (priv->kernel != NULL) {
+ UFO_RESOURCES_CHECK_CLERR (clRetainKernel (priv->kernel));
+ }
+}
+
+static void
+ufo_cone_beam_projection_weight_task_get_requisition (UfoTask *task,
+ UfoBuffer **inputs,
+ UfoRequisition *requisition,
+ GError **error)
+{
+ ufo_buffer_get_requisition (inputs[0], requisition);
+}
+
+static guint
+ufo_cone_beam_projection_weight_task_get_num_inputs (UfoTask *task)
+{
+ return 1;
+}
+
+static guint
+ufo_cone_beam_projection_weight_task_get_num_dimensions (UfoTask *task,
+ guint input)
+{
+ g_return_val_if_fail (input == 0, 0);
+
+ return 2;
+}
+
+static UfoTaskMode
+ufo_cone_beam_projection_weight_task_get_mode (UfoTask *task)
+{
+ return UFO_TASK_MODE_PROCESSOR | UFO_TASK_MODE_GPU;
+}
+
+static gboolean
+ufo_cone_beam_projection_weight_task_process (UfoTask *task,
+ UfoBuffer **inputs,
+ UfoBuffer *output,
+ UfoRequisition *requisition)
+{
+ UfoConeBeamProjectionWeightTaskPrivate *priv;
+ UfoGpuNode *node;
+ UfoProfiler *profiler;
+ cl_command_queue cmd_queue;
+ cl_mem in_mem;
+ cl_mem out_mem;
+ gfloat source_distance, detector_distance, overall_distance, magnification_recip, cos_angle;
+ gfloat center[2];
+
+ priv = UFO_CONE_BEAM_PROJECTION_WEIGHT_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));
+ in_mem = ufo_buffer_get_device_array (inputs[0], cmd_queue);
+ out_mem = ufo_buffer_get_device_array (output, cmd_queue);
+ cos_angle = cos (ufo_scarray_get_float (priv->axis_angle_x, priv->count));
+ center[0] = ufo_scarray_get_float (priv->center_position_x, priv->count);
+ center[1] = ufo_scarray_get_float (priv->center_position_z, priv->count);
+ source_distance = ufo_scarray_get_float (priv->source_distance, priv->count);
+ detector_distance = ufo_scarray_get_float (priv->detector_distance, priv->count);
+ overall_distance = source_distance + detector_distance;
+ magnification_recip = source_distance / overall_distance;
+ if (cos_angle > 0.9999999f) {
+ overall_distance = source_distance;
+ }
+
+ 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_float2), &center));
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 3, sizeof(cl_float), &source_distance));
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 4, sizeof(cl_float), &overall_distance));
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 5, sizeof(cl_float), &magnification_recip));
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (priv->kernel, 6, sizeof(cl_float), &cos_angle));
+ ufo_profiler_call (profiler, cmd_queue, priv->kernel, 2, requisition->dims, NULL);
+ priv->count++;
+
+ return TRUE;
+}
+
+
+static void
+ufo_cone_beam_projection_weight_task_set_property (GObject *object,
+ guint property_id,
+ const GValue *value,
+ GParamSpec *pspec)
+{
+ UfoConeBeamProjectionWeightTaskPrivate *priv = UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_PRIVATE (object);
+
+ switch (property_id) {
+ case PROP_CENTER_POSITION_X:
+ ufo_scarray_get_value (priv->center_position_x, value);
+ break;
+ case PROP_CENTER_POSITION_Z:
+ ufo_scarray_get_value (priv->center_position_z, value);
+ break;
+ case PROP_SOURCE_DISTANCE:
+ ufo_scarray_get_value (priv->source_distance, value);
+ break;
+ case PROP_DETECTOR_DISTANCE:
+ ufo_scarray_get_value (priv->detector_distance, value);
+ break;
+ case PROP_AXIS_ANGLE_X:
+ ufo_scarray_get_value (priv->axis_angle_x, value);
+ break;
+ default:
+ G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
+ break;
+ }
+}
+
+static void
+ufo_cone_beam_projection_weight_task_get_property (GObject *object,
+ guint property_id,
+ GValue *value,
+ GParamSpec *pspec)
+{
+ UfoConeBeamProjectionWeightTaskPrivate *priv = UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_PRIVATE (object);
+
+ switch (property_id) {
+ case PROP_CENTER_POSITION_X:
+ ufo_scarray_set_value (priv->center_position_x, value);
+ break;
+ case PROP_CENTER_POSITION_Z:
+ ufo_scarray_set_value (priv->center_position_z, value);
+ break;
+ case PROP_SOURCE_DISTANCE:
+ ufo_scarray_set_value (priv->source_distance, value);
+ break;
+ case PROP_DETECTOR_DISTANCE:
+ ufo_scarray_set_value (priv->detector_distance, value);
+ break;
+ case PROP_AXIS_ANGLE_X:
+ ufo_scarray_set_value (priv->axis_angle_x, value);
+ break;
+ default:
+ G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
+ break;
+ }
+}
+
+static void
+ufo_cone_beam_projection_weight_task_finalize (GObject *object)
+{
+ UfoConeBeamProjectionWeightTaskPrivate *priv;
+
+ priv = UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_PRIVATE (object);
+ if (priv->kernel) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->kernel));
+ priv->kernel = NULL;
+ }
+
+ if (priv->context) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseContext (priv->context));
+ priv->context = NULL;
+ }
+
+ G_OBJECT_CLASS (ufo_cone_beam_projection_weight_task_parent_class)->finalize (object);
+}
+
+static void
+ufo_task_interface_init (UfoTaskIface *iface)
+{
+ iface->setup = ufo_cone_beam_projection_weight_task_setup;
+ iface->get_num_inputs = ufo_cone_beam_projection_weight_task_get_num_inputs;
+ iface->get_num_dimensions = ufo_cone_beam_projection_weight_task_get_num_dimensions;
+ iface->get_mode = ufo_cone_beam_projection_weight_task_get_mode;
+ iface->get_requisition = ufo_cone_beam_projection_weight_task_get_requisition;
+ iface->process = ufo_cone_beam_projection_weight_task_process;
+}
+
+static void
+ufo_cone_beam_projection_weight_task_class_init (UfoConeBeamProjectionWeightTaskClass *klass)
+{
+ GObjectClass *oclass = G_OBJECT_CLASS (klass);
+
+ oclass->set_property = ufo_cone_beam_projection_weight_task_set_property;
+ oclass->get_property = ufo_cone_beam_projection_weight_task_get_property;
+ oclass->finalize = ufo_cone_beam_projection_weight_task_finalize;
+
+ GParamSpec *float_region_vals = g_param_spec_float ("float-region-values",
+ "Float Region values",
+ "Elements in float regions",
+ -INFINITY,
+ INFINITY,
+ 0.0f,
+ G_PARAM_READWRITE);
+
+ properties[PROP_CENTER_POSITION_X] =
+ g_param_spec_value_array ("center-position-x",
+ "Global x center (horizontal in a projection) of the volume with respect to projections",
+ "Global x center (horizontal in a projection) of the volume with respect to projections",
+ float_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_CENTER_POSITION_Z] =
+ g_param_spec_value_array ("center-position-z",
+ "Global z center (vertical in a projection) of the volume with respect to projections",
+ "Global z center (vertical in a projection) of the volume with respect to projections",
+ float_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_SOURCE_DISTANCE] =
+ g_param_spec_value_array ("source-distance",
+ "Distance from source to the volume center",
+ "Distance from source to the volume center",
+ float_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_DISTANCE] =
+ g_param_spec_value_array ("detector-distance",
+ "Distance from detector to the volume center",
+ "Distance from detector to the volume center",
+ float_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_AXIS_ANGLE_X] =
+ g_param_spec_value_array("axis-angle-x",
+ "Rotation axis rotation around the x-axis (laminographic angle [rad], 0 = tomography)",
+ "Rotation axis rotation around the x-axis (laminographic angle [rad], 0 = tomography)",
+ float_region_vals,
+ 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 (oclass, sizeof(UfoConeBeamProjectionWeightTaskPrivate));
+}
+
+static void
+ufo_cone_beam_projection_weight_task_init(UfoConeBeamProjectionWeightTask *self)
+{
+ self->priv = UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_PRIVATE(self);
+
+ self->priv->center_position_x = ufo_scarray_new (0, G_TYPE_FLOAT, NULL);
+ self->priv->center_position_z = ufo_scarray_new (0, G_TYPE_FLOAT, NULL);
+ self->priv->source_distance = ufo_scarray_new (0, G_TYPE_FLOAT, NULL);
+ self->priv->detector_distance = ufo_scarray_new (0, G_TYPE_FLOAT, NULL);
+ self->priv->axis_angle_x = ufo_scarray_new (1, G_TYPE_FLOAT, NULL);
+ self->priv->count = 0;
+}
diff --git a/src/ufo-cone-beam-projection-weight-task.h b/src/ufo-cone-beam-projection-weight-task.h
new file mode 100644
index 0000000..a4ec3bc
--- /dev/null
+++ b/src/ufo-cone-beam-projection-weight-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_CONE_BEAM_PROJECTION_WEIGHT_TASK_H
+#define __UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_H
+
+#include <ufo/ufo.h>
+
+G_BEGIN_DECLS
+
+#define UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK (ufo_cone_beam_projection_weight_task_get_type())
+#define UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK(obj) (G_TYPE_CHECK_INSTANCE_CAST((obj), UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK, UfoConeBeamProjectionWeightTask))
+#define UFO_IS_CONE_BEAM_PROJECTION_WEIGHT_TASK(obj) (G_TYPE_CHECK_INSTANCE_TYPE((obj), UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK))
+#define UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_CLASS(klass) (G_TYPE_CHECK_CLASS_CAST((klass), UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK, UfoConeBeamProjectionWeightTaskClass))
+#define UFO_IS_CONE_BEAM_PROJECTION_WEIGHT_TASK_CLASS(klass) (G_TYPE_CHECK_CLASS_TYPE((klass), UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK))
+#define UFO_CONE_BEAM_PROJECTION_WEIGHT_TASK_GET_CLASS(obj) (G_TYPE_INSTANCE_GET_CLASS((obj), UFO_TYPE_CONE_BEAM_PROJECTION_WEIGHT_TASK, UfoConeBeamProjectionWeightTaskClass))
+
+typedef struct _UfoConeBeamProjectionWeightTask UfoConeBeamProjectionWeightTask;
+typedef struct _UfoConeBeamProjectionWeightTaskClass UfoConeBeamProjectionWeightTaskClass;
+typedef struct _UfoConeBeamProjectionWeightTaskPrivate UfoConeBeamProjectionWeightTaskPrivate;
+
+struct _UfoConeBeamProjectionWeightTask {
+ UfoTaskNode parent_instance;
+
+ UfoConeBeamProjectionWeightTaskPrivate *priv;
+};
+
+struct _UfoConeBeamProjectionWeightTaskClass {
+ UfoTaskNodeClass parent_class;
+};
+
+UfoNode *ufo_cone_beam_projection_weight_task_new (void);
+GType ufo_cone_beam_projection_weight_task_get_type (void);
+
+G_END_DECLS
+
+#endif
diff --git a/src/ufo-general-backproject-task.c b/src/ufo-general-backproject-task.c
new file mode 100644
index 0000000..ac4547c
--- /dev/null
+++ b/src/ufo-general-backproject-task.c
@@ -0,0 +1,2310 @@
+/*
+ * 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/>.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <glib.h>
+#include <glib/gprintf.h>
+#ifdef __APPLE__
+#include <OpenCL/cl.h>
+#else
+#include <CL/cl.h>
+#endif
+
+#include <config.h>
+#include <common/ufo-math.h>
+#include "common/ufo-conebeam.h"
+#include "common/ufo-scarray.h"
+#include "common/ufo-ctgeometry.h"
+#include "common/ufo-addressing.h"
+#include "ufo-general-backproject-task.h"
+
+#define NUM_VECTOR_ARGUMENTS 11
+#define REAL_SIZE_ARG_INDEX 1
+#define STATIC_ARG_OFFSET 18
+#define G_LOG_LEVEL_DOMAIN "gbp"
+#define REGION_SIZE(region) ((ufo_scarray_get_int ((region), 2)) == 0) ? 0 : \
+ ((ufo_scarray_get_int ((region), 1) - ufo_scarray_get_int ((region), 0) - 1) /\
+ ufo_scarray_get_int ((region), 2) + 1)
+#define NEXT_DIVISOR(dividend, divisor) ((dividend) + (divisor) - (dividend) % (divisor))
+#define DEFINE_FILL_SINCOS(type) \
+static void \
+fill_sincos_##type (type *array, const gdouble angle) \
+{ \
+ array[0] = (type) sin (angle); \
+ array[1] = (type) cos (angle); \
+}
+
+#define DEFINE_CREATE_REGIONS(type) \
+static void \
+create_regions_##type (UfoGeneralBackprojectTaskPrivate *priv, \
+ const cl_command_queue cmd_queue, \
+ const gdouble start, \
+ const gdouble step) \
+{ \
+ guint i, j; \
+ gsize region_size; \
+ type *region_values; \
+ cl_int cl_error; \
+ gdouble value; \
+ gboolean is_angular = is_parameter_angular (priv->parameter); \
+ \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Start, step: %g %g", start, step); \
+ \
+ region_size = priv->num_slices_per_chunk * 2 * sizeof (type); \
+ region_values = (type *) g_malloc0 (region_size); \
+ \
+ for (i = 0; i < priv->num_chunks; i++) { \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Chunk %d region:", i); \
+ for (j = 0; j < priv->num_slices_per_chunk; j++) { \
+ value = start + i * priv->num_slices_per_chunk * step + j * step; \
+ if (is_angular) { \
+ region_values[2 * j] = (type) sin (value); \
+ region_values[2 * j + 1] = (type) cos (value); \
+ } else { \
+ region_values[2 * j] = (type) value; \
+ } \
+ } \
+ priv->cl_regions[i] = clCreateBuffer (priv->context, \
+ CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,\
+ region_size, \
+ region_values, \
+ &cl_error); \
+ UFO_RESOURCES_CHECK_CLERR (cl_error); \
+ } \
+ \
+ g_free (region_values); \
+}
+
+#define DEFINE_TRANSFER_ANGULAR_ARGUMENT(type) \
+static cl_mem \
+transfer_angular_argument_##type (UfoGeneralBackprojectTaskPrivate *priv, \
+ UfoScarray *source) \
+{ \
+ gsize size; \
+ guint i; \
+ cl_mem device_array; \
+ type *host_array; \
+ \
+ size = 2 * priv->num_projections * sizeof (type); \
+ if (!(host_array = (type *) g_malloc (size))) { \
+ g_warning ("Error allocating vectorized parameter host memory"); \
+ return NULL; \
+ } \
+ \
+ for (i = 0; i < priv->num_projections; i++) { \
+ fill_sincos_##type (host_array + 2 * i, ufo_scarray_get_double (source, i)); \
+ } \
+ \
+ device_array = transfer_host_to_device (priv->context, host_array, size); \
+ g_free (host_array); \
+ \
+ return device_array; \
+} \
+
+#define DEFINE_TRANSFER_POSITINAL_ARGUMENT(type) \
+static cl_mem \
+transfer_positional_argument_##type (UfoGeneralBackprojectTaskPrivate *priv, \
+ UfoScpoint *source) \
+{ \
+ gsize size; \
+ guint i; \
+ type *host_array; \
+ cl_mem device_array; \
+ \
+ size = 4 * priv->num_projections * sizeof (type); \
+ if (!(host_array = (type *) g_malloc (size))) { \
+ g_warning ("Error allocating vectorized parameter host memory"); \
+ return NULL; \
+ } \
+ \
+ for (i = 0; i < priv->num_projections; i++) { \
+ host_array[4 * i] = (type) ufo_scarray_get_double (source->x, i); \
+ host_array[4 * i + 1] = (type) ufo_scarray_get_double (source->y, i); \
+ host_array[4 * i + 2] = (type) ufo_scarray_get_double (source->z, i); \
+ } \
+ \
+ device_array = transfer_host_to_device (priv->context, host_array, size); \
+ g_free (host_array); \
+ \
+ return device_array; \
+}
+
+#define DEFINE_SET_ANGULAR_VECTOR_KERNEL_ARGUMENT(type) \
+static void \
+set_angular_vector_kernel_argument_##type (UfoGeneralBackprojectTaskPrivate *priv, \
+ cl_kernel kernel, \
+ UfoScpoint *point, \
+ gint mem_index, \
+ gint arg_index) \
+{ \
+ priv->vector_arguments[mem_index] = transfer_angular_argument_##type (priv, point->x); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[mem_index++]));\
+ priv->vector_arguments[mem_index] = transfer_angular_argument_##type (priv, point->y); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[mem_index++]));\
+ priv->vector_arguments[mem_index] = transfer_angular_argument_##type (priv, point->z); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index, sizeof (cl_mem), &priv->vector_arguments[mem_index])); \
+}
+
+#define DEFINE_SET_STATIC_VECTOR_ARGUMENTS(type) \
+static gint \
+set_static_vector_arguments_##type (UfoTask *task, \
+ cl_kernel kernel, \
+ gint arg_index) \
+{ \
+ UfoGeneralBackprojectTaskPrivate *priv; \
+ gint mem_index; \
+ \
+ priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (task); \
+ priv->vector_arguments = (cl_mem *) g_malloc (NUM_VECTOR_ARGUMENTS * sizeof (cl_mem)); \
+ \
+ /* Axis angle has only two vector components, the z one is the tomographic angle with offset set in process () */ \
+ priv->vector_arguments[0] = transfer_angular_argument_##type (priv, priv->geometry->axis->angle->x); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[0])); \
+ priv->vector_arguments[1] = transfer_angular_argument_##type (priv, priv->geometry->axis->angle->y); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[1])); \
+ \
+ set_angular_vector_kernel_argument_##type (priv, kernel, priv->geometry->volume_angle, 2, arg_index); \
+ set_angular_vector_kernel_argument_##type (priv, kernel, priv->geometry->detector->angle, 5, arg_index + 3); \
+ mem_index = 8; \
+ arg_index += 6; \
+ priv->vector_arguments[mem_index] = transfer_positional_argument_##type (priv, priv->geometry->axis->position); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[mem_index++]));\
+ priv->vector_arguments[mem_index] = transfer_positional_argument_##type (priv, priv->geometry->source_position); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[mem_index++]));\
+ priv->vector_arguments[mem_index] = transfer_positional_argument_##type (priv, priv->geometry->detector->position); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (cl_mem), &priv->vector_arguments[mem_index++]));\
+ \
+ return arg_index; \
+}
+
+#define DEFINE_SET_STATIC_SCALAR_ARGUMENTS(type) \
+static gint \
+set_static_scalar_arguments_##type (UfoTask *task, \
+ cl_kernel kernel, \
+ gint arg_index) \
+{ \
+ UfoGeneralBackprojectTaskPrivate *priv; \
+ type axis_angle_x[2], axis_angle_y[2], axis_angle_z[2], volume_angle_x[2], volume_angle_y[2], volume_angle_z[2], \
+ detector_angle_x[2], detector_angle_y[2], detector_angle_z[2], center_position[4], source_position[4], \
+ detector_position[4]; \
+ guint count; \
+ priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (task); \
+ g_object_get (task, "num_processed", &count, NULL); \
+ \
+ fill_sincos_##type (axis_angle_x, ufo_scarray_get_double (priv->geometry->axis->angle->x, count)); \
+ fill_sincos_##type (axis_angle_y, ufo_scarray_get_double (priv->geometry->axis->angle->y, count)); \
+ fill_sincos_##type (axis_angle_z, ufo_scarray_get_double (priv->geometry->axis->angle->z, count)); \
+ fill_sincos_##type (volume_angle_x, ufo_scarray_get_double (priv->geometry->volume_angle->x, count)); \
+ fill_sincos_##type (volume_angle_y, ufo_scarray_get_double (priv->geometry->volume_angle->y, count)); \
+ fill_sincos_##type (volume_angle_z, ufo_scarray_get_double (priv->geometry->volume_angle->z, count)); \
+ fill_sincos_##type (detector_angle_x, ufo_scarray_get_double (priv->geometry->detector->angle->x, count)); \
+ fill_sincos_##type (detector_angle_y, ufo_scarray_get_double (priv->geometry->detector->angle->y, count)); \
+ fill_sincos_##type (detector_angle_z, ufo_scarray_get_double (priv->geometry->detector->angle->z, count)); \
+ center_position[0] = (type) ufo_scarray_get_double (priv->geometry->axis->position->x, count); \
+ center_position[2] = (type) ufo_scarray_get_double (priv->geometry->axis->position->z, count); \
+ /* TODO: use only 2D center in the kernel */ \
+ center_position[1] = 0.0f; \
+ source_position[0] = (type) ufo_scarray_get_double (priv->geometry->source_position->x, count); \
+ source_position[1] = (type) ufo_scarray_get_double (priv->geometry->source_position->y, count); \
+ source_position[2] = (type) ufo_scarray_get_double (priv->geometry->source_position->z, count); \
+ detector_position[0] = (type) ufo_scarray_get_double (priv->geometry->detector->position->x, count); \
+ detector_position[1] = (type) ufo_scarray_get_double (priv->geometry->detector->position->y, count); \
+ detector_position[2] = (type) ufo_scarray_get_double (priv->geometry->detector->position->z, count); \
+ \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), axis_angle_x)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), axis_angle_y)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), volume_angle_x)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), volume_angle_y)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), volume_angle_z)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), detector_angle_x)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), detector_angle_y)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##2), detector_angle_z)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##3), center_position)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##3), source_position)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, arg_index++, sizeof (type##3), detector_position)); \
+ \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "axis: %g %g, %g %g, %g %g", \
+ axis_angle_x[0], axis_angle_x[1], axis_angle_y[0], axis_angle_y[1], axis_angle_z[0], axis_angle_z[1]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "volume: %g %g, %g %g, %g %g", \
+ volume_angle_x[0], volume_angle_x[1], volume_angle_y[0], volume_angle_y[1], volume_angle_z[0], \
+ volume_angle_z[1]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "detector_x: %g %g, %g %g, %g %g", \
+ detector_angle_x[0], detector_angle_x[1], detector_angle_y[0], detector_angle_y[1], detector_angle_z[0], \
+ detector_angle_z[1]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "center_position: %g %g %g", \
+ center_position[0], center_position[1], center_position[2]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "source_position: %g %g %g", \
+ source_position[0], source_position[1], source_position[2]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "detector_position: %g %g %g", \
+ detector_position[0], detector_position[1], detector_position[2]); \
+ return arg_index; \
+}
+
+#define DEFINE_COMPUTE_SLICE_REGION(type) \
+static void \
+compute_slice_region_##type (UfoGeneralBackprojectTaskPrivate *priv, \
+ gsize length, \
+ UfoScarray *region, \
+ type region_for_opencl[2]) \
+{ \
+ if (ufo_scarray_get_int (region, 2)) { \
+ region_for_opencl[0] = (type) ufo_scarray_get_int (region, 0); \
+ region_for_opencl[1] = (type) ufo_scarray_get_int (region, 2); \
+ } else { \
+ region_for_opencl[0] = (type) -(length / 2); \
+ region_for_opencl[1] = (type) 1; \
+ } \
+}
+
+#define DEFINE_SET_STATIC_ARGS(type) \
+static void \
+set_static_args_##type (UfoTask *task, \
+ UfoRequisition *requisition, \
+ const cl_kernel kernel) \
+{ \
+ UfoGeneralBackprojectTaskPrivate *priv; \
+ gdouble gray_delta_recip; \
+ type slice_z_position, region_x[2], region_y[2], gray_limit[2], norm_factor; \
+ /* 0 = sampler, 1 = real size set in process (), thus we start with 2 */ \
+ guint burst, j, i = 2; \
+ priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (task); \
+ gray_delta_recip = (gdouble) get_integer_maximum (st_values[priv->store_type].value_nick) / \
+ (priv->gray_map_max - priv->gray_map_min); \
+ norm_factor = fabs (priv->overall_angle) / priv->num_projections; \
+ burst = kernel == priv->kernel ? priv->burst : priv->num_projections % priv->burst; \
+ \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, 0, sizeof (cl_sampler), &priv->sampler)); \
+ \
+ compute_slice_region_##type (priv, requisition->dims[0], priv->region_x, region_x); \
+ compute_slice_region_##type (priv, requisition->dims[1], priv->region_y, region_y); \
+ slice_z_position = (type) priv->z; \
+ norm_factor = (type) norm_factor; \
+ gray_limit[0] = (type) priv->gray_map_min; \
+ gray_limit[1] = (type) gray_delta_recip; \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (type##2), region_x)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (type##2), region_y)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (type), &slice_z_position)); \
+ if (priv->vectorized) { \
+ i = set_static_vector_arguments_##type (task, kernel, i); \
+ } else { \
+ i = set_static_scalar_arguments_##type (task, kernel, i); \
+ } \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (type), &norm_factor)); \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (type##2), gray_limit)); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "region_x: %g %g", region_x[0], region_x[1]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "region_y: %g %g", region_y[0], region_y[1]); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "slice_z_position: %g", slice_z_position); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "norm_factor: %g", norm_factor); \
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "gray_limit: %g %g", gray_limit[0], gray_limit[1]); \
+ \
+ for (j = 0; j < burst; j++) { \
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, i++, sizeof (cl_mem), &priv->projections[j])); \
+ } \
+}
+
+/*{{{ Enumerations */
+typedef enum {
+ FT_HALF,
+ FT_FLOAT,
+ FT_DOUBLE
+} FloatType;
+
+typedef enum {
+ CT_FLOAT,
+ CT_DOUBLE
+} ComputeType;
+
+typedef enum {
+ ST_HALF,
+ ST_FLOAT,
+ ST_DOUBLE,
+ ST_UCHAR,
+ ST_USHORT,
+ ST_UINT
+} StoreType;
+
+static const GEnumValue parameter_values[] = {
+ { UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_X, "UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_X", "axis-angle-x" },
+ { UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Y, "UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Y", "axis-angle-y" },
+ { UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Z, "UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Z", "axis-angle-z" },
+ { UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_X, "UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_X", "volume-angle-x" },
+ { UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Y, "UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Y", "volume-angle-y" },
+ { UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Z, "UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Z", "volume-angle-z" },
+ { UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_X, "UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_X", "detector-angle-x" },
+ { UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Y, "UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Y", "detector-angle-y" },
+ { UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Z, "UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Z", "detector-angle-z" },
+ { UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_X, "UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_X", "detector-position-x" },
+ { UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Y, "UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Y", "detector-position-y" },
+ { UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Z, "UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Z", "detector-position-z" },
+ { UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_X, "UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_X", "source-position-x" },
+ { UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Y, "UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Y", "source-position-y" },
+ { UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Z, "UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Z", "source-position-z" },
+ { UFO_UNI_RECO_PARAMETER_CENTER_POSITION_X, "UFO_UNI_RECO_PARAMETER_CENTER_POSITION_X", "center-position-x" },
+ { UFO_UNI_RECO_PARAMETER_CENTER_POSITION_Z, "UFO_UNI_RECO_PARAMETER_CENTER_POSITION_Z", "center-position-z" },
+ { UFO_UNI_RECO_PARAMETER_Z, "UFO_UNI_RECO_PARAMETER_Z", "z" },
+ { 0, NULL, NULL}
+};
+
+static GEnumValue compute_type_values[] = {
+ {CT_FLOAT, "FT_FLOAT", "float"},
+ {CT_DOUBLE, "FT_DOUBLE", "double"},
+ { 0, NULL, NULL}
+};
+
+static GEnumValue ft_values[] = {
+ {FT_HALF, "FT_HALF", "half"},
+ {FT_FLOAT, "FT_FLOAT", "float"},
+ {FT_DOUBLE, "FT_DOUBLE", "double"},
+ { 0, NULL, NULL}
+};
+
+static GEnumValue st_values[] = {
+ {ST_HALF, "ST_HALF", "half"},
+ {ST_FLOAT, "ST_FLOAT", "float"},
+ {ST_DOUBLE, "ST_DOUBLE", "double"},
+ {ST_UCHAR, "ST_UCHAR", "uchar"},
+ {ST_USHORT, "ST_USHORT", "ushort"},
+ {ST_UINT, "ST_UINT", "uint"},
+ { 0, NULL, NULL}
+};
+/*}}}*/
+
+struct _UfoGeneralBackprojectTaskPrivate {
+ /* Properties */
+ guint burst;
+ gdouble z;
+ UfoScarray *region, *region_x, *region_y;
+ UfoCTGeometry *geometry;
+ ComputeType compute_type, result_type;
+ StoreType store_type;
+ UfoUniRecoParameter parameter;
+ gdouble gray_map_min, gray_map_max;
+ /* Private */
+ gboolean vectorized;
+ guint generated;
+ UfoResources *resources;
+ cl_mem *projections;
+ cl_mem *chunks;
+ cl_mem *cl_regions, *vector_arguments;
+ guint num_slices, num_slices_per_chunk, num_chunks;
+ guint num_projections;
+ gdouble overall_angle;
+ AddressingMode addressing_mode;
+ GHashTable *node_props_table;
+ /* OpenCL */
+ cl_context context;
+ cl_kernel kernel, rest_kernel;
+ cl_sampler sampler;
+};
+
+static void ufo_task_interface_init (UfoTaskIface *iface);
+
+G_DEFINE_TYPE_WITH_CODE (UfoGeneralBackprojectTask, ufo_general_backproject_task, UFO_TYPE_TASK_NODE,
+ G_IMPLEMENT_INTERFACE (UFO_TYPE_TASK,
+ ufo_task_interface_init))
+
+#define UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE(obj) (G_TYPE_INSTANCE_GET_PRIVATE((obj), UFO_TYPE_GENERAL_BACKPROJECT_TASK, UfoGeneralBackprojectTaskPrivate))
+
+enum {
+ PROP_0,
+ PROP_BURST,
+ PROP_PARAMETER,
+ PROP_Z,
+ PROP_REGION,
+ PROP_REGION_X,
+ PROP_REGION_Y,
+ PROP_CENTER_POSITION_X,
+ PROP_CENTER_POSITION_Z,
+ PROP_SOURCE_POSITION_X,
+ PROP_SOURCE_POSITION_Y,
+ PROP_SOURCE_POSITION_Z,
+ PROP_DETECTOR_POSITION_X,
+ PROP_DETECTOR_POSITION_Y,
+ PROP_DETECTOR_POSITION_Z,
+ PROP_DETECTOR_ANGLE_X,
+ PROP_DETECTOR_ANGLE_Y,
+ PROP_DETECTOR_ANGLE_Z,
+ PROP_AXIS_ANGLE_X,
+ PROP_AXIS_ANGLE_Y,
+ PROP_AXIS_ANGLE_Z,
+ PROP_VOLUME_ANGLE_X,
+ PROP_VOLUME_ANGLE_Y,
+ PROP_VOLUME_ANGLE_Z,
+ PROP_NUM_PROJECTIONS,
+ PROP_COMPUTE_TYPE,
+ PROP_RESULT_TYPE,
+ PROP_STORE_TYPE,
+ PROP_OVERALL_ANGLE,
+ PROP_ADDRESSING_MODE,
+ PROP_GRAY_MAP_MIN,
+ PROP_GRAY_MAP_MAX,
+ N_PROPERTIES
+};
+
+static GParamSpec *properties[N_PROPERTIES] = { NULL, };
+
+/*{{{ General helper functions*/
+DEFINE_FILL_SINCOS (cl_float)
+DEFINE_FILL_SINCOS (cl_double)
+
+static gsize
+get_type_size (StoreType type)
+{
+ gsize size;
+
+ switch (type) {
+ case ST_HALF:
+ size = sizeof (cl_half);
+ break;
+ case ST_FLOAT:
+ size = sizeof (cl_float);
+ break;
+ case ST_DOUBLE:
+ size = sizeof (cl_double);
+ break;
+ case ST_UCHAR:
+ size = sizeof (cl_uchar);
+ break;
+ case ST_USHORT:
+ size = sizeof (cl_ushort);
+ break;
+ case ST_UINT:
+ size = sizeof (cl_uint);
+ break;
+ default:
+ g_warning ("Uknown store type");
+ size = 0;
+ break;
+ }
+
+ return size;
+}
+
+static gulong
+get_integer_maximum (const gchar *type_name)
+{
+ if (!g_strcmp0 (type_name, "uchar"))
+ return 0xFF;
+
+ if (!g_strcmp0 (type_name, "ushort"))
+ return 0xFFFF;
+
+ if (!g_strcmp0 (type_name, "uint"))
+ return 0xFFFFFFFF;
+
+ return 0;
+}
+
+static gboolean
+is_axis_parameter (UfoUniRecoParameter parameter)
+{
+ return parameter == UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_X ||
+ parameter == UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Y ||
+ parameter == UFO_UNI_RECO_PARAMETER_AXIS_ROTATION_Z;
+}
+
+static gboolean
+is_volume_parameter (UfoUniRecoParameter parameter)
+{
+ return parameter == UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_X ||
+ parameter == UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Y ||
+ parameter == UFO_UNI_RECO_PARAMETER_VOLUME_ROTATION_Z;
+}
+
+static gboolean
+is_detector_rotation_parameter (UfoUniRecoParameter parameter)
+{
+ return parameter == UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_X ||
+ parameter == UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Y ||
+ parameter == UFO_UNI_RECO_PARAMETER_DETECTOR_ROTATION_Z;
+}
+
+static gboolean
+is_center_position_parameter (UfoUniRecoParameter parameter)
+{
+ return parameter == UFO_UNI_RECO_PARAMETER_CENTER_POSITION_X ||
+ parameter == UFO_UNI_RECO_PARAMETER_CENTER_POSITION_Z;
+}
+
+static gboolean
+is_source_position_parameter (UfoUniRecoParameter parameter)
+{
+ return parameter == UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_X ||
+ parameter == UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Y ||
+ parameter == UFO_UNI_RECO_PARAMETER_SOURCE_POSITION_Z;
+}
+
+static gboolean
+is_detector_position_parameter (UfoUniRecoParameter parameter)
+{
+ return parameter == UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_X ||
+ parameter == UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Y ||
+ parameter == UFO_UNI_RECO_PARAMETER_DETECTOR_POSITION_Z;
+}
+
+static gboolean
+is_parameter_positional (UfoUniRecoParameter parameter)
+{
+ return is_center_position_parameter (parameter) ||
+ is_source_position_parameter (parameter) ||
+ is_detector_position_parameter (parameter);
+}
+
+static gboolean
+is_parameter_angular (UfoUniRecoParameter parameter)
+{
+ return is_axis_parameter (parameter) ||
+ is_volume_parameter (parameter) ||
+ is_detector_position_parameter (parameter);
+}
+/*}}}*/
+
+/*{{{ String Helper functions*/
+static gchar *
+replace_substring (const gchar *haystack, const gchar *needle, const gchar *replacement)
+{
+ GRegex *regex;
+ gchar *result;
+
+ regex = g_regex_new (needle, 0, 0, NULL);
+ result = g_regex_replace_literal (regex, haystack, -1, 0, replacement, 0, NULL);
+ g_regex_unref (regex);
+ return result;
+}
+/*}}}*/
+
+/*{{{ Kernel creation*/
+static gchar *
+replace_parameter_dashes (UfoUniRecoParameter parameter)
+{
+ return replace_substring (parameter_values[parameter].value_nick, "-", "_");
+}
+
+static gchar *
+get_kernel_parameter_name (UfoUniRecoParameter parameter)
+{
+ gchar **entries, *result, *param_kernel_name;
+ param_kernel_name = replace_parameter_dashes (parameter);
+
+ if (is_parameter_positional (parameter)) {
+ entries = g_strsplit (param_kernel_name, "_", 3);
+ if (g_strcmp0 (entries[0], param_kernel_name)) {
+ result = g_strconcat (entries[0], "_", entries[1], NULL);
+ }
+ g_strfreev (entries);
+ } else {
+ result = g_strdup (param_kernel_name);
+ }
+ g_free (param_kernel_name);
+
+ return result;
+}
+
+static gchar *
+make_template (UfoGeneralBackprojectTaskPrivate *priv)
+{
+ gchar *template, *definitions, *header, *header_1, *body, *kernel_parameter_name, *tmp;
+
+ if (!(definitions = ufo_resources_get_kernel_source (priv->resources, "general_bp_definitions.in", NULL))) {
+ g_warning ("Error obtaining general backprojection kernel header template");
+ return NULL;
+ }
+ if (!(body = ufo_resources_get_kernel_source (priv->resources, "general_bp_body.in", NULL))) {
+ g_warning ("Error obtaining general backprojection kernel body template");
+ return NULL;
+ }
+ if (priv->vectorized) {
+ if (!(header = ufo_resources_get_kernel_source (priv->resources, "general_bp_header_vector.in", NULL))) {
+ g_warning ("Error obtaining general backprojection kernel header template");
+ return NULL;
+ }
+ if (priv->parameter != UFO_UNI_RECO_PARAMETER_Z) {
+ kernel_parameter_name = get_kernel_parameter_name (priv->parameter);
+ tmp = g_strconcat (kernel_parameter_name, "_global", NULL);
+ header_1 = replace_substring (header, kernel_parameter_name, tmp);
+ g_free (tmp);
+ g_free (header);
+ header = header_1;
+ }
+ header_1 = replace_substring (header, "%memspace%", "global ");
+ g_free (header);
+ header = header_1;
+ } else {
+ if (!(header = ufo_resources_get_kernel_source (priv->resources, "general_bp_header_scalar.in", NULL))) {
+ g_warning ("Error obtaining general backprojection kernel header template");
+ return NULL;
+ }
+ }
+ template = g_strconcat (definitions, header, body, NULL);
+ g_free (definitions);
+ g_free (header);
+ g_free (body);
+
+ return template;
+}
+
+/**
+ * make_args:
+ * @burst: (in): number of processed projections in the kernel
+ * @fmt: (in): format string which will be transformed to the projection index
+ *
+ * Make kernel arguments.
+ */
+static gchar *
+make_args (gint burst, const gchar *fmt)
+{
+ gint i;
+ gulong size, written;
+ gchar *one, *str, *ptr;
+
+ size = strlen (fmt) + 1;
+ one = g_strnfill (size, 0);
+ str = g_strnfill (burst * size, 0);
+ ptr = str;
+
+ for (i = 0; i < burst; i++) {
+ written = g_snprintf (one, size, fmt, i);
+ if (written > size) {
+ g_free (one);
+ g_free (str);
+ return NULL;
+ }
+ ptr = g_stpcpy (ptr, one);
+ }
+ g_free (one);
+
+ return str;
+}
+
+/**
+ * make_type_conversion:
+ * @compute_type: (in): data type for internal computations
+ * @store_type: (in): output volume data type
+ *
+ * Make conversions necessary for computation and output data types.
+ */
+static gchar *
+make_type_conversion (const gchar *compute_type, const gchar *store_type)
+{
+ gchar *code;
+ gulong maxval = get_integer_maximum (store_type);
+
+ if (maxval) {
+ code = g_strdup_printf ("(%s) clamp ((%s)(gray_limit.y * (norm_factor * result - gray_limit.x)), (%s) 0.0, (%s) %lu.0)",
+ store_type, compute_type, compute_type, compute_type, maxval);
+ } else {
+ code = g_strdup_printf ("(%s) (norm_factor * result)", store_type);
+ }
+
+ return code;
+}
+
+/**
+ * make_parameter_initial_assignment:
+ * @parameter: (in): parameter which represents the third reconstruction axis
+ *
+ * Make initial parameter assignment for vectorized kernels which need to first
+ * copy the global values to a private variable.
+ */
+static gchar *
+make_parameter_initial_assignment (UfoUniRecoParameter parameter)
+{
+ gchar *code = NULL, *kernel_parameter_name;
+
+ if (parameter == UFO_UNI_RECO_PARAMETER_Z) {
+ code = g_strdup ("");
+ } else {
+ kernel_parameter_name = get_kernel_parameter_name (parameter);
+ if (is_parameter_positional (parameter)) {
+ code = g_strconcat ("cfloat3 ", kernel_parameter_name, ";\n", NULL);
+ } else if (is_parameter_angular (parameter)) {
+ code = g_strconcat ("cfloat2 ", kernel_parameter_name, ";\n", NULL);
+ }
+ g_free (kernel_parameter_name);
+ }
+
+ return code;
+}
+
+/**
+ * make_parameter_assignment:
+ * @parameter: (in): parameter which represents the third reconstruction axis
+ *
+ * Make parameter assignment.
+ */
+static gchar *
+make_parameter_assignment (UfoUniRecoParameter parameter)
+{
+ gchar **entries, *param_kernel_name;
+ gchar *code = NULL;
+ param_kernel_name = replace_parameter_dashes (parameter);
+
+ if (parameter == UFO_UNI_RECO_PARAMETER_Z) {
+ code = g_strdup ("\tvoxel_0.z = region[idz].x;\n");
+ } else if (is_parameter_positional (parameter)) {
+ entries = g_strsplit (param_kernel_name, "_", 3);
+ if (g_strcmp0 (entries[0], param_kernel_name)) {
+ code = g_strconcat ("\t", entries[0], "_", entries[1], ".", entries[2], " = region[idz].x;\n", NULL);
+ }
+ g_strfreev (entries);
+ } else if (is_parameter_angular (parameter)) {
+ code = g_strconcat ("\t", param_kernel_name, " = region[idz];\n", NULL);
+ }
+ g_free (param_kernel_name);
+
+ return code;
+}
+
+/**
+ * make_volume_transformation:
+ * @values: (in): sine and cosine angle values
+ * @point: 3D point which will be rotated
+ * @suffix: suffix which comes after values
+ *
+ * Inplace point rotation about the three coordinate axes.
+ */
+static gchar *
+make_volume_transformation (const gchar *values, const gchar *point, const gchar *suffix)
+{
+ return g_strdup_printf ("\t%s = rotate_z (%s_z%s, %s);"
+ "\n\t%s = rotate_y (%s_y%s, %s);"
+ "\n\t%s = rotate_x (%s_x%s, %s);\n",
+ point, values, suffix, point,
+ point, values, suffix, point,
+ point, values, suffix, point);
+}
+
+/**
+ * make_static_transformations:
+ * @vectorized (in): geometry parameters are vectors
+ * @with_volume: (in): rotate reconstructed volume
+ * @perpendicular_detector: (in): is the detector perpendicular to the beam
+ * @parallel_beam: (in): is the beam parallel
+ *
+ * Make static transformations independent from the tomographic rotation angle.
+ */
+static gchar *
+make_static_transformations (gboolean vectorized, gboolean with_volume,
+ gboolean perpendicular_detector, gboolean parallel_beam)
+{
+ gchar *code = g_strnfill (8192, 0);
+ gchar *current = code;
+ gchar *detector_transformation, *volume_transformation;
+ const gchar *voxel_0 = vectorized ? "voxel" : "voxel_0";
+
+ if (parallel_beam) {
+ if (vectorized) {
+ current = g_stpcpy (current,"\tvoxel = voxel_0;\n");
+ }
+ } else {
+ if (vectorized) {
+ current = g_stpcpy (current, "\t// Magnification\n\tvoxel = voxel_0 * -native_divide(source_position[%d].y, "
+ "(detector_position[%d].y - source_position[%d].y));\n");
+ } else {
+ current = g_stpcpy (current, "// Magnification\n\tvoxel_0 *= -native_divide(source_position.y, "
+ "(detector_position.y - source_position.y));\n");
+ }
+ }
+ if (!perpendicular_detector) {
+ detector_transformation = make_volume_transformation ("detector_angle", "detector_normal", "[%d]");
+ current = g_stpcpy (current, "\tdetector_normal = (float3)(0.0, -1.0, 0.0);\n");
+ current = g_stpcpy (current, detector_transformation);
+ current = g_stpcpy (current, "\n\tdetector_offset = -dot (detector_position[%d], detector_normal);\n");
+ g_free (detector_transformation);
+ } else if (!parallel_beam) {
+ current = g_stpcpy (current, "\n\tproject_tmp = detector_position[%d].y - source_position[%d].y;\n");
+ }
+ if (with_volume) {
+ volume_transformation = make_volume_transformation ("volume_angle", voxel_0, "[%d]");
+ current = g_stpcpy (current, volume_transformation);
+ g_free (volume_transformation);
+ }
+ if (!(perpendicular_detector || parallel_beam)) {
+ current = g_stpcpy (current,
+ "\n\ttmp_transformation = "
+ "- (detector_offset + dot (source_position[%d], detector_normal));\n");
+ }
+
+ return code;
+}
+
+/**
+ * make_projection_computation:
+ * @perpendicular_detector: (in): is the detector perpendicular to the beam
+ * @parallel_beam: (in): is the beam parallel
+ *
+ * Make voxel projection calculation with the least possible operations based on
+ * geometry settings.
+ */
+static const gchar *
+make_projection_computation (gboolean perpendicular_detector, gboolean parallel_beam)
+{
+ const gchar *code;
+
+ if (perpendicular_detector) {
+ if (parallel_beam) {
+ code = "\t// Perpendicular detector in combination with parallel beam geometry, i.e.\n"
+ "\t// voxel.xz is directly the detector coordinate, no transformation necessary\n";
+ } else {
+ code = "\tvoxel = mad (native_divide (project_tmp, (voxel.y - source_position[%d].y)), voxel, source_position[%d]);\n";
+ }
+ } else {
+ if (parallel_beam) {
+ code = "\tvoxel.y = -native_divide (mad (voxel.z, detector_normal.z, "
+ "mad (voxel.x, detector_normal.x, detector_offset)), detector_normal.y);\n";
+ } else {
+ code = "\tvoxel -= source_position[%d];\n"
+ "\tvoxel = mad (native_divide (tmp_transformation, dot (voxel, detector_normal)), voxel, source_position[%d]);\n";
+ }
+ }
+
+ return code;
+}
+
+/**
+ * make_transformations:
+ * @parameter: parameter reconstructed along the third dimension
+ * @vectorized (in): geometry parameters are vectors
+ * @burst (in): number of projections processed by the kernel
+ * @with_axis: (in): do computations related with rotation axis
+ * @perpendicular_detector: (in): is the detector perpendicular to the the mean
+ * beam direction
+ * @parallel_beam: (in): is the beam parallel
+ * @compute_type: (in): data type for internal computations
+ *
+ * Make voxel projection calculation with the least possible operations based on
+ * geometry settings.
+ */
+static gchar *
+make_transformations (UfoUniRecoParameter parameter, gboolean vectorized, gint burst, gboolean with_volume,
+ gboolean with_axis, gboolean perpendicular_detector, gboolean parallel_beam,
+ const gchar *compute_type)
+{
+ gint i;
+ size_t written = 0;
+ gchar *code_fmt, *code, *current, *volume_transformation, *pretransformation, *pixel_reading,
+ *str_iteration, *parameter_assignment, *kernel_parameter_name, *tmp, *tmp_2, *str_index;
+ const guint snippet_size = 16384;
+ const guint size = burst * snippet_size;
+ /* Based on eq. 30 from "Direct cone beam SPECT reconstruction with camera tilt" */
+ const gchar *slice_coefficient =
+ "\t// Get the value and weigh it (source_position is negative, so -voxel.y\n"
+ "\tcoeff = native_divide (source_position[%d].y - detector_position[%d].y, (source_position[%d].y - voxel.y));\n";
+ const gchar *detector_transformation =
+ "\tvoxel -= detector_position[%d];\n"
+ "\tvoxel = rotate_x ((cfloat2)(-detector_angle_x[%d].x, detector_angle_x[%d].y), voxel);\n"
+ "\tvoxel = rotate_y ((cfloat2)(-detector_angle_y[%d].x, detector_angle_y[%d].y), voxel);\n"
+ "\tvoxel = rotate_z ((cfloat2)(-detector_angle_z[%d].x, detector_angle_z[%d].y), voxel);\n";
+
+ code_fmt = g_strnfill (snippet_size, 0);
+ code = g_strnfill (size, 0);
+ kernel_parameter_name = get_kernel_parameter_name (parameter);
+
+ current = g_stpcpy (code_fmt, "\t/* Tomographic rotation angle %02d */\n");
+
+ if (vectorized) {
+ if (is_parameter_positional (parameter)) {
+ /* If the parameter is positional, first load the global 3-tuple to
+ * which it belongs (e.g. if it is source_position_x we need to load the
+ * source_position) and change only the part specified by parameter.
+ * This way the two components can still be governed by the tomographic
+ * angle and the third can be optimized for. */
+ tmp = g_strconcat ("\t", kernel_parameter_name, " = ", kernel_parameter_name, "_global", "[%d];\n", NULL);
+ parameter_assignment = make_parameter_assignment (parameter);
+ current = g_stpcpy (current, tmp);
+ current = g_stpcpy (current, parameter_assignment);
+ g_free (parameter_assignment);
+ g_free (tmp);
+ }
+ /* For vectorized kernel all static transformations become non-static */
+ if ((pretransformation = make_static_transformations (vectorized, with_volume, perpendicular_detector,
+ parallel_beam)) == NULL) {
+ g_warning ("Error making static transformations");
+ return NULL;
+ }
+ current = g_stpcpy (current, pretransformation);
+ g_free (pretransformation);
+ }
+
+ current = g_stpcpy (current, vectorized ? "\tvoxel = rotate_z (tomo_%02d, voxel);\n" :
+ "\tvoxel = rotate_z (tomo_%02d, voxel_0);\n");
+
+ if (with_axis) {
+ /* Tilted axis of rotation */
+ volume_transformation = g_strdup_printf ("\t%s = rotate_y (%s_y%s, %s);\n"
+ "\t%s = rotate_x (%s_x%s, %s);\n",
+ "voxel", "axis_angle", "[%d]", "voxel",
+ "voxel", "axis_angle", "[%d]", "voxel");
+ current = g_stpcpy (current, volume_transformation);
+ current = g_stpcpy (current, "\n");
+ g_free (volume_transformation);
+ }
+ if (!parallel_beam) {
+ /* FDK normalization computation */
+ current = g_stpcpy (current, slice_coefficient);
+ }
+
+ /* Voxel projection on the detector */
+ current = g_stpcpy (current,
+ "\t// Compute the voxel projected on the detector plane in the global coordinates\n"
+ "\t// V = S + u * (V - S)\n");
+ current = g_stpcpy (current, make_projection_computation (perpendicular_detector, parallel_beam));
+
+ if (!perpendicular_detector) {
+ /* Transform global coordinates to detector coordinates */
+ current = g_stpcpy (current,
+ "\t// Transform the projected coordinates to the detector coordinates, i.e. rotate the\n"
+ "\t// projected voxel to the detector plane\n");
+ current = g_stpcpy (current, detector_transformation);
+ }
+
+ /* Computational data type adjustment */
+ pixel_reading = g_strconcat ("\tresult += read_imagef (projection_%02d, sampler, ",
+ !g_strcmp0 (compute_type, "float") ? "(" : "convert_float2(",
+ "voxel.xz + center_position[%d].xz)).x", NULL);
+ current = g_stpcpy (current, pixel_reading);
+ g_free (pixel_reading);
+
+ /* FDK normalization application */
+ if (parallel_beam) {
+ current = g_stpcpy (current, ";\n\n");
+ } else {
+ current = g_stpcpy (current, " * coeff * coeff;\n\n");
+ }
+
+ if (vectorized) {
+ if (parameter != UFO_UNI_RECO_PARAMETER_Z) {
+ tmp_2 = g_strconcat (kernel_parameter_name, "\\[%d\\]", NULL);
+ tmp = replace_substring (code_fmt, tmp_2, kernel_parameter_name);
+ g_free (code_fmt);
+ g_free (tmp_2);
+ code_fmt = tmp;
+ }
+ } else {
+ tmp = replace_substring (code_fmt, "\\[%d\\]", "");
+ g_free (code_fmt);
+ code_fmt = tmp;
+ }
+
+ current = code;
+ for (i = 0; i < burst; i++) {
+ /* %02d would result in octa-based indexing which would crash the kernel for burst > 7 */
+ str_index = g_strdup_printf ("iteration + %d", i);
+ tmp = replace_substring (code_fmt, "%d", str_index);
+ g_free (str_index);
+ str_index = g_strdup_printf ("%02d", i);
+ str_iteration = replace_substring (tmp, "%02d", str_index);
+ g_free (str_index);
+ g_free (tmp);
+ written += strlen (str_iteration);
+ if (written > size) {
+ g_free (code_fmt);
+ g_free (code);
+ g_free (str_iteration);
+ return NULL;
+ }
+ current = g_stpcpy (current, str_iteration);
+ g_free (str_iteration);
+ }
+ g_free (code_fmt);
+ g_free (kernel_parameter_name);
+
+ return code;
+}
+
+/**
+ * make_kernel:
+ * @template (in): kernel template string
+ * @vectorized (in): geometry parameters are vectors
+ * @burst (in): how many projections to process in one kernel invocation
+ * @with_axis (in): rotate the rotation axis
+ * @with_volume: (in): rotate reconstructed volume
+ * @perpendicular_detector: (in): is the detector perpendicular to the beam
+ * @parallel_beam: (in): is the beam parallel
+ * @compute_type (in): data type for calculations (one of "half", "float", "double")
+ * @result_type (in): data type for storing the intermediate result (one of "half", "float", "double")
+ * @store_type (in): data type of the output volume (one of "half", "float", "double",
+ * "uchar", "ushort", "uint")
+ * @parameter: (in): parameter which represents the third reconstruction axis
+ *
+ * Make backprojection kernel.
+ */
+static gchar *
+make_kernel (gchar *template, gboolean vectorized, gint burst, gboolean with_axis, gboolean with_volume,
+ gboolean perpendicular_detector, gboolean parallel_beam, const gchar *compute_type,
+ const gchar *result_type, const gchar *store_type, UfoUniRecoParameter parameter)
+{
+ const gchar *double_pragma_def, *double_pragma, *half_pragma_def, *half_pragma,
+ *image_args_fmt, *trigonomoerty_args_fmt;
+ gchar *image_args, *trigonometry_args, *type_conversion, *parameter_assignment, *local_assignment,
+ *static_transformations, *transformations, *code_tmp, *code, *tmp, **parts;
+ gboolean positional_param = is_parameter_positional (parameter);
+
+ double_pragma_def = "#pragma OPENCL EXTENSION cl_khr_fp64 : enable\n";
+ half_pragma_def = "#pragma OPENCL EXTENSION cl_khr_fp16 : enable\n\n";
+ image_args_fmt = "\t\t\t read_only image2d_t projection_%02d,\n";
+ trigonomoerty_args_fmt = "\t\t\t const cfloat2 tomo_%02d,\n";
+ parts = g_strsplit (template, "%tmpl%", 9);
+
+ if ((image_args = make_args (burst, image_args_fmt)) == NULL) {
+ g_warning ("Error making image arguments");
+ return NULL;
+ }
+ if ((trigonometry_args = make_args (burst, trigonomoerty_args_fmt)) == NULL) {
+ g_warning ("Error making trigonometric arguments");
+ return NULL;
+ }
+ if ((type_conversion = make_type_conversion (compute_type, store_type)) == NULL) {
+ g_warning ("Error making type conversion");
+ return NULL;
+ }
+ if (vectorized) {
+ /* First make a private variable with the same name as the kernel
+ * argument name in case of scalar kernel. The buffer value will be
+ * loaded for each tomographic angle. */
+ local_assignment = make_parameter_initial_assignment (parameter);
+ if (local_assignment == NULL) {
+ g_warning ("Wrong parameter name");
+ return NULL;
+ }
+ if (positional_param) {
+ /* In case of positional parameter two components can still be
+ * controlled by the tomographic angle and one can be the actual
+ * parameter specified by region, thus the variable needs to be set
+ * for every tomographic angle separately */
+ parameter_assignment = g_strdup ("");
+ }
+ }
+ if (!vectorized || (vectorized && !positional_param)) {
+ parameter_assignment = make_parameter_assignment (parameter);
+ if (parameter_assignment == NULL) {
+ g_warning ("Wrong parameter name");
+ return NULL;
+ }
+ }
+
+ if (!vectorized) {
+ if ((static_transformations = make_static_transformations(FALSE, with_volume, perpendicular_detector,
+ parallel_beam)) == NULL) {
+ g_warning ("Error making static transformations");
+ return NULL;
+ }
+ tmp = replace_substring (static_transformations, "\\[%d\\]", "");
+ g_free (static_transformations);
+ static_transformations = tmp;
+ }
+ if ((transformations = make_transformations (parameter, vectorized, burst, with_volume, with_axis,
+ perpendicular_detector, parallel_beam, compute_type)) == NULL) {
+ g_warning ("Error making tomographic-angle-based transformations");
+ return NULL;
+ }
+ if (!(g_strcmp0 (compute_type, "double") && g_strcmp0 (result_type, "double"))) {
+ double_pragma = double_pragma_def;
+ } else {
+ double_pragma = "";
+ }
+ if (!(g_strcmp0 (compute_type, "half") && g_strcmp0 (result_type, "half"))) {
+ half_pragma = half_pragma_def;
+ } else {
+ half_pragma = "";
+ }
+ code_tmp = g_strconcat (double_pragma, half_pragma,
+ parts[0], image_args,
+ parts[1], trigonometry_args,
+ parts[2], vectorized ? local_assignment : "",
+ parts[3], parameter_assignment,
+ parts[4], vectorized ? "" : static_transformations,
+ parts[5], transformations,
+ parts[6], type_conversion,
+ parts[7], type_conversion,
+ parts[8], NULL);
+ code = replace_substring (code_tmp, "cfloat", compute_type);
+ g_free (code_tmp);
+ code_tmp = replace_substring (code, "rtype", result_type);
+ g_free (code);
+ code = replace_substring (code_tmp, "stype", store_type);
+
+ g_free (image_args);
+ g_free (trigonometry_args);
+ g_free (type_conversion);
+ g_free (parameter_assignment);
+ if (vectorized) {
+ g_free (local_assignment);
+ } else {
+ g_free (static_transformations);
+ }
+ g_free (transformations);
+ g_free (code_tmp);
+ g_strfreev (parts);
+
+ return code;
+}
+/*}}}*/
+
+/*{{{ OpenCL helper functions */
+DEFINE_CREATE_REGIONS (cl_float)
+DEFINE_CREATE_REGIONS (cl_double)
+
+static void
+create_images (UfoGeneralBackprojectTaskPrivate *priv, gsize width, gsize height)
+{
+ cl_image_format image_fmt;
+ cl_int cl_error;
+ guint i;
+
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Creating images %lu x %lu", width, height);
+
+ /* TODO: dangerous, don't rely on the ufo-buffer */
+ image_fmt.image_channel_order = CL_INTENSITY;
+ image_fmt.image_channel_data_type = CL_FLOAT;
+
+ for (i = 0; i < priv->burst; i++) {
+ /* TODO: what about the "other" API? */
+ priv->projections[i] = clCreateImage2D (priv->context,
+ CL_MEM_READ_ONLY,
+ &image_fmt,
+ width,
+ height,
+ 0,
+ NULL,
+ &cl_error);
+ UFO_RESOURCES_CHECK_CLERR (cl_error);
+ }
+}
+
+static cl_mem
+transfer_host_to_device (cl_context context, gpointer host_array, gsize num_bytes)
+{
+ cl_int cl_error;
+ cl_mem device_array;
+
+ device_array = clCreateBuffer (context,
+ CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
+ num_bytes,
+ host_array,
+ &cl_error);
+ UFO_RESOURCES_CHECK_CLERR (cl_error);
+
+ return device_array;
+}
+
+DEFINE_TRANSFER_ANGULAR_ARGUMENT (cl_float)
+DEFINE_TRANSFER_ANGULAR_ARGUMENT (cl_double)
+DEFINE_TRANSFER_POSITINAL_ARGUMENT (cl_float)
+DEFINE_TRANSFER_POSITINAL_ARGUMENT (cl_double)
+DEFINE_SET_ANGULAR_VECTOR_KERNEL_ARGUMENT (cl_float)
+DEFINE_SET_ANGULAR_VECTOR_KERNEL_ARGUMENT (cl_double)
+DEFINE_SET_STATIC_SCALAR_ARGUMENTS (cl_float)
+DEFINE_SET_STATIC_SCALAR_ARGUMENTS (cl_double)
+DEFINE_SET_STATIC_VECTOR_ARGUMENTS (cl_float)
+DEFINE_SET_STATIC_VECTOR_ARGUMENTS (cl_double)
+DEFINE_COMPUTE_SLICE_REGION (cl_float)
+DEFINE_COMPUTE_SLICE_REGION (cl_double)
+DEFINE_SET_STATIC_ARGS (cl_float)
+DEFINE_SET_STATIC_ARGS (cl_double)
+
+/**
+ * Copy buffer to OpenCL image.
+ */
+static void
+copy_to_image (const cl_command_queue cmd_queue,
+ UfoBuffer *input,
+ cl_mem output,
+ gsize width,
+ gsize height)
+{
+ cl_event event;
+ cl_int errcode;
+ cl_mem input_array;
+ const size_t origin[] = {0, 0, 0};
+ const size_t region[] = {width, height, 1};
+
+ input_array = ufo_buffer_get_device_array (input, cmd_queue);
+ errcode = clEnqueueCopyBufferToImage (cmd_queue,
+ input_array,
+ output,
+ 0, origin, region,
+ 0, NULL, &event);
+
+ UFO_RESOURCES_CHECK_CLERR (errcode);
+ UFO_RESOURCES_CHECK_CLERR (clWaitForEvents (1, &event));
+ UFO_RESOURCES_CHECK_CLERR (clReleaseEvent (event));
+}
+
+static void
+node_setup (UfoGeneralBackprojectTaskPrivate *priv,
+ UfoGpuNode *node)
+{
+ guint i;
+ gboolean with_axis, with_volume, parallel_beam, perpendicular_detector;
+ gchar *template, *kernel_code, *compiler_options = NULL;
+ const gchar *node_name;
+ GValue *node_name_gvalue;
+ const gchar compiler_options_tmpl[] = "-cl-nv-maxrregcount=%u";
+ UfoUniRecoNodeProps *node_props;
+
+ /* GPU type specific settings */
+ node_name_gvalue = ufo_gpu_node_get_info (node, UFO_GPU_NODE_INFO_NAME);
+ node_name = g_value_get_string (node_name_gvalue);
+ if (!(node_props = g_hash_table_lookup (priv->node_props_table, node_name))) {
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "GPU with name %s not in database", node_name);
+ node_props = g_hash_table_lookup (priv->node_props_table, "GENERIC");
+ }
+ if (!priv->burst) {
+ priv->burst = node_props->burst;
+ }
+ if (node_props->max_regcount) {
+ compiler_options = g_strdup_printf (compiler_options_tmpl, node_props->max_regcount);
+ }
+ g_log ("gbp", G_LOG_LEVEL_DEBUG,
+ "GPU node %s properties: burst: %u, compiler options: '%s'",
+ node_name, priv->burst, compiler_options);
+ g_value_unset (node_name_gvalue);
+
+ /* Initialization */
+ /* Assume the most efficient geometry, change if necessary */
+ with_axis = is_axis_parameter (priv->parameter) ||
+ !(ufo_scarray_is_almost_zero (priv->geometry->axis->angle->x) &&
+ ufo_scarray_is_almost_zero (priv->geometry->axis->angle->y));
+ with_volume = is_volume_parameter (priv->parameter) || !ufo_scpoint_are_almost_zero (priv->geometry->volume_angle);
+ perpendicular_detector = !is_detector_rotation_parameter (priv->parameter) &&
+ !is_detector_position_parameter (priv->parameter) &&
+ ufo_scpoint_are_almost_zero (priv->geometry->detector->angle);
+ parallel_beam = TRUE;
+ /* Actual parameter setup */
+ for (i = 0; i < priv->num_projections; i++) {
+ parallel_beam = parallel_beam && isinf (ufo_scarray_get_double (priv->geometry->source_position->y, i));
+ }
+ priv->vectorized = (ufo_scarray_has_n_values (priv->geometry->axis->angle->x, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->axis->angle->y, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->volume_angle->x, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->volume_angle->y, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->volume_angle->z, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->detector->angle->x, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->detector->angle->y, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->detector->angle->z, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->detector->position->x, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->detector->position->y, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->detector->position->z, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->source_position->x, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->source_position->y, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->source_position->z, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->axis->position->x, priv->num_projections) ||
+ ufo_scarray_has_n_values (priv->geometry->axis->position->z, priv->num_projections));
+
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "vectorized: %d, parameter: %s with axis: %d, with volume: %d, "
+ "perpendicular detector: %d, parallel beam: %d, "
+ "compute type: %s, result type: %s, store type: %s",
+ priv->vectorized, parameter_values[priv->parameter].value_nick, with_axis, with_volume,
+ perpendicular_detector, parallel_beam,
+ compute_type_values[priv->compute_type].value_nick,
+ ft_values[priv->result_type].value_nick,
+ st_values[priv->store_type].value_nick);
+
+ if ((template = make_template (priv)) == NULL) {
+ return;
+ }
+
+ /* Create kernel source code based on geometry settings */
+ kernel_code = make_kernel (template, priv->vectorized, priv->burst, with_axis, with_volume,
+ perpendicular_detector, parallel_beam,
+ compute_type_values[priv->compute_type].value_nick,
+ ft_values[priv->result_type].value_nick,
+ st_values[priv->store_type].value_nick,
+ priv->parameter);
+ if (!kernel_code) {
+ g_free (template);
+ g_free (compiler_options);
+ return;
+ }
+ priv->kernel = ufo_resources_get_kernel_from_source (priv->resources,
+ kernel_code,
+ "backproject",
+ compiler_options,
+ NULL);
+ g_free (kernel_code);
+
+ if (priv->num_projections % priv->burst) {
+ kernel_code = make_kernel (template, priv->vectorized, priv->num_projections % priv->burst,
+ with_axis, with_volume,
+ perpendicular_detector, parallel_beam,
+ compute_type_values[priv->compute_type].value_nick,
+ ft_values[priv->result_type].value_nick,
+ st_values[priv->store_type].value_nick,
+ priv->parameter);
+ if (!kernel_code) {
+ g_free (template);
+ g_free (compiler_options);
+ return;
+ }
+
+ /* If num_projections % priv->burst != 0 we need one more kernel to process the remaining projections */
+ priv->rest_kernel = ufo_resources_get_kernel_from_source (priv->resources,
+ kernel_code,
+ "backproject",
+ compiler_options,
+ NULL);
+ /* g_printf ("%s", kernel_code); */
+ g_free (kernel_code);
+ }
+ g_free (template);
+ g_free (compiler_options);
+
+ if (priv->kernel) {
+ UFO_RESOURCES_CHECK_CLERR (clRetainKernel (priv->kernel));
+ }
+ if (priv->rest_kernel) {
+ UFO_RESOURCES_CHECK_CLERR (clRetainKernel (priv->rest_kernel));
+ }
+}
+/*}}}*/
+
+UfoNode *
+ufo_general_backproject_task_new (void)
+{
+ return UFO_NODE (g_object_new (UFO_TYPE_GENERAL_BACKPROJECT_TASK, NULL));
+}
+
+static void
+ufo_general_backproject_task_setup (UfoTask *task,
+ UfoResources *resources,
+ GError **error)
+{
+ cl_int cl_error;
+ guint i;
+ GValue tomo_angle = G_VALUE_INIT;
+ UfoGeneralBackprojectTaskPrivate *priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (task);
+ priv->resources = g_object_ref (resources);
+ priv->kernel = NULL;
+ priv->rest_kernel = NULL;
+ priv->projections = NULL;
+ priv->chunks = NULL;
+ priv->cl_regions = NULL;
+ priv->vector_arguments = NULL;
+
+ /* Check parameter values */
+ if (!priv->num_projections) {
+ g_set_error (error, UFO_TASK_ERROR, UFO_TASK_ERROR_SETUP,
+ "Number of projections not set");
+ return;
+ }
+
+ if (!ufo_scarray_has_n_values (priv->geometry->axis->angle->z, priv->num_projections)) {
+ /* Create equidistant tomographic angles and treat the one specified in
+ * the current priv->geometry->axis->angle->z as offset */
+ ufo_scarray_free (priv->geometry->axis->angle->z);
+ priv->geometry->axis->angle->z = ufo_scarray_new (priv->num_projections, G_TYPE_DOUBLE, NULL);
+ g_value_init (&tomo_angle, G_TYPE_DOUBLE);
+ for (i = 0; i < priv->num_projections; i++) {
+ g_value_set_double (&tomo_angle, ((gdouble) i) / priv->num_projections * priv->overall_angle);
+ ufo_scarray_insert (priv->geometry->axis->angle->z, i, &tomo_angle);
+ }
+ g_value_unset (&tomo_angle);
+ }
+
+ if (priv->gray_map_min >= priv->gray_map_max &&
+ (priv->store_type == ST_UCHAR || priv->store_type == ST_USHORT || priv->store_type == ST_UINT)) {
+ g_set_error (error, UFO_TASK_ERROR, UFO_TASK_ERROR_SETUP,
+ "Gray mapping minimum must be less then the maximum");
+ return;
+ }
+
+ priv->node_props_table = ufo_get_node_props_table ();
+
+ /* Set OpenCL variables */
+ priv->context = ufo_resources_get_context (resources);
+ UFO_RESOURCES_CHECK_CLERR (clRetainContext (priv->context));
+ priv->sampler = clCreateSampler (priv->context, (cl_bool) FALSE, priv->addressing_mode, CL_FILTER_LINEAR, &cl_error);
+ UFO_RESOURCES_CHECK_CLERR (cl_error);
+}
+
+static void
+ufo_general_backproject_task_get_requisition (UfoTask *task,
+ UfoBuffer **inputs,
+ UfoRequisition *requisition,
+ GError **error)
+{
+ UfoGeneralBackprojectTaskPrivate *priv;
+ UfoGpuNode *node;
+ cl_command_queue cmd_queue;
+ UfoRequisition in_req;
+ gdouble region_start, region_stop, region_step;
+ gsize slice_size, chunk_size, volume_size, projections_size;
+ GValue *max_global_mem_size_gvalue, *max_mem_alloc_size_gvalue;
+ cl_ulong max_global_mem_size, max_mem_alloc_size;
+ cl_int cl_error;
+ guint i;
+ typedef void (*CreateRegionFunc) (UfoGeneralBackprojectTaskPrivate *, const cl_command_queue,
+ const gdouble, const gdouble);
+ typedef void (*SetStaticArgsFunc) (UfoTask *, UfoRequisition *, const cl_kernel);
+ CreateRegionFunc create_regions[2] = {create_regions_cl_float, create_regions_cl_double};
+ SetStaticArgsFunc set_static_args[2] = {set_static_args_cl_float, set_static_args_cl_double};
+ GValue g_value_int = G_VALUE_INIT;
+ g_value_init (&g_value_int, G_TYPE_INT);
+
+ priv = UFO_GENERAL_BACKPROJECT_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);
+ g_assert_true (ufo_scarray_has_n_values (priv->region, 3));
+ requisition->n_dims = 2;
+ ufo_buffer_get_requisition (inputs[0], &in_req);
+
+ if (ufo_scarray_get_int (priv->region_x, 2) == 0) {
+ /* If the slice width is not set, reconstruct full width */
+ requisition->dims[0] = in_req.dims[0];
+ } else {
+ requisition->dims[0] = REGION_SIZE (priv->region_x);
+ }
+ if (ufo_scarray_get_int (priv->region_y, 2) == 0) {
+ /* If the slice height is not set, reconstruct full height, which is the
+ * same as width */
+ requisition->dims[1] = in_req.dims[0];
+ } else {
+ requisition->dims[1] = REGION_SIZE (priv->region_y);
+ }
+
+ if (!priv->kernel) {
+ /* First iteration, setup kernels */
+ node_setup (priv, node);
+ if (!priv->kernel) {
+ g_set_error_literal (error, UFO_TASK_ERROR, UFO_TASK_ERROR_GET_REQUISITION,
+ "Error creating backprojection kernels");
+ return;
+ }
+
+ if (UFO_MATH_ARE_ALMOST_EQUAL (ufo_scarray_get_double (priv->region, 2), 0)) {
+ /* Conservative approach, reconstruct just one slice */
+ region_start = 0.0f;
+ region_stop = 1.0f;
+ region_step = 1.0f;
+ } else {
+ region_start = ufo_scarray_get_double (priv->region, 0);
+ region_stop = ufo_scarray_get_double (priv->region, 1);
+ region_step = ufo_scarray_get_double (priv->region, 2);
+ }
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "region: %g %g %g", region_start, region_stop, region_step);
+ priv->num_slices = (gsize) ceil ((region_stop - region_start) / region_step);
+ max_global_mem_size_gvalue = ufo_gpu_node_get_info (node, UFO_GPU_NODE_INFO_GLOBAL_MEM_SIZE);
+ max_global_mem_size = g_value_get_ulong (max_global_mem_size_gvalue);
+ g_value_unset (max_global_mem_size_gvalue);
+ projections_size = priv->burst * in_req.dims[0] * in_req.dims[1] * sizeof (cl_float);
+ slice_size = requisition->dims[0] * requisition->dims[1] * get_type_size (priv->store_type);
+ volume_size = slice_size * priv->num_slices;
+ max_mem_alloc_size_gvalue = ufo_gpu_node_get_info (node, UFO_GPU_NODE_INFO_MAX_MEM_ALLOC_SIZE);
+ max_mem_alloc_size = g_value_get_ulong (max_mem_alloc_size_gvalue);
+ g_value_unset (max_mem_alloc_size_gvalue);
+ priv->num_slices_per_chunk = (guint) floor ((gdouble) MIN (max_mem_alloc_size, volume_size) / ((gdouble) slice_size));
+ if (projections_size + volume_size > max_global_mem_size) {
+ g_set_error_literal (error, UFO_TASK_ERROR, UFO_TASK_ERROR_GET_REQUISITION,
+ "Volume size doesn't fit to memory");
+ return;
+ }
+
+ priv->projections = (cl_mem *) g_malloc (priv->burst * sizeof (cl_mem));
+ /* Create subvolumes (because one large volume might be larger than the maximum allocatable memory chunk */
+ priv->num_chunks = (priv->num_slices - 1) / priv->num_slices_per_chunk + 1;
+ chunk_size = priv->num_slices_per_chunk * slice_size;
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Max alloc size: %lu, max global size: %lu", max_mem_alloc_size, max_global_mem_size);
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Num chunks: %d, chunk size: %lu, num slices per chunk: %u",
+ priv->num_chunks, chunk_size, priv->num_slices_per_chunk);
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Volume size: %lu, num slices: %u", volume_size, priv->num_slices);
+ priv->chunks = (cl_mem *) g_malloc (priv->num_chunks * sizeof (cl_mem));
+ if (!priv->chunks) {
+ g_set_error_literal (error, UFO_TASK_ERROR, UFO_TASK_ERROR_GET_REQUISITION,
+ "Error allocating volume chunks");
+ return;
+ }
+ priv->cl_regions = (cl_mem *) g_malloc (priv->num_chunks * sizeof (cl_mem));
+ if (!priv->cl_regions) {
+ g_set_error_literal (error, UFO_TASK_ERROR, UFO_TASK_ERROR_GET_REQUISITION,
+ "Error allocating volume chunks");
+ return;
+ }
+ for (i = 0; i < priv->num_chunks; i++) {
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Creating chunk %d with size %lu",
+ i, MIN (volume_size, (i + 1) * chunk_size) - i * chunk_size);
+ priv->chunks[i] = clCreateBuffer (priv->context,
+ CL_MEM_WRITE_ONLY,
+ MIN (volume_size, (i + 1) * chunk_size) - i * chunk_size,
+ NULL,
+ &cl_error);
+ UFO_RESOURCES_CHECK_CLERR (cl_error);
+ }
+ create_images (priv, in_req.dims[0], in_req.dims[1]);
+ create_regions[priv->compute_type] (priv, cmd_queue, region_start, region_step);
+ set_static_args[priv->compute_type] (task, requisition, priv->kernel);
+ if (priv->rest_kernel) {
+ set_static_args[priv->compute_type] (task, requisition, priv->rest_kernel);
+ }
+ }
+
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "requisition (x, y, z): %lu %lu %d", requisition->dims[0], requisition->dims[1], 1);
+}
+
+static guint
+ufo_general_backproject_task_get_num_inputs (UfoTask *task)
+{
+ return 1;
+}
+
+static guint
+ufo_general_backproject_task_get_num_dimensions (UfoTask *task,
+ guint input)
+{
+ g_return_val_if_fail (input == 0, 0);
+
+ return 2;
+}
+
+static UfoTaskMode
+ufo_general_backproject_task_get_mode (UfoTask *task)
+{
+ return UFO_TASK_MODE_REDUCTOR | UFO_TASK_MODE_GPU;
+}
+
+static gboolean
+ufo_general_backproject_task_process (UfoTask *task,
+ UfoBuffer **inputs,
+ UfoBuffer *output,
+ UfoRequisition *requisition)
+{
+ UfoGeneralBackprojectTaskPrivate *priv;
+ UfoRequisition in_req;
+ UfoGpuNode *node;
+ UfoProfiler *profiler;
+ guint i, index, ki;
+ guint count, burst, num_slices_current_chunk;
+ cl_kernel kernel;
+ cl_command_queue cmd_queue;
+ gdouble rot_angle;
+ cl_float f_tomo_angle[2];
+ cl_double d_tomo_angle[2];
+ cl_int iteration;
+ const gsize local_work_size[3] = {16, 8, 8};
+ gsize global_work_size[3];
+ gint real_size[4];
+
+ priv = UFO_GENERAL_BACKPROJECT_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);
+ ufo_buffer_get_requisition (inputs[0], &in_req);
+ g_object_get (task, "num_processed", &count, NULL);
+
+ if (count >= priv->num_projections / priv->burst * priv->burst) {
+ kernel = priv->rest_kernel;
+ burst = priv->num_projections % priv->burst;
+ index = (count - priv->num_projections / priv->burst * priv->burst) % burst;
+ } else {
+ kernel = priv->kernel;
+ burst = priv->burst;
+ index = count % burst;
+ }
+
+ global_work_size[0] = requisition->dims[0] % local_work_size[0] ?
+ NEXT_DIVISOR (requisition->dims[0], local_work_size[0]) :
+ requisition->dims[0];
+ global_work_size[1] = requisition->dims[1] % local_work_size[1] ?
+ NEXT_DIVISOR (requisition->dims[1], local_work_size[1]) :
+ requisition->dims[1];
+ global_work_size[2] = priv->num_slices_per_chunk % local_work_size[2] ?
+ NEXT_DIVISOR (priv->num_slices_per_chunk, local_work_size[2]) :
+ priv->num_slices_per_chunk;
+ real_size[0] = requisition->dims[0];
+ real_size[1] = requisition->dims[1];
+ real_size[3] = 0;
+ if (!count) {
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Global work size: %lu %lu %lu, local: %lu %lu %lu",
+ global_work_size[0], global_work_size[1], global_work_size[2],
+ local_work_size[0], local_work_size[1], local_work_size[2]);
+ }
+
+ /* Setup tomographic rotation angle dependent arguments */
+ ki = STATIC_ARG_OFFSET + burst;
+ rot_angle = ufo_scarray_get_double (priv->geometry->axis->angle->z, count);
+ if (priv->compute_type == CT_FLOAT) {
+ fill_sincos_cl_float (f_tomo_angle, rot_angle);
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, ki + index, sizeof (cl_float2), f_tomo_angle));
+ } else {
+ fill_sincos_cl_double (d_tomo_angle, rot_angle);
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, ki + index, sizeof (cl_double2), d_tomo_angle));
+ }
+ copy_to_image (cmd_queue, inputs[0], priv->projections[index], in_req.dims[0], in_req.dims[1]);
+
+ if (index + 1 == burst) {
+ profiler = ufo_task_node_get_profiler (UFO_TASK_NODE (task));
+ ki += index + 1;
+ iteration = (cl_int) (count + 1 - burst);
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, ki++, sizeof (cl_int), &iteration));
+ for (i = 0; i < priv->num_chunks; i++) {
+ /* The last chunk might be smaller */
+ num_slices_current_chunk = MIN (priv->num_slices, (i + 1) * priv->num_slices_per_chunk) - i * priv->num_slices_per_chunk;
+ real_size[2] = (gint) num_slices_current_chunk;
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, REAL_SIZE_ARG_INDEX, sizeof (cl_int3), real_size));
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, ki, sizeof (cl_mem), &priv->chunks[i]));
+ UFO_RESOURCES_CHECK_CLERR (clSetKernelArg (kernel, ki + 1, sizeof (cl_mem), &priv->cl_regions[i]));
+ ufo_profiler_call_blocking (profiler, cmd_queue, kernel, 3, global_work_size, local_work_size);
+ }
+ }
+
+ return TRUE;
+}
+
+static gboolean
+ufo_general_backproject_task_generate (UfoTask *task,
+ UfoBuffer *output,
+ UfoRequisition *requisition)
+{
+ UfoGeneralBackprojectTaskPrivate *priv;
+ UfoGpuNode *node;
+ cl_command_queue cmd_queue;
+ cl_mem out_mem;
+ guint count, chunk_index;
+ /* TODO: handle other data types */
+ size_t bpp;
+ size_t src_row_pitch, src_slice_pitch;
+ size_t src_origin[3] = {0, 0, 0};
+ size_t dst_origin[3] = {0, 0, 0};
+ size_t region[3] = {0, 0, 1};
+
+ priv = UFO_GENERAL_BACKPROJECT_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);
+ out_mem = ufo_buffer_get_device_array (output, cmd_queue);
+ chunk_index = priv->generated / priv->num_slices_per_chunk;
+ bpp = get_type_size (priv->store_type);
+ g_object_get (task, "num_processed", &count, NULL);
+
+ if (count != priv->num_projections || priv->generated >= priv->num_slices) {
+ /* Don't send volume if not enough projections came */
+ return FALSE;
+ }
+
+ src_row_pitch = requisition->dims[0] * bpp;
+ src_slice_pitch = src_row_pitch * requisition->dims[1];
+ src_origin[2] = priv->generated % priv->num_slices_per_chunk;
+ region[0] = src_row_pitch;
+ region[1] = requisition->dims[1];
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "Generating slice %u from chunk %u", priv->generated + 1, chunk_index);
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "src_origin: %lu %lu %lu", src_origin[0], src_origin[1], src_origin[2]);
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "region: %lu %lu %lu", region[0], region[1], region[2]);
+ g_log ("gbp", G_LOG_LEVEL_DEBUG, "row pitch %lu, slice pitch %lu", src_row_pitch, src_slice_pitch);
+
+ UFO_RESOURCES_CHECK_CLERR (clEnqueueCopyBufferRect (cmd_queue,
+ priv->chunks[chunk_index], out_mem,
+ src_origin, dst_origin, region,
+ src_row_pitch, src_slice_pitch,
+ src_row_pitch, 0,
+ 0, NULL, NULL));
+
+ priv->generated++;
+
+ return TRUE;
+}
+
+/*{{{ Setters and getters and properties initialization */
+static void
+ufo_general_backproject_task_set_property (GObject *object,
+ guint property_id,
+ const GValue *value,
+ GParamSpec *pspec)
+{
+ UfoGeneralBackprojectTaskPrivate *priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (object);
+
+ switch (property_id) {
+ case PROP_BURST:
+ priv->burst = g_value_get_uint (value);
+ break;
+ case PROP_PARAMETER:
+ priv->parameter = g_value_get_enum (value);
+ break;
+ case PROP_Z:
+ priv->z = g_value_get_double (value);
+ break;
+ case PROP_REGION:
+ ufo_scarray_get_value (priv->region, value);
+ break;
+ case PROP_REGION_X:
+ ufo_scarray_get_value (priv->region_x, value);
+ break;
+ case PROP_REGION_Y:
+ ufo_scarray_get_value (priv->region_y, value);
+ break;
+ case PROP_CENTER_POSITION_X:
+ ufo_scarray_get_value (priv->geometry->axis->position->x, value);
+ break;
+ case PROP_CENTER_POSITION_Z:
+ ufo_scarray_get_value (priv->geometry->axis->position->z, value);
+ break;
+ case PROP_SOURCE_POSITION_X:
+ ufo_scarray_get_value (priv->geometry->source_position->x, value);
+ break;
+ case PROP_SOURCE_POSITION_Y:
+ ufo_scarray_get_value (priv->geometry->source_position->y, value);
+ break;
+ case PROP_SOURCE_POSITION_Z:
+ ufo_scarray_get_value (priv->geometry->source_position->z, value);
+ break;
+ case PROP_DETECTOR_POSITION_X:
+ ufo_scarray_get_value (priv->geometry->detector->position->x, value);
+ break;
+ case PROP_DETECTOR_POSITION_Y:
+ ufo_scarray_get_value (priv->geometry->detector->position->y, value);
+ break;
+ case PROP_DETECTOR_POSITION_Z:
+ ufo_scarray_get_value (priv->geometry->detector->position->z, value);
+ break;
+ case PROP_DETECTOR_ANGLE_X:
+ ufo_scarray_get_value (priv->geometry->detector->angle->x, value);
+ break;
+ case PROP_DETECTOR_ANGLE_Y:
+ ufo_scarray_get_value (priv->geometry->detector->angle->y, value);
+ break;
+ case PROP_DETECTOR_ANGLE_Z:
+ ufo_scarray_get_value (priv->geometry->detector->angle->z, value);
+ break;
+ case PROP_AXIS_ANGLE_X:
+ ufo_scarray_get_value (priv->geometry->axis->angle->x, value);
+ break;
+ case PROP_AXIS_ANGLE_Y:
+ ufo_scarray_get_value (priv->geometry->axis->angle->y, value);
+ break;
+ case PROP_AXIS_ANGLE_Z:
+ ufo_scarray_get_value (priv->geometry->axis->angle->z, value);
+ break;
+ case PROP_VOLUME_ANGLE_X:
+ ufo_scarray_get_value (priv->geometry->volume_angle->x, value);
+ break;
+ case PROP_VOLUME_ANGLE_Y:
+ ufo_scarray_get_value (priv->geometry->volume_angle->y, value);
+ break;
+ case PROP_VOLUME_ANGLE_Z:
+ ufo_scarray_get_value (priv->geometry->volume_angle->z, value);
+ break;
+ case PROP_NUM_PROJECTIONS:
+ priv->num_projections = g_value_get_uint (value);
+ break;
+ case PROP_COMPUTE_TYPE:
+ priv->compute_type = g_value_get_enum (value);
+ break;
+ case PROP_RESULT_TYPE:
+ priv->result_type = g_value_get_enum (value);
+ break;
+ case PROP_STORE_TYPE:
+ priv->store_type = g_value_get_enum (value);
+ break;
+ case PROP_OVERALL_ANGLE:
+ priv->overall_angle = g_value_get_double (value);
+ break;
+ case PROP_ADDRESSING_MODE:
+ priv->addressing_mode = g_value_get_enum (value);
+ break;
+ case PROP_GRAY_MAP_MIN:
+ priv->gray_map_min = g_value_get_double (value);
+ break;
+ case PROP_GRAY_MAP_MAX:
+ priv->gray_map_max = g_value_get_double (value);
+ break;
+ default:
+ G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
+ break;
+ }
+}
+
+static void
+ufo_general_backproject_task_get_property (GObject *object,
+ guint property_id,
+ GValue *value,
+ GParamSpec *pspec)
+{
+ UfoGeneralBackprojectTaskPrivate *priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (object);
+
+ switch (property_id) {
+ case PROP_BURST:
+ g_value_set_uint (value, priv->burst);
+ break;
+ case PROP_PARAMETER:
+ g_value_set_enum (value, priv->parameter);
+ break;
+ case PROP_Z:
+ g_value_set_double (value, priv->z);
+ break;
+ case PROP_REGION:
+ ufo_scarray_set_value (priv->region, value);
+ break;
+ case PROP_REGION_X:
+ ufo_scarray_set_value (priv->region_x, value);
+ break;
+ case PROP_REGION_Y:
+ ufo_scarray_set_value (priv->region_y, value);
+ break;
+ case PROP_CENTER_POSITION_X:
+ ufo_scarray_set_value (priv->geometry->axis->position->x, value);
+ break;
+ case PROP_CENTER_POSITION_Z:
+ ufo_scarray_set_value (priv->geometry->axis->position->z, value);
+ break;
+ case PROP_SOURCE_POSITION_X:
+ ufo_scarray_set_value (priv->geometry->source_position->x, value);
+ break;
+ case PROP_SOURCE_POSITION_Y:
+ ufo_scarray_set_value (priv->geometry->source_position->y, value);
+ break;
+ case PROP_SOURCE_POSITION_Z:
+ ufo_scarray_set_value (priv->geometry->source_position->z, value);
+ break;
+ case PROP_DETECTOR_POSITION_X:
+ ufo_scarray_set_value (priv->geometry->detector->position->x, value);
+ break;
+ case PROP_DETECTOR_POSITION_Y:
+ ufo_scarray_set_value (priv->geometry->detector->position->y, value);
+ break;
+ case PROP_DETECTOR_POSITION_Z:
+ ufo_scarray_set_value (priv->geometry->detector->position->z, value);
+ break;
+ case PROP_DETECTOR_ANGLE_X:
+ ufo_scarray_set_value (priv->geometry->detector->angle->x, value);
+ break;
+ case PROP_DETECTOR_ANGLE_Y:
+ ufo_scarray_set_value (priv->geometry->detector->angle->y, value);
+ break;
+ case PROP_DETECTOR_ANGLE_Z:
+ ufo_scarray_set_value (priv->geometry->detector->angle->z, value);
+ break;
+ case PROP_AXIS_ANGLE_X:
+ ufo_scarray_set_value (priv->geometry->axis->angle->x, value);
+ break;
+ case PROP_AXIS_ANGLE_Y:
+ ufo_scarray_set_value (priv->geometry->axis->angle->y, value);
+ break;
+ case PROP_AXIS_ANGLE_Z:
+ ufo_scarray_set_value (priv->geometry->axis->angle->z, value);
+ break;
+ case PROP_VOLUME_ANGLE_X:
+ ufo_scarray_set_value (priv->geometry->volume_angle->x, value);
+ break;
+ case PROP_VOLUME_ANGLE_Y:
+ ufo_scarray_set_value (priv->geometry->volume_angle->y, value);
+ break;
+ case PROP_VOLUME_ANGLE_Z:
+ ufo_scarray_set_value (priv->geometry->volume_angle->z, value);
+ break;
+ case PROP_NUM_PROJECTIONS:
+ g_value_set_uint (value, priv->num_projections);
+ break;
+ case PROP_COMPUTE_TYPE:
+ g_value_set_enum (value, priv->compute_type);
+ break;
+ case PROP_RESULT_TYPE:
+ g_value_set_enum (value, priv->result_type);
+ break;
+ case PROP_STORE_TYPE:
+ g_value_set_enum (value, priv->store_type);
+ break;
+ case PROP_OVERALL_ANGLE:
+ g_value_set_double (value, priv->overall_angle);
+ break;
+ case PROP_GRAY_MAP_MIN:
+ g_value_set_double (value, priv->gray_map_min);
+ break;
+ case PROP_GRAY_MAP_MAX:
+ g_value_set_double (value, priv->gray_map_max);
+ break;
+ case PROP_ADDRESSING_MODE:
+ g_value_set_enum (value, priv->addressing_mode);
+ break;
+ default:
+ G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
+ break;
+ }
+}
+
+static void
+ufo_general_backproject_task_finalize (GObject *object)
+{
+ guint i;
+ UfoGeneralBackprojectTaskPrivate *priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE (object);
+ g_object_unref (priv->resources);
+ priv->resources = NULL;
+
+ ufo_scarray_free (priv->region);
+ ufo_scarray_free (priv->region_x);
+ ufo_scarray_free (priv->region_y);
+ ufo_ctgeometry_free (priv->geometry);
+ g_hash_table_destroy (priv->node_props_table);
+
+ if (priv->projections) {
+ for (i = 0; i < priv->burst; i++) {
+ if (priv->projections[i] != NULL) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (priv->projections[i]));
+ priv->projections[i] = NULL;
+ }
+ }
+ g_free (priv->projections);
+ priv->projections = NULL;
+ }
+
+
+ if (priv->chunks) {
+ for (i = 0; i < priv->num_chunks; i++) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (priv->chunks[i]));
+ }
+ g_free (priv->chunks);
+ priv->chunks = NULL;
+ }
+
+ if (priv->cl_regions) {
+ for (i = 0; i < priv->num_chunks; i++) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (priv->cl_regions[i]));
+ }
+ g_free (priv->cl_regions);
+ priv->cl_regions = NULL;
+ }
+
+ if (priv->vector_arguments) {
+ for (i = 0; i < NUM_VECTOR_ARGUMENTS; i++) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseMemObject (priv->vector_arguments[i]));
+ }
+ g_free (priv->vector_arguments);
+ priv->vector_arguments = NULL;
+ }
+
+ if (priv->kernel) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->kernel));
+ priv->kernel = NULL;
+ }
+ if (priv->rest_kernel) {
+ UFO_RESOURCES_CHECK_CLERR (clReleaseKernel (priv->rest_kernel));
+ priv->rest_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;
+ }
+
+ G_OBJECT_CLASS (ufo_general_backproject_task_parent_class)->finalize (object);
+}
+
+static void
+ufo_task_interface_init (UfoTaskIface *iface)
+{
+ iface->setup = ufo_general_backproject_task_setup;
+ iface->get_num_inputs = ufo_general_backproject_task_get_num_inputs;
+ iface->get_num_dimensions = ufo_general_backproject_task_get_num_dimensions;
+ iface->get_mode = ufo_general_backproject_task_get_mode;
+ iface->get_requisition = ufo_general_backproject_task_get_requisition;
+ iface->process = ufo_general_backproject_task_process;
+ iface->generate = ufo_general_backproject_task_generate;
+}
+
+static void
+ufo_general_backproject_task_class_init (UfoGeneralBackprojectTaskClass *klass)
+{
+ GObjectClass *oclass = G_OBJECT_CLASS (klass);
+
+ oclass->set_property = ufo_general_backproject_task_set_property;
+ oclass->get_property = ufo_general_backproject_task_get_property;
+ oclass->finalize = ufo_general_backproject_task_finalize;
+
+ GParamSpec *region_vals = g_param_spec_int ("region-values",
+ "Region values",
+ "Elements in regions",
+ G_MININT,
+ G_MAXINT,
+ (gint) 0,
+ G_PARAM_READWRITE);
+
+ 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_BURST] =
+ g_param_spec_uint ("burst",
+ "Number of projections processed per one kernel invocation",
+ "Number of projections processed per one kernel invocation",
+ 0, 128, 0,
+ G_PARAM_READWRITE);
+
+ properties[PROP_PARAMETER] =
+ g_param_spec_enum ("parameter",
+ "Which parameter will be varied along the z-axis",
+ "Which parameter will be varied along the z-axis",
+ g_enum_register_static ("GBPParameter", parameter_values),
+ UFO_UNI_RECO_PARAMETER_Z,
+ G_PARAM_READWRITE);
+
+ properties[PROP_Z] =
+ g_param_spec_double ("z",
+ "Z coordinate of the reconstructed slice",
+ "Z coordinate of the reconstructed slice",
+ -G_MAXDOUBLE, G_MAXDOUBLE, 0.0,
+ 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)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_REGION_X] =
+ g_param_spec_value_array ("x-region",
+ "X region for reconstruction (horizontal axis) as (from, to, step)",
+ "X region for reconstruction (horizontal axis) as (from, to, step)",
+ region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_REGION_Y] =
+ g_param_spec_value_array ("y-region",
+ "Y region for reconstruction (beam direction axis) as (from, to, step)",
+ "Y region for reconstruction (beam direction axis) as (from, to, step)",
+ region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_CENTER_POSITION_X] =
+ g_param_spec_value_array ("center-position-x",
+ "Global x center (horizontal in a projection) of the volume with respect to projections",
+ "Global x center (horizontal in a projection) of the volume with respect to projections",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_CENTER_POSITION_Z] =
+ g_param_spec_value_array ("center-position-z",
+ "Global z center (vertical in a projection) of the volume with respect to projections",
+ "Global z center (vertical in a projection) of the volume with respect to projections",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_SOURCE_POSITION_X] =
+ g_param_spec_value_array ("source-position-x",
+ "X source position (horizontal) in global coordinates [pixels]",
+ "X source position (horizontal) in global coordinates [pixels]",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_SOURCE_POSITION_Y] =
+ g_param_spec_value_array ("source-position-y",
+ "Y source position (beam direction) in global coordinates [pixels]",
+ "Y source position (beam direction) in global coordinates [pixels]",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_SOURCE_POSITION_Z] =
+ g_param_spec_value_array ("source-position-z",
+ "Z source position (vertical) in global coordinates [pixels]",
+ "Z source position (vertical) in global coordinates [pixels]",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_POSITION_X] =
+ g_param_spec_value_array ("detector-position-x",
+ "X detector position (horizontal) in global coordinates [pixels]",
+ "X detector position (horizontal) in global coordinates [pixels]",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_POSITION_Y] =
+ g_param_spec_value_array ("detector-position-y",
+ "Y detector position (along beam direction) in global coordinates [pixels]",
+ "Y detector position (along beam direction) in global coordinates [pixels]",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_POSITION_Z] =
+ g_param_spec_value_array ("detector-position-z",
+ "Z detector position (vertical) in global coordinates [pixels]",
+ "Z detector position (vertical) in global coordinates [pixels]",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_ANGLE_X] =
+ g_param_spec_value_array("detector-angle-x",
+ "Detector rotation around the x axis [rad] (horizontal)",
+ "Detector rotation around the x axis [rad] (horizontal)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_ANGLE_Y] =
+ g_param_spec_value_array("detector-angle-y",
+ "Detector rotation around the y axis [rad] (along beam direction)",
+ "Detector rotation around the y axis [rad] (balong eam direction)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_DETECTOR_ANGLE_Z] =
+ g_param_spec_value_array("detector-angle-z",
+ "Detector rotation around the z axis [rad] (vertical)",
+ "Detector rotation around the z axis [rad] (vertical)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_AXIS_ANGLE_X] =
+ g_param_spec_value_array("axis-angle-x",
+ "Rotation axis rotation around the x axis [rad] (laminographic angle, 0 = tomography)",
+ "Rotation axis rotation around the x axis [rad] (laminographic angle, 0 = tomography)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_AXIS_ANGLE_Y] =
+ g_param_spec_value_array("axis-angle-y",
+ "Rotation axis rotation around the y axis [rad] (along beam direction)",
+ "Rotation axis rotation around the y axis [rad] (along beam direction)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_AXIS_ANGLE_Z] =
+ g_param_spec_value_array("axis-angle-z",
+ "Rotation axis rotation around the z axis [rad] (vertical)",
+ "Rotation axis rotation around the z axis [rad] (vertical)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_VOLUME_ANGLE_X] =
+ g_param_spec_value_array("volume-angle-x",
+ "Volume rotation around the x axis [rad] (horizontal)",
+ "Volume rotation around the x axis [rad] (horizontal)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_VOLUME_ANGLE_Y] =
+ g_param_spec_value_array("volume-angle-y",
+ "Volume rotation around the y axis [rad] (along beam direction)",
+ "Volume rotation around the y axis [rad] (along beam direction)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_VOLUME_ANGLE_Z] =
+ g_param_spec_value_array("volume-angle-z",
+ "Volume rotation around the z axis [rad] (vertical)",
+ "Volume rotation around the z axis [rad] (vertical)",
+ double_region_vals,
+ G_PARAM_READWRITE);
+
+ properties[PROP_COMPUTE_TYPE] =
+ g_param_spec_enum ("compute-type",
+ "Data type for performing kernel math operations",
+ "Data type for performing kernel math operations "
+ "(\"half\", \"float\", \"double\")",
+ g_enum_register_static ("compute-type", compute_type_values),
+ FT_FLOAT,
+ G_PARAM_READWRITE);
+
+ properties[PROP_ADDRESSING_MODE] =
+ g_param_spec_enum ("addressing-mode",
+ "Outlier treatment (\"none\", \"clamp\", \"clamp_to_edge\", \"repeat\")",
+ "Outlier treatment (\"none\", \"clamp\", \"clamp_to_edge\", \"repeat\")",
+ g_enum_register_static ("bp_addressing_mode", addressing_values),
+ CL_ADDRESS_CLAMP,
+ G_PARAM_READWRITE);
+
+ properties[PROP_RESULT_TYPE] =
+ g_param_spec_enum ("result-type",
+ "Data type for storing the intermediate gray value for a voxel from various rotation angles",
+ "Data type for storing the intermediate gray value for a voxel from various rotation angles "
+ "(\"half\", \"float\", \"double\")",
+ g_enum_register_static ("result-type", ft_values),
+ FT_FLOAT,
+ G_PARAM_READWRITE);
+
+ properties[PROP_STORE_TYPE] =
+ g_param_spec_enum ("store-type",
+ "Data type of the output volume",
+ "Data type of the output volume "
+ "(\"half\", \"float\", \"double\", \"uchar\", \"ushort\", \"uint\")",
+ g_enum_register_static ("store-type", st_values),
+ ST_FLOAT,
+ G_PARAM_READWRITE);
+
+ properties[PROP_OVERALL_ANGLE] =
+ g_param_spec_double ("overall-angle",
+ "Angle covered by all projections [rad]",
+ "Angle covered by all projections [rad] (can be negative for negative steps "
+ "in case only num-projections is specified",
+ -G_MAXDOUBLE, G_MAXDOUBLE, 2 * G_PI,
+ G_PARAM_READWRITE);
+
+ properties[PROP_GRAY_MAP_MIN] =
+ g_param_spec_double ("gray-map-min",
+ "Gray valye which maps to 0 in case of integer store type",
+ "Gray valye which maps to 0 in case of integer store type",
+ -G_MAXDOUBLE, G_MAXDOUBLE, 0,
+ G_PARAM_READWRITE);
+
+ properties[PROP_GRAY_MAP_MAX] =
+ g_param_spec_double ("gray-map-max",
+ "Gray valye which maps to maximum of the chosen integer type in case of integer store type",
+ "Gray valye which maps to maximum of the chosen integer type in case of integer store type",
+ -G_MAXDOUBLE, G_MAXDOUBLE, 0,
+ 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);
+
+ for (guint i = PROP_0 + 1; i < N_PROPERTIES; i++)
+ g_object_class_install_property (oclass, i, properties[i]);
+
+ g_type_class_add_private (oclass, sizeof(UfoGeneralBackprojectTaskPrivate));
+}
+
+static void
+ufo_general_backproject_task_init(UfoGeneralBackprojectTask *self)
+{
+ self->priv = UFO_GENERAL_BACKPROJECT_TASK_GET_PRIVATE(self);
+
+ /* Scalars */
+ self->priv->burst = 0;
+ self->priv->parameter = UFO_UNI_RECO_PARAMETER_Z;
+ self->priv->z = 0.0;
+ self->priv->num_projections = 0;
+ self->priv->compute_type = CT_FLOAT;
+ self->priv->result_type = FT_FLOAT;
+ self->priv->store_type = ST_FLOAT;
+ self->priv->overall_angle = 2 * G_PI;
+ self->priv->addressing_mode = CL_ADDRESS_CLAMP;
+ self->priv->gray_map_min = 0.0;
+ self->priv->gray_map_max = 0.0;
+
+ /* Value arrays */
+ self->priv->region = ufo_scarray_new (3, G_TYPE_DOUBLE, NULL);
+ self->priv->region_x = ufo_scarray_new (3, G_TYPE_INT, NULL);
+ self->priv->region_y = ufo_scarray_new (3, G_TYPE_INT, NULL);
+ self->priv->geometry = ufo_ctgeometry_new ();
+
+ /* Private */
+ self->priv->vectorized = FALSE;
+ self->priv->num_slices = 0;
+ self->priv->num_slices_per_chunk = 0;
+ self->priv->generated = 0;
+}
+/*}}}*/
diff --git a/src/ufo-general-backproject-task.h b/src/ufo-general-backproject-task.h
new file mode 100644
index 0000000..42a8a71
--- /dev/null
+++ b/src/ufo-general-backproject-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_GENERAL_BACKPROJECT_TASK_H
+#define __UFO_GENERAL_BACKPROJECT_TASK_H
+
+#include <ufo/ufo.h>
+
+G_BEGIN_DECLS
+
+#define UFO_TYPE_GENERAL_BACKPROJECT_TASK (ufo_general_backproject_task_get_type())
+#define UFO_GENERAL_BACKPROJECT_TASK(obj) (G_TYPE_CHECK_INSTANCE_CAST((obj), UFO_TYPE_GENERAL_BACKPROJECT_TASK, UfoGeneralBackprojectTask))
+#define UFO_IS_GENERAL_BACKPROJECT_TASK(obj) (G_TYPE_CHECK_INSTANCE_TYPE((obj), UFO_TYPE_GENERAL_BACKPROJECT_TASK))
+#define UFO_GENERAL_BACKPROJECT_TASK_CLASS(klass) (G_TYPE_CHECK_CLASS_CAST((klass), UFO_TYPE_GENERAL_BACKPROJECT_TASK, UfoGeneralBackprojectTaskClass))
+#define UFO_IS_GENERAL_BACKPROJECT_TASK_CLASS(klass) (G_TYPE_CHECK_CLASS_TYPE((klass), UFO_TYPE_GENERAL_BACKPROJECT_TASK))
+#define UFO_GENERAL_BACKPROJECT_TASK_GET_CLASS(obj) (G_TYPE_INSTANCE_GET_CLASS((obj), UFO_TYPE_GENERAL_BACKPROJECT_TASK, UfoGeneralBackprojectTaskClass))
+
+typedef struct _UfoGeneralBackprojectTask UfoGeneralBackprojectTask;
+typedef struct _UfoGeneralBackprojectTaskClass UfoGeneralBackprojectTaskClass;
+typedef struct _UfoGeneralBackprojectTaskPrivate UfoGeneralBackprojectTaskPrivate;
+
+struct _UfoGeneralBackprojectTask {
+ UfoTaskNode parent_instance;
+
+ UfoGeneralBackprojectTaskPrivate *priv;
+};
+
+struct _UfoGeneralBackprojectTaskClass {
+ UfoTaskNodeClass parent_class;
+};
+
+UfoNode *ufo_general_backproject_task_new (void);
+GType ufo_general_backproject_task_get_type (void);
+
+G_END_DECLS
+
+#endif
generated by cgit v1.2.1 (git 2.18.0) at 2024-05-20 15:10:05 +0000