1 files changed, 20 insertions, 29 deletions

diff --git a/src/kernels/nlm.cl b/src/kernels/nlm.cl
index 8b1bcb8..f904ac7 100644
--- a/src/kernels/nlm.cl
+++ b/src/kernels/nlm.cl
@@ -17,35 +17,34 @@
  * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#define flatten(x,y,r,w) ((y-r)*w + (x-r))
-
-/* Compute the distance of two neighbourhood vectors _starting_ from index i
-   and j and edge length radius */
 float
-dist (global float *input,
-      int i,
-      int j,
+dist (read_only image2d_t input,
+      sampler_t sampler,
+      float2 p,
+      float2 q,
       int radius,
-      int image_width)
+      int width,
+      int height)
 {
     float dist = 0.0f, tmp;
     float wsize = (2.0f * radius + 1.0f);
     wsize *= wsize;
 
-    const int nb_width = 2 * radius + 1;
-    const int stride = image_width - nb_width;
-    for (int k = 0; k < nb_width; k++, i += stride, j += stride) {
-        for (int l = 0; l < nb_width; l++, i++, j++) {
-            tmp = input[i] - input[j];
+    for (int i = -radius; i < radius + 1; i++) {
+        for (int j = -radius; j < radius + 1; j++) {
+            tmp = read_imagef (input, sampler, (float2) ((p.x + i) / width, (p.y + j) / height)).x -
+                    read_imagef (input, sampler, (float2) ((q.x + i) / width, (q.y + j) / height)).x;
             dist += tmp * tmp;
         }
     }
+
     return dist / wsize;
 }
 
 kernel void
-nlm_noise_reduction (global float *input,
+nlm_noise_reduction (read_only image2d_t input,
                      global float *output,
+                     sampler_t sampler,
                      const int search_radius,
                      const int patch_radius,
                      const float sigma)
@@ -55,25 +54,17 @@ nlm_noise_reduction (global float *input,
     const int width = get_global_size (0);
     const int height = get_global_size (1);
     const float sigma_2 = sigma * sigma;
+    float d, weight;
 
     float total_weight = 0.0f;
     float pixel_value = 0.0f;
 
-    /* 
-     * Compute the upper left (sx,sy) and lower right (tx, ty) corner points of
-     * our search window.
-     */
-    int r = min (patch_radius, min(width - 1 - x, min (height - 1 - y, min (x, y))));
-    int sx = max (x - search_radius, r);
-    int sy = max (y - search_radius, r);
-    int tx = min (x + search_radius, width - 1 - r);
-    int ty = min (y + search_radius, height - 1 - r);
-
-    for (int i = sx; i < tx; i++) {
-        for (int j = sy; j < ty; j++) {
-            float d = dist (input, flatten(x, y, r, width), flatten (i,j,r,width), r, width);
-            float weight = exp (- sigma_2 * d);
-            pixel_value += weight * input[j * width + i];
+    for (int i = x - search_radius; i < x + search_radius; i++) {
+        for (int j = y - search_radius; j < y + search_radius; j++) {
+            d = dist (input, sampler, (float2) (x + 0.5f, y + 0.5f), (float2) (i + 0.5f, j + 0.5f),
+                      patch_radius, width, height);
+            weight = exp (- sigma_2 * d);
+            pixel_value += weight * read_imagef (input, sampler, (float2) ((i + 0.5f) / width, (j + 0.5f) / height)).x;
             total_weight += weight;
         }
     }