[Scipy-svn] r3499 - in trunk/scipy/ndimage: . segment tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Nov 7 12:54:00 EST 2007
Author: tom.waite
Date: 2007-11-07 11:53:37 -0600 (Wed, 07 Nov 2007)
New Revision: 3499
Added:
trunk/scipy/ndimage/segment/
trunk/scipy/ndimage/segment/Segmenter_EXT.c
trunk/scipy/ndimage/segment/Segmenter_IMPL.c
trunk/scipy/ndimage/segment/__init__.py
trunk/scipy/ndimage/segment/ndImage_Segmenter_structs.h
trunk/scipy/ndimage/segment/setup.py
trunk/scipy/ndimage/tests/slice112.raw
trunk/scipy/ndimage/tests/test_segmenter.py
Log:
Add beginning of Tom Waite's ndimage.segment package.
Added: trunk/scipy/ndimage/segment/Segmenter_EXT.c
===================================================================
--- trunk/scipy/ndimage/segment/Segmenter_EXT.c 2007-11-05 23:13:01 UTC (rev 3498)
+++ trunk/scipy/ndimage/segment/Segmenter_EXT.c 2007-11-07 17:53:37 UTC (rev 3499)
@@ -0,0 +1,467 @@
+#include "ndImage_Segmenter_structs.h"
+#include "Python.h"
+#include "numpy/arrayobject.h"
+
+static PyObject *NDI_Segmenter_CannyEdges(PyObject *self, PyObject *args)
+{
+
+ double sigma;
+ double cannyLow;
+ double cannyHigh;
+ double BPHigh;
+ int lowThreshold;
+ int highThreshold;
+ int apearture;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int mode;
+ int groups;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_Parse(args, "(dddiiidiO)", &sigma, &cannyLow, &cannyHigh, &mode, &lowThreshold, &highThreshold,
+ &BPHigh, &apearture, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ itype = 4;
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+ if(!NI_CannyEdges(num, (int)dims[0], (int)dims[1], sigma, cannyLow, cannyHigh, mode, lowThreshold,
+ highThreshold, BPHigh, apearture, fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups);
+
+}
+
+static PyObject *NDI_Segmenter_SobelEdges(PyObject *self, PyObject *args)
+{
+
+ double sobelLow;
+ double BPHigh;
+ int lowThreshold;
+ int highThreshold;
+ int apearture;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int groups;
+ int mode;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_Parse(args, "(diiidiO)", &sobelLow, &mode, &lowThreshold, &highThreshold, &BPHigh, &apearture, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // this is int type and hard-wirred. pass this in from Python code
+ itype = 4; // unsigned short
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+
+ if(!NI_SobelEdges(num, (int)dims[0], (int)dims[1], sobelLow, mode, lowThreshold, highThreshold, BPHigh, apearture,
+ fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups-1);
+
+}
+
+
+
+static PyObject *NDI_Segmenter_ShenCastanEdges(PyObject *self, PyObject *args)
+{
+ int window;
+ int lowThreshold;
+ int highThreshold;
+ double ShenCastanLow;
+ double b;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ int groups;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ if(!PyArg_Parse(args, "(ddiiiO)", &ShenCastanLow, &b, &window, &lowThreshold, &highThreshold, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // this is int type and hard-wirred. pass this in from Python code
+ itype = 4; // unsigned short
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+ if(!NI_ShenCastanEdges(num, (int)dims[0], (int)dims[1], b, ShenCastanLow, window, lowThreshold, highThreshold,
+ fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups-1);
+
+}
+
+static PyObject *NDI_Segmenter_GetObjectStats(PyObject *self, PyObject *args)
+{
+
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ unsigned short *fP1;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ objStruct *myData;
+
+ if(!PyArg_Parse(args, "(OO)", &iArray, &nArray))
+ goto exit;
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(nArray))
+ goto exit;
+
+ //
+ // PyArray_ContiguousFromObject or PyArray_ContiguousFromAny to be explored
+ // for non-contiguous
+ //
+
+
+ // pointer to the edge-labeled image
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+ fP1 = (unsigned short *)PyArray_DATA(iArray);
+
+ // the object descriptor array that was allocated from numpy
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ myData = (objStruct*)PyArray_DATA(nArray);
+
+ if(!NI_GetObjectStats((int)dims[0], (int)dims[1], (int)objNumber[0], fP1, myData))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_MorphoThinFilt(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ unsigned short *fP1;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OO)", &iArray, &nArray))
+ goto exit;
+
+ fP1 = (unsigned short *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ if(!NI_ThinFilter(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_BuildBoundary(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ unsigned short *fP1;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OO)", &iArray, &nArray))
+ goto exit;
+
+ fP1 = (unsigned short *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+ //
+ // this is int type and hard-wirred. pass this in from Python code
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ //
+ // pass in ROI list and labeled edges
+ // return an augmented ROI list
+ // replace the edgeImage with maskImage
+ //
+ if(!NI_BuildBoundary(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+
+static PyObject *NDI_Segmenter_VoxelMeasures(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ PyObject *eArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OOO)", &iArray, &eArray, &nArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // eArray and iArray are same dims
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ //
+ // pass in ROI list and labeled edges
+ // return an augmented ROI list
+ // replace the edgeImage with maskImage
+ //
+
+ if(!NI_VoxelMeasures(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, fP2, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_TextureMeasures(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ PyObject *eArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OOO)", &iArray, &eArray, &nArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // eArray and iArray are same dims
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ //
+ // pass in ROI list and labeled edges
+ // return an augmented ROI list
+ // replace the edgeImage with maskImage
+ //
+
+ if(!NI_TextureMeasures(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, fP2, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_RegionGrow(PyObject *self, PyObject *args)
+{
+
+ int lowThreshold;
+ int highThreshold;
+ int closeWindow;
+ int openWindow;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int groups;
+ int mode;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_Parse(args, "(iiiiO)", &lowThreshold, &highThreshold, &closeWindow, &openWindow, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // this is int type and hard-wirred. pass this in from Python code
+ itype = 4; // unsigned short
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+
+ if(!NI_RegionGrow(num, (int)dims[0], (int)dims[1], lowThreshold, highThreshold, closeWindow, openWindow,
+ fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups-1);
+
+}
+
+static PyMethodDef NDI_SegmenterMethods[] =
+{
+ { "canny_edges", NDI_Segmenter_CannyEdges, METH_VARARGS },
+ { "shen_castan_edges", NDI_Segmenter_ShenCastanEdges, METH_VARARGS },
+ { "sobel_edges", NDI_Segmenter_SobelEdges, METH_VARARGS },
+ { "get_object_stats", NDI_Segmenter_GetObjectStats, METH_VARARGS },
+ { "morpho_thin_filt", NDI_Segmenter_MorphoThinFilt, METH_VARARGS },
+ { "build_boundary", NDI_Segmenter_BuildBoundary, METH_VARARGS },
+ { "voxel_measures", NDI_Segmenter_VoxelMeasures, METH_VARARGS },
+ { "texture_measures", NDI_Segmenter_TextureMeasures, METH_VARARGS },
+ { "region_grow", NDI_Segmenter_RegionGrow, METH_VARARGS },
+ { NULL, NULL },
+};
+
+void init_segmenter()
+{
+ Py_InitModule("_segmenter", NDI_SegmenterMethods);
+ import_array();
+}
+
+/*
+static PyMethodDef NDI_SegmenterMethods[] =
+{
+ { "CannyEdges", NDI_Segmenter_CannyEdges, METH_VARARGS },
+ { "ShenCastanEdges", NDI_Segmenter_ShenCastanEdges, METH_VARARGS },
+ { "SobelEdges", NDI_Segmenter_SobelEdges, METH_VARARGS },
+ { "GetObjectStats", NDI_Segmenter_GetObjectStats, METH_VARARGS },
+ { "MorphoThinFilt", NDI_Segmenter_MorphoThinFilt, METH_VARARGS },
+ { "BuildBoundary", NDI_Segmenter_BuildBoundary, METH_VARARGS },
+ { "VoxelMeasures", NDI_Segmenter_VoxelMeasures, METH_VARARGS },
+ { "TextureMeasures", NDI_Segmenter_TextureMeasures, METH_VARARGS },
+ { "RegionGrow", NDI_Segmenter_RegionGrow, METH_VARARGS },
+ { NULL, NULL },
+};
+
+void initNDI_Segmenter()
+{
+ Py_InitModule("NDI_Segmenter", NDI_SegmenterMethods);
+ import_array();
+}
+*/
+
Added: trunk/scipy/ndimage/segment/Segmenter_IMPL.c
===================================================================
--- trunk/scipy/ndimage/segment/Segmenter_IMPL.c 2007-11-05 23:13:01 UTC (rev 3498)
+++ trunk/scipy/ndimage/segment/Segmenter_IMPL.c 2007-11-07 17:53:37 UTC (rev 3499)
@@ -0,0 +1,2976 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "ndImage_Segmenter_structs.h"
+
+// these are for this standalone and come out with the full build
+//
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+#define FALSE 0
+#define TRUE 1
+
+int NI_GetObjectStats(int rows, int cols, int numberObjects, unsigned short *labeledEdges, objStruct objectMetrics[]){
+
+ int i, j, k, m;
+ int offset;
+ int count;
+ int LowX;
+ int LowY;
+ int HighX;
+ int HighY;
+ int status;
+ float centerX;
+ float centerY;
+
+ for(k = 1; k < numberObjects; ++k){
+ offset = cols;
+ LowX = 32767;
+ LowY = 32767;
+ HighX = 0;
+ HighY = 0;
+ count = 0;
+ centerX = (float)0.0;
+ centerY = (float)0.0;
+ for(i = 1; i < (rows-1); ++i){
+ for(j = 1; j < (cols-1); ++j){
+ m = labeledEdges[offset+j];
+ if(k == m){
+ if(i < LowY) LowY = i;
+ if(j < LowX) LowX = j;
+ if(i > HighY) HighY = i;
+ if(j > HighX) HighX = j;
+ centerX += (float)j;
+ centerY += (float)i;
+ ++count;
+ }
+ }
+ offset += cols;
+ }
+ // the bounding box for the 2D blob
+ objectMetrics[k-1].L = LowX;
+ objectMetrics[k-1].R = HighX;
+ objectMetrics[k-1].B = LowY;
+ objectMetrics[k-1].T = HighY;
+ objectMetrics[k-1].Area = count;
+ objectMetrics[k-1].cX = centerX/(float)count;
+ objectMetrics[k-1].cY = centerY/(float)count;
+ objectMetrics[k-1].Label = k;
+ }
+
+ status = numberObjects;
+ return status;
+
+}
+
+
+void buildKernel(double BPHigh, int HalfFilterTaps, int apearture, float *kernel){
+
+ int i, j;
+ float r, t1, t2, t3, t4;
+ float LC, HC, tLOW, tHIGH;
+ float pi = (float)3.14159, rad = (float)0.01745;
+
+ LC = (float)0.0;
+ HC = BPHigh * rad;
+ t2 = (float)2.0*pi;
+ t1 = (float)2.0*HalfFilterTaps + (float)1.0;
+ //
+ // build the Filter Kernel
+ // the kernel starts at 1 only because it is linked to the internal filter2D routine
+ // the code is not a Fortran code
+ //
+ j = 1;
+ for(i = -HalfFilterTaps; i <= HalfFilterTaps; ++i){
+ r = (float)i;
+ if(r == (float)0.0){
+ tLOW = LC;
+ tHIGH = HC;
+ }
+ else{
+ tLOW = (float)(sin(r*LC))/r;
+ tHIGH = (float)(sin(r*HC))/r;
+ }
+ t3 = (float)0.54 + (float)0.46*((float)cos(r*t2/t1));
+ t4 = t3*(tHIGH-tLOW);
+ kernel[j++] = t4;
+ }
+
+ // normalize the kernel so unity gain (as is LP filter this is easy)
+ t1 = (float)0.0;
+ for(j = 1; j <= apearture; ++j){
+ t1 += kernel[j];
+ }
+ for(j = 1; j <= apearture; ++j){
+ kernel[j] /= t1;
+ }
+
+ t1 = (float)0.0;
+ for(j = 1; j <= apearture; ++j){
+ t1 += kernel[j];
+ }
+ return;
+}
+
+void filter2D(int HalfFilterTaps, int rows, int cols, int lowThreshold, int highThreshold, float *kernel, double *Image){
+
+ int i, j, k, n, num1;
+ int offset;
+ float sum, value;
+ float buffer[1024];
+
+ num1 = HalfFilterTaps + 1;
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ // copy image row to local buffer
+ for(j = 0; j < cols; ++j){
+ buffer[num1+j] = Image[offset+j];
+ }
+ // constant pad the ends of the buffer
+ for(j = 0; j < num1; ++j){
+ buffer[j] = buffer[num1];
+ }
+ for(j = cols+num1; j < cols+2*num1; ++j){
+ buffer[j] = buffer[cols-1+num1];
+ }
+
+ // Perform Symmetric Convolution in the X dimension.
+ for(n = 0, j = num1; j < (cols+num1); ++j, ++n){
+ sum = buffer[j] * kernel[num1];
+ for(k = 1; k < num1; ++k){
+ sum += kernel[num1-k] * (buffer[j+k] + buffer[j-k]);
+ }
+ Image[offset+n] = sum;
+ }
+ offset += cols;
+ }
+
+ offset = 0;
+ for(i = 0; i < cols; ++i){
+ // copy image column to local buffer
+ offset = 0;
+ for(j = 0; j < rows; ++j){
+ buffer[num1+j] = Image[offset+i];
+ offset += cols;
+ }
+ // constant pad the ends of the buffer
+ for(j = 0; j < num1; ++j){
+ buffer[j] = buffer[num1];
+ }
+ for(j = rows+num1; j < rows+2*num1; ++j){
+ buffer[j] = buffer[rows-1+num1];
+ }
+
+ // Perform Symmetric Convolution in the Y dimension.
+ offset = 0;
+ for(j = num1; j < (rows+num1); ++j){
+ sum = buffer[j] * kernel[num1];
+ for(k = 1; k < num1; ++k){
+ sum += kernel[num1-k] * (buffer[j+k] + buffer[j-k]);
+ }
+ Image[offset+i] = sum;
+ offset += cols;
+ }
+ }
+
+ // threshold the image
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = Image[offset+j];
+ if(value < (float)lowThreshold) value = (float)0.0;
+ if(value > (float)highThreshold) value = (float)0.0;
+ Image[offset+j] = value;
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+void doPreProcess(int samples, int rows, int cols, double *rawImage, double BPHigh, int apearture, int lowThreshold, int highThreshold){
+
+ //
+ // 2D low pass filter using bisinc and threshold
+ // this specific example is on cardiac CT and focuses on segmenting the
+ // aorta and blood-filled chambers. for MRI the threshold will be different
+ //
+
+ float *kernel;
+ int HalfFilterTaps = (apearture-1)/2;
+ kernel = calloc(apearture+16, sizeof(float));
+
+ buildKernel(BPHigh, HalfFilterTaps, apearture, kernel);
+ filter2D(HalfFilterTaps, rows, cols, lowThreshold, highThreshold, kernel, rawImage);
+
+ free(kernel);
+
+ return;
+
+}
+
+
+int ConnectedEdgePoints(int rows, int cols, unsigned short *connectedEdges){
+
+ int i, j, k, l, m;
+ int offset;
+ int Label;
+ int Classes[4096];
+ bool NewLabel;
+ bool Change;
+ unsigned short T[12];
+
+ //
+ // connected components labeling. pixels touch within 3x3 mask for edge connectedness.
+ //
+ Label = 1;
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(connectedEdges[offset+j] == 1){
+ connectedEdges[offset+j] = Label++;
+ }
+ }
+ offset += cols;
+ }
+
+ while(1){
+ Change = FALSE;
+ //
+ // TOP-DOWN Pass for labeling
+ //
+ offset = cols;
+ for(i = 1; i < rows-1; ++i){
+ for(j = 1; j < cols-1; ++j){
+ if(connectedEdges[offset+j] != 0){
+ T[0] = connectedEdges[offset+j];
+ T[1] = connectedEdges[offset+j+1];
+ T[2] = connectedEdges[offset-cols+j+1];
+ T[3] = connectedEdges[offset-cols+j];
+ T[4] = connectedEdges[offset-cols+j-1];
+ T[5] = connectedEdges[offset+j-1];
+ T[6] = connectedEdges[offset+cols+j-1];
+ T[7] = connectedEdges[offset+cols+j];
+ T[8] = connectedEdges[offset+cols+j+1];
+ m = T[0];
+ for(l = 1; l < 9; ++l){
+ if(T[l] != 0){
+ if(T[l] < m) m = T[l];
+ }
+ }
+ if(m != connectedEdges[offset+j]){
+ Change = TRUE;
+ connectedEdges[offset+j] = m;
+ }
+ }
+ }
+ offset += cols;
+ }
+ //
+ // BOTTOM-UP Pass for labeling
+ //
+ offset = (rows-1)*cols;
+ for(i = (rows-1); i > 1; --i){
+ for(j = (cols-1); j > 1; --j){
+ if(connectedEdges[offset+j] != 0){
+ T[0] = connectedEdges[offset+j];
+ T[1] = connectedEdges[offset+j+1];
+ T[2] = connectedEdges[offset-cols+j+1];
+ T[3] = connectedEdges[offset-cols+j];
+ T[4] = connectedEdges[offset-cols+j-1];
+ T[5] = connectedEdges[offset+j-1];
+ T[6] = connectedEdges[offset+cols+j-1];
+ T[7] = connectedEdges[offset+cols+j];
+ T[8] = connectedEdges[offset+cols+j+1];
+ m = T[0];
+ for(l = 1; l < 9; ++l){
+ if(T[l] != 0){
+ if(T[l] < m) m = T[l];
+ }
+ }
+ if(m != connectedEdges[offset+j]){
+ Change = TRUE;
+ connectedEdges[offset+j] = m;
+ }
+ }
+ }
+ offset -= cols;
+ }
+ if(!Change) break;
+ } // end while loop
+
+ Classes[0] = 0;
+ Label = 1;
+ offset = cols;
+ for(i = 1; i < (rows-1); ++i){
+ for(j = 1; j < (cols-1); ++j){
+ m = connectedEdges[offset+j];
+ if(m > 0){
+ NewLabel = TRUE;
+ for(k = 1; k < Label; ++k){
+ if(Classes[k] == m) NewLabel = FALSE;
+ }
+ if(NewLabel){
+ Classes[Label++] = m;
+ if(Label > 4000){
+ return 0; // too many labeled regions. this is a pathology in the image slice
+ }
+ }
+ }
+ }
+ offset += cols;
+ }
+
+ //
+ // re-label the connected blobs in continuous label order
+ //
+ offset = cols;
+ for(i = 1; i < (rows-1); ++i){
+ for(j = 1; j < (cols-1); ++j){
+ m = connectedEdges[offset+j];
+ if(m > 0){
+ for(k = 1; k < Label; ++k){
+ if(Classes[k] == m){
+ connectedEdges[offset+j] = (unsigned short)k;
+ break;
+ }
+ }
+ }
+ }
+ offset += cols;
+ }
+
+ return Label;
+}
+
+float magnitude(float X, float Y){
+
+ return (float)sqrt(X*X + Y*Y);
+}
+
+int traceEdge(int i, int j, int rows, int cols, double cannyLow, float *magImage, float *HYSImage){
+
+ int n, m;
+ int ptr;
+ int flag;
+
+ ptr = i * cols;
+ if(HYSImage[ptr+j] == (float)0.0){
+ //
+ // this point is above high threshold
+ //
+ HYSImage[ptr+j] = (float)1.0;
+ flag = 0;
+ for(n = -1; n <= 1; ++n){
+ for(m = -1; m <= 1; ++m){
+ if(n == 0 && m == 0) continue;
+ if(((i+n) > 0) && ((j+m) > 0) && ((i+n) < rows) && ((j+m) < cols)){
+ ptr = (i+n) * cols;
+ if(magImage[ptr+j+m] > cannyLow){
+ //
+ // this point is above low threshold
+ //
+ if(traceEdge(i+n, j+m, rows, cols, cannyLow, magImage, HYSImage)){
+ flag = 1;
+ break;
+ }
+ }
+ }
+ }
+ if(flag) break;
+ }
+ return(1);
+ }
+
+ return(0);
+
+}
+
+
+void edgeThreshold(int rows, int cols, double cannyLow, float *magImage, float *HYSImage){
+
+ int i, j;
+ int ptr;
+
+ for(i = 0; i < rows; ++i){
+ ptr = i * cols;
+ for(j = 0; j < cols; ++j){
+ if(magImage[ptr+j] > cannyLow){
+ HYSImage[ptr+j] = (float)1.0;
+ }
+ }
+ }
+
+ return;
+
+}
+
+void edgeHysteresis(int rows, int cols, double cannyLow, double cannyHigh, float *magImage, float *HYSImage){
+
+ int i, j;
+ int ptr;
+
+ for(i = 0; i < rows; ++i){
+ ptr = i * cols;
+ for(j = 0; j < cols; ++j){
+ if(magImage[ptr+j] > cannyHigh){
+ traceEdge(i, j, rows, cols, cannyLow, magImage, HYSImage);
+ }
+ }
+ }
+
+ return;
+
+}
+
+void nonMaxSupress(int rows, int cols, float aveXValue, float aveYValue, double *cannyLow, double *cannyHigh,
+ int mode, float *hDGImage, float *vDGImage, float *magImage){
+
+ int i, j;
+ int ptr, ptr_m1, ptr_p1;
+ float xSlope, ySlope, G1, G2, G3, G4, G, xC, yC;
+ float scale;
+ float maxValue = (float)0.0;
+ float minValue = (float)-1.0;
+ int histogram[256];
+ int value;
+ int mValue;
+ int mIndex;
+ int count;
+ double step;
+ double tAve;
+
+ for(i = 1; i < rows-1; ++i){
+ ptr = i * cols;
+ ptr_m1 = ptr - cols;
+ ptr_p1 = ptr + cols;
+ for(j = 1; j < cols; ++j){
+ magImage[ptr+j] = (float)0.0;
+ xC = hDGImage[ptr+j];
+ yC = vDGImage[ptr+j];
+ if((fabs(xC) < aveXValue) && (fabs(yC) < aveYValue)) continue;
+ G = magnitude(xC, yC);
+ if(fabs(yC) > fabs(xC)){
+ // vertical gradient
+ xSlope = (float)(fabs(xC) / fabs(yC));
+ ySlope = (float)1.0;
+ G2 = magnitude(hDGImage[ptr_m1+j], vDGImage[ptr_m1+j]);
+ G4 = magnitude(hDGImage[ptr_p1+j], vDGImage[ptr_p1+j]);
+ if((xC*yC) > (float)0.0){
+ G1 = magnitude(hDGImage[ptr_m1+j-1], vDGImage[ptr_m1+j-1]);
+ G3 = magnitude(hDGImage[ptr_p1+j+1], vDGImage[ptr_p1+j+1]);
+ }
+ else{
+ G1 = magnitude(hDGImage[ptr_m1+j+1], vDGImage[ptr_m1+j+1]);
+ G3 = magnitude(hDGImage[ptr_p1+j-1], vDGImage[ptr_p1+j-1]);
+ }
+ }
+ else{
+ // horizontal gradient
+ xSlope = (float)(fabs(yC) / fabs(xC));
+ ySlope = (float)1.0;
+ G2 = magnitude(hDGImage[ptr+j+1], vDGImage[ptr+j+1]);
+ G4 = magnitude(hDGImage[ptr+j-1], vDGImage[ptr+j-1]);
+ if((xC*yC) > (float)0.0){
+ G1 = magnitude(hDGImage[ptr_p1+j+1], vDGImage[ptr_p1+j+1]);
+ G3 = magnitude(hDGImage[ptr_m1+j-1], vDGImage[ptr_m1+j-1]);
+ }
+ else{
+ G1 = magnitude(hDGImage[ptr_m1+j+1], vDGImage[ptr_m1+j+1]);
+ G3 = magnitude(hDGImage[ptr_p1+j-1], vDGImage[ptr_p1+j-1]);
+ }
+ }
+ if( (G > (xSlope*G1 + (ySlope-xSlope)*G2)) && (G > (xSlope*G3 + (ySlope-xSlope)*G4)) ){
+ magImage[ptr+j] = G;
+ }
+ if(magImage[ptr+j] > maxValue) maxValue = magImage[ptr+j];
+ if(magImage[ptr+j] < minValue) minValue = magImage[ptr+j];
+ }
+ }
+
+ scale = (float)1.0 / (maxValue-minValue);
+ ptr = 0;
+ count = 0;
+ tAve = 0.0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ magImage[ptr] = scale * (magImage[ptr]-minValue);
+ if(magImage[ptr] > 0.0){
+ tAve += magImage[ptr];
+ ++count;
+ }
+ ++ptr;
+ }
+ }
+ tAve /= (float)count;
+
+ step = 255.0;
+ for(i = 0; i < 256; ++i){
+ histogram[i] = 0;
+ }
+ ptr = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = (int)(step*(magImage[ptr]));
+ ++histogram[value];
+ ++ptr;
+ }
+ }
+ //
+ // now get the max after skipping the low values
+ //
+ mValue = -1;
+ mIndex = 0;
+ for(i = 10; i < 256; ++i){
+ if(histogram[i] > mValue){
+ mValue = histogram[i];
+ mIndex = i;
+ }
+ }
+
+ if(mode == 1){
+ // based on the mean value of edge energy
+ *cannyLow = ((*cannyLow) * tAve);
+ *cannyHigh = ((*cannyHigh) * tAve);
+ }
+ else{
+ // based on the mode value of edge energy
+ *cannyLow = ((*cannyLow) * ((float)mIndex/step));
+ *cannyHigh = ((*cannyHigh) * ((float)mIndex/step));
+ }
+
+ return;
+
+}
+
+void DGFilters(int samples, int rows, int cols, double cannySigma, int gWidth,
+ float *aveXValue, float *aveYValue, double *rawImage,
+ double *dgKernel, float *hDGImage, float *vDGImage){
+
+ //
+ // implements the derivative of Gaussian filter. kernel set by CannyEdges
+ //
+ int i, j, k;
+ int ptr;
+ int mLength;
+ int count;
+ float *tBuffer = NULL;
+ double sum;
+
+ *aveXValue = (float)0.0;
+ *aveYValue = (float)0.0;
+
+ mLength = MAX(rows, cols) + 64;
+ tBuffer = calloc(mLength, sizeof(float));
+
+ //
+ // filter X
+ //
+ count = 0;
+ for(i = 0; i < rows; ++i){
+ ptr = i * cols;
+ for(j = gWidth; j < cols-gWidth; ++j){
+ sum = dgKernel[0] * rawImage[ptr+j];
+ for(k = 1; k < gWidth; ++k){
+ sum += dgKernel[k] * (-rawImage[ptr+j+k] + rawImage[ptr+j-k]);
+ }
+ hDGImage[ptr+j] = (float)sum;
+ if(sum != (float)0.0){
+ ++count;
+ *aveXValue += (float)fabs(sum);
+ }
+ }
+ }
+ if(count){
+ *aveXValue /= (float)count;
+ *aveXValue = (float)0.5 * (*aveXValue);
+ // this is 50% of the max, hardwirred for now, and is part of the threshold
+ }
+ //
+ // filter Y
+ //
+ count = 0;
+ for(i = 0; i < cols; ++i){
+ for(j = 0; j < rows; ++j){
+ ptr = j * cols;
+ tBuffer[j] = rawImage[ptr+i];
+ }
+ for(j = gWidth; j < rows-gWidth; ++j){
+ ptr = j * cols;
+ sum = dgKernel[0] * tBuffer[j];
+ for(k = 1; k < gWidth; ++k){
+ sum += dgKernel[k] * (-tBuffer[j+k] + tBuffer[j-k]);
+ }
+ vDGImage[ptr+i] = sum;
+ if(sum != (float)0.0){
+ ++count;
+ *aveYValue += (float)fabs(sum);
+ }
+ }
+ }
+ if(count){
+ *aveYValue /= (float)count;
+ *aveYValue = (float)0.5 * (*aveYValue);
+ // this is 50% of the max, hardwirred for now, and is part of the threshold
+ }
+
+ free(tBuffer);
+
+ return;
+
+}
+
+
+int NI_CannyEdges(int samples, int rows, int cols, double cannySigma, double cannyLow, double cannyHigh, int mode,
+ int lowThreshold, int highThreshold, double BPHigh, int apearture, double *rawImage,
+ unsigned short *edgeImage, int *groups){
+
+ int i, j;
+ int offset;
+ int doHysteresis = 0;
+ int gWidth;
+ int mLength;
+ int status;
+ float aveXValue;
+ float aveYValue;
+ double t;
+ double dgKernel[20];
+ float *HYSImage = NULL;
+ float *hDGImage = NULL;
+ float *vDGImage = NULL;
+ float *magImage = NULL;
+ float *tBuffer = NULL;
+
+ // filter
+ printf("do preProcess\n");
+ doPreProcess(samples, rows, cols, rawImage, BPHigh, apearture, lowThreshold, highThreshold);
+ printf("do Canny\n");
+
+ //
+ // memory for magnitude, horizontal and vertical derivative of Gaussian filter
+ //
+ mLength = MAX(rows, cols) + 64;
+ HYSImage = calloc(samples, sizeof(float));
+ hDGImage = calloc(samples, sizeof(float));
+ vDGImage = calloc(samples, sizeof(float));
+ magImage = calloc(samples, sizeof(float));
+ tBuffer = calloc(mLength, sizeof(float));
+
+ //
+ // build derivative of Gaussian filter kernel
+ // kernel is anti-symmetric so convolution is k[j]*(v[i+j] - v[i-j])
+ //
+ gWidth = 20;
+ for(i = 0; i < gWidth; ++i){
+ t = (float)i;
+ dgKernel[i] = (float)exp((double)((-t*t)/((float)2.0 * cannySigma * cannySigma)));
+ dgKernel[i] *= -(t / (cannySigma * cannySigma));
+ }
+ for(i = 0; i < samples; ++i){
+ HYSImage[i] = (float)0.0;
+ }
+
+ DGFilters(samples, rows, cols, cannySigma, gWidth, &aveXValue, &aveYValue, rawImage, dgKernel, hDGImage, vDGImage);
+ nonMaxSupress(rows, cols, aveXValue, aveYValue, &cannyLow, &cannyHigh, mode, hDGImage, vDGImage, magImage);
+ if(doHysteresis){
+ edgeHysteresis(rows, cols, cannyLow, cannyHigh, magImage, HYSImage);
+ }
+ else{
+ edgeThreshold(rows, cols, cannyLow, magImage, HYSImage);
+ }
+
+ //
+ // edge image
+ //
+ for(i = 0; i < samples; ++i){
+ edgeImage[i] = (unsigned short)HYSImage[i];
+ }
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+
+ free(tBuffer);
+ free(hDGImage);
+ free(vDGImage);
+ free(magImage);
+ free(HYSImage);
+
+ status = *groups;
+ return status;
+
+}
+
+void doSobel(int samples, int rows, int cols, double sobelLow, int mode, double *rawImage, unsigned short *edgeImage){
+
+ int i, j;
+ int p, m, n;
+ int offset;
+ int offsetM1;
+ int offsetP1;
+ int minValue, maxValue;
+ int pAve = 0;
+ int count = 0;
+ int histogram[256];
+ int value;
+ int maxIndex;
+ float pThreshold;
+ double scale;
+ double step;
+ float *filteredImage = NULL;
+
+ filteredImage = calloc(samples, sizeof(float));
+
+ minValue = 10000;
+ maxValue = -10000;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ filteredImage[offset+j] = 0;
+ edgeImage[offset+j] = 0;
+ }
+ offset += cols;
+ }
+
+ //
+ // Sobel
+ //
+ offset = cols;
+ for(i = 1; i < rows-1; ++i){
+ offsetM1 = offset - cols;
+ offsetP1 = offset + cols;
+ for(j = 1; j < cols-1; ++j){
+ n = 2*rawImage[offsetM1+j] + rawImage[offsetM1+j-1] + rawImage[offsetM1+j+1] -
+ 2*rawImage[offsetP1+j] - rawImage[offsetP1+j-1] - rawImage[offsetP1+j+1];
+ m = 2*rawImage[offset+j-1] + rawImage[offsetM1+j-1] + rawImage[offsetP1+j-1] -
+ 2*rawImage[offset+j+1] - rawImage[offsetM1+j+1] - rawImage[offsetP1+j+1];
+ p = (int)sqrt((float)(m*m) + (float)(n*n));
+ if(p > 0){
+ pAve += p;
+ ++count;
+ if(p > maxValue) maxValue = p;
+ if(p < minValue) minValue = p;
+ }
+ filteredImage[offset+j] = p;
+ }
+ offset += cols;
+ }
+
+ // threshold based on ave
+ pAve /= count;
+ scale = 1.0 / maxValue;
+
+ step = 255.0/(maxValue-minValue);
+ for(i = 0; i < 256; ++i){
+ histogram[i] = 0;
+ }
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = (int)(step*(filteredImage[offset+j]-minValue));
+ ++histogram[value];
+ }
+ offset += cols;
+ }
+ //
+ // now get the max after skipping the low values
+ //
+ maxValue = -1;
+ maxIndex = 0;
+ for(i = 10; i < 256; ++i){
+ if(histogram[i] > maxValue){
+ maxValue = histogram[i];
+ maxIndex = i;
+ }
+ }
+
+ if(mode == 1){
+ // based on the mean value of edge energy
+ pThreshold = (int)(sobelLow * (float)pAve);
+ }
+ else{
+ // based on the mode value of edge energy
+ pThreshold = (sobelLow * (minValue + ((float)maxIndex/step)));
+ }
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(filteredImage[offset+j] > pThreshold){
+ edgeImage[offset+j] = 1;
+ }
+ else{
+ edgeImage[offset+j] = 0;
+ }
+ filteredImage[offset+j] *= scale;
+ }
+ offset += cols;
+ }
+
+ free(filteredImage);
+
+ return;
+
+
+}
+
+void estimateThreshold(float *lowThreshold, float *highThreshold, float ShenCastanLow, int rows, int cols, float *SourceImage){
+
+ int i, j;
+ int offset;
+ int value;
+ int mIndex;
+ int histogram[256];
+ float low, high;
+ float scale;
+
+ low = (float)1000.0;
+ high = (float)-1000.0;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(fabs(SourceImage[offset+j]) > high) high = fabs(SourceImage[offset+j]);
+ if(fabs(SourceImage[offset+j]) < low) low = fabs(SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ scale = (float)255.0 / (high-low);
+ for(i = 0; i < 256; ++i){
+ histogram[i] = 0;
+ }
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = (int)(scale*(fabs(SourceImage[offset+j]) - low));
+ ++histogram[value];
+ }
+ offset += cols;
+ }
+
+ //
+ // now get the edge energy mode
+ //
+ value = 0;
+ mIndex = 10;
+ for(i = 10; i < 256; ++i){
+ if(histogram[i] > value){
+ value = histogram[i];
+ mIndex = i;
+ }
+ }
+
+ *highThreshold = ((float)mIndex / scale) + low;
+ *lowThreshold = ((float)mIndex / scale) + low;
+
+ *highThreshold *= ShenCastanLow;
+ *lowThreshold *= ShenCastanLow;
+
+ return;
+
+}
+
+void thresholdEdges(float *SourceImage, unsigned short *EdgeImage, double ShenCastanLow, int rows, int cols){
+
+ int i, j;
+ int offset;
+ float tLow, tHigh;
+
+ //
+ // SourceImage contains the adaptive gradient
+ // get threshold from the mode of the edge energy
+ //
+ estimateThreshold(&tLow, &tHigh, ShenCastanLow, rows, cols, SourceImage);
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(SourceImage[offset+j] > tLow){
+ EdgeImage[offset+j] = 1;
+ }
+ else{
+ EdgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+float adaptiveGradient(float *BLImage, float *FilterImage, int nrow, int ncol, int cols, int window){
+
+ int i, j;
+ int offset;
+ int numOn, numOff;
+ int hWindow = window/2;
+ float sumOn, sumOff;
+ float aveOn, aveOff;
+
+ numOn = 0;
+ numOff = 0;
+
+ sumOn = (float)0.0;
+ sumOff = (float)0.0;
+
+ aveOn = (float)0.0;
+ aveOff = (float)0.0;
+
+ offset = nrow * cols;
+ for(i = -hWindow; i < hWindow; ++i){
+ for(j = -hWindow; j < hWindow; ++j){
+ if(BLImage[offset+(i*cols)+(j+ncol)] == 1){
+ sumOn += FilterImage[offset+(i*cols)+(j+ncol)];
+ ++numOn;
+ }
+ else{
+ sumOff += FilterImage[offset+(i*cols)+(j+ncol)];
+ ++numOff;
+ }
+ }
+ }
+
+ if(numOn){
+ aveOn = sumOn / numOn;
+ }
+
+ if(numOff){
+ aveOff = sumOff / numOff;
+ }
+
+ return (aveOff-aveOn);
+
+}
+
+void getZeroCrossings(float *SourceImage, float *FilterImage, float *BLImage, int rows, int cols, int window){
+
+ int i, j;
+ int offset;
+ bool validEdge;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ SourceImage[offset+j] = 0.0;
+ }
+ offset += cols;
+ }
+
+ offset = window*cols;
+ for(i = window; i < rows-window; ++i){
+ for(j = window; j < cols-window; ++j){
+ validEdge = FALSE;
+ if((BLImage[offset+j] == 1) && (BLImage[offset+cols+j] == 0)){
+ if((FilterImage[offset+cols+j] - FilterImage[offset-cols+j]) > 0.0){
+ validEdge = TRUE;
+ }
+ }
+ else if((BLImage[offset+j] == 1) && (BLImage[offset+j+1] == 0)){
+ if((FilterImage[offset+j+1] - FilterImage[offset+j-1]) > 0.0){
+ validEdge = TRUE;
+ }
+ }
+ else if((BLImage[offset+j] == 1) && (BLImage[offset-cols+j] == 0)){
+ if((FilterImage[offset+cols+j] - FilterImage[offset-cols+j]) < 0.0){
+ validEdge = TRUE;
+ }
+ }
+ else if((BLImage[offset+j] == 1) && (BLImage[offset+j-1] == 0)){
+ if((FilterImage[offset+j+1] - FilterImage[offset+j-1]) < 0.0){
+ validEdge = TRUE;
+ }
+ }
+ if(validEdge){
+ // adaptive gradeint is signed
+ SourceImage[offset+j] = (float)fabs(adaptiveGradient(BLImage, FilterImage, i, j, cols, window));
+ }
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+
+void computeBandedLaplacian(float *image1, float *image2, float *BLImage, int rows, int cols){
+
+ int i, j;
+ int offset;
+ float t;
+
+ //
+ // like an unsharp mask
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ t = image1[offset+j] - image2[offset+j];
+ if(t < (float)0.0){
+ t = (float)0.0;
+ }
+ else{
+ t = (float)1.0;
+ }
+ BLImage[offset+j] = t;
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+void thresholdImage(float *Raw, float *Filtered, int rows, int cols, int tLow, int tHigh){
+
+ int i, j;
+ int ptr;
+
+ ptr = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(Raw[ptr] > tHigh){
+ Raw[ptr] = 0.0;
+ Filtered[ptr] = 0.0;
+ }
+ if(Raw[ptr] < tLow){
+ Raw[ptr] = 0.0;
+ Filtered[ptr] = 0.0;
+ }
+ ++ptr;
+ }
+ }
+
+ return;
+
+}
+
+void ISEF_Vertical(float *SourceImage, float *FilterImage, float *A, float *B, int rows, int cols, double b){
+
+
+ int i, j;
+ int offset;
+ float b1, b2;
+
+ b1 = ((float)1.0 - b)/((float)1.0 + b);
+ b2 = b * b1;
+
+ //
+ // set the boundaries
+ //
+ offset = (rows-1)*cols;
+ for(i = 0; i < cols; ++i){
+ // process row 0
+ A[i] = b1 * SourceImage[i];
+ // process row N-1
+ B[offset+i] = b2 * SourceImage[offset+i];
+ }
+
+ //
+ // causal component of IIR filter
+ //
+ offset = cols;
+ for(i = 1; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ //
+ // IIR ISEF filter applied across rows
+ //
+ A[offset+j] = (b * A[offset-cols+j]) + (b1 * SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ //
+ // anti-causal component of IIR filter
+ //
+ offset = (rows-2)*cols;
+ for(i = rows-2; i >= 0; --i){
+ for(j = 0; j < cols; ++j){
+ //
+ // IIR ISEF filter applied across rows
+ //
+ B[offset+j] = (b * B[offset+cols+j]) + (b2 * SourceImage[offset+j]);
+ }
+ offset -= cols;
+ }
+
+ offset = (rows-1)*cols;
+ for(j = 0; j < cols-1; ++j){
+ FilterImage[offset+j] = A[offset+j];
+ }
+
+ //
+ // add causal and anti-causal IIR parts
+ //
+ offset = 0;
+ for(i = 1; i < rows-2; ++i){
+ for(j = 0; j < cols-1; ++j){
+ FilterImage[offset+j] = A[offset+j] + B[offset+cols+j];
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+void ISEF_Horizontal(float *SourceImage, float *FilterImage, float *A, float *B, int rows, int cols, double b){
+
+
+ //
+ // source and smooth are the same in this pass of the 2D IIR
+ //
+
+ int i, j;
+ int offset;
+ float b1, b2;
+
+ b1 = ((float)1.0 - b)/((float)1.0 + b);
+ b2 = b * b1;
+
+ //
+ // columns boundaries
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ // col 0
+ A[offset] = b1 * SourceImage[offset];
+ // col N-1
+ B[offset+cols-1] = b2 * SourceImage[offset+cols-1];
+ }
+
+ //
+ // causal IIR part
+ //
+ offset = 0;
+ for(j = 1; j < cols; ++j){
+ for(i = 0; i < rows; ++i){
+ A[offset+j] = (b * A[offset+j-1]) + (b1 * SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ //
+ // anti-causal IIR part
+ //
+ offset = 0;
+ for(j = cols-2; j > 0; --j){
+ for(i = 0; i < rows; ++i){
+ B[offset+j] = (b * B[offset+j+1]) + (b2 * SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ //
+ // filtered output. this is 2-pass IIR and pass 1 is vertical
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ FilterImage[offset+cols-1] = A[offset+cols-1];
+ }
+
+ //
+ // add causal and anti-causal IIR parts
+ //
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols-1; ++j){
+ FilterImage[offset+j] = A[offset+j] + B[offset+j+1];
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+
+void computeISEF(float *SourceImage, float *FilterImage, int rows, int cols, double b){
+
+ int imageSize = rows*cols;
+ float *A;
+ float *B;
+
+ A = calloc(imageSize, sizeof(float));
+ B = calloc(imageSize, sizeof(float));
+
+ ISEF_Vertical(SourceImage, FilterImage, A, B, rows, cols, b);
+ ISEF_Horizontal(FilterImage, FilterImage, A, B, rows, cols, b);
+
+ free(A);
+ free(B);
+
+ return;
+
+}
+
+void Shen_Castan(double b, double ShenCastanLow, int rows, int cols, int window, int lowThreshold, int highThreshold,
+ double *RawImage, unsigned short *EdgeImage){
+
+ int i;
+ int imageSize = rows*cols;
+ float *FilterImage;
+ float *BinaryLaplacianImage;
+ float *SourceImage;
+
+ FilterImage = calloc(imageSize, sizeof(float));
+ BinaryLaplacianImage = calloc(imageSize, sizeof(float));
+ SourceImage = calloc(imageSize, sizeof(float));
+
+ for(i = 0; i < imageSize; ++i){
+ SourceImage[i] = RawImage[i];
+ }
+ computeISEF(SourceImage, FilterImage, rows, cols, b);
+ // optional thresholding based on low, high
+ thresholdImage(SourceImage, FilterImage, rows, cols, lowThreshold, highThreshold);
+ computeBandedLaplacian(FilterImage, SourceImage, BinaryLaplacianImage, rows, cols);
+ // the new source image is now the adaptive gradient
+ getZeroCrossings(SourceImage, FilterImage, BinaryLaplacianImage, rows, cols, window);
+ thresholdEdges(SourceImage, EdgeImage, ShenCastanLow, rows, cols);
+
+ free(FilterImage);
+ free(BinaryLaplacianImage);
+ free(SourceImage);
+
+ return;
+
+}
+
+int NI_ShenCastanEdges(int samples, int rows, int cols, double b, double ShenCastanLow, int window, int lowThreshold,
+ int highThreshold, double *rawImage, unsigned short *edgeImage, int *groups){
+
+
+ int i, j;
+ int offset;
+ int status = 0;
+
+ Shen_Castan(b, ShenCastanLow, rows, cols, window, lowThreshold, highThreshold, rawImage, edgeImage);
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ status = *groups;
+
+ return status;
+
+}
+
+void buildBinaryImage(int rows, int cols, double *rawImage, unsigned short *edgeImage, int lowThreshold, int highThreshold){
+
+ int i, j;
+ int offset;
+ double value;
+ int maskValue;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = rawImage[offset+j];
+ maskValue = 1;
+ if(value < (double)lowThreshold) maskValue = 0;
+ if(value > (double)highThreshold) maskValue = 0;
+ edgeImage[offset+j] = maskValue;
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+
+
+void morphoFilterBinaryImage(int rows, int cols, unsigned short *edgeImage, int CloseSize, int OpenSize){
+
+
+ int i, j;
+ int offset, offset2;
+ unsigned short cmask[11][11];
+ unsigned short omask[11][11];
+ int olapValuesC[4];
+ int olapValuesO[4];
+ int CloseMaskSize;
+ int OpenMaskSize;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+ int spadSize;
+ unsigned char *ImageE;
+ unsigned char *ImageC;
+
+ spadSize = MAX(rows, cols);
+
+ ImageE = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageC = calloc(spadSize*spadSize, sizeof(unsigned char));
+
+ //
+ // Close filter
+ //
+ if(CloseSize){
+ CloseMaskSize = (CloseSize-1)/2;
+ for(i = 0; i < 2*CloseMaskSize+1; ++i){
+ for(j = 0; j < 2*CloseMaskSize+1; ++j){
+ cmask[i][j] = 1;
+ }
+ }
+ LowValue1 = 0;
+ HighValue1 = 1;
+ LowValue2 = 1;
+ HighValue2 = 0;
+ olapValuesC[0] = LowValue1;
+ olapValuesC[1] = HighValue1;
+ olapValuesC[2] = LowValue2;
+ olapValuesC[3] = HighValue2;
+ }
+
+ //
+ // Open filter
+ //
+ if(OpenSize){
+ OpenMaskSize = (OpenSize-1)/2;
+ for(i = 0; i < 2*OpenMaskSize+1; ++i){
+ for(j = 0; j < 2*OpenMaskSize+1; ++j){
+ omask[i][j] = 1;
+ }
+ }
+ LowValue1 = 1;
+ HighValue1 = 0;
+ LowValue2 = 0;
+ HighValue2 = 1;
+ olapValuesO[0] = LowValue1;
+ olapValuesO[1] = HighValue1;
+ olapValuesO[2] = LowValue2;
+ olapValuesO[3] = HighValue2;
+ }
+
+ offset = 0;
+ offset2 = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ ImageE[offset2+j] = (unsigned char)edgeImage[offset+j];
+ }
+ offset2 += spadSize;
+ offset += cols;
+ }
+
+ if(OpenSize){
+ OpenCloseFilter(olapValuesO, OpenMaskSize, rows, cols, spadSize, ImageE, ImageC, omask);
+ }
+
+ if(CloseSize){
+ OpenCloseFilter(olapValuesC, CloseMaskSize, rows, cols, spadSize, ImageE, ImageC, cmask);
+ }
+
+ offset = 0;
+ offset2 = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(ImageE[offset2+j] == 1){
+ // this will activate some original off-pixels
+ edgeImage[offset+j] = 1;
+ }
+ else{
+ // this will zero some original on-pixels
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset2 += spadSize;
+ offset += cols;
+ }
+
+ free(ImageE);
+ free(ImageC);
+
+ return;
+
+}
+
+void doRegionGrow(int samples, int rows, int cols, double *rawImage, unsigned short *edgeImage, int lowThreshold,
+ int highThreshold, int closeWindow, int openWindow){
+
+ buildBinaryImage(rows, cols, rawImage, edgeImage, lowThreshold, highThreshold);
+ morphoFilterBinaryImage(rows, cols, edgeImage, closeWindow, openWindow);
+
+ return;
+
+}
+
+int NI_RegionGrow(int samples, int rows, int cols, int lowThreshold, int highThreshold, int closeWindow,
+ int openWindow, double *rawImage, unsigned short *edgeImage, int *groups){
+
+ int i, j;
+ int offset;
+ int status;
+
+ doRegionGrow(samples, rows, cols, rawImage, edgeImage, lowThreshold, highThreshold, closeWindow, openWindow);
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ status = *groups;
+ return status;
+
+}
+
+int NI_SobelEdges(int samples, int rows, int cols, double sobelLow, int mode, int lowThreshold, int highThreshold, double BPHigh,
+ int apearture, double *rawImage, unsigned short *edgeImage, int *groups){
+
+
+ int i, j;
+ int offset;
+ int status;
+
+ doPreProcess(samples, rows, cols, rawImage, BPHigh, apearture, lowThreshold, highThreshold);
+ doSobel(samples, rows, cols, sobelLow, mode, rawImage, edgeImage);
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ status = *groups;
+ return status;
+
+}
+
+void initThinFilter(int J_mask[3][30], int K_mask[3][30]){
+
+ int i, j;
+ int Column;
+
+ for(i = 0; i < 3; ++i){
+ for(j = 0; j < 30; ++j){
+ J_mask[i][j] = 0;
+ K_mask[i][j] = 0;
+ }
+ }
+
+ Column = 0;
+ J_mask[0][Column+0] = 1;
+ J_mask[0][Column+1] = 1;
+ J_mask[0][Column+2] = 1;
+ J_mask[1][Column+1] = 1;
+
+ Column += 3;
+ J_mask[0][Column+1] = 1;
+ J_mask[1][Column+1] = 1;
+ J_mask[1][Column+2] = 1;
+
+ Column += 3;
+ J_mask[0][Column+0] = 1;
+ J_mask[1][Column+0] = 1;
+ J_mask[2][Column+0] = 1;
+ J_mask[1][Column+1] = 1;
+
+ Column += 3;
+ J_mask[0][Column+1] = 1;
+ J_mask[1][Column+0] = 1;
+ J_mask[1][Column+1] = 1;
+
+ Column += 3;
+ J_mask[0][Column+2] = 1;
+ J_mask[1][Column+1] = 1;
+ J_mask[1][Column+2] = 1;
+ J_mask[2][Column+2] = 1;
+
+ Column += 3;
+ J_mask[1][Column+0] = 1;
+ J_mask[1][Column+1] = 1;
+ J_mask[2][Column+1] = 1;
+
+ Column += 3;
+ J_mask[1][Column+1] = 1;
+ J_mask[2][Column+0] = 1;
+ J_mask[2][Column+1] = 1;
+ J_mask[2][Column+2] = 1;
+
+ Column += 3;
+ J_mask[1][Column+1] = 1;
+ J_mask[1][Column+2] = 1;
+ J_mask[2][Column+1] = 1;
+
+ Column = 0;
+ K_mask[2][Column+0] = 1;
+ K_mask[2][Column+1] = 1;
+ K_mask[2][Column+2] = 1;
+
+ Column += 3;
+ K_mask[1][Column+0] = 1;
+ K_mask[2][Column+0] = 1;
+ K_mask[2][Column+1] = 1;
+
+ Column += 3;
+ K_mask[0][Column+2] = 1;
+ K_mask[1][Column+2] = 1;
+ K_mask[2][Column+2] = 1;
+
+ Column += 3;
+ K_mask[1][Column+2] = 1;
+ K_mask[2][Column+1] = 1;
+ K_mask[2][Column+2] = 1;
+
+ Column += 3;
+ K_mask[0][Column+0] = 1;
+ K_mask[1][Column+0] = 1;
+ K_mask[2][Column+0] = 1;
+
+ Column += 3;
+ K_mask[0][Column+1] = 1;
+ K_mask[0][Column+2] = 1;
+ K_mask[1][Column+2] = 1;
+
+ Column += 3;
+ K_mask[0][Column+0] = 1;
+ K_mask[0][Column+1] = 1;
+ K_mask[0][Column+2] = 1;
+
+ Column += 3;
+ K_mask[0][Column+0] = 1;
+ K_mask[0][Column+1] = 1;
+ K_mask[1][Column+0] = 1;
+
+ return;
+
+}
+
+void ThinningFilter(int regRows, int regColumns, int spadSize, int J_mask[3][30], int K_mask[3][30],
+ unsigned char *Input, unsigned char *CInput, unsigned char *ErosionStage,
+ unsigned char *DialationStage, unsigned char *HMT, unsigned char *Copy){
+
+ int i, j, k, l, m, n, overlap, hit;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+ int Column, T, nloop;
+ int Offset;
+ int N, M;
+ int j_mask[3][3];
+ int k_mask[3][3];
+
+ N = regRows;
+ M = regColumns;
+
+ LowValue1 = 1;
+ HighValue1 = 0;
+
+ LowValue2 = 0;
+ HighValue2 = 1;
+
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ Copy[Offset+j] = Input[Offset+j];
+ }
+ Offset += spadSize;
+ }
+
+ nloop = 0;
+ while(1){
+ // erode
+ Column = 0;
+ for(n = 0; n < 8; ++n){
+ for(i = 0; i < 3; ++i){
+ for(j = 0; j < 3; ++j){
+ j_mask[i][j] = J_mask[i][Column+j];
+ }
+ }
+ for(i = 0; i < 3; ++i){
+ for(j = 0; j < 3; ++j){
+ k_mask[i][j] = K_mask[i][Column+j];
+ }
+ }
+ Column += 3;
+
+ Offset = spadSize;
+ for(i = 1; i < N-1; ++i){
+ for(j = 1; j < M-1; ++j){
+ hit = LowValue1;
+ for(k = -1; k < 2; ++k){
+ for(l = -1; l < 2; ++l){
+ T = j_mask[k+1][l+1];
+ if(T == 1){
+ overlap = T*Input[Offset+(k*spadSize)+j+l];
+ if(overlap == HighValue1) hit = HighValue1;
+ }
+ }
+ }
+ ErosionStage[Offset+j] = hit;
+ }
+ Offset += spadSize;
+ }
+
+ // dialate
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ CInput[Offset+j] = (~Input[Offset+j]) & 0x1;
+ }
+ Offset += spadSize;
+ }
+
+ Offset = spadSize;
+ for(i = 1; i < N-1; ++i){
+ for(j = 1; j < M-1; ++j){
+ hit = LowValue1;
+ for(k = -1; k < 2; ++k){
+ for(l = -1; l < 2; ++l){
+ T = k_mask[k+1][l+1];
+ if(T == 1){
+ overlap = T*CInput[Offset+(k*spadSize)+j+l];
+ if(overlap == HighValue1) hit = HighValue1;
+ }
+ }
+ }
+ DialationStage[Offset+j] = hit;
+ }
+ Offset += spadSize;
+ }
+
+ // form the HMT
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ m = (ErosionStage[Offset+j]*DialationStage[Offset+j]);
+ HMT[Offset+j] = m;
+ }
+ Offset += spadSize;
+ }
+
+ // Thin for stage n
+
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ HMT[Offset+j] = (~HMT[Offset+j]) & 0x1;
+ }
+ Offset += spadSize;
+ }
+
+ Offset = 0;
+ for (i = 0; i < N; ++i){
+ for (j = 0; j < M; ++j){
+ m = (Input[Offset+j]*HMT[Offset+j]);
+ Input[Offset+j] = m;
+ }
+ Offset += spadSize;
+ }
+ }
+
+ // check for the NULL set
+ hit = 0;
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ hit += abs(Copy[Offset+j]-Input[Offset+j]);
+ }
+ Offset += spadSize;
+ }
+ if(!hit) break;
+
+ hit = 0;
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ Copy[Offset+j] = Input[Offset+j];
+ if(Input[Offset+j]) ++hit;
+ }
+ Offset += spadSize;
+ }
+ ++nloop;
+ }
+
+
+ return;
+
+}
+
+
+int NI_ThinFilter(int samples, int rows, int cols, int numberObjects, unsigned short *edgeImage, objStruct objectMetrics[]){
+
+ int i, j;
+ int loop;
+ int label;
+ int left, right, top, bottom;
+ int roiRows, roiCols;
+ int srcOffset;
+ int dstOffset;
+ int status;
+ int inflate = 1;
+ int J_mask[3][30];
+ int K_mask[3][30];
+
+ unsigned char *Input;
+ unsigned char *CInput;
+ unsigned char *ErosionStage;
+ unsigned char *DialationStage;
+ unsigned char *HMT;
+ unsigned char *Copy;
+ unsigned short *thinEdgeImage;
+
+ //
+ // scratch pad (spad) memory
+ //
+ Input = calloc(samples, sizeof(unsigned char));
+ CInput = calloc(samples, sizeof(unsigned char));
+ ErosionStage = calloc(samples, sizeof(unsigned char));
+ DialationStage = calloc(samples, sizeof(unsigned char));
+ HMT = calloc(samples, sizeof(unsigned char));
+ Copy = calloc(samples, sizeof(unsigned char));
+ thinEdgeImage = calloc(samples, sizeof(unsigned short));
+
+ initThinFilter(J_mask, K_mask);
+ for(loop = 0; loop < numberObjects; ++loop){
+ label = objectMetrics[loop].Label;
+ left = objectMetrics[loop].L;
+ right = objectMetrics[loop].R;
+ top = objectMetrics[loop].T;
+ bottom = objectMetrics[loop].B;
+ roiRows = top-bottom+2*inflate;
+ roiCols = right-left+2*inflate;
+
+ //
+ // clear the scratch pad
+ //
+ srcOffset = 0;
+ for(i = 0; i < roiRows; ++i){
+ for(j = 0; j < roiCols; ++j){
+ Input[srcOffset+j] = 0;
+ }
+ srcOffset += cols;
+ }
+
+ //
+ // copy the ROI for MAT (medial axis transformation) filter
+ //
+ dstOffset = inflate*rows;
+ for(i = bottom; i < top; ++i){
+ srcOffset = i*cols;
+ for(j = left; j < right; ++j){
+ if(edgeImage[srcOffset+j] == label){
+ Input[dstOffset+j-left+inflate] = 1;
+ }
+ }
+ dstOffset += cols;
+ }
+ ThinningFilter(roiRows, roiCols, cols, J_mask, K_mask, Input, CInput, ErosionStage, DialationStage, HMT, Copy);
+
+ //
+ // copy the MAT roi to the new edgeImage (clip the inflate border)
+ //
+ dstOffset = inflate*rows;
+ for(i = bottom; i < top; ++i){
+ srcOffset = i*cols;
+ for(j = left; j < right; ++j){
+ if(Input[dstOffset+j-left+inflate]){
+ thinEdgeImage[srcOffset+j] = label;
+ }
+ }
+ dstOffset += cols;
+ }
+ }
+
+ //
+ // copy the MAT edges and return the thinned edges
+ // this will prune the isolated edge points from the edgeImage source
+ //
+ for(i = 0; i < rows*cols; ++i){
+ edgeImage[i] = thinEdgeImage[i];
+ }
+
+ free(Input);
+ free(CInput);
+ free(ErosionStage);
+ free(DialationStage);
+ free(HMT);
+ free(Copy);
+ free(thinEdgeImage);
+
+ status = 1;
+
+ return status;
+
+}
+
+
+void generateMask(unsigned char *ImageH, bPOINT *boundary, int newSamples, int label, int cols){
+
+ //
+ // get the boundary point pairs (left, right) for each line
+ // if there is no pair, then the boundary is open
+ // then fill the image in with the current label
+ //
+
+ int i, j, k, m;
+ int list[2048];
+ int distance;
+ int neighbor = 4;
+ int index;
+ int offset;
+ int maxDistance = 1024;
+ int x, y;
+ int low, high;
+
+ for(i = 0; i < newSamples; ++i){
+ boundary[i].haveLink = FALSE;
+ boundary[i].linkIndex = -1;
+ }
+
+ for(i = 0; i < newSamples; ++i){
+ if(!boundary[i].haveLink){
+ boundary[i].haveLink = TRUE;
+ x = boundary[i].x;
+ y = boundary[i].y;
+ for(k = 0, j = 0; j < newSamples; ++j){
+ if((j != i)){
+ if(boundary[j].y == y){
+ list[k] = j;
+ ++k;
+ }
+ }
+ }
+ // now get the closest boundary
+ if(k){
+ distance = maxDistance;
+ index = -1;
+ for(j = 0; j < k; ++j){
+ m = abs(x - boundary[list[j]].x);
+ if((m < distance) && (m > neighbor)){
+ distance = m;
+ index = list[j];
+ }
+ else if(m <= neighbor){
+ boundary[list[j]].haveLink = TRUE;
+ }
+ }
+ if(index != -1){
+ boundary[i].linkIndex = index;
+ boundary[index].linkIndex = i;
+ boundary[index].haveLink = TRUE;
+ if(boundary[i].x < boundary[index].x){
+ low = boundary[i].x;
+ high = boundary[index].x;
+ }
+ else{
+ low = boundary[index].x;
+ high = boundary[i].x;
+ }
+ //
+ // do the fill
+ //
+ offset = y * cols;
+ for(j = low; j <= high; ++j){
+ ImageH[offset+j] = label;
+ }
+ }
+ }
+ else{
+ // boundary point is isolated
+ boundary[i].linkIndex = i;
+ }
+ }
+ }
+
+ return;
+
+}
+
+void getBoundaryMetrics(bPOINT *boundary, float *length, float *minRadius, float *maxRadius, float *aveRadius,
+ float Xcenter, float Ycenter, int newSamples){
+
+ int j;
+ float dX, dY;
+ float distance;
+
+ if(newSamples < 2){
+ *length = (float)0.0;
+ *minRadius = (float)0.0;
+ *maxRadius = (float)0.0;
+ *aveRadius = (float)0.0;
+ return;
+ }
+
+ *length = (float)0.0;
+ for(j = 1; j < newSamples; ++j){
+ dX = (float)(boundary[j].x - boundary[j-1].x);
+ dY = (float)(boundary[j].y - boundary[j-1].y);
+ distance = (float)sqrt(dX*dX + dY*dY);
+ *length += distance;
+ }
+
+ *minRadius = (float)10000.0;
+ *maxRadius = (float)-10000.0;
+ *aveRadius = (float)0.0;
+ for(j = 0; j < newSamples; ++j){
+ dX = (float)(boundary[j].x - Xcenter);
+ dY = (float)(boundary[j].y - Ycenter);
+ distance = (float)sqrt(dX*dX + dY*dY);
+ *aveRadius += distance;
+ if(distance < *minRadius) *minRadius = distance;
+ if(distance > *maxRadius) *maxRadius = distance;
+ }
+
+ if(newSamples){
+ *aveRadius /= (float)newSamples;
+ }
+
+ return;
+
+}
+
+void trackBoundary(unsigned char *Input, blobBoundary lBoundary[], int mcount, int spadSize,
+ blobBoundary seedValue, int searchWindow){
+
+
+ int i, j, k, m, p;
+ int offset;
+ int CurI;
+ int CurJ;
+ int StrI;
+ int StrJ;
+ int NewI;
+ int NewJ;
+ int MinD;
+ int inflate = searchWindow;
+
+ CurI = seedValue.xy.x;
+ CurJ = seedValue.xy.y;
+ StrI = CurI;
+ StrJ = CurJ;
+
+ p = 0;
+ lBoundary[p].xy.x = StrI;
+ lBoundary[p].xy.y = StrJ;
+ offset = StrI * spadSize;
+
+ p = 1;
+ while(p < mcount){
+ offset = (CurI-inflate)*spadSize;
+ if(offset < 0){
+ printf("offset < 0 ");
+ printf("CurI [%d]. p [%d]. mcount [%d]\n", CurI, p, mcount);
+ getchar();
+ }
+ MinD = 1024;
+ NewI = -1;
+ NewJ = -1;
+ for(i = CurI-inflate; i < CurI+inflate; ++i){
+ for(j = CurJ-inflate; j < CurJ+inflate; ++j){
+ m = Input[offset+j];
+ if(m == 1){
+ // city block distance
+ k = abs(i-CurI) + abs(j-CurJ);
+ if(k < MinD){
+ MinD = k;
+ NewI = i;
+ NewJ = j;
+ }
+ }
+ }
+ offset += spadSize;
+ }
+ if(NewI != -1) CurI = NewI;
+ if(NewJ != -1) CurJ = NewJ;
+ offset = CurI * spadSize;
+ Input[offset+CurJ] = 0;
+ lBoundary[p].xy.x = CurJ;
+ lBoundary[p].xy.y = CurI;
+ ++p;
+ }
+
+ return;
+
+}
+
+
+void OpenCloseFilter(int olapValues[], int maskSize, int rows, int columns, int spadSize,
+ unsigned char *input, unsigned char *output, unsigned short mask[11][11]){
+
+
+ //
+ // do morphological open/close image filtering. the olapValues array determines
+ // if the filter is Open or Close.
+ //
+ int i, j, k, l, m, overlap, hit;
+ int offset;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+
+ LowValue1 = olapValues[0];
+ HighValue1 = olapValues[1];
+ LowValue2 = olapValues[2];
+ HighValue2 = olapValues[3];
+
+ // close - step 1 is dialate
+ // open - step 1 is erode
+ offset = maskSize*spadSize;
+ for(i = maskSize; i < rows-maskSize; ++i){
+ for(j = maskSize; j < columns-maskSize; ++j){
+ hit = LowValue1;
+ for(k = -maskSize; k < maskSize; ++k){
+ m = k*spadSize;
+ for(l = -maskSize; l < maskSize; ++l){
+ overlap = mask[k+maskSize][l+maskSize]*input[offset+m+j+l];
+ if(overlap == HighValue1){
+ hit = HighValue1;
+ }
+ }
+ }
+ output[offset+j] = hit;
+ }
+ offset += spadSize;
+ }
+
+ // close - step 2 is erode
+ // open - step 2 is dialate
+ offset = maskSize*spadSize;
+ for(i = maskSize; i < rows-maskSize; ++i){
+ for(j = maskSize; j < columns-maskSize; ++j){
+ hit = LowValue2;
+ for(k = -maskSize; k < maskSize; ++k){
+ m = k*spadSize;
+ for(l = -maskSize; l < maskSize; ++l){
+ overlap = mask[k+maskSize][l+maskSize]*output[offset+m+j+l];
+ if(overlap == HighValue2){
+ hit = HighValue2;
+ }
+ }
+ }
+ input[offset+j] = hit;
+ }
+ offset += spadSize;
+ }
+
+ return;
+}
+
+void getCompactness(unsigned char *Input, RECT roi, int label, int spadSize, float *vCompactness, float length){
+
+ int i, j;
+ int maskOffset;
+ int area;
+ static float fpi = (float)(4.0 * 3.14159);
+
+ area = 0;
+ for(i = roi.bottom; i < roi.top; ++i){
+ maskOffset = i*spadSize;
+ for(j = roi.left; j < roi.right; ++j){
+ if(Input[maskOffset+j] == label){
+ ++area;
+ }
+ }
+ }
+ if(area && (length != (float)0.0)){
+ *vCompactness = (fpi * (float)area) / (length*length);
+ }
+ else{
+ *vCompactness = (float)0.0;
+ }
+
+ return;
+}
+
+
+void doMorphology(unsigned char *Input, unsigned char *ImageE, unsigned char *ImageC, unsigned char *ImageH,
+ int olapValuesC[],int olapValuesO[], unsigned short cmask[11][11], unsigned short omask[11][11],
+ RECT roi, int label, int CloseMaskSize, int OpenMaskSize, int spadSize){
+
+ int i, j;
+ int rows, cols;
+ int srcOffset;
+ int dstOffset;
+ int maskSize;
+
+ cols = roi.right - roi.left;
+ rows = roi.top - roi.bottom;
+
+ for(i = 0; i < spadSize*spadSize; ++i){
+ ImageE[i] = 0;
+ ImageC[i] = 0;
+ }
+
+ //
+ // put the ROI in the ImageE array centered in ULC
+ //
+ dstOffset = 0;
+ for(i = roi.bottom; i < roi.top; ++i){
+ srcOffset = i*spadSize;
+ for(j = roi.left; j < roi.right; ++j){
+ if(ImageH[srcOffset+j] == label){
+ ImageE[dstOffset+j-roi.left] = 1;
+ }
+ }
+ dstOffset += spadSize;
+ }
+
+ //
+ // open
+ //
+ maskSize = OpenMaskSize;
+ OpenCloseFilter(olapValuesO, maskSize, rows, cols, spadSize, ImageE, ImageC, omask);
+ //
+ // close
+ //
+ maskSize = CloseMaskSize;
+ OpenCloseFilter(olapValuesC, maskSize, rows, cols, spadSize, ImageE, ImageC, cmask);
+
+ //
+ // put the closed ROI (in ImageE) back in its roi space
+ //
+
+ srcOffset = 0;
+ for(i = roi.bottom; i < roi.top+2*maskSize+1; ++i){
+ dstOffset = (i-(2*maskSize+1))*spadSize;
+ for(j = roi.left-maskSize-1; j < roi.right+maskSize+1; ++j){
+ if(ImageE[srcOffset+j-roi.left] == 1){
+ Input[dstOffset+j-maskSize+1] = label;
+ }
+ }
+ srcOffset += spadSize;
+ }
+
+ return;
+
+}
+
+
+void getBoundary(unsigned short *ThinEdgeImage, unsigned char *Input, blobBoundary *pBoundary, blobBoundary *lBoundary,
+ boundaryIndex *pBoundaryIndex, RECT boundBox, int label, int bBox, int nextSlot, int memOffset,
+ int spadSize, int searchWindow){
+
+ int i, j;
+ int dstOffset;
+ int srcOffset;
+ int mcount;
+ int rows;
+ int columns;
+ bool first;
+ blobBoundary value;
+ int inflate = searchWindow+1;
+ int count;
+
+ pBoundaryIndex[bBox+1].rectangle.left = boundBox.left;
+ pBoundaryIndex[bBox+1].rectangle.right = boundBox.right;
+ pBoundaryIndex[bBox+1].rectangle.top = boundBox.top;
+ pBoundaryIndex[bBox+1].rectangle.bottom = boundBox.bottom;
+
+ for(i = 0; i < spadSize*spadSize; ++i){
+ Input[i] = 0;
+ }
+
+ //copy to spad
+
+ count = 0;
+ rows = boundBox.top-boundBox.bottom+2*inflate;
+ columns = boundBox.right-boundBox.left+2*inflate;
+ dstOffset = inflate*spadSize;
+ for(i = boundBox.bottom; i < boundBox.top; ++i){
+ srcOffset = i*spadSize;
+ for(j = boundBox.left; j < boundBox.right; ++j){
+ if(ThinEdgeImage[srcOffset+j] == label){
+ Input[dstOffset+j-boundBox.left+inflate] = 1;
+ ++count;
+ }
+ }
+ dstOffset += spadSize;
+ }
+
+ mcount = 0;
+ first = TRUE;
+ srcOffset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < columns; ++j){
+ if(Input[srcOffset+j]){
+ if(first){
+ first = FALSE;
+ // index of the seed sample
+ value.xy.x = i;
+ value.xy.y = j;
+ }
+ ++mcount;
+ }
+ }
+ srcOffset += spadSize;
+ }
+
+ trackBoundary(Input, lBoundary, mcount, spadSize, value, searchWindow);
+
+ pBoundaryIndex[nextSlot].numberPoints = mcount;
+ for(i = 0; i < mcount; ++i){
+ value.xy.x = lBoundary[i].xy.x + boundBox.left - inflate;
+ value.xy.y = lBoundary[i].xy.y + boundBox.bottom - inflate + 1;
+ //printf("[%d, %d]\n", value.xy.x, value.xy.y);
+ pBoundary[memOffset].xy.x = value.xy.x;
+ pBoundary[memOffset].xy.y = value.xy.y;
+ ++memOffset;
+ }
+ //getchar();
+
+ return;
+
+}
+
+
+void buildBoundary(objStruct objectMetrics[], int searchWindow, unsigned short *ThinEdgeImage,
+ int numberObjects, int srcRows, int srcCols){
+
+ int i, j, k;
+ int count;
+ int numBoundaries;
+ int numSamples;
+ int offset;
+ int offset2;
+ int end;
+ int label;
+ int distance;
+ // these should be setup parameters
+ int closureDistance = 12;
+ int CloseSize = 5;
+ int OpenSize = 5;
+ int threshold = 3;
+ int newSamples;
+ int spadSize;
+ POINT rectPoint[4];
+ int in[4];
+ float length;
+ float minRadius;
+ float maxRadius;
+ float aveRadius;
+ float vCompactness;
+ // for morphological close of mask. max structuring element is 11x11
+ unsigned short cmask[11][11];
+ unsigned short omask[11][11];
+ int olapValuesC[4];
+ int olapValuesO[4];
+ int CloseMaskSize;
+ int OpenMaskSize;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+ RECT bBox;
+
+ boundaryIndex *pBoundaryIndex;
+ blobBoundary *pBoundary;
+ blobBoundary *lBoundary;
+ bPOINT *boundary;
+ unsigned char *Input;
+ unsigned char *ImageE;
+ unsigned char *ImageC;
+ unsigned char *ImageH;
+
+ //
+ // Close filter
+ //
+ CloseMaskSize = (CloseSize-1)/2;
+ for(i = 0; i < 2*CloseMaskSize+1; ++i){
+ for(j = 0; j < 2*CloseMaskSize+1; ++j){
+ cmask[i][j] = 1;
+ }
+ }
+ LowValue1 = 0;
+ HighValue1 = 1;
+ LowValue2 = 1;
+ HighValue2 = 0;
+ olapValuesC[0] = LowValue1;
+ olapValuesC[1] = HighValue1;
+ olapValuesC[2] = LowValue2;
+ olapValuesC[3] = HighValue2;
+
+ //
+ // Open filter
+ //
+ OpenMaskSize = (OpenSize-1)/2;
+ for(i = 0; i < 2*OpenMaskSize+1; ++i){
+ for(j = 0; j < 2*OpenMaskSize+1; ++j){
+ omask[i][j] = 1;
+ }
+ }
+ LowValue1 = 1;
+ HighValue1 = 0;
+ LowValue2 = 0;
+ HighValue2 = 1;
+ olapValuesO[0] = LowValue1;
+ olapValuesO[1] = HighValue1;
+ olapValuesO[2] = LowValue2;
+ olapValuesO[3] = HighValue2;
+
+
+ //spadSize = MAX(srcRows, srcCols);
+ spadSize = srcCols;
+
+ pBoundaryIndex = calloc(srcRows+srcCols, sizeof(boundaryIndex));
+ Input = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageE = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageC = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageH = calloc(spadSize*spadSize, sizeof(unsigned char));
+ pBoundary = calloc(srcRows*srcCols, sizeof(blobBoundary));
+ lBoundary = calloc(32767, sizeof(blobBoundary));
+ boundary = calloc(32767, sizeof(POINT));
+
+ for(i = 0; i < (srcRows+srcCols); ++i){
+ pBoundaryIndex[i].numberPoints = 0;
+ pBoundaryIndex[i].curveClose = 0;
+ pBoundaryIndex[i].isWithin = FALSE;
+ pBoundaryIndex[i].criticalSize = FALSE;
+ pBoundaryIndex[i].closedCurve = FALSE;
+ }
+
+
+ for(i = 0; i < numberObjects; ++i){
+ //for(i = 0; i < 1; ++i){
+ ++pBoundaryIndex[0].numberPoints;
+ count = 0;
+ j = 1;
+ while(pBoundaryIndex[j].numberPoints){
+ count += pBoundaryIndex[j++].numberPoints;
+ }
+ bBox.left = objectMetrics[i].L;
+ bBox.right = objectMetrics[i].R;
+ bBox.top = objectMetrics[i].T;
+ bBox.bottom = objectMetrics[i].B;
+ label = objectMetrics[i].Label;
+ pBoundaryIndex[i+1].Label = label;
+ //printf("(%d, %d, %d, %d [%d])\n", bBox.left, bBox.right, bBox.top, bBox.bottom, label);
+ getBoundary(ThinEdgeImage, Input, pBoundary, lBoundary, pBoundaryIndex, bBox, label,
+ i, pBoundaryIndex[0].numberPoints, count, spadSize, searchWindow);
+ }
+
+ //
+ // Input will now be used in the fill. Copy the labeled edge image
+ //
+
+ //
+ // numBoundaries = numberObjects
+ //
+ offset = 0;
+ numBoundaries = pBoundaryIndex[0].numberPoints;
+ //printf("numBoundaries [%d]\n", numBoundaries);
+ for(i = 0; i < numBoundaries; ++i){
+ numSamples = pBoundaryIndex[i+1].numberPoints;
+ end = numSamples-2;
+ newSamples = numSamples-1;
+ for(j = 0; j < numSamples; ++j){
+ boundary[j].x = pBoundary[offset+j+1].xy.x;
+ boundary[j].y = pBoundary[offset+j+1].xy.y;
+ }
+
+ //
+ // clip off the ends where stray boundary pixels were left over
+ //
+ while(1){
+ distance = abs(boundary[end].x-boundary[end-1].x) + abs(boundary[end].y-boundary[end-1].y);
+ if(distance > threshold){
+ --end;
+ --newSamples;
+ }
+ else{
+ break;
+ }
+ }
+ //printf("[%d] newSamples [%d]\n", i, newSamples);
+
+ distance = abs(boundary[0].x-boundary[end-2].x) + abs(boundary[0].y-boundary[end-2].y);
+ pBoundaryIndex[i+1].curveClose = distance;
+
+ if(pBoundaryIndex[i+1].curveClose < closureDistance){
+ pBoundaryIndex[i+1].closedCurve = TRUE;
+ }
+ pBoundaryIndex[i+1].centroid.x = 0;
+ pBoundaryIndex[i+1].centroid.y = 0;
+ for(j = 0; j < newSamples; ++j){
+ pBoundaryIndex[i+1].centroid.x += boundary[j].x;
+ pBoundaryIndex[i+1].centroid.y += boundary[j].y;
+ }
+ if(newSamples){
+ pBoundaryIndex[i+1].centroid.x /= newSamples;
+ pBoundaryIndex[i+1].centroid.y /= newSamples;
+ }
+ getBoundaryMetrics(boundary, &length, &minRadius, &maxRadius, &aveRadius,
+ (float)pBoundaryIndex[i+1].centroid.x, (float)pBoundaryIndex[i+1].centroid.y, newSamples);
+ pBoundaryIndex[i+1].boundaryLength = length;
+ pBoundaryIndex[i+1].minRadius = minRadius;
+ pBoundaryIndex[i+1].maxRadius = maxRadius;
+ pBoundaryIndex[i+1].aveRadius = aveRadius;
+ if(minRadius != 0.0){
+ pBoundaryIndex[i+1].ratio = maxRadius / minRadius;
+ }
+ else{
+ pBoundaryIndex[i+1].ratio = -1.0;
+ }
+
+ //
+ // augment the ROI boundary
+ //
+ pBoundaryIndex[i+1].rectangle.left -= 2*CloseMaskSize;
+ pBoundaryIndex[i+1].rectangle.right += 2*CloseMaskSize;
+ pBoundaryIndex[i+1].rectangle.bottom -= 2*CloseMaskSize;
+ pBoundaryIndex[i+1].rectangle.top += 2*CloseMaskSize;
+ label = pBoundaryIndex[i+1].Label;
+
+ //
+ // mask goes in ImageH. morpho filter the mask first
+ //
+ generateMask(ImageH, boundary, newSamples, label, spadSize);
+
+ //
+ // open-close the mask
+ //
+ doMorphology(Input, ImageE, ImageC, ImageH, olapValuesC, olapValuesO, cmask, omask,
+ pBoundaryIndex[i+1].rectangle, label, CloseMaskSize, OpenMaskSize, spadSize);
+
+ //
+ // now get the compactness metrics
+ //
+ getCompactness(Input, pBoundaryIndex[i+1].rectangle, label, spadSize, &vCompactness, length);
+ pBoundaryIndex[i+1].compactness = vCompactness;
+
+ //
+ // reset the ROI boundary
+ //
+ pBoundaryIndex[i+1].rectangle.left += 2*CloseMaskSize;
+ pBoundaryIndex[i+1].rectangle.right -= 2*CloseMaskSize;
+ pBoundaryIndex[i+1].rectangle.bottom += 2*CloseMaskSize;
+ pBoundaryIndex[i+1].rectangle.top -= 2*CloseMaskSize;
+ offset += numSamples;
+ }
+
+
+ for(i = 0; i < numBoundaries; ++i){
+ for(j = 0; j < numBoundaries; ++j){
+ if(j != i){
+ rectPoint[0].x = pBoundaryIndex[j+1].rectangle.left;
+ rectPoint[0].y = pBoundaryIndex[j+1].rectangle.bottom;
+ rectPoint[1].x = pBoundaryIndex[j+1].rectangle.left;
+ rectPoint[1].y = pBoundaryIndex[j+1].rectangle.top;
+ rectPoint[2].x = pBoundaryIndex[j+1].rectangle.right;
+ rectPoint[2].y = pBoundaryIndex[j+1].rectangle.bottom;
+ rectPoint[3].x = pBoundaryIndex[j+1].rectangle.right;
+ rectPoint[3].y = pBoundaryIndex[j+1].rectangle.top;
+ in[0] = 0;
+ in[1] = 0;
+ in[2] = 0;
+ in[3] = 0;
+ for(k = 0; k < 4; ++k){
+ if((rectPoint[k].x > pBoundaryIndex[i+1].rectangle.left) &&
+ (rectPoint[k].x < pBoundaryIndex[i+1].rectangle.right)){
+ if((rectPoint[k].y > pBoundaryIndex[i+1].rectangle.bottom) &&
+ (rectPoint[k].y < pBoundaryIndex[i+1].rectangle.top)){
+ in[k] = 1;
+ }
+ }
+ }
+ if(in[0] && in[1] && in[2] && in[3]){
+ pBoundaryIndex[j+1].isWithin = TRUE;
+ }
+ }
+ }
+ }
+
+ //
+ // fill in the Python features
+ //
+ for(i = 0; i < numBoundaries; ++i){
+ objectMetrics[i].curveClose = pBoundaryIndex[i+1].curveClose;
+ objectMetrics[i].cXBoundary = pBoundaryIndex[i+1].centroid.x;
+ objectMetrics[i].cYBoundary = pBoundaryIndex[i+1].centroid.y;
+ objectMetrics[i].boundaryLength = pBoundaryIndex[i+1].boundaryLength;
+ objectMetrics[i].minRadius = pBoundaryIndex[i+1].minRadius;
+ objectMetrics[i].maxRadius = pBoundaryIndex[i+1].maxRadius;
+ objectMetrics[i].aveRadius = pBoundaryIndex[i+1].aveRadius;
+ objectMetrics[i].ratio = pBoundaryIndex[i+1].ratio;
+ objectMetrics[i].compactness = pBoundaryIndex[i+1].compactness;
+ }
+
+ // debug only
+ if(1){
+ for(i = 0; i < numBoundaries; ++i){
+ if(pBoundaryIndex[i+1].boundaryLength != (float)0.0){
+ printf("boundary %d:\n", i);
+ printf("\t\tRect (%d, %d, %d, %d)\n", pBoundaryIndex[i+1].rectangle.left,
+ pBoundaryIndex[i+1].rectangle.right,
+ pBoundaryIndex[i+1].rectangle.top,
+ pBoundaryIndex[i+1].rectangle.bottom);
+ printf("\t\tCentroid (%d, %d)\n", pBoundaryIndex[i+1].centroid.x, pBoundaryIndex[i+1].centroid.y);
+ printf("\t\tLength (%f)\n", pBoundaryIndex[i+1].boundaryLength);
+ printf("\t\tRatio (%f)\n", pBoundaryIndex[i+1].ratio);
+ printf("\t\taveRadius (%f)\n", pBoundaryIndex[i+1].aveRadius);
+ printf("\t\tLabel (%d)\n", pBoundaryIndex[i+1].Label);
+ printf("\t\tCompactness (%f)\n", pBoundaryIndex[i+1].compactness);
+ printf("\t\tCurveClose (%d)\n", pBoundaryIndex[i+1].curveClose);
+ if(pBoundaryIndex[i+1].isWithin){
+ printf("\t\tContained (T)\n");
+ }
+ else{
+ printf("\t\tContained (F)\n");
+ }
+ if(pBoundaryIndex[i+1].closedCurve){
+ printf("\t\tclosedCurve (T)\n");
+ }
+ else{
+ printf("\t\tclosedCurve (F)\n");
+ }
+ }
+ }
+ }
+
+ //
+ // need to return input which is now mask image
+ //
+
+ offset = 0;
+ offset2 = 0;
+ for(i = 0; i < srcRows; ++i){
+ for(j = 0; j < srcCols; ++j){
+ ThinEdgeImage[offset+j] = (unsigned short)Input[offset2+j];
+ }
+ offset += srcCols;
+ offset2 += spadSize;
+ }
+
+ free(pBoundaryIndex);
+ free(Input);
+ free(ImageE);
+ free(ImageC);
+ free(ImageH);
+ free(pBoundary);
+ free(lBoundary);
+ free(boundary);
+
+ return;
+
+}
+
+
+void initLaws(LawsFilter7 *lawsFilter){
+
+ int i;
+ float sum;
+ float L7[7] = { 1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0};
+ float E7[7] = {-1.0, -4.0, -5.0, 0.0, 5.0, 4.0, 1.0};
+ float S7[7] = {-1.0, -2.0, 1.0, 4.0, 1.0, -2.0, -1.0};
+ float W7[7] = {-1.0, 0.0, 3.0, 0.0, -3.0, 0.0, 1.0};
+ float R7[7] = { 1.0, -2.0, -1.0, 4.0, -1.0, -2.0, 1.0};
+ float O7[7] = {-1.0, 6.0, -15.0, 20.0, -15.0, 6.0, -1.0};
+
+ lawsFilter->numberKernels = 6;
+ lawsFilter->kernelLength = 7;
+ lawsFilter->numberFilterLayers = 21;
+ lawsFilter->name[0] = 'L';
+ lawsFilter->name[1] = 'E';
+ lawsFilter->name[2] = 'S';
+ lawsFilter->name[3] = 'W';
+ lawsFilter->name[4] = 'R';
+ lawsFilter->name[5] = 'O';
+ for(i = 0; i < 7; ++i){
+ lawsFilter->lawsKernel[0][i] = L7[i];
+ lawsFilter->lawsKernel[1][i] = E7[i];
+ lawsFilter->lawsKernel[2][i] = S7[i];
+ lawsFilter->lawsKernel[3][i] = W7[i];
+ lawsFilter->lawsKernel[4][i] = R7[i];
+ lawsFilter->lawsKernel[5][i] = O7[i];
+ }
+
+ // DC filter is unity gain
+ sum = (float)0.0;
+ for(i = 0; i < 7; ++i){
+ sum += lawsFilter->lawsKernel[0][i];
+ }
+ for(i = 0; i < 7; ++i){
+ lawsFilter->lawsKernel[0][i] /= sum;
+ }
+
+ return;
+
+}
+
+float lawsConvolution(float *image, float *rowFilter, float *colFilter, int kernelSize){
+
+ int i, j;
+ int offset;
+ float result[7];
+ float sum;
+
+ // filter rows
+ for(i = 0; i < kernelSize; ++i){
+ sum = (float)0.0;
+ offset = i * kernelSize;
+ for(j = 0; j < kernelSize; ++j){
+ sum += (rowFilter[j]*image[offset+j]);
+ }
+ result[i] = sum;
+ }
+
+ //filter columns
+ sum = (float)0.0;
+ for(j = 0; j < kernelSize; ++j){
+ sum += (rowFilter[j]*result[j]);
+ }
+
+ return(sum);
+
+}
+
+
+void getLawsTexture(LawsFilter7 lawsFilter, tTEM LawsFeatures[], objStruct objectMetrics[], double *sourceImage,
+ unsigned short *MaskImage, int numberObjects, int srcRows, int srcCols){
+
+ int i, j;
+ int label;
+ RECT bBox;
+ int aperature = (lawsFilter.kernelLength-1)/2;
+ unsigned char *ImageH;
+ float *ImageT;
+ float *lawsImage;
+
+ ImageH = calloc(srcRows*srcCols, sizeof(unsigned char));
+ ImageT = calloc(srcRows*srcCols, sizeof(float));
+ lawsImage = calloc(lawsFilter.numberFilterLayers*srcRows*srcCols, sizeof(float));
+
+ for(i = 0; i < numberObjects; ++i){
+ bBox.left = objectMetrics[i].L;
+ bBox.right = objectMetrics[i].R;
+ bBox.top = objectMetrics[i].T;
+ bBox.bottom = objectMetrics[i].B;
+ label = objectMetrics[i].Label;
+ if(objectMetrics[i].voxelMean != (float)0.0){
+ //
+ // valid size region
+ //
+ computeLaws(lawsFilter, LawsFeatures, bBox, label, aperature, srcRows, srcCols, ImageH, ImageT,
+ MaskImage, lawsImage, sourceImage);
+ for(j = 1; j < lawsFilter.numberFilterLayers; ++j){
+ objectMetrics[i].TEM[j-1] = LawsFeatures[j].Variance;
+ }
+ /*
+ int index;
+ int offset;
+ int layerStep = srcRows*srcCols;
+ if(label == debugBlob){
+ index = 0;
+ for(j = 1; j < lawsFilter.numberFilterLayers; ++j){
+ if(LawsFeatures[j].Variance == (float)1.0) index = j;
+ }
+ // overwrite the raw image
+ offset = index * layerStep;
+ for(j = 0; j < layerStep; ++j){
+ sourceImage[j] = lawsImage[offset+j];
+ }
+ }
+ */
+ }
+ }
+
+ free(ImageH);
+ free(ImageT);
+ free(lawsImage);
+
+ return;
+
+}
+
+void computeLaws(LawsFilter7 lawsFilter, tTEM LawsFeatures[], RECT roi, int label, int aperature, int srcRows, int srcCols,
+ unsigned char *ImageH, float *ImageT, unsigned short *MaskImage, float *lawsImage, double *sourceImage){
+
+ //
+ // hard-wirred to Law's 7 kernels
+ //
+ int i, j, k;
+ int lawsLayer;
+ int column, row;
+ int offset;
+ int maskOffset[7];
+ int dataOffset[7];
+ float myImage[49];
+ int count;
+ int outerKernelNumber;
+ int innerKernelNumber;
+ int rowNumber;
+ int kernelSize = lawsFilter.kernelLength;
+ int fullMask = kernelSize*kernelSize;
+ int layerStep = srcRows*srcCols;
+ float *rowFilter;
+ float *colFilter;
+ float filterResult1;
+ float filterResult2;
+ float lawsLL;
+ float t;
+ float maxValue;
+ float scale;
+ char I, J;
+ char combo[24];
+ char dual[24];
+
+
+ // zero the laws mask memory first
+ for(i = 0; i < srcRows*srcCols; ++i){
+ ImageH[i] = 0;
+ }
+ for(j = 0; j < lawsFilter.numberFilterLayers; ++j){
+ LawsFeatures[j].Mean = (float)0.0;
+ LawsFeatures[j].Variance = (float)0.0;
+ }
+
+ for(i = roi.bottom+aperature; i < roi.top-aperature; ++i){
+ // get the row array offset for mask and data source.
+ for(row = -aperature; row <= aperature; ++row){
+ maskOffset[row+aperature] = (i+row)*srcCols;
+ dataOffset[row+aperature] = maskOffset[row+aperature];
+ }
+ for(j = roi.left+aperature; j < roi.right-aperature; ++j){
+ //
+ // get 7x7 segment and make sure have 100% mask coverage
+ //
+ count = 0;
+ for(row = -aperature; row <= aperature; ++row){
+ rowNumber = (row+aperature)*kernelSize;
+ for(column = -aperature; column <= aperature; ++column){
+ if(MaskImage[maskOffset[row+aperature]+j+column] == label){
+ myImage[rowNumber+column+aperature] = sourceImage[dataOffset[row+aperature]+j+column];
+ ++count;
+ }
+ }
+ }
+ if(count == fullMask){
+ //
+ // full mask. 100% coverage. now do the Law's texture filters
+ //
+ ImageH[i*srcCols+j] = 1;
+ lawsLayer = 0;
+ for(outerKernelNumber = 0; outerKernelNumber < lawsFilter.numberKernels; ++outerKernelNumber){
+ //
+ // outer loop pulls the i'th kernel. kernel 0 is the LP kernel
+ // the outer loop is the iso-kernel
+ //
+ I = lawsFilter.name[outerKernelNumber];
+ sprintf(dual, "%c_%c", I, I);
+ rowFilter = &lawsFilter.lawsKernel[outerKernelNumber][0];
+ colFilter = &lawsFilter.lawsKernel[outerKernelNumber][0];
+ filterResult1 = lawsConvolution(myImage, rowFilter, colFilter, kernelSize);
+ // lawsLayer 0 is the LP and needs to be used to scale.
+ if(outerKernelNumber){
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] = (float)2.0 * filterResult1 / lawsLL;
+ }
+ else{
+ lawsLL = (float)2.0 * filterResult1;
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] = (float)2.0 * filterResult1;
+ }
+ strcpy(&LawsFeatures[lawsLayer].filterName[0], dual);
+ ++lawsLayer;
+ //
+ // now do the inner loop and get the column filters for the other laws kernels
+ //
+ for(innerKernelNumber = outerKernelNumber+1; innerKernelNumber < lawsFilter.numberKernels; ++innerKernelNumber){
+ J = lawsFilter.name[innerKernelNumber];
+ sprintf(combo, "%c_%c", I, J);
+ strcpy(&LawsFeatures[lawsLayer].filterName[0], combo);
+ colFilter = &lawsFilter.lawsKernel[innerKernelNumber][0];
+ filterResult1 = lawsConvolution(myImage, rowFilter, colFilter, kernelSize);
+ filterResult2 = lawsConvolution(myImage, colFilter, rowFilter, kernelSize);
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] = (filterResult1 / lawsLL) + (filterResult2 / lawsLL);
+ ++lawsLayer;
+ }
+ }
+ }
+ }
+ }
+
+ for(i = 0; i < lawsFilter.numberFilterLayers; ++i){
+ LawsFeatures[i].Mean = (float)0.0;
+ LawsFeatures[i].Variance = (float)0.0;
+ }
+
+ count = 0;
+ for(i = roi.bottom+aperature; i < roi.top-aperature; ++i){
+ row = i * srcCols;
+ for(j = roi.left+aperature; j < roi.right-aperature; ++j){
+ if(ImageH[row+j]){
+ ++count;
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ offset = k * layerStep + row;
+ LawsFeatures[k].Mean += lawsImage[offset+j];
+ }
+ }
+ }
+ }
+
+ if(count == 0){
+ printf("no samples for texture\n");
+ return;
+ }
+
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ LawsFeatures[k].Mean /= (float)count;
+ }
+ for(i = roi.bottom+aperature; i < roi.top-aperature; ++i){
+ row = i * srcCols;
+ for(j = roi.left+aperature; j < roi.right-aperature; ++j){
+ if(ImageH[row+j]){
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ offset = k * layerStep + row;
+ t = lawsImage[offset+j] - LawsFeatures[k].Mean;
+ LawsFeatures[k].Variance += (t * t);
+ }
+ }
+ }
+ }
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ LawsFeatures[k].Variance /= (float)count;
+ LawsFeatures[k].Variance = (float)(sqrt(LawsFeatures[k].Variance));
+ }
+
+ //
+ // now normalize the variance feature (TEM)
+ //
+ maxValue = (float)0.0;
+ for(i = 1; i < lawsFilter.numberFilterLayers; ++i){
+ if((LawsFeatures[i].Variance) > maxValue) maxValue = LawsFeatures[i].Variance;
+ }
+ scale = (float)1.0 / maxValue;
+ for(i = 1; i < lawsFilter.numberFilterLayers; ++i){
+ LawsFeatures[i].Variance = scale * LawsFeatures[i].Variance;
+ }
+
+
+ return;
+
+}
+
+void getVoxelMeasures(objStruct objectMetrics[], double *sourceImage, unsigned short *MaskImage,
+ int numberObjects, int srcRows, int srcCols){
+
+ int i, j, k;
+ int label;
+ int offset;
+ int count;
+ float mean, std, t;
+ RECT bBox;
+
+ for(i = 0; i < numberObjects; ++i){
+ bBox.left = objectMetrics[i].L;
+ bBox.right = objectMetrics[i].R;
+ bBox.top = objectMetrics[i].T;
+ bBox.bottom = objectMetrics[i].B;
+ label = objectMetrics[i].Label;
+ count = 0;
+ mean = (float)0.0;
+ for(j = bBox.bottom; j < bBox.top; ++j){
+ offset = j * srcCols;
+ for(k = bBox.left; k < bBox.right; ++k){
+ if(MaskImage[offset+k] == label){
+ mean += sourceImage[offset+k];
+ ++count;
+ }
+ }
+ }
+ if(count){
+ mean /= (float)count;
+ std = (float)0.0;
+ for(j = bBox.bottom; j < bBox.top; ++j){
+ offset = j * srcCols;
+ for(k = bBox.left; k < bBox.right; ++k){
+ if(MaskImage[offset+k] == label){
+ t = (sourceImage[offset+k]-mean);
+ std += (t * t);
+ }
+ }
+ }
+ }
+ if(count){
+ std /= (float)count;
+ std = sqrt(std);
+ objectMetrics[i].voxelMean = mean - 2048.0; // the 2048 is only for the cardiac CT volume
+ objectMetrics[i].voxelVar = std;
+ }
+ else{
+ objectMetrics[i].voxelMean = 0.0;
+ objectMetrics[i].voxelVar = 0.0;
+ }
+ if(0) printf("(%d) mean %f, std %f\n", label, objectMetrics[i].voxelMean, objectMetrics[i].voxelVar);
+ }
+
+ return;
+
+}
+
+int NI_BuildBoundary(int samples, int rows, int cols, int numberObjects, unsigned short *edgeImage,
+ objStruct objectMetrics[]){
+
+ int searchWindow = 5; // 5 is good value for Sobel - (should be 13 for Canny edges)
+ int status = 1;
+
+ buildBoundary(objectMetrics, searchWindow, edgeImage, numberObjects, rows, cols);
+
+ return status;
+
+}
+
+int NI_VoxelMeasures(int samples, int rows, int cols, int numberObjects, double *sourceImage,
+ unsigned short *maskImage, objStruct objectMetrics[]){
+
+ int status = 1;
+ getVoxelMeasures(objectMetrics, sourceImage, maskImage, numberObjects, rows, cols);
+
+ return status;
+
+}
+
+
+int NI_TextureMeasures(int samples, int rows, int cols, int numberObjects, double *sourceImage,
+ unsigned short *maskImage, objStruct objectMetrics[]){
+
+ int status = 1;
+ LawsFilter7 lawsFilter;
+ tTEM LawsFeatures[21];
+
+ initLaws(&lawsFilter);
+ getLawsTexture(lawsFilter, LawsFeatures, objectMetrics, sourceImage, maskImage, numberObjects, rows, cols);
+
+ return status;
+
+}
+
+
+
Added: trunk/scipy/ndimage/segment/__init__.py
===================================================================
--- trunk/scipy/ndimage/segment/__init__.py 2007-11-05 23:13:01 UTC (rev 3498)
+++ trunk/scipy/ndimage/segment/__init__.py 2007-11-07 17:53:37 UTC (rev 3499)
@@ -0,0 +1,3 @@
+
+ver = '0.0.1'
+
Added: trunk/scipy/ndimage/segment/ndImage_Segmenter_structs.h
===================================================================
--- trunk/scipy/ndimage/segment/ndImage_Segmenter_structs.h 2007-11-05 23:13:01 UTC (rev 3498)
+++ trunk/scipy/ndimage/segment/ndImage_Segmenter_structs.h 2007-11-07 17:53:37 UTC (rev 3499)
@@ -0,0 +1,155 @@
+#ifndef V1_STRUCTSH
+#define V1_STRUCTSH
+
+#define bool unsigned char
+
+typedef struct{
+ int x;
+ int y;
+}POINT;
+
+typedef struct{
+ int x;
+ int y;
+ int linkIndex;
+ bool haveLink;
+}bPOINT;
+
+typedef struct{
+ int left;
+ int right;
+ int top;
+ int bottom;
+}RECT;
+
+typedef struct{
+ char filterName[20];
+ float Mean;
+ float Variance;
+}tTEM;
+
+typedef struct{
+ int numberKernels;
+ int kernelLength;
+ int numberFilterLayers;
+ float lawsKernel[6][7];
+ char name[7];
+}LawsFilter7;
+
+typedef struct{
+ // filled in GetObjectStats
+ int L;
+ int R;
+ int T;
+ int B;
+ int Label;
+ int Area;
+ float cX;
+ float cY;
+ // filled in BuildBoundary
+ int curveClose;
+ float cXBoundary;
+ float cYBoundary;
+ float boundaryLength;
+ float minRadius;
+ float maxRadius;
+ float aveRadius;
+ float ratio;
+ float compactness;
+ // filled in VoxelMeasures
+ float voxelMean;
+ float voxelVar;
+ // filled in TextureMeasures
+ float TEM[20];
+}objStruct;
+
+typedef struct{
+ int numberPoints;
+ int curveClose;
+ int classify;
+ float boundaryLength;
+ float minRadius;
+ float maxRadius;
+ float aveRadius;
+ float ratio;
+ float compactness;
+ float voxelMean;
+ float voxelVar;
+ RECT rectangle;
+ POINT centroid;
+ bool isWithin;
+ bool closedCurve;
+ bool criticalSize;
+ int Label;
+}boundaryIndex;
+
+
+typedef struct{
+ POINT xy;
+}blobBoundary;
+
+
+//
+// prototypes
+//
+int NI_RegionGrow(int, int, int, int, int, int, int, double *, unsigned short *, int *);
+int NI_TextureMeasures(int, int, int, int, double *, unsigned short *, objStruct objectMetrics[]);
+int NI_VoxelMeasures(int, int, int, int, double *, unsigned short *, objStruct objectMetrics[]);
+int NI_BuildBoundary(int, int, int, int, unsigned short *, objStruct objectMetrics[]);
+int NI_GetObjectStats(int, int, int, unsigned short *, objStruct objectMetrics[]);
+int NI_ThinFilter(int, int, int, int, unsigned short *, objStruct objectMetrics[]);
+int NI_SobelEdges(int, int, int, double, int, int, int, double, int, double *, unsigned short *, int *);
+int NI_ShenCastanEdges(int, int, int, double, double, int, int, int, double *, unsigned short *, int *);
+int NI_CannyEdges(int, int, int, double, double, double, int, int, int, double, int,
+ double *, unsigned short *, int *);
+
+void computeLaws(LawsFilter7, tTEM LawsFeatures[], RECT, int, int, int, int, unsigned char *, float *,
+ unsigned short *, float *, double *);
+float lawsConvolution(float *, float *, float *, int);
+void initLaws(LawsFilter7*);
+void getVoxelMeasures(objStruct objectMetrics[], double *, unsigned short *, int, int, int);
+void getLawsTexture(LawsFilter7, tTEM LawsFeatures[], objStruct objectMetrics[], double *, unsigned short *, int, int, int);
+
+void morphoFilterBinaryImage(int, int, unsigned short *, int, int);
+void buildBinaryImage(int, int, double *, unsigned short *, int, int);
+void doRegionGrow(int, int, int, double *, unsigned short *, int, int, int, int);
+void buildBoundary(objStruct objectMetrics[], int, unsigned short *, int, int, int);
+void getBoundary(unsigned short *, unsigned char *, blobBoundary *, blobBoundary *,
+ boundaryIndex *, RECT, int, int, int, int, int, int);
+void doMorphology(unsigned char *, unsigned char *, unsigned char *, unsigned char *, int olapValuesC[],
+ int olapValuesO[], unsigned short cmask[11][11], unsigned short omask[11][11],
+ RECT, int, int, int, int);
+void getCompactness(unsigned char *, RECT, int, int, float *, float);
+void OpenCloseFilter(int olapValues[], int, int, int, int, unsigned char *,
+ unsigned char *, unsigned short mask[11][11]);
+void trackBoundary(unsigned char *, blobBoundary lBoundary[], int, int, blobBoundary, int);
+void getBoundaryMetrics(bPOINT *, float *, float *, float *, float *, float, float, int);
+void generateMask(unsigned char *, bPOINT *, int, int, int);
+void ThinningFilter(int, int, int, int J_mask[3][30], int K_mask[3][30], unsigned char *,
+ unsigned char *, unsigned char *, unsigned char *, unsigned char *, unsigned char *);
+void initThinFilter(int J_mask[3][30], int K_mask[3][30]);
+void Shen_Castan(double, double, int, int, int, int, int, double *, unsigned short *);
+void computeISEF(float *, float *, int, int, double);
+void ISEF_Horizontal(float *, float *, float *, float *, int, int, double);
+void ISEF_Vertical(float *, float *, float *, float *, int, int, double);
+void thresholdImage(float *, float *, int, int, int, int);
+void computeBandedLaplacian(float *, float *, float *, int, int);
+void getZeroCrossings(float *, float *, float *, int, int, int);
+float adaptiveGradient(float *, float *, int, int, int, int);
+void thresholdEdges(float *, unsigned short *, double, int, int);
+void estimateThreshold(float *, float *, float, int, int, float *);
+void doSobel(int, int, int, double, int, double *, unsigned short *);
+void DGFilters(int, int, int, double, int, float *, float *, double *, double *, float *, float *);
+void nonMaxSupress(int, int, float, float, double *, double *, int, float *, float *, float *);
+void edgeHysteresis(int, int, double, double, float *, float *);
+void edgeThreshold(int, int, double, float *, float *);
+int traceEdge(int, int, int, int, double, float *, float *);
+float magnitude(float, float);
+int ConnectedEdgePoints(int, int, unsigned short *);
+void doPreProcess(int, int, int, double *, double, int, int, int);
+void filter2D(int, int, int, int, int, float *, double *);
+void buildKernel(double, int, int, float *);
+
+
+
+#endif
Added: trunk/scipy/ndimage/segment/setup.py
===================================================================
--- trunk/scipy/ndimage/segment/setup.py 2007-11-05 23:13:01 UTC (rev 3498)
+++ trunk/scipy/ndimage/segment/setup.py 2007-11-07 17:53:37 UTC (rev 3499)
@@ -0,0 +1,21 @@
+
+#!/usr/bin/env python
+
+def configuration(parent_package='',top_path=None):
+ from numpy.distutils.misc_util import Configuration
+
+ config = Configuration('segment', parent_package, top_path)
+
+ config.add_extension('_segmenter',
+ sources=['Segmenter_EXT.c',
+ 'Segmenter_IMPL.c'],
+ depends = ['ndImage_Segmenter_structs.h']
+ )
+
+ return config
+
+if __name__ == '__main__':
+ from numpy.distutils.core import setup
+ setup(**configuration(top_path='').todict())
+
+
Added: trunk/scipy/ndimage/tests/slice112.raw
===================================================================
(Binary files differ)
Property changes on: trunk/scipy/ndimage/tests/slice112.raw
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: trunk/scipy/ndimage/tests/test_segmenter.py
===================================================================
--- trunk/scipy/ndimage/tests/test_segmenter.py 2007-11-05 23:13:01 UTC (rev 3498)
+++ trunk/scipy/ndimage/tests/test_segmenter.py 2007-11-07 17:53:37 UTC (rev 3499)
@@ -0,0 +1,161 @@
+
+import numpy as N
+import scipy.ndimage.segment as S
+
+def shen_castan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048, highThreshold=600+2048, dust=16):
+ datatype = [('L', 'i', 1), ('R', 'i', 1), ('T', 'i', 1), ('B', 'i', 1),
+ ('Label', 'i', 1), ('Area', 'i', 1), ('cX', 'f', 1), ('cY', 'f', 1),
+ ('curveClose', 'i', 1), ('cXB', 'f', 1), ('cYB', 'f', 1), ('bLength', 'f', 1),
+ ('minRadius', 'f', 1), ('maxRadius', 'f', 1), ('aveRadius', 'f', 1), ('ratio', 'f', 1),
+ ('compactness', 'f', 1), ('voxelMean', 'f', 1), ('voxelVar', 'f', 1), ('TEM', 'f', 20)]
+ labeledEdges, numberObjects = S.shen_castan_edges(scLow, IIRFilter, window, lowThreshold, highThreshold, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=datatype)
+ # return the bounding box for each connected edge
+ S.get_object_stats(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def sobel(image, sLow=0.3, tMode=1, lowThreshold=220+2048, highThreshold=600+2048, BPHigh=10.0, apearture=21, dust=16):
+ datatype = [('L', 'i', 1), ('R', 'i', 1), ('T', 'i', 1), ('B', 'i', 1),
+ ('Label', 'i', 1), ('Area', 'i', 1), ('cX', 'f', 1), ('cY', 'f', 1),
+ ('curveClose', 'i', 1), ('cXB', 'f', 1), ('cYB', 'f', 1), ('bLength', 'f', 1),
+ ('minRadius', 'f', 1), ('maxRadius', 'f', 1), ('aveRadius', 'f', 1), ('ratio', 'f', 1),
+ ('compactness', 'f', 1), ('voxelMean', 'f', 1), ('voxelVar', 'f', 1), ('TEM', 'f', 20)]
+ # get sobel edge points. return edges that are labeled (1..numberObjects)
+ labeledEdges, numberObjects = S.sobel_edges(sLow, tMode, lowThreshold, highThreshold, BPHigh, apearture, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=datatype)
+ # return the bounding box for each connected edge
+ S.get_object_stats(labeledEdges, ROIList)
+ # thin (medial axis transform) of the sobel edges as the sobel produces a 'band edge'
+ S.morpho_thin_filt(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def canny(image, cSigma=1.0, cLow=0.5, cHigh=0.8, tMode=1, lowThreshold=220+2048, highThreshold=600+2048,
+ BPHigh=10.0, apearture=21, dust=16):
+ datatype = [('L', 'i', 1), ('R', 'i', 1), ('T', 'i', 1), ('B', 'i', 1),
+ ('Label', 'i', 1), ('Area', 'i', 1), ('cX', 'f', 1), ('cY', 'f', 1),
+ ('curveClose', 'i', 1), ('cXB', 'f', 1), ('cYB', 'f', 1), ('bLength', 'f', 1),
+ ('minRadius', 'f', 1), ('maxRadius', 'f', 1), ('aveRadius', 'f', 1), ('ratio', 'f', 1),
+ ('compactness', 'f', 1), ('voxelMean', 'f', 1), ('voxelVar', 'f', 1), ('TEM', 'f', 20)]
+ # get canny edge points. return edges that are labeled (1..numberObjects)
+ labeledEdges, numberObjects = S.canny_edges(cSigma, cLow, cHigh, tMode, lowThreshold, highThreshold,
+ BPHigh, apearture, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=datatype)
+ # return the bounding box for each connected edge
+ S.get_object_stats(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def get_shape_mask(labeledEdges, ROIList):
+ # pass in Sobel morph-thinned labeled edge image (LEI) and ROIList
+ # GetShapeMask will augment the ROI list
+ # labeledEdges is the original edge image and overwritten as mask image
+ # maskImage is the mask that is used for blob texture / pixel features
+ S.build_boundary(labeledEdges, ROIList)
+ return
+
+def get_voxel_measures(rawImage, labeledEdges, ROIList):
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the voxel features in the list
+ #
+ S.voxel_measures(rawImage, labeledEdges, ROIList)
+ return
+
+def get_texture_measures(rawImage, labeledEdges, ROIList):
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the texture (Law's, co-occurence, Gabor) features in the list
+ #
+ S.texture_measures(rawImage, labeledEdges, ROIList)
+ return
+
+def segment_regions():
+ # get slice from the CT volume
+ image = get_slice('slice112.raw')
+ # need a copy of original image as filtering will occur on the extracted slice
+ sourceImage = image.copy()
+ # Sobel is the first level segmenter. Sobel magnitude and MAT (medial axis transform)
+ # followed by connected component analysis. What is returned is labeled edges and the object list
+ labeledMask, ROIList = sobel(image)
+ # From the labeled edges and the object list get the labeled mask for each blob object
+ get_shape_mask(labeledMask, ROIList)
+ # Use the labeled mask and source image (raw) to get voxel features
+ get_voxel_measures(sourceImage, labeledMask, ROIList)
+ # Use the labeled mask and source image (raw) to get texture features
+ get_texture_measures(sourceImage, labeledMask, ROIList)
+ return sourceImage, labeledMask, ROIList
+
+def grow_regions():
+ # get slice from the CT volume
+ image = get_slice('slice112.raw')
+ regionMask, numberRegions = region_grow(image)
+ return regionMask, numberRegions
+
+
+def region_grow(image, lowThreshold=220+2048, highThreshold=600+2048, open=7, close=7):
+ # morphology filters need to be clipped to 11 max and be odd
+ regionMask, numberRegions = S.region_grow(lowThreshold, highThreshold, close, open, image)
+ return regionMask, numberRegions
+
+
+def get_slice(imageName='junk.raw', bytes=2, rows=512, columns=512):
+ # get a slice alrady extracted from the CT volume
+ #image = open(imageName, 'rb')
+ #slice = image.read(rows*columns*bytes)
+ #values = struct.unpack('h'*rows*columns, slice)
+ #ImageSlice = N.array(values, dtype=float).reshape(rows, columns)
+
+ ImageSlice = N.fromfile(imageName, dtype=N.uint16).reshape(rows, columns);
+
+ # clip the ends for this test CT image file as the spine runs off the end of the image
+ ImageSlice[505:512, :] = 0
+ return (ImageSlice).astype(float)
+
+def get_slice2(image_name='junk.raw', bytes=2, shape=(512,512)):
+ import mmap
+ file = open(image_name, 'rb')
+ mm = mmap.mmap(file.fileno(), 0, access=mmap.ACCESS_READ)
+ slice = N.frombuffer(mm, dtype='u%d' % bytes).reshape(shape)
+ slice = slice.astype(float)
+ slice[505:512,:] = 0
+ return slice
+
+def save_slice(mySlice, filename='junk.raw', bytes=4):
+ # just save the slice to a fixed file
+ slice = mySlice.astype('u%d' % bytes)
+ slice.tofile(filename)
+
+
+# test 1
+def test1():
+ import Segmenter as S
+ image = S.get_slice('slice112.raw')
+ sourceImage = image.copy()
+ edges, objects = S.sobel(image)
+ S.get_shape_mask(edges, objects)
+ S.get_voxel_measures(sourceImage, edges, objects)
+ S.get_texture_measures(sourceImage, edges, objects)
+
+# test 2
+def test2():
+ import volumeInput as V
+ import Segmenter as S
+ sourceImage, labeledMask, ROIList = S.segment_regions()
+
+
+# test 3
+def test3():
+ import Segmenter as S
+ regionMask, numberRegions = S.grow_regions()
+ regionMask.max()
+ S.save_slice(regionMask, 'regionMask.raw')
+
+
+
+class TestSegmenter(NumpyTestCase):
+ def test_1(self):
+
+
+
More information about the Scipy-svn
mailing list