summaryrefslogtreecommitdiff
path: root/src/process/im_analyze.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/process/im_analyze.cpp')
-rw-r--r--src/process/im_analyze.cpp1262
1 files changed, 1262 insertions, 0 deletions
diff --git a/src/process/im_analyze.cpp b/src/process/im_analyze.cpp
new file mode 100644
index 0000000..50bcbcd
--- /dev/null
+++ b/src/process/im_analyze.cpp
@@ -0,0 +1,1262 @@
+/** \file
+ * \brief Image Analysis
+ *
+ * See Copyright Notice in im_lib.h
+ * $Id: im_analyze.cpp,v 1.1 2008/10/17 06:16:33 scuri Exp $
+ */
+
+
+#include <im.h>
+#include <im_util.h>
+#include <im_math.h>
+
+#include "im_process_ana.h"
+#include "im_process_pon.h"
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <memory.h>
+#include <string.h>
+
+#define MAX_COUNT 65536 // maximum number of regions
+
+/* ajust the alias table to be a remap table (final step) */
+static void alias_update(imushort* alias_table, int &region_count)
+{
+ int i, real_count = region_count;
+
+ for (i = 0; i < region_count; i++)
+ {
+ if (alias_table[i])
+ {
+ // search for the first alias
+ imushort prev = alias_table[i];
+ while (alias_table[prev])
+ prev = alias_table[prev];
+
+ alias_table[i] = prev;
+ real_count--; // decrement aliases from the region count
+ }
+ }
+
+ // now all the aliases in the same group point to only one alias
+ // transform the alias table into a remap table
+
+ alias_table[0] = 0;
+ alias_table[1] = 0; // border is mapped to background
+
+ int r = 1;
+ for (i = 2; i < region_count; i++)
+ {
+ if (!alias_table[i])
+ {
+ alias_table[i] = (imushort)r; // only non alias get real values
+ r++;
+ }
+ else
+ alias_table[i] = (imushort)(alias_table[alias_table[i]]);
+ }
+
+ region_count = real_count-2; // remove the regions (background,border) from the count
+}
+
+/* find the smallest region number to be set as alias. */
+static void alias_getmin(imushort* alias_table, imushort region, imushort &min)
+{
+ while (alias_table[region])
+ {
+ if (min > alias_table[region])
+ min = alias_table[region];
+
+ region = alias_table[region];
+ }
+}
+
+/* replace all the aliases of a region by its smallest value. */
+static void alias_setmin(imushort* alias_table, imushort region, imushort min)
+{
+ while (alias_table[region])
+ {
+ imushort next_region = alias_table[region];
+ alias_table[region] = min;
+ region = next_region;
+ }
+
+ if (region != min)
+ alias_table[region] = min;
+}
+
+/* set a region number to be an alias of another */
+static void alias_set(imushort* alias_table, imushort region1, imushort region2)
+{
+ if (region1 == region2)
+ return;
+
+ imushort min = region1<region2? region1: region2;
+
+ alias_getmin(alias_table, region1, min);
+ alias_getmin(alias_table, region2, min);
+
+ if (region1 != min && alias_table[region1] != min)
+ alias_setmin(alias_table, region1, min);
+ if (region2 != min && alias_table[region2] != min)
+ alias_setmin(alias_table, region2, min);
+}
+
+static int DoAnalyzeFindRegions(int width, int height, imbyte* map, imushort* new_map, int connect)
+{
+ int i, j;
+
+ // mark the pixels that touch the border
+ // if a region touch the border, is the invalid region 1
+
+ imbyte* pmap = map;
+ imushort* new_pmap = new_map;
+ for (j = 0; j < width; j++) // first line
+ {
+ if (pmap[j])
+ new_pmap[j] = 1;
+ }
+ pmap += width;
+ new_pmap += width;
+
+ for (i = 1; i < height-1; i++) // first column
+ {
+ if (pmap[0])
+ new_pmap[0] = 1;
+
+ pmap += width;
+ new_pmap += width;
+ }
+
+ // find and connect the regions
+
+ imbyte* pmap1 = map; // previous line (line 0)
+ imushort* new_pmap1 = new_map;
+
+ pmap = map + width; // current line (line 1)
+ new_pmap = new_map + width;
+
+ int region_count = 2; // 0- background, 1-border
+ imushort* alias_table = new imushort [MAX_COUNT];
+ memset(alias_table, 0, MAX_COUNT); // aliases are all zero at start (not used)
+
+ for (i = 1; i < height; i++)
+ {
+ for (j = 1; j < width; j++)
+ {
+ int has_j1 = j < width-1? 1: 0;
+ if (pmap[j])
+ {
+ if (pmap[j-1] || pmap1[j] ||
+ (connect == 8 && (pmap1[j-1] || (has_j1&&pmap1[j+1])))) // 4 or 8 connected to the previous neighbors
+ {
+ imushort region = 0;
+ if (i == height-1 || j == width-1)
+ {
+ region = new_pmap[j] = 1;
+ }
+
+ if (pmap[j-1])
+ {
+ if (!region)
+ region = new_pmap[j-1]; // horizontal neighbor -00
+ else // X1
+ {
+ // this is a right border pixel that connects to an horizontal neighbor
+
+ // this pixel can connect two different regions
+ alias_set(alias_table, region, new_pmap[j-1]);
+ }
+ }
+
+ if (pmap1[j]) // vertical neighbor
+ {
+ if (!region)
+ region = new_pmap1[j]; // isolated vertical neighbor -X-
+ else // 01
+ {
+ // an horizontal neighbor connects to a vertical neighbor -X-
+ // X1
+
+ // this pixel can connect two different regions
+ alias_set(alias_table, region, new_pmap1[j]);
+ }
+ }
+ else if (region && connect==8 && (has_j1&&pmap1[j+1]))
+ {
+ // an horizontal neighbor connects to a right corner neighbor 00X
+ // X1
+
+ // this pixel can connect two different regions
+ alias_set(alias_table, region, new_pmap1[j+1]);
+ }
+
+ if (connect == 8 && (pmap1[j-1] || (has_j1&&pmap1[j+1])) && !region) // isolated corner
+ {
+ // a left corner neighbor or a right corner neighbor X0X
+ // 01
+
+ if (pmap1[j-1]) // left corner
+ region = new_pmap1[j-1];
+
+ if (pmap1[j+1]) // right corner
+ {
+ if (!region) // isolated right corner
+ region = new_pmap1[j+1];
+ else
+ {
+ // this pixel can connect two different regions
+ alias_set(alias_table, new_pmap1[j-1], new_pmap1[j+1]);
+ }
+ }
+ }
+
+ new_pmap[j] = region;
+ }
+ else
+ {
+ // this pixel touches no pixels
+
+ if (i == height-1 || j == width-1)
+ new_pmap[j] = 1;
+ else
+ {
+ // create a new region 000
+ // 01
+ new_pmap[j] = (imushort)region_count;
+ region_count++;
+
+ if (region_count > MAX_COUNT)
+ {
+ delete [] alias_table;
+ return -1;
+ }
+ }
+ }
+ }
+ }
+
+ pmap1 = pmap;
+ new_pmap1 = new_pmap;
+ pmap += width;
+ new_pmap += width;
+ }
+
+ // now all pixels are marked,
+ // but some marks are aliases to others
+
+ // ajust the alias table to be a remap table
+ // and return the real region count
+ alias_update(alias_table, region_count);
+
+ int count = width*height;
+ for (i = 0; i < count; i++)
+ {
+ new_map[i] = alias_table[new_map[i]];
+ }
+
+ delete [] alias_table;
+
+ return region_count;
+}
+
+static int DoAnalyzeFindRegionsBorder(int width, int height, imbyte* map, imushort* new_map, int connect)
+{
+ int i, j;
+
+ imbyte* pmap1 = map - width; // previous line (line -1 = invalid)
+ imushort* new_pmap1 = new_map - width;
+
+ imbyte* pmap = map; // current line (line 0)
+ imushort* new_pmap = new_map;
+
+ int region_count = 2; // still consider: 0- background, 1-border
+ imushort* alias_table = new imushort [MAX_COUNT];
+ memset(alias_table, 0, MAX_COUNT); // aliases are all zero at start (not used)
+
+ for (i = 0; i < height; i++)
+ {
+ for (j = 0; j < width; j++)
+ {
+ if (pmap[j])
+ {
+ int b01 = j > 0? 1: 0; // valid for pmap[j-1]
+ int b10 = i > 0? 1: 0; // valid for pmap1[j]
+ int b11 = i > 0 && j > 0? 1: 0; // valid for pmap1[j-1]
+ int b12 = i > 0 && j < width-1? 1: 0; // valid for pmap1[j+1]
+
+ if ((b01&&pmap[j-1]) || (b10&&pmap1[j]) ||
+ (connect == 8 && ((b11&&pmap1[j-1]) || (b12&&pmap1[j+1])))) // 4 or 8 connected to the previous neighbors
+ {
+ imushort region = 0;
+
+ if (b01&&pmap[j-1])
+ {
+ if (!region)
+ region = new_pmap[j-1]; // horizontal neighbor -00
+ else // X1
+ {
+ // this is a right border pixel that connects to an horizontal neighbor
+
+ // this pixel can connect two different regions
+ alias_set(alias_table, region, new_pmap[j-1]);
+ }
+ }
+
+ if (b10&&pmap1[j]) // vertical neighbor
+ {
+ if (!region)
+ region = new_pmap1[j]; // isolated vertical neighbor -X-
+ else // 01
+ {
+ // an horizontal neighbor connects to a vertical neighbor -X-
+ // X1
+
+ // this pixel can connect two different regions
+ alias_set(alias_table, region, new_pmap1[j]);
+ }
+ }
+ else if (region && connect == 8 && (b12&&pmap1[j+1]))
+ {
+ // an horizontal neighbor connects to a right corner neighbor 00X
+ // X1
+
+ // this pixel can connect two different regions
+ alias_set(alias_table, region, new_pmap1[j+1]);
+ }
+
+ if (connect == 8 && ((b11&&pmap1[j-1]) || (b12&&pmap1[j+1])) && !region) // isolated corner
+ {
+ // a left corner neighbor or a right corner neighbor X0X
+ // 01
+
+ if (b11&&pmap1[j-1]) // left corner
+ region = new_pmap1[j-1];
+
+ if (b12&&pmap1[j+1]) // right corner
+ {
+ if (!region) // isolated right corner
+ region = new_pmap1[j+1];
+ else
+ {
+ // this pixel can connect two different regions
+ alias_set(alias_table, new_pmap1[j-1], new_pmap1[j+1]);
+ }
+ }
+ }
+
+ new_pmap[j] = region;
+ }
+ else
+ {
+ // this pixel touches no pixels
+
+ // create a new region 000
+ // 01
+ new_pmap[j] = (imushort)region_count;
+ region_count++;
+
+ if (region_count > MAX_COUNT)
+ {
+ delete [] alias_table;
+ return -1;
+ }
+ }
+ }
+ }
+
+ pmap1 = pmap;
+ new_pmap1 = new_pmap;
+ pmap += width;
+ new_pmap += width;
+ }
+
+ // now all pixels are marked,
+ // but some marks are aliases to others
+
+ // ajust the alias table to be a remap table
+ // and return the real region count
+ alias_update(alias_table, region_count);
+
+ int count = width*height;
+ for (i = 0; i < count; i++)
+ {
+ new_map[i] = alias_table[new_map[i]];
+ }
+
+ delete [] alias_table;
+
+ return region_count;
+}
+
+int imAnalyzeFindRegions(const imImage* image, imImage* NewImage, int connect, int touch_border)
+{
+ imImageSetAttribute(NewImage, "REGION_CONNECT", IM_BYTE, 1, connect==4?"4":"8");
+ if (touch_border)
+ return DoAnalyzeFindRegionsBorder(image->width, image->height, (imbyte*)image->data[0], (imushort*)NewImage->data[0], connect);
+ else
+ return DoAnalyzeFindRegions(image->width, image->height, (imbyte*)image->data[0], (imushort*)NewImage->data[0], connect);
+}
+
+void imAnalyzeMeasureArea(const imImage* image, int* data_area, int region_count)
+{
+ imushort* img_data = (imushort*)image->data[0];
+
+ memset(data_area, 0, region_count*sizeof(int));
+
+ for (int i = 0; i < image->count; i++)
+ {
+ if (*img_data)
+ data_area[(*img_data) - 1]++;
+ img_data++;
+ }
+}
+
+void imAnalyzeMeasureCentroid(const imImage* image, const int* data_area, int region_count, float* data_cx, float* data_cy)
+{
+ imushort* img_data = (imushort*)image->data[0];
+ int* local_data_area = 0;
+
+ if (!data_area)
+ {
+ local_data_area = (int*)malloc(region_count*sizeof(int));
+ imAnalyzeMeasureArea(image, local_data_area, region_count);
+ data_area = (const int*)local_data_area;
+ }
+
+ if (data_cx) memset(data_cx, 0, region_count*sizeof(float));
+ if (data_cy) memset(data_cy, 0, region_count*sizeof(float));
+
+ for (int y = 0; y < image->height; y++)
+ {
+ int offset = y*image->width;
+
+ for (int x = 0; x < image->width; x++)
+ {
+ int region_index = img_data[offset+x];
+ if (region_index)
+ {
+ if (data_cx) data_cx[region_index-1] += (float)x;
+ if (data_cy) data_cy[region_index-1] += (float)y;
+ }
+ }
+ }
+
+ for (int i = 0; i < region_count; i++)
+ {
+ if (data_cx) data_cx[i] /= (float)data_area[i];
+ if (data_cy) data_cy[i] /= (float)data_area[i];
+ }
+
+ if (local_data_area)
+ free(local_data_area);
+}
+
+static inline double ipow(double x, int j)
+{
+ double r = 1.0;
+ for (int i = 0; i < j; i++)
+ r *= x;
+ return r;
+}
+
+static void iCalcMoment(double* cm, int px, int py, const imImage* image, const float* cx, const float* cy, int region_count)
+{
+ imushort* img_data = (imushort*)image->data[0];
+
+ memset(cm, 0, region_count*sizeof(double));
+
+ for (int y = 0; y < image->height; y++)
+ {
+ int offset = y*image->width;
+
+ for (int x = 0; x < image->width; x++)
+ {
+ int region_index = img_data[offset+x];
+ if (region_index)
+ {
+ int i = region_index-1;
+
+ if (px == 0)
+ cm[i] += ipow(y-cy[i],py);
+ else if (py == 0)
+ cm[i] += ipow(x-cx[i],px);
+ else
+ cm[i] += ipow(x-cx[i],px)*ipow(y-cy[i],py);
+ }
+ }
+ }
+}
+
+template<class T>
+static inline int IsPerimeterPoint(T* map, int width, int height, int x, int y)
+{
+ // map here points to the start of the line, even if its an invalid line.
+
+ // if outside the image, then is not a perimeter line.
+ if (x == -1 || x == width ||
+ y == -1 || y == height)
+ return 0;
+
+ T v = map[x]; // here v is image(x,y)
+ if (!v)
+ return 0;
+
+ // if touches the border, then is a perimeter line.
+ if (x == 0 || x == width-1 ||
+ y == 0 || y == height-1)
+ return 1;
+
+ // if has 4 connected neighbors, then is a perimeter line.
+ if (map[width+x] != v ||
+ map[x+1] != v ||
+ map[x-1] != v ||
+ map[-width+x] != v)
+ return 1;
+
+ return 0;
+}
+
+void imAnalyzeMeasurePrincipalAxis(const imImage* image, const int* data_area, const float* data_cx, const float* data_cy,
+ const int region_count, float* major_slope, float* major_length,
+ float* minor_slope, float* minor_length)
+{
+ int i;
+ int *local_data_area = 0;
+ float *local_data_cx = 0, *local_data_cy = 0;
+
+ if (!data_area)
+ {
+ local_data_area = (int*)malloc(region_count*sizeof(int));
+ imAnalyzeMeasureArea(image, local_data_area, region_count);
+ data_area = (const int*)local_data_area;
+ }
+
+ if (!data_cx || !data_cy)
+ {
+ if (!data_cx)
+ {
+ local_data_cx = (float*)malloc(region_count*sizeof(float));
+ data_cx = (const float*)local_data_cx;
+ }
+
+ if (!data_cy)
+ {
+ local_data_cy = (float*)malloc(region_count*sizeof(float));
+ data_cy = (const float*)local_data_cy;
+ }
+
+ if (local_data_cx && local_data_cy)
+ imAnalyzeMeasureCentroid(image, data_area, region_count, local_data_cx, local_data_cy);
+ else if (local_data_cx)
+ imAnalyzeMeasureCentroid(image, data_area, region_count, local_data_cx, NULL);
+ else if (local_data_cy)
+ imAnalyzeMeasureCentroid(image, data_area, region_count, NULL, local_data_cy);
+ }
+
+ // additional moments
+ double* cm20 = (double*)malloc(region_count*sizeof(double));
+ double* cm02 = (double*)malloc(region_count*sizeof(double));
+ double* cm11 = (double*)malloc(region_count*sizeof(double));
+
+ iCalcMoment(cm20, 2, 0, image, data_cx, data_cy, region_count);
+ iCalcMoment(cm02, 0, 2, image, data_cx, data_cy, region_count);
+ iCalcMoment(cm11, 1, 1, image, data_cx, data_cy, region_count);
+
+ float *local_major_slope = 0, *local_minor_slope = 0;
+ if (!major_slope)
+ {
+ local_major_slope = (float*)malloc(region_count*sizeof(float));
+ major_slope = local_major_slope;
+ }
+ if (!minor_slope)
+ {
+ local_minor_slope = (float*)malloc(region_count*sizeof(float));
+ minor_slope = local_minor_slope;
+ }
+
+#define RAD2DEG 57.296
+
+ // We are going to find 2 axis parameters.
+ // Axis 1 are located in quadrants 1-3
+ // Axis 2 are located in quadrants 2-4
+
+ // Quadrants
+ // 2 | 1
+ // -----
+ // 3 | 4
+
+ // line coeficients for lines that belongs to axis 1 and 2
+ float* A1 = (float*)malloc(region_count*sizeof(float));
+ float* A2 = (float*)malloc(region_count*sizeof(float));
+ float* C1 = (float*)malloc(region_count*sizeof(float));
+ float* C2 = (float*)malloc(region_count*sizeof(float));
+
+ float *slope1 = major_slope; // Use major_slope as a storage place,
+ float *slope2 = minor_slope; // and create an alias to make code clear.
+
+ for (i = 0; i < region_count; i++)
+ {
+ if (cm11[i] == 0)
+ {
+ slope1[i] = 0;
+ slope2[i] = 90;
+
+ // These should not be used
+ A1[i] = 0;
+ A2[i] = 0; // infinite
+ C1[i] = 0; // data_cy[i]
+ C2[i] = 0;
+ }
+ else
+ {
+ double b = (cm20[i] - cm02[i])/cm11[i];
+ double delta = sqrt(b*b + 4.0);
+ double r1 = (-b-delta)/2.0;
+ double r2 = (-b+delta)/2.0;
+ float a1 = (float)(atan(r1)*RAD2DEG + 90); // to avoid negative results
+ float a2 = (float)(atan(r2)*RAD2DEG + 90);
+
+ if (a1 == 180) a1 = 0;
+ if (a2 == 180) a2 = 0;
+
+ if (a1 < 90) // a1 is quadrants q1-q3
+ {
+ slope1[i] = a1;
+ slope2[i] = a2;
+ A1[i] = (float)r1;
+ A2[i] = (float)r2;
+ }
+ else // a2 is quadrants q1-q3
+ {
+ slope1[i] = a2;
+ slope2[i] = a1;
+ A1[i] = (float)r2;
+ A2[i] = (float)r1;
+ }
+
+ C1[i] = data_cy[i] - A1[i] * data_cx[i];
+ C2[i] = data_cy[i] - A2[i] * data_cx[i];
+ }
+ }
+
+ // moments are not necessary anymore
+ free(cm20); free(cm02); free(cm11);
+ cm20 = 0; cm02 = 0; cm11 = 0;
+
+ // maximum distance from a point in the perimeter to an axis in each side of the axis
+ // D1 is distance to axis 1, a and b are sides
+ float* D1a = (float*)malloc(region_count*sizeof(float));
+ float* D1b = (float*)malloc(region_count*sizeof(float));
+ float* D2a = (float*)malloc(region_count*sizeof(float));
+ float* D2b = (float*)malloc(region_count*sizeof(float));
+ memset(D1a, 0, region_count*sizeof(float));
+ memset(D1b, 0, region_count*sizeof(float));
+ memset(D2a, 0, region_count*sizeof(float));
+ memset(D2b, 0, region_count*sizeof(float));
+
+ imushort* img_data = (imushort*)image->data[0];
+ int width = image->width;
+ int height = image->height;
+ for (int y = 0; y < height; y++)
+ {
+ int offset = y*width;
+
+ for (int x = 0; x < width; x++)
+ {
+ if (IsPerimeterPoint(img_data+offset, width, height, x, y))
+ {
+ i = img_data[offset+x] - 1;
+
+ float d1, d2;
+ if (slope2[i] == 90)
+ {
+ d2 = y - data_cy[i]; // I ckecked this many times, looks odd but it is correct.
+ d1 = x - data_cx[i];
+ }
+ else
+ {
+ d1 = A1[i]*x - y + C1[i];
+ d2 = A2[i]*x - y + C2[i];
+ }
+
+ if (d1 < 0)
+ {
+ d1 = (float)fabs(d1);
+ if (d1 > D1a[i])
+ D1a[i] = d1;
+ }
+ else
+ {
+ if (d1 > D1b[i])
+ D1b[i] = d1;
+ }
+
+ if (d2 < 0)
+ {
+ d2 = (float)fabs(d2);
+ if (d2 > D2a[i])
+ D2a[i] = d2;
+ }
+ else
+ {
+ if (d2 > D2b[i])
+ D2b[i] = d2;
+ }
+ }
+ }
+ }
+
+ for (i = 0; i < region_count; i++)
+ {
+ float AB1 = (float)sqrt(A1[i]*A1[i] + 1);
+ float AB2 = (float)sqrt(A2[i]*A2[i] + 1);
+
+ float D1 = (D1a[i] + D1b[i]) / AB1;
+ float D2 = (D2a[i] + D2b[i]) / AB2;
+
+ if (D1 < D2) // Major Axis in 2-4 quadrants
+ {
+ // now remember that we did an alias before
+ // slope1 -> major_slope
+ // slope2 -> minor_slope
+
+ float tmp = major_slope[i];
+ major_slope[i] = minor_slope[i];
+ minor_slope[i] = tmp;
+
+ if (minor_length) minor_length[i] = D1;
+ if (major_length) major_length[i] = D2;
+ }
+ else
+ {
+ if (minor_length) minor_length[i] = D2;
+ if (major_length) major_length[i] = D1;
+ }
+ }
+
+ if (local_major_slope) free(local_major_slope);
+ if (local_minor_slope) free(local_minor_slope);
+ if (local_data_area) free(local_data_area);
+ if (local_data_cx) free(local_data_cx);
+ if (local_data_cy) free(local_data_cy);
+
+ free(A1);
+ free(A2);
+ free(C1);
+ free(C2);
+
+ free(D1b);
+ free(D2b);
+ free(D1a);
+ free(D2a);
+}
+
+void imAnalyzeMeasureHoles(const imImage* image, int connect, int* count_data, int* area_data, float* perim_data)
+{
+ int i;
+ imImage *inv_image = imImageCreate(image->width, image->height, IM_BINARY, IM_BYTE);
+ imbyte* inv_data = (imbyte*)inv_image->data[0];
+ imushort* img_data = (imushort*)image->data[0];
+
+ // finds the holes in the inverted image
+ for (i = 0; i < image->count; i++)
+ {
+ if (*img_data)
+ *inv_data = 0;
+ else
+ *inv_data = 1;
+
+ img_data++;
+ inv_data++;
+ }
+
+ imImage *holes_image = imImageClone(image);
+ if (!holes_image)
+ return;
+
+ int holes_count = imAnalyzeFindRegions(inv_image, holes_image, connect, 0);
+ imImageDestroy(inv_image);
+
+ if (!holes_count)
+ {
+ imImageDestroy(holes_image);
+ return;
+ }
+
+ // measure the holes area
+ int* holes_area = (int*)malloc(holes_count*sizeof(int));
+ imAnalyzeMeasureArea(holes_image, holes_area, holes_count);
+
+ float* holes_perim = 0;
+ if (perim_data)
+ {
+ holes_perim = (float*)malloc(holes_count*sizeof(int));
+ imAnalyzeMeasurePerimeter(holes_image, holes_perim, holes_count);
+ }
+
+ imushort* holes_data = (imushort*)holes_image->data[0];
+ img_data = (imushort*)image->data[0];
+
+ // holes do not touch the border
+ for (int y = 1; y < image->height-1; y++)
+ {
+ int offset_up = (y+1)*image->width;
+ int offset = y*image->width;
+ int offset_dw = (y-1)*image->width;
+
+ for (int x = 1; x < image->width-1; x++)
+ {
+ int hole_index = holes_data[offset+x];
+
+ if (hole_index && holes_area[hole_index-1]) // a hole not yet used
+ {
+ // if the hole has not been used,
+ // it is the first time we encounter a pixel of this hole.
+ // then it is a pixel from the hole border.
+ // now find which region this hole is inside.
+ // a 4 connected neighbour is necessarilly a valid region or 0.
+
+ int region_index = 0;
+ if (img_data[offset_up + x]) region_index = img_data[offset_up + x];
+ else if (img_data[offset + x+1]) region_index = img_data[offset + x+1];
+ else if (img_data[offset + x-1]) region_index = img_data[offset + x-1];
+ else if (img_data[offset_dw+x]) region_index = img_data[offset_dw+x];
+
+ if (!region_index) continue;
+
+ if (count_data) count_data[region_index-1]++;
+ if (area_data) area_data[region_index-1] += holes_area[hole_index-1];
+ if (perim_data) perim_data[region_index-1] += holes_perim[hole_index-1];
+ holes_area[hole_index-1] = 0; // mark hole as used
+ }
+ }
+ }
+
+ if (holes_perim) free(holes_perim);
+ free(holes_area);
+ imImageDestroy(holes_image);
+}
+
+template<class T>
+static void DoPerimeterLine(T* map, T* new_map, int width, int height)
+{
+ int x, y, offset;
+
+ for (y = 0; y < height; y++)
+ {
+ offset = y*width;
+
+ for (x = 0; x < width; x++)
+ {
+ if (IsPerimeterPoint(map+offset, width, height, x, y))
+ new_map[offset+x] = map[offset+x];
+ else
+ new_map[offset+x] = 0;
+ }
+ }
+}
+
+void imProcessPerimeterLine(const imImage* src_image, imImage* dst_image)
+{
+ switch(src_image->data_type)
+ {
+ case IM_BYTE:
+ DoPerimeterLine((imbyte*)src_image->data[0], (imbyte*)dst_image->data[0], src_image->width, src_image->height);
+ break;
+ case IM_USHORT:
+ DoPerimeterLine((imushort*)src_image->data[0], (imushort*)dst_image->data[0], src_image->width, src_image->height);
+ break;
+ case IM_INT:
+ DoPerimeterLine((int*)src_image->data[0], (int*)dst_image->data[0], src_image->width, src_image->height);
+ break;
+ }
+}
+
+/* Perimeter Templates idea based in
+ Parker, Pratical Computer Vision Using C
+
+For 1.414 (sqrt(2)/2 + sqrt(2)/2) [1]:
+ 1 0 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 1
+ 0 x 0 0 x 0 0 x 0 0 x 0 0 x 0 0 x 0
+ 0 0 1 1 0 0 1 0 0 0 0 1 1 0 1 0 0 0
+ 129 36 132 33 5 160
+
+For 1.207 (sqrt(2)/2 + 1.0/2) [2]:
+ 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0
+ 1 x 0 1 x 0 0 x 0 0 x 0 0 x 0 0 x 0 0 x 1 0 x 1
+ 0 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0
+ 17 48 68 65 130 34 12 136
+
+ 0 0 0 1 0 0 1 1 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0
+ 1 x 0 1 x 0 0 x 0 0 x 0 0 x 1 0 x 1 0 x 0 0 x 0
+ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 0
+ 20 144 192 96 40 9 3 6
+
+For 1.0 (1.0/2 + 1.0/2) [0]:
+ 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0
+ 1 x 1 0 x 0 1 x 0 0 x 1 1 x 0 0 x 1
+ 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0
+ 24 66 18 10 80 72
+
+For 0.707 (sqrt(2)/2) [3]:
+ 1 0 0 0 0 1 0 0 0 0 0 0
+ 0 x 0 0 x 0 0 x 0 0 x 0 (For Line Length)
+ 0 0 0 0 0 0 0 0 1 1 0 0
+ 128 32 1 4
+
+For 0.5 (1.0/2) [4]:
+ 0 1 0 0 0 0 0 0 0 0 0 0
+ 0 x 0 0 x 1 0 x 0 1 x 0 (For Line Length)
+ 0 0 0 0 0 0 0 1 0 0 0 0
+ 64 8 2 16
+
+*/
+static void iInitPerimTemplate(imbyte *templ, float *v)
+{
+ memset(templ, 0, 256);
+
+ templ[129] = 1;
+ templ[36] = 1;
+ templ[132] = 1;
+ templ[33] = 1;
+ templ[5] = 1;
+ templ[160] = 1;
+
+ templ[17] = 2;
+ templ[48] = 2;
+ templ[68] = 2;
+ templ[65] = 2;
+ templ[130] = 2;
+ templ[34] = 2;
+ templ[12] = 2;
+ templ[136] = 2;
+ templ[20] = 2;
+ templ[144] = 2;
+ templ[192] = 2;
+ templ[96] = 2;
+ templ[40] = 2;
+ templ[9] = 2;
+ templ[3] = 2;
+ templ[6] = 2;
+
+ templ[24] = 0;
+ templ[66] = 0;
+ templ[18] = 0;
+ templ[10] = 0;
+ templ[80] = 0;
+ templ[72] = 0;
+
+ templ[128] = 3;
+ templ[32] = 3;
+ templ[1] = 3;
+ templ[4] = 3;
+
+ templ[64] = 4;
+ templ[8] = 4;
+ templ[2] = 4;
+ templ[16] = 4;
+
+const float DT_SQRT2 = 1.414213562373f;
+const float DT_SQRT2D2 = 0.707106781187f;
+
+ v[1] = DT_SQRT2;
+ v[2] = DT_SQRT2D2 + 0.5f;
+ v[0] = 1.0f;
+ v[3] = DT_SQRT2D2;
+ v[4] = 0.5f;
+}
+
+void imAnalyzeMeasurePerimeter(const imImage* image, float* perim_data, int region_count)
+{
+ static imbyte templ[256];
+ static float vt[5];
+ static int first = 1;
+ if (first)
+ {
+ iInitPerimTemplate(templ, vt);
+ first = 0;
+ }
+
+ imushort* map = (imushort*)image->data[0];
+
+ memset(perim_data, 0, region_count*sizeof(int));
+
+ int width = image->width;
+ int height = image->height;
+ for (int y = 0; y < height; y++)
+ {
+ int offset = y*image->width;
+
+ for (int x = 0; x < width; x++)
+ {
+ if (IsPerimeterPoint(map+offset, width, height, x, y))
+ {
+ int T = 0;
+
+ // check the 8 neighboors if they belong to the perimeter
+ if (IsPerimeterPoint(map+offset+width, width, height, x-1, y+1))
+ T |= 0x01;
+ if (IsPerimeterPoint(map+offset+width, width, height, x, y+1))
+ T |= 0x02;
+ if (IsPerimeterPoint(map+offset+width, width, height, x+1, y+1))
+ T |= 0x04;
+
+ if (IsPerimeterPoint(map+offset, width, height, x-1, y))
+ T |= 0x08;
+ if (IsPerimeterPoint(map+offset, width, height, x+1, y))
+ T |= 0x10;
+
+ if (IsPerimeterPoint(map+offset-width, width, height, x-1, y-1))
+ T |= 0x20;
+ if (IsPerimeterPoint(map+offset-width, width, height, x, y-1))
+ T |= 0x40;
+ if (IsPerimeterPoint(map+offset-width, width, height, x+1, y-1))
+ T |= 0x80;
+
+ if (T)
+ perim_data[map[offset+x] - 1] += vt[templ[T]];
+ }
+ }
+ }
+}
+
+/* Perimeter Area Templates
+
+For "1.0" (0):
+
+ 1 1 1
+ 1 x 1
+ 1 1 1
+ 255
+
+For "0.75" (1):
+
+ 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1
+ 1 x 1 1 x 1 1 x 1 1 x 1 0 x 1 1 x 0 1 x 1 1 x 1
+ 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1
+ 251 254 127 223 239 247 253 191
+
+For "0.625" (2):
+
+ 1 1 1 0 0 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 0
+ 1 x 1 1 x 1 0 x 1 1 x 0 0 x 1 1 x 0 1 x 1 1 x 1
+ 0 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 0 0 1 1 1
+ 249 63 111 215 235 246 252 159
+
+For "0.5" (3):
+
+ 0 0 0 0 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 0 0 1 1 1
+ 1 x 1 0 x 1 1 x 1 1 x 0 0 x 1 0 x 1 1 x 0 1 x 0
+ 1 1 1 0 1 1 0 0 0 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0
+ 31 107 248 214 233 47 151 244
+
+For "0.375" (4):
+
+ 0 0 0 1 1 1 1 1 0 0 1 1 1 0 0 0 0 1 0 0 0 1 1 1
+ 1 x 0 1 x 0 1 x 0 0 x 1 1 x 0 0 x 1 0 x 1 0 x 1
+ 1 1 1 0 0 0 1 0 0 0 0 1 1 1 0 0 1 1 1 1 1 0 0 0
+ 23 240 212 105 150 43 15 232
+
+For "0.25" (5):
+
+ 0 0 0 0 0 0 1 1 0 0 1 1 1 0 0 0 0 1 0 0 0 1 1 1
+ 1 x 0 0 x 1 1 x 0 0 x 1 1 x 0 0 x 1 0 x 0 0 x 0
+ 1 1 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 1 1 1 1 0 0 0
+ 22 11 208 104 148 41 7 224
+
+For "0.125" (6):
+
+ 0 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 1
+ 1 x 0 0 x 0 0 x 0 0 x 1 1 x 0 0 x 1 0 x 0 0 x 0
+ 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0
+ 20 3 192 40 144 9 6 96
+
+*/
+static void iInitPerimAreaTemplate(imbyte *templ, float *v)
+{
+ memset(templ, 0, 256);
+
+ templ[255] = 0;
+
+ templ[251] = 1;
+ templ[254] = 1;
+ templ[127] = 1;
+ templ[223] = 1;
+ templ[239] = 1;
+ templ[247] = 1;
+ templ[253] = 1;
+ templ[191] = 1;
+
+ templ[249] = 2;
+ templ[63] = 2;
+ templ[111] = 2;
+ templ[215] = 2;
+ templ[235] = 2;
+ templ[246] = 2;
+ templ[252] = 2;
+ templ[159] = 2;
+
+ templ[31] = 3;
+ templ[107] = 3;
+ templ[248] = 3;
+ templ[214] = 3;
+ templ[233] = 3;
+ templ[47] = 3;
+ templ[151] = 3;
+ templ[244] = 3;
+
+ templ[23] = 4;
+ templ[240] = 4;
+ templ[212] = 4;
+ templ[105] = 4;
+ templ[150] = 4;
+ templ[43] = 4;
+ templ[15] = 4;
+ templ[232] = 4;
+
+ templ[22] = 5;
+ templ[11] = 5;
+ templ[208] = 5;
+ templ[104] = 5;
+ templ[148] = 5;
+ templ[41] = 5;
+ templ[7] = 5;
+ templ[224] = 5;
+
+ templ[20] = 6;
+ templ[3] = 6;
+ templ[192] = 6;
+ templ[40] = 6;
+ templ[144] = 6;
+ templ[9] = 6;
+ templ[6] = 6;
+ templ[96] = 6;
+
+ v[0] = 1.0f;
+ v[1] = 0.75f;
+ v[2] = 0.625f;
+ v[3] = 0.5f;
+ v[4] = 0.375f;
+ v[5] = 0.25f;
+ v[6] = 0.125f;
+}
+
+void imAnalyzeMeasurePerimArea(const imImage* image, float* area_data)
+{
+ static imbyte templ[256];
+ static float vt[7];
+ static int first = 1;
+ if (first)
+ {
+ iInitPerimAreaTemplate(templ, vt);
+ first = 0;
+ }
+
+ imushort* map = (imushort*)image->data[0];
+
+ int width = image->width;
+ int height = image->height;
+ for (int y = 0; y < height; y++)
+ {
+ int offset_up = (y+1)*width;
+ int offset = y*width;
+ int offset_dw = (y-1)*width;
+
+ for (int x = 0; x < width; x++)
+ {
+ imushort v = map[offset+x];
+ if (v)
+ {
+ int T = 0;
+ if (x>0 && y<height-1 && map[offset_up + x-1] == v) T |= 0x01;
+ if (y<height-1 && map[offset_up + x ] == v) T |= 0x02;
+ if (x<width-1 && y<height-1 && map[offset_up + x+1] == v) T |= 0x04;
+ if (x>0 && map[offset + x-1] == v) T |= 0x08;
+ if (x<width-1 && map[offset + x+1] == v) T |= 0x10;
+ if (x>0 && y>0 && map[offset_dw + x-1] == v) T |= 0x20;
+ if (y>0 && map[offset_dw + x ] == v) T |= 0x40;
+ if (x<width-1 && y>0 && map[offset_dw + x+1] == v) T |= 0x80;
+
+ if (T)
+ area_data[v-1] += vt[templ[T]];
+ }
+ }
+ }
+}
+
+void imProcessPrune(const imImage* image, imImage* NewImage, int connect, int start_size, int end_size)
+{
+ imImage *region_image = imImageCreate(image->width, image->height, IM_GRAY, IM_USHORT);
+ if (!region_image)
+ return;
+
+ int region_count = imAnalyzeFindRegions(image, region_image, connect, 1);
+ if (!region_count)
+ {
+ imImageClear(NewImage);
+ imImageDestroy(region_image);
+ return;
+ }
+
+ int* area_data = (int*)malloc(region_count*sizeof(int));
+ imAnalyzeMeasureArea(region_image, area_data, region_count);
+
+ imushort* region_data = (imushort*)region_image->data[0];
+ imbyte* img_data = (imbyte*)NewImage->data[0];
+
+ for (int i = 0; i < image->count; i++)
+ {
+ if (*region_data)
+ {
+ int area = area_data[(*region_data) - 1];
+ if (area < start_size || (end_size && area > end_size))
+ *img_data = 0;
+ else
+ *img_data = 1;
+ }
+ else
+ *img_data = 0;
+
+ region_data++;
+ img_data++;
+ }
+
+ free(area_data);
+ imImageDestroy(region_image);
+}
+
+void imProcessFillHoles(const imImage* image, imImage* NewImage, int connect)
+{
+ // finding regions in the inverted image will isolate only the holes.
+ imProcessNegative(image, NewImage);
+
+ imImage *region_image = imImageCreate(image->width, image->height, IM_GRAY, IM_USHORT);
+ if (!region_image)
+ return;
+
+ int holes_count = imAnalyzeFindRegions(NewImage, region_image, connect, 0);
+ if (!holes_count)
+ {
+ imImageCopy(image, NewImage);
+ imImageDestroy(region_image);
+ return;
+ }
+
+ imushort* region_data = (imushort*)region_image->data[0];
+ imbyte* dst_data = (imbyte*)NewImage->data[0];
+
+ for (int i = 0; i < image->count; i++)
+ {
+ if (*region_data)
+ *dst_data = 1;
+ else
+ *dst_data = !(*dst_data); // Fix negative data.
+
+ region_data++;
+ dst_data++;
+ }
+
+ imImageDestroy(region_image);
+}