summaryrefslogtreecommitdiff
path: root/src/process/im_convolve.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/process/im_convolve.cpp')
-rw-r--r--src/process/im_convolve.cpp1512
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;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-21 22:46:17 +0000