First commit - moving from LuaForge to SourceForge

1 files changed, 1512 insertions, 0 deletions

diff --git a/src/process/im_convolve.cpp b/src/process/im_convolve.cpp
new file mode 100644
index 0000000..bca2dcd
--- /dev/null
+++ b/src/process/im_convolve.cpp
@@ -0,0 +1,1512 @@
+/** \file
+ * \brief Convolution Operations
+ *
+ * See Copyright Notice in im_lib.h
+ * $Id: im_convolve.cpp,v 1.1 2008/10/17 06:16:33 scuri Exp $
+ */
+
+
+#include <im.h>
+#include <im_util.h>
+#include <im_counter.h>
+#include <im_complex.h>
+#include <im_math_op.h>
+#include <im_image.h>
+#include <im_kernel.h>
+
+#include "im_process_loc.h"
+#include "im_process_pon.h"
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <memory.h>
+#include <string.h>
+#include <math.h>
+
+/* Rotating Kernels
+3x3
+  6 7 8   7 8 5
+  3 4 5   6 4 2
+  0 1 2   3 0 1
+
+5x5
+  20 21 22 23 24   22 23 24 19 14
+  15 16 17 18 19   21 17 18 13  9
+  10 11 12 13 14   20 16 12  8  4
+   5  6  7  8  9   15 11  6  7  3
+   0  1  2  3  4   10  5  0  1  2
+
+7x7
+  42 43 44 45 46 47 48     45 46 47 48 41 34 27
+  35 36 37 38 39 40 41     44 38 39 40 33 26 20
+  28 29 30 31 32 33 34     43 37 31 32 25 19 13
+  21 22 23 24 25 26 27     42 36 30 24 18 12  6
+  14 15 16 17 18 19 20     35 29 23 16 17 11  5
+   7  8  9 10 11 12 13     28 22 15  8  9 10  4
+   0  1  2  3  4  5  6     21 14  7  0  1  2  3
+
+ TO DO: a generic odd rotation function...
+*/
+
+template <class KT> 
+static void iRotateKernel(KT* kernel_map, int kernel_size)
+{
+  KT temp;
+
+  switch (kernel_size)
+  {
+  case 3:
+    {
+      temp = kernel_map[0];
+      kernel_map[0] = kernel_map[3];
+      kernel_map[3] = kernel_map[6];
+      kernel_map[6] = kernel_map[7];
+      kernel_map[7] = kernel_map[8];
+      kernel_map[8] = kernel_map[5];
+      kernel_map[5] = kernel_map[2];
+      kernel_map[2] = kernel_map[1];
+      kernel_map[1] = temp;
+    }
+    break;
+  case 5:
+    {
+      temp = kernel_map[0];
+      kernel_map[0] = kernel_map[10];
+      kernel_map[10] = kernel_map[20];
+      kernel_map[20] = kernel_map[22];
+      kernel_map[22] = kernel_map[24];
+      kernel_map[24] = kernel_map[14];
+      kernel_map[14] = kernel_map[4];
+      kernel_map[4] = kernel_map[2];
+      kernel_map[2] = temp;
+
+      temp = kernel_map[5];
+      kernel_map[5] = kernel_map[15];
+      kernel_map[15] = kernel_map[21];
+      kernel_map[21] = kernel_map[23];
+      kernel_map[23] = kernel_map[19];
+      kernel_map[19] = kernel_map[9];
+      kernel_map[9] = kernel_map[3];
+      kernel_map[3] = kernel_map[1];
+      kernel_map[1] = temp;
+
+      temp = kernel_map[6];
+      kernel_map[6] = kernel_map[11];
+      kernel_map[11] = kernel_map[16];
+      kernel_map[16] = kernel_map[17];
+      kernel_map[17] = kernel_map[18];
+      kernel_map[18] = kernel_map[13];
+      kernel_map[13] = kernel_map[8];
+      kernel_map[8] = kernel_map[7];
+      kernel_map[7] = temp;
+    }
+    break;
+  case 7:
+    {
+      temp = kernel_map[2];
+      kernel_map[2] = kernel_map[7];
+      kernel_map[7] = kernel_map[28];
+      kernel_map[28] = kernel_map[43];
+      kernel_map[43] = kernel_map[46];
+      kernel_map[46] = kernel_map[41];
+      kernel_map[41] = kernel_map[20];
+      kernel_map[20] = kernel_map[5];
+      kernel_map[5] = temp;
+
+      temp = kernel_map[1];
+      kernel_map[1] = kernel_map[14];
+      kernel_map[14] = kernel_map[35];
+      kernel_map[35] = kernel_map[44];
+      kernel_map[44] = kernel_map[47];
+      kernel_map[47] = kernel_map[34];
+      kernel_map[34] = kernel_map[13];
+      kernel_map[13] = kernel_map[4];
+      kernel_map[4] = temp;
+
+      temp = kernel_map[0];
+      kernel_map[0] = kernel_map[21];
+      kernel_map[21] = kernel_map[42];
+      kernel_map[42] = kernel_map[45];
+      kernel_map[45] = kernel_map[48];
+      kernel_map[48] = kernel_map[27];
+      kernel_map[27] = kernel_map[6];
+      kernel_map[6] = kernel_map[3];
+      kernel_map[3] = temp;
+
+      temp = kernel_map[9];
+      kernel_map[9] = kernel_map[15];
+      kernel_map[15] = kernel_map[29];
+      kernel_map[29] = kernel_map[37];
+      kernel_map[37] = kernel_map[39];
+      kernel_map[39] = kernel_map[33];
+      kernel_map[33] = kernel_map[19];
+      kernel_map[19] = kernel_map[11];
+      kernel_map[11] = temp;
+
+      temp = kernel_map[8];
+      kernel_map[8] = kernel_map[22];
+      kernel_map[22] = kernel_map[36];
+      kernel_map[36] = kernel_map[38];
+      kernel_map[38] = kernel_map[40];
+      kernel_map[40] = kernel_map[26];
+      kernel_map[26] = kernel_map[12];
+      kernel_map[12] = kernel_map[10];
+      kernel_map[10] = temp;
+
+      temp = kernel_map[16];
+      kernel_map[16] = kernel_map[23];
+      kernel_map[23] = kernel_map[30];
+      kernel_map[30] = kernel_map[31];
+      kernel_map[31] = kernel_map[32];
+      kernel_map[32] = kernel_map[25];
+      kernel_map[25] = kernel_map[18];
+      kernel_map[18] = kernel_map[17];
+      kernel_map[17] = temp;
+    }
+    break;
+  }
+}
+
+void imProcessRotateKernel(imImage* kernel)
+{
+  if (kernel->data_type == IM_INT)
+    iRotateKernel((int*)kernel->data[0], kernel->width);
+  else
+    iRotateKernel((float*)kernel->data[0], kernel->width);
+}
+
+template <class T, class KT, class CT> 
+static int DoCompassConvolve(T* map, T* new_map, int width, int height, KT* orig_kernel_map, int kernel_size, int counter, CT)
+{
+  CT value;
+  KT total, *kernel_line, kvalue;
+  int offset, new_offset, i, j, x, y, kcount;
+
+  // duplicate the kernel data so we can rotate it
+  kcount = kernel_size*kernel_size;
+  KT* kernel_map = (KT*)malloc(kcount*sizeof(KT));
+
+  int ks2 = kernel_size/2;
+
+  total = 0;
+  for(j = 0; j < kcount; j++) 
+  {
+    kvalue = orig_kernel_map[j];
+    kernel_map[j] = kvalue;
+    total += kvalue;
+  }
+
+  if (total == 0)
+    total = 1;
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      CT max_value = 0;
+
+      for(int k = 0; k < 8; k++) // Rotate 8 times
+      {
+        value = 0;
+      
+        for(y = -ks2; y <= ks2; y++)
+        {
+          kernel_line = kernel_map + (y+ks2)*kernel_size;
+
+          if (j + y < 0)             // pass the bottom border
+            offset = -(y + j + 1) * width;
+          else if (j + y >= height)  // pass the top border
+            offset = (2*height - 1 - (j + y)) * width;
+          else
+            offset = (j + y) * width;
+
+          for(x = -ks2; x <= ks2; x++)
+          {
+            if (i + x < 0)            // pass the left border
+              value += kernel_line[x+ks2] * map[offset - (i + x + 1)];
+            else if (i + x >= width)  // pass the right border
+              value += kernel_line[x+ks2] * map[offset + 2*width - 1 - (i + x)];
+            else if (offset != -1)
+              value += kernel_line[x+ks2] * map[offset + (i + x)];
+          }
+        }
+
+        if (abs_op(value) > max_value)
+          max_value = abs_op(value);
+
+        iRotateKernel(kernel_map, kernel_size);
+      }  
+
+      max_value /= total;
+
+      int size_of = sizeof(imbyte);
+      if (sizeof(T) == size_of)
+        new_map[new_offset + i] = (T)IM_BYTECROP(max_value);
+      else
+        new_map[new_offset + i] = (T)max_value;
+    }    
+
+    if (!imCounterInc(counter))
+    {
+      free(kernel_map);
+      return 0;
+    }
+  }
+
+  free(kernel_map);
+  return 1;
+}
+
+int imProcessCompassConvolve(const imImage* src_image, imImage* dst_image, imImage *kernel)
+{
+  int ret = 0;
+
+  int counter = imCounterBegin("Compass Convolution");
+  const char* msg = (const char*)imImageGetAttribute(kernel, "Description", NULL, NULL);
+  if (!msg) msg = "Filtering...";
+  imCounterTotal(counter, src_image->depth*src_image->height, msg);
+
+  for (int i = 0; i < src_image->depth; i++)
+  {
+    switch(src_image->data_type)
+    {
+    case IM_BYTE:
+      if (kernel->data_type == IM_INT)
+        ret = DoCompassConvolve((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, counter, (int)0);
+      else
+        ret = DoCompassConvolve((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, counter, (float)0);
+      break;                                                                                
+    case IM_USHORT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoCompassConvolve((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, counter, (int)0);
+      else
+        ret = DoCompassConvolve((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, counter, (float)0);
+      break;                                                                                
+    case IM_INT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoCompassConvolve((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, counter, (int)0);
+      else
+        ret = DoCompassConvolve((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, counter, (float)0);
+      break;                                                                                
+    case IM_FLOAT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoCompassConvolve((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, counter, (float)0);
+      else
+        ret = DoCompassConvolve((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, counter, (float)0);
+      break;                                                                                
+    }
+    
+    if (!ret) 
+      break;
+  }
+
+  imCounterEnd(counter);
+
+  return ret;
+}
+
+template <class T, class KT, class CT> 
+static int DoConvolveDual(T* map, T* new_map, int width, int height, KT* kernel_map1, KT* kernel_map2, int kernel_width, int kernel_height, int counter, CT)
+{
+  CT value1, value2, value;
+  KT total1, total2, *kernel_line;
+  int offset, new_offset, i, j, x, y;
+
+  int kh2 = kernel_height/2;
+  int kw2 = kernel_width/2;
+
+  if (kernel_height % 2 == 0) kh2--;
+  if (kernel_width % 2 == 0) kw2--;
+
+  total1 = 0;
+  for(j = 0; j < kernel_height; j++) 
+  {
+    for(i = 0; i < kernel_width; i++)
+      total1 += kernel_map1[j*kernel_width + i];
+  }
+
+  if (total1 == 0)
+    total1 = 1;
+
+  total2 = 0;
+  for(j = 0; j < kernel_height; j++) 
+  {
+    for(i = 0; i < kernel_width; i++)
+      total2 += kernel_map2[j*kernel_width + i];
+  }
+
+  if (total2 == 0)
+    total2 = 1;
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      value1 = 0;
+      value2 = 0;
+    
+      for(y = -kh2; y <= kh2; y++)
+      {
+        if (j + y < 0)             // pass the bottom border
+          offset = -(y + j + 1) * width;
+        else if (j + y >= height)  // pass the top border
+          offset = (2*height - 1 - (j + y)) * width;
+        else
+          offset = (j + y) * width;
+
+        kernel_line = kernel_map1 + (y+kh2)*kernel_width;
+        for(x = -kw2; x <= kw2; x++)
+        {
+          if (i + x < 0)            // pass the left border
+            value1 += kernel_line[x+kw2] * map[offset - (i + x + 1)];
+          else if (i + x >= width)  // pass the right border
+            value1 += kernel_line[x+kw2] * map[offset + 2*width - 1 - (i + x)];
+          else if (offset != -1)
+            value1 += kernel_line[x+kw2] * map[offset + (i + x)];
+        }
+
+        kernel_line = kernel_map2 + (y+kh2)*kernel_width;
+        for(x = -kw2; x <= kw2; x++)
+        {
+          if (i + x < 0)            // pass the left border
+            value2 += kernel_line[x+kw2] * map[offset - (i + x + 1)];
+          else if (i + x >= width)  // pass the right border
+            value2 += kernel_line[x+kw2] * map[offset + 2*width - 1 - (i + x)];
+          else if (offset != -1)
+            value2 += kernel_line[x+kw2] * map[offset + (i + x)];
+        }
+      }
+      
+      value1 /= total1;
+      value2 /= total2;
+
+      value = (CT)sqrt((double)(value1*value1 + value2*value2));
+
+      int size_of = sizeof(imbyte);
+      if (sizeof(T) == size_of)
+        new_map[new_offset + i] = (T)IM_BYTECROP(value);
+      else
+        new_map[new_offset + i] = (T)value;
+    }    
+
+    if (!imCounterInc(counter))
+      return 0;
+  }
+
+  return 1;
+}
+
+template <class KT> 
+static int DoConvolveDualCpx(imcfloat* map, imcfloat* new_map, int width, int height, KT* kernel_map1, KT* kernel_map2, int kernel_width, int kernel_height, int counter)
+{
+  imcfloat value1, value2;
+  KT total1, total2, *kernel_line;
+  int offset, new_offset, i, j, x, y;
+
+  int kh2 = kernel_height/2;
+  int kw2 = kernel_width/2;
+
+  if (kernel_height % 2 == 0) kh2--;
+  if (kernel_width % 2 == 0) kw2--;
+
+  total1 = 0;
+  for(j = 0; j < kernel_height; j++) 
+  {
+    for(i = 0; i < kernel_width; i++)
+      total1 += kernel_map1[j*kernel_width + i];
+  }
+
+  if (total1 == 0)
+    total1 = 1;
+
+  total2 = 0;
+  for(j = 0; j < kernel_height; j++) 
+  {
+    for(i = 0; i < kernel_width; i++)
+      total2 += kernel_map1[j*kernel_width + i];
+  }
+
+  if (total2 == 0)
+    total2 = 1;
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      value1 = 0;
+      value2 = 0;
+    
+      for(y = -kh2; y <= kh2; y++)
+      {
+        if (j + y < 0)             // pass the bottom border
+          offset = -(y + j + 1) * width;
+        else if (j + y >= height)  // pass the top border
+          offset = (2*height - 1 - (j + y)) * width;
+        else
+          offset = (j + y) * width;
+
+        kernel_line = kernel_map1 + (y+kh2)*kernel_width;
+        for(x = -kw2; x <= kw2; x++)
+        {
+          if (i + x < 0)            // pass the left border
+            value1 += map[offset - (i + x + 1)] * (float)kernel_line[x+kw2];
+          else if (i + x >= width)  // pass the right border
+            value1 += map[offset + 2*width - 1 - (i + x)] * (float)kernel_line[x+kw2];
+          else if (offset != -1)
+            value1 += map[offset + (i + x)] * (float)kernel_line[x+kw2];
+        }
+
+        kernel_line = kernel_map2 + (y+kh2)*kernel_width;
+        for(x = -kw2; x <= kw2; x++)
+        {
+          if (i + x < 0)            // pass the left border
+            value2 += map[offset - (i + x + 1)] * (float)kernel_line[x+kw2];
+          else if (i + x >= width)  // pass the right border
+            value2 += map[offset + 2*width - 1 - (i + x)] * (float)kernel_line[x+kw2];
+          else if (offset != -1)
+            value2 += map[offset + (i + x)] * (float)kernel_line[x+kw2];
+        }
+      }
+      
+      value1 /= (float)total1;
+      value2 /= (float)total2;
+
+      new_map[new_offset + i] = sqrt(value1*value1 + value2*value2);
+    }    
+
+    if (!imCounterInc(counter))
+      return 0;
+  }
+
+  return 1;
+}
+
+int imProcessConvolveDual(const imImage* src_image, imImage* dst_image, const imImage *kernel1, const imImage *kernel2)
+{
+  int counter = imCounterBegin("Convolution");
+  const char* msg = (const char*)imImageGetAttribute(kernel1, "Description", NULL, NULL);
+  if (!msg) msg = "Filtering...";
+  imCounterTotal(counter, src_image->depth*src_image->height, msg);
+
+  int ret = 0;
+
+  for (int i = 0; i < src_image->depth; i++)
+  {
+    switch(src_image->data_type)
+    {
+    case IM_BYTE:
+      if (kernel1->data_type == IM_INT)
+        ret = DoConvolveDual((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel1->data[0], (int*)kernel2->data[0], kernel1->width, kernel1->height, counter, (int)0);
+      else
+        ret = DoConvolveDual((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel1->data[0], (float*)kernel2->data[0], kernel1->width, kernel1->height, counter, (float)0);
+      break;                                                                                
+    case IM_USHORT:                                                                           
+      if (kernel1->data_type == IM_INT)
+        ret = DoConvolveDual((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel1->data[0], (int*)kernel2->data[0], kernel1->width, kernel1->height, counter, (int)0);
+      else
+        ret = DoConvolveDual((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel1->data[0], (float*)kernel2->data[0], kernel1->width, kernel1->height, counter, (float)0);
+      break;                                                                                
+    case IM_INT:                                                                           
+      if (kernel1->data_type == IM_INT)
+        ret = DoConvolveDual((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel1->data[0], (int*)kernel2->data[0], kernel1->width, kernel1->height, counter, (int)0);
+      else
+        ret = DoConvolveDual((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel1->data[0], (float*)kernel2->data[0], kernel1->width, kernel1->height, counter, (float)0);
+      break;                                                                                
+    case IM_FLOAT:                                                                           
+      if (kernel1->data_type == IM_INT)
+        ret = DoConvolveDual((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel1->data[0], (int*)kernel2->data[0], kernel1->width, kernel1->height, counter, (float)0);
+      else
+        ret = DoConvolveDual((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel1->data[0], (float*)kernel2->data[0], kernel1->width, kernel1->height, counter, (float)0);
+      break;                                                                                
+    case IM_CFLOAT:            
+      if (kernel1->data_type == IM_INT)
+        ret = DoConvolveDualCpx((imcfloat*)src_image->data[i], (imcfloat*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel1->data[0], (int*)kernel2->data[0], kernel1->width, kernel1->height, counter);
+      else
+        ret = DoConvolveDualCpx((imcfloat*)src_image->data[i], (imcfloat*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel1->data[0], (float*)kernel2->data[0], kernel1->width, kernel1->height, counter);
+      break;
+    }
+    
+    if (!ret) 
+      break;
+  }
+
+  imCounterEnd(counter);
+
+  return ret;
+}
+
+template <class T, class KT, class CT> 
+static int DoConvolve(T* map, T* new_map, int width, int height, KT* kernel_map, int kernel_width, int kernel_height, int counter, CT)
+{
+  CT value;
+  KT total, *kernel_line;
+  int offset, new_offset, i, j, x, y;
+
+  int kh2 = kernel_height/2;
+  int kw2 = kernel_width/2;
+
+  if (kernel_height % 2 == 0) kh2--;
+  if (kernel_width % 2 == 0) kw2--;
+
+  total = 0;
+  for(j = 0; j < kernel_height; j++) 
+  {
+    for(i = 0; i < kernel_width; i++)
+      total += kernel_map[j*kernel_width + i];
+  }
+
+  if (total == 0)
+    total = 1;
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      value = 0;
+    
+      for(y = -kh2; y <= kh2; y++)
+      {
+        kernel_line = kernel_map + (y+kh2)*kernel_width;
+
+        if (j + y < 0)             // pass the bottom border
+          offset = -(y + j + 1) * width;
+        else if (j + y >= height)  // pass the top border
+          offset = (2*height - 1 - (j + y)) * width;
+        else
+          offset = (j + y) * width;
+
+        for(x = -kw2; x <= kw2; x++)
+        {
+          if (i + x < 0)            // pass the left border
+            value += kernel_line[x+kw2] * map[offset - (i + x + 1)];
+          else if (i + x >= width)  // pass the right border
+            value += kernel_line[x+kw2] * map[offset + 2*width - 1 - (i + x)];
+          else if (offset != -1)
+            value += kernel_line[x+kw2] * map[offset + (i + x)];
+        }
+      }
+      
+      value /= total;
+
+      int size_of = sizeof(imbyte);
+      if (sizeof(T) == size_of)
+        new_map[new_offset + i] = (T)IM_BYTECROP(value);
+      else
+        new_map[new_offset + i] = (T)value;
+    }    
+
+    if (!imCounterInc(counter))
+      return 0;
+  }
+
+  return 1;
+}
+
+template <class KT> 
+static int DoConvolveCpx(imcfloat* map, imcfloat* new_map, int width, int height, KT* kernel_map, int kernel_width, int kernel_height, int counter)
+{
+  imcfloat value;
+  KT total, *kernel_line;
+  int offset, new_offset, i, j, x, y;
+
+  int kh2 = kernel_height/2;
+  int kw2 = kernel_width/2;
+
+  if (kernel_height % 2 == 0) kh2--;
+  if (kernel_width % 2 == 0) kw2--;
+
+  total = 0;
+  for(j = 0; j < kernel_height; j++) 
+  {
+    for(i = 0; i < kernel_width; i++)
+      total += kernel_map[j*kernel_width + i];
+  }
+
+  if (total == 0)
+    total = 1;
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      value = 0;
+    
+      for(y = -kh2; y <= kh2; y++)
+      {
+        kernel_line = kernel_map + (y+kh2)*kernel_width;
+
+        if (j + y < 0)             // pass the bottom border
+          offset = -(y + j + 1) * width;
+        else if (j + y >= height)  // pass the top border
+          offset = (2*height - 1 - (j + y)) * width;
+        else
+          offset = (j + y) * width;
+
+        for(x = -kw2; x <= kw2; x++)
+        {
+          if (i + x < 0)            // pass the left border
+            value += map[offset - (i + x + 1)] * (float)kernel_line[x+kw2];
+          else if (i + x >= width)  // pass the right border
+            value += map[offset + 2*width - 1 - (i + x)] * (float)kernel_line[x+kw2];
+          else if (offset != -1)
+            value += map[offset + (i + x)] * (float)kernel_line[x+kw2];
+        }
+      }
+      
+      value /= (float)total;
+
+      new_map[new_offset + i] = value;
+    }    
+
+    if (!imCounterInc(counter))
+      return 0;
+  }
+
+  return 1;
+}
+
+static int DoConvolveStep(const imImage* src_image, imImage* dst_image, const imImage *kernel, int counter)
+{
+  int ret = 0;
+
+  for (int i = 0; i < src_image->depth; i++)
+  {
+    switch(src_image->data_type)
+    {
+    case IM_BYTE:
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolve((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (int)0);
+      else
+        ret = DoConvolve((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_USHORT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolve((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (int)0);
+      else
+        ret = DoConvolve((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_INT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolve((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (int)0);
+      else
+        ret = DoConvolve((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_FLOAT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolve((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      else
+        ret = DoConvolve((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_CFLOAT:            
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolveCpx((imcfloat*)src_image->data[i], (imcfloat*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter);
+      else
+        ret = DoConvolveCpx((imcfloat*)src_image->data[i], (imcfloat*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter);
+      break;
+    }
+    
+    if (!ret) 
+      break;
+  }
+
+  return ret;
+}
+
+int imProcessConvolve(const imImage* src_image, imImage* dst_image, const imImage *kernel)
+{
+  int counter = imCounterBegin("Convolution");
+  const char* msg = (const char*)imImageGetAttribute(kernel, "Description", NULL, NULL);
+  if (!msg) msg = "Filtering...";
+  imCounterTotal(counter, src_image->depth*src_image->height, msg);
+
+  int ret = DoConvolveStep(src_image, dst_image, kernel, counter);
+
+  imCounterEnd(counter);
+
+  return ret;
+}
+
+int imProcessConvolveRep(const imImage* src_image, imImage* dst_image, const imImage *kernel, int ntimes)
+{
+  imImage *AuxImage = imImageClone(dst_image);
+  if (!AuxImage)
+    return 0;
+
+  int counter = imCounterBegin("Repeated Convolution");
+  const char* msg = (const char*)imImageGetAttribute(kernel, "Description", NULL, NULL);
+  if (!msg) msg = "Filtering...";
+  imCounterTotal(counter, src_image->depth*src_image->height*ntimes, msg);
+
+  const imImage *image1 = src_image;
+  imImage *image2 = dst_image;
+
+  for (int i = 0; i < ntimes; i++)
+  {
+    if (!DoConvolveStep(image1, image2, kernel, counter))
+    {
+      imCounterEnd(counter);
+      imImageDestroy(AuxImage);
+      return 0;
+    }
+    
+    image1 = image2;
+
+    if (image1 == dst_image)
+      image2 = AuxImage;
+    else
+      image2 = dst_image;
+  }
+
+  // The result is in image1, if in the Aux swap the data
+  if (image1 == AuxImage)
+  {
+    void** temp = (void**)dst_image->data;
+    dst_image->data = AuxImage->data;
+    AuxImage->data = (void**)temp;
+  }
+
+  imCounterEnd(counter);
+  imImageDestroy(AuxImage);
+
+  return 1;
+}
+
+template <class T, class KT, class CT> 
+static int DoConvolveSep(T* map, T* new_map, int width, int height, KT* kernel_map, int kernel_width, int kernel_height, int counter, CT)
+{
+  CT value;
+  KT totalV, totalH, *kernel_line;
+  T* aux_line;
+  int offset, new_offset, i, j;
+
+  int kh2 = kernel_height/2;
+  int kw2 = kernel_width/2;
+
+  if (kernel_height % 2 == 0) kh2--;
+  if (kernel_width % 2 == 0) kw2--;
+
+  // use only the first line and the first column of the kernel
+
+  totalV = 0;
+  for(j = 0; j < kernel_height; j++) 
+    totalV += kernel_map[j*kernel_width];
+
+  if (totalV == 0)
+    totalV = 1;
+
+  totalH = 0;
+  for(i = 0; i < kernel_width; i++)
+    totalH += kernel_map[i];
+
+  if (totalH == 0)
+    totalH = 1;
+
+  aux_line = (T*)malloc(width*sizeof(T));
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      int y;
+      value = 0;
+
+      // first pass, only for columns
+    
+      for(y = -kh2; y <= kh2; y++)
+      {
+        kernel_line = kernel_map + (y+kh2)*kernel_width;
+
+        if (j + y < 0)             // pass the bottom border
+          offset = -(y + j + 1) * width;
+        else if (j + y >= height)  // pass the top border
+          offset = (2*height - 1 - (j + y)) * width;
+        else
+          offset = (j + y) * width;
+
+        if (offset != -1)
+          value += kernel_line[0] * map[offset + i];
+      }
+      
+      value /= totalV;
+
+      int size_of = sizeof(imbyte);
+      if (sizeof(T) == size_of)
+        new_map[new_offset + i] = (T)IM_BYTECROP(value);
+      else
+        new_map[new_offset + i] = (T)value;
+    }    
+
+    if (!imCounterInc(counter))
+    {
+      free(aux_line);
+      return 0;
+    }
+  }
+
+  for(j = 0; j < height; j++)
+  {
+    offset = new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      int x;
+      value = 0;
+
+      // second pass, only for lines, but has to use an auxiliar buffer
+    
+      kernel_line = kernel_map;
+
+      for(x = -kw2; x <= kw2; x++)
+      {
+        if (i + x < 0)            // pass the left border
+          value += kernel_line[x+kw2] * new_map[offset - (i + x + 1)];
+        else if (i + x >= width)  // pass the right border
+          value += kernel_line[x+kw2] * new_map[offset + 2*width - 1 - (i + x)];
+        else
+          value += kernel_line[x+kw2] * new_map[offset + (i + x)];
+      }
+      
+      value /= totalH;
+
+      int size_of = sizeof(imbyte);
+      if (sizeof(T) == size_of)
+        aux_line[i] = (T)IM_BYTECROP(value);
+      else
+        aux_line[i] = (T)value;
+    }    
+
+    memcpy(new_map + new_offset, aux_line, width*sizeof(T));
+
+    if (!imCounterInc(counter))
+    {
+      free(aux_line);
+      return 0;
+    }
+  }
+
+  free(aux_line);
+  return 1;
+}
+
+
+template <class KT> 
+static int DoConvolveSepCpx(imcfloat* map, imcfloat* new_map, int width, int height, KT* kernel_map, int kernel_width, int kernel_height, int counter)
+{
+  imcfloat value;
+  KT totalV, totalH, *kernel_line;
+  imcfloat* aux_line;
+  int offset, new_offset, i, j;
+
+  int kh2 = kernel_height/2;
+  int kw2 = kernel_width/2;
+
+  if (kernel_height % 2 == 0) kh2--;
+  if (kernel_width % 2 == 0) kw2--;
+
+  // use only the first line and the first column of the kernel
+
+  totalV = 0;
+  for(j = 0; j < kernel_height; j++) 
+    totalV += kernel_map[j*kernel_width];
+
+  if (totalV == 0)
+    totalV = 1;
+
+  totalH = 0;
+  for(i = 0; i < kernel_width; i++)
+    totalH += kernel_map[i];
+
+  if (totalH == 0)
+    totalH = 1;
+
+  aux_line = (imcfloat*)malloc(width*sizeof(imcfloat));
+
+  for(j = 0; j < height; j++)
+  {
+    new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      int y;
+      value = 0;
+
+      // first pass, only for columns
+    
+      for(y = -kh2; y <= kh2; y++)
+      {
+        kernel_line = kernel_map + (y+kh2)*kernel_width;
+
+        if (j + y < 0)             // pass the bottom border
+          offset = -(y + j + 1) * width;
+        else if (j + y >= height)  // pass the top border
+          offset = (2*height - 1 - (j + y)) * width;
+        else
+          offset = (j + y) * width;
+
+        if (offset != -1)
+          value += map[offset + i] * (float)kernel_line[0];
+      }
+      
+      value /= (float)totalV;
+
+      new_map[new_offset + i] = value;
+    }    
+
+    if (!imCounterInc(counter))
+    {
+      free(aux_line);
+      return 0;
+    }
+  }
+
+  for(j = 0; j < height; j++)
+  {
+    offset = new_offset = j * width;
+
+    for(i = 0; i < width; i++)
+    {
+      int x;
+      value = 0;
+
+      // second pass, only for lines, but has to use an auxiliar buffer
+    
+      kernel_line = kernel_map;
+    
+      for(x = -kw2; x <= kw2; x++)
+      {
+        if (i + x < 0)            // pass the left border
+          value += new_map[offset - (i + x + 1)] * (float)kernel_line[x+kw2];
+        else if (i + x >= width)  // pass the right border
+          value += new_map[offset + 2*width - 1 - (i + x)] * (float)kernel_line[x+kw2];
+        else if (offset != -1)
+          value += new_map[offset + (i + x)] * (float)kernel_line[x+kw2];
+      }
+      
+      value /= (float)totalH;
+
+      aux_line[i] = value;
+    }    
+
+    memcpy(new_map + new_offset, aux_line, width*sizeof(imcfloat));
+
+    if (!imCounterInc(counter))
+    {
+      free(aux_line);
+      return 0;
+    }
+  }
+
+  free(aux_line);
+  return 1;
+}
+
+int imProcessConvolveSep(const imImage* src_image, imImage* dst_image, const imImage *kernel)
+{
+  int counter = imCounterBegin("Separable Convolution");
+  const char* msg = (const char*)imImageGetAttribute(kernel, "Description", NULL, NULL);
+  if (!msg) msg = "Filtering...";
+  imCounterTotal(counter, 2*src_image->depth*src_image->height, msg);
+
+  int ret = 0;
+
+  for (int i = 0; i < src_image->depth; i++)
+  {
+    switch(src_image->data_type)
+    {
+    case IM_BYTE:
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolveSep((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (int)0);
+      else
+        ret = DoConvolveSep((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_USHORT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolveSep((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (int)0);
+      else
+        ret = DoConvolveSep((imushort*)src_image->data[i], (imushort*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_INT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolveSep((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (int)0);
+      else
+        ret = DoConvolveSep((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_FLOAT:                                                                           
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolveSep((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      else
+        ret = DoConvolveSep((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter, (float)0);
+      break;                                                                                
+    case IM_CFLOAT:            
+      if (kernel->data_type == IM_INT)
+        ret = DoConvolveSepCpx((imcfloat*)src_image->data[i], (imcfloat*)dst_image->data[i], src_image->width, src_image->height, (int*)kernel->data[0], kernel->width, kernel->height, counter);
+      else
+        ret = DoConvolveSepCpx((imcfloat*)src_image->data[i], (imcfloat*)dst_image->data[i], src_image->width, src_image->height, (float*)kernel->data[0], kernel->width, kernel->height, counter);
+      break;
+    }
+    
+    if (!ret) 
+      break;
+  }
+
+  imCounterEnd(counter);
+
+  return ret;
+}
+
+/*
+Description:	
+    Can be used to find zero crossing of second derivative,
+		laplace. Can also be used to determine any other kind
+		of crossing. Pixels below or equal to 't' are set if the pixel
+		to the right or below is above 't', pixels above 't' are
+		set if the pixel to the right or below is below or equal to
+		't'. Pixels that are "set" are set to the maximum absolute
+		difference of the two neighbours, to indicate the strength
+		of the edge.
+
+		| IF (crossing t)
+		|   out(x,y) = MAX(ABS(in(x,y)-in(x+1,y)), ABS(in(x,y)-in(x,y+1)))
+		| ELSE
+		|   out(x,y) = 0
+
+Author:		Tor Lønnestad, BLAB, Ifi, UiO
+
+Copyright 1991, Blab, UiO
+Image processing lab, Department of Informatics
+University of Oslo
+*/
+template <class T> 
+static void do_crossing(T* iband, T* oband, int width, int height, T t)
+{
+  int x, y, offset00 = 0, offset10 = 0, offset01 = 0;
+  T v, diff;
+
+  for (y=0; y < height-1; y++)
+  {
+    offset00 = y*width;
+    offset10 = (y+1)*width;
+    offset01 = offset00 + 1;
+
+    for (x=0; x < width-1; x++)
+    {
+      v = 0;
+
+      if (iband[offset00] <= t)
+      {
+        if (iband[offset10] > t) 
+          v = iband[offset10]-iband[offset00];
+
+	      if (iband[offset01] > t) 
+        {
+          diff = iband[offset01]-iband[offset00];
+          if (diff > v) v = diff;
+        }
+      }
+      else
+      {
+	      if (iband[offset10] <= t) 
+          v = iband[offset00]-iband[offset10];
+
+	      if (iband[offset01] <= t) 
+        {
+          diff = iband[offset00]-iband[offset01];
+          if (diff > v) v = diff;
+        }
+      }
+
+      oband[offset00] = v;
+
+      offset00++;
+      offset10++;
+      offset01++;
+    }
+
+    /* last pixel on line */
+    offset00++;
+    offset10++;
+
+    v = 0;
+
+    if (iband[offset00] <= t)
+    {
+      if (iband[offset10] > t)
+        v = iband[offset10]-iband[offset00];
+    }
+    else
+    {
+      if (iband[offset10] <= t)
+        v = iband[offset00]-iband[offset10];
+    }
+
+    oband[offset00] = v;
+  }
+
+  /* last line */
+  offset00 = y*width;
+  offset01 = offset00 + 1;
+
+  for (x=0; x < width-1; x++)
+  {
+    v = 0;
+
+    if (iband[offset00] <= t)
+    {
+      if (iband[offset01] > t)
+        v = iband[offset01]-iband[offset00];
+    }
+    else
+    {
+      if (iband[offset01] <= t)
+        v = iband[offset00]-iband[offset01];
+    }
+
+    oband[offset00] = v;
+
+    offset00++;
+    offset01++;
+  }
+
+  offset00++;
+
+  /* last pixel */
+  oband[offset00] = 0;
+}
+
+void imProcessZeroCrossing(const imImage* src_image, imImage* dst_image)
+{
+  for (int i = 0; i < src_image->depth; i++)
+  {
+    switch(src_image->data_type)
+    {
+    case IM_INT:                                                                           
+      do_crossing((int*)src_image->data[i], (int*)dst_image->data[i], src_image->width, src_image->height, 0);
+      break;                                                                                
+    case IM_FLOAT:                                                                           
+      do_crossing((float*)src_image->data[i], (float*)dst_image->data[i], src_image->width, src_image->height, 0.0f);
+      break;                                                                                
+    }
+  }
+}
+
+int imProcessBarlettConvolve(const imImage* src_image, imImage* dst_image, int kernel_size)
+{
+  imImage* kernel = imImageCreate(kernel_size, kernel_size, IM_GRAY, IM_INT);
+  if (!kernel)
+    return 0;
+
+  imImageSetAttribute(kernel, "Description", IM_BYTE, -1, (void*)"Barlett");
+
+  int* kernel_data = (int*)kernel->data[0];
+  int half = kernel_size / 2;
+  for (int i = 0; i < kernel_size; i++)
+  {
+    if (i <= half)
+      kernel_data[i] = i+1;
+    else
+      kernel_data[i] = kernel_size-i;
+  }
+  for (int j = 0; j < kernel_size; j++)
+  {
+    if (j <= half)
+      kernel_data[j*kernel_size] = j+1;
+    else
+      kernel_data[j*kernel_size] = kernel_size-j;
+  }
+
+  int ret = imProcessConvolveSep(src_image, dst_image, kernel);
+
+  imImageDestroy(kernel);
+
+  return ret;
+}
+
+int imProcessSobelConvolve(const imImage* src_image, imImage* dst_image)
+{
+	int ret = 0;
+
+  imImage* kernel1 = imKernelSobel();
+  imImage* kernel2 = imImageCreate(3, 3, IM_GRAY, IM_INT);
+  imProcessRotate90(kernel1, kernel2, 1);
+
+  ret = imProcessConvolveDual(src_image, dst_image, kernel1, kernel2);
+
+  imImageDestroy(kernel1);
+  imImageDestroy(kernel2);
+
+  return ret;
+}
+
+int imProcessPrewittConvolve(const imImage* src_image, imImage* dst_image)
+{
+	int ret = 0;
+
+  imImage* kernel1 = imKernelPrewitt();
+  imImage* kernel2 = imImageClone(kernel1);
+  imProcessRotate90(kernel1, kernel2, 1);
+
+  ret = imProcessConvolveDual(src_image, dst_image, kernel1, kernel2);
+
+  imImageDestroy(kernel1);
+  imImageDestroy(kernel2);
+
+  return ret;
+}
+
+int imProcessSplineEdgeConvolve(const imImage* src_image, imImage* dst_image)
+{
+	int ret = 0;
+
+  imImage* tmp_image = imImageClone(src_image);
+  if (!tmp_image) return 0;
+
+  imImage* kernel1 = imImageCreate(5, 5, IM_GRAY, IM_INT);
+  imImageSetAttribute(kernel1, "Description", IM_BYTE, -1, (void*)"SplineEdge");
+
+  int* kernel_data = (int*)kernel1->data[0];
+  kernel_data[10] = -1;
+  kernel_data[11] = 8;
+  kernel_data[12] = 0;
+  kernel_data[13] = -8;
+  kernel_data[14] = 1;
+
+  imImage* kernel2 = imImageClone(kernel1);
+  imProcessRotate90(kernel1, kernel2, 1);
+
+  imImage* kernel3 = imImageClone(kernel1);
+  imProcessRotateKernel(kernel3);
+
+  imImage* kernel4 = imImageClone(kernel1);
+  imProcessRotate90(kernel3, kernel4, 1);
+
+  ret = imProcessConvolveDual(src_image, tmp_image, kernel1, kernel2);
+  ret = imProcessConvolveDual(src_image, dst_image, kernel3, kernel4);
+
+  imProcessArithmeticConstOp(tmp_image, (float)sqrt(2.0), tmp_image, IM_BIN_MUL);
+  imProcessArithmeticOp(tmp_image, dst_image, dst_image, IM_BIN_ADD);
+
+  imImageDestroy(kernel1);
+  imImageDestroy(kernel2);
+  imImageDestroy(kernel3);
+  imImageDestroy(kernel4);
+  imImageDestroy(tmp_image);
+
+  return ret;
+}
+
+int imGaussianStdDev2KernelSize(float stddev)
+{
+  if (stddev < 0)
+    return (int)-stddev;
+  else
+  {
+	  int width = (int)(3.35*stddev + 0.3333);
+    return 2*width + 1;
+  }
+}
+
+float imGaussianKernelSize2StdDev(int kernel_size)
+{
+  int width = (kernel_size - 1)/2;
+	return (width - 0.3333f)/3.35f;
+}
+
+int imProcessGaussianConvolve(const imImage* src_image, imImage* dst_image, float stddev)
+{
+  int kernel_size = imGaussianStdDev2KernelSize(stddev);
+
+  imImage* kernel = imImageCreate(kernel_size, kernel_size, IM_GRAY, IM_FLOAT);
+  if (!kernel)
+    return 0;
+
+  imImageSetAttribute(kernel, "Description", IM_BYTE, -1, (void*)"Gaussian");
+  imProcessRenderGaussian(kernel, stddev);
+
+  int ret = imProcessConvolveSep(src_image, dst_image, kernel);
+
+  imImageDestroy(kernel);
+
+  return ret;
+}
+
+int imProcessLapOfGaussianConvolve(const imImage* src_image, imImage* dst_image, float stddev)
+{
+  int kernel_size = imGaussianStdDev2KernelSize(stddev);
+
+  imImage* kernel = imImageCreate(kernel_size, kernel_size, IM_GRAY, IM_FLOAT);
+  if (!kernel)
+    return 0;
+
+  imImageSetAttribute(kernel, "Description", IM_BYTE, -1, (void*)"Laplacian Of Gaussian");
+  imProcessRenderLapOfGaussian(kernel, stddev);
+
+  int ret;
+  if (src_image->data_type == IM_BYTE || src_image->data_type == IM_USHORT)
+  {
+    imImage* aux_image = imImageClone(dst_image);
+    if (!aux_image)
+    {
+      imImageDestroy(kernel);
+      return 0;
+    }
+
+    imProcessUnArithmeticOp(src_image, aux_image, IM_UN_EQL);  // Convert to IM_INT
+    ret = imProcessConvolve(aux_image, dst_image, kernel);
+    imImageDestroy(aux_image);
+  }
+  else
+    ret = imProcessConvolve(src_image, dst_image, kernel);
+
+  imImageDestroy(kernel);
+
+  return ret;
+}
+
+int imProcessDiffOfGaussianConvolve(const imImage* src_image, imImage* dst_image, float stddev1, float stddev2)
+{
+  imImage* aux_image1 = imImageClone(src_image);
+  imImage* aux_image2 = imImageClone(src_image);
+  if (!aux_image1 || !aux_image2)
+  {
+    if (aux_image1) imImageDestroy(aux_image1);
+    return 0;
+  }
+
+  int kernel_size1 = imGaussianStdDev2KernelSize(stddev1);
+  int kernel_size2 = imGaussianStdDev2KernelSize(stddev2);
+  int size = kernel_size1;
+  if (kernel_size1 < kernel_size2) size = kernel_size2;
+
+  imImage* kernel1 = imImageCreate(size, size, IM_GRAY, IM_FLOAT);
+  imImage* kernel2 = imImageCreate(size, size, IM_GRAY, IM_FLOAT);
+  if (!kernel1 || !kernel2)
+  {
+    if (kernel1) imImageDestroy(kernel1);
+    if (kernel2) imImageDestroy(kernel2);
+    imImageDestroy(aux_image1);
+    imImageDestroy(aux_image2);
+    return 0;
+  }
+
+  imImageSetAttribute(kernel1, "Description", IM_BYTE, -1, (void*)"Gaussian1");
+  imImageSetAttribute(kernel2, "Description", IM_BYTE, -1, (void*)"Gaussian2");
+
+  imProcessRenderGaussian(kernel1, stddev1);
+  imProcessRenderGaussian(kernel2, stddev2);
+
+  if (!imProcessConvolve(src_image, aux_image1, kernel1) ||
+      !imProcessConvolve(src_image, aux_image2, kernel2))
+  {
+    imImageDestroy(kernel1);
+    imImageDestroy(kernel2);
+    imImageDestroy(aux_image1);
+    imImageDestroy(aux_image2);
+    return 0;
+  }
+
+  imProcessArithmeticOp(aux_image1, aux_image2, dst_image, IM_BIN_SUB);
+
+  imImageDestroy(kernel1);
+  imImageDestroy(kernel2);
+  imImageDestroy(aux_image1);
+  imImageDestroy(aux_image2);
+
+  return 1;
+}
+
+#ifdef _TEST_CODE_
+int imProcessDiffOfGaussianConvolveTEST(const imImage* src_image, imImage* dst_image, float stddev1, float stddev2)
+{
+  int kernel_size1 = imGaussianStdDev2KernelSize(stddev1);
+  int kernel_size2 = imGaussianStdDev2KernelSize(stddev2);
+  int size = kernel_size1;
+  if (kernel_size1 < kernel_size2) size = kernel_size2;
+
+  imImage* kernel1 = imImageCreate(size, size, IM_GRAY, IM_FLOAT);
+  imImage* kernel2 = imImageCreate(size, size, IM_GRAY, IM_FLOAT);
+  if (!kernel1 || !kernel2)
+  {
+    if (kernel1) imImageDestroy(kernel1);
+    if (kernel2) imImageDestroy(kernel2);
+    return 0;
+  }
+
+  imImageSetAttribute(kernel1, "Description", IM_BYTE, -1, (void*)"Gaussian");
+  imImageSetAttribute(kernel2, "Description", IM_BYTE, -1, (void*)"Gaussian");
+
+  imProcessRenderGaussian(kernel1, stddev1);
+  imProcessRenderGaussian(kernel2, stddev2);
+
+  // ERROR: kernel 1 should be multiplied by a factor to improve the difference.
+
+  imProcessArithmeticOp(kernel1, kernel2, kernel1, IM_BIN_SUB);
+  imImageSetAttribute(kernel1, "Description", IM_BYTE, -1, (void*)"Difference of Gaussian");
+
+  int ret = 0;
+  if (src_image->data_type == IM_BYTE || src_image->data_type == IM_USHORT)
+  {
+    imImage* aux_image = imImageClone(dst_image);
+    if (!aux_image)
+    {
+      imImageDestroy(kernel1);
+      imImageDestroy(kernel2);
+      return 0;
+    }
+
+    imProcessUnArithmeticOp(src_image, aux_image, IM_UN_EQL);  // Convert to IM_INT
+    ret = imProcessConvolve(aux_image, dst_image, kernel1);
+    imImageDestroy(aux_image);
+  }
+  else
+    ret = imProcessConvolve(src_image, dst_image, kernel1);
+
+  imImageDestroy(kernel1);
+  imImageDestroy(kernel2);
+
+  return ret;
+}
+#endif
+
+int imProcessMeanConvolve(const imImage* src_image, imImage* dst_image, int ks)
+{
+  int counter = imCounterBegin("Mean Convolve");
+  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");
+
+  imImage* kernel = imImageCreate(ks, ks, IM_GRAY, IM_INT);
+
+  int* kernel_data = (int*)kernel->data[0];
+
+  int ks2 = ks/2;
+  for(int ky = 0; ky < ks; ky++)
+  {
+    int ky2 = ky-ks2;
+    ky2 = ky2*ky2;
+    for(int kx = 0; kx < ks; kx++) 
+    {
+      int kx2 = kx-ks2;
+      kx2 = kx2*kx2;
+      int radius = imRound(sqrt(double(kx2 + ky2)));
+      if (radius <= ks2)
+        kernel_data[ky*ks + kx] = 1;
+    }
+  }
+
+  int ret = DoConvolveStep(src_image, dst_image, kernel, counter);
+
+  imImageDestroy(kernel);
+  imCounterEnd(counter);
+
+  return ret;
+}