Flat-field correction (without detector plane handling)

2 files changed, 53 insertions, 0 deletions

diff --git a/src/kernels/meson.build b/src/kernels/meson.build
index 9803fb9..6472953 100644
--- a/src/kernels/meson.build
+++ b/src/kernels/meson.build
@@ -1,4 +1,5 @@
 kernel_files = [
+    'roof-ffc.cl'
 ]
 
 install_data(kernel_files,
diff --git a/src/kernels/roof-ffc.cl b/src/kernels/roof-ffc.cl
new file mode 100644
index 0000000..8752990
--- /dev/null
+++ b/src/kernels/roof-ffc.cl
@@ -0,0 +1,52 @@
+/*
+ * 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/>.
+ */
+
+kernel void
+flat_correct (global float *corrected,
+              global float *data,
+              global const float *dark,
+              global const float *flat,
+              const int sinogram_input,
+              const int absorptivity,
+              const int fix_abnormal,
+              const float dark_scale,
+              const float flat_scale)
+{
+    const float eps = 1E-5;
+    const int gid = get_global_id(1) * get_global_size(0) + get_global_id(0);
+    const int corr_idx = sinogram_input ? get_global_id(0) : gid;
+    const float cdark = dark[corr_idx] * dark_scale;
+    const float cflat = flat[corr_idx] * flat_scale;
+    float result;
+
+    float nom = data[gid] - cdark;
+    float denom = cflat - cdark;
+
+    if (nom < eps) nom = eps;
+//    if (denom < eps) denom = eps;
+    if (denom < eps) denom = eps;
+    
+    result = -log(nom / denom);
+
+    if (fix_abnormal && (isnan (result) || isinf (result))) {
+        result = 0.0f;
+    }
+
+    corrected[gid] = result;
+}