[Scipy-svn] r3916 - in trunk/scipy/ndimage: . src src/segment src/segment/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Feb 11 18:32:21 EST 2008
Author: chris.burns
Date: 2008-02-11 17:32:11 -0600 (Mon, 11 Feb 2008)
New Revision: 3916
Added:
trunk/scipy/ndimage/src/segment/
trunk/scipy/ndimage/src/segment/Segmenter_EXT.c
trunk/scipy/ndimage/src/segment/Segmenter_IMPL.c
trunk/scipy/ndimage/src/segment/ndImage_Segmenter_structs.h
trunk/scipy/ndimage/src/segment/objectdata.py
trunk/scipy/ndimage/src/segment/tests/
Removed:
trunk/scipy/ndimage/segment/
trunk/scipy/ndimage/src/segment/Segmenter_EXT.c
trunk/scipy/ndimage/src/segment/Segmenter_IMPL.c
trunk/scipy/ndimage/src/segment/__init__.py
trunk/scipy/ndimage/src/segment/ndImage_Segmenter_structs.h
trunk/scipy/ndimage/src/segment/objectdata.py
trunk/scipy/ndimage/src/segment/setup.py
trunk/scipy/ndimage/src/segment/tests/
Modified:
trunk/scipy/ndimage/segmenter.py
trunk/scipy/ndimage/setup.py
trunk/scipy/ndimage/src/segment/tests/test_segment.py
Log:
Reorg package structure for ndimage.segment.
Modified: trunk/scipy/ndimage/segmenter.py
===================================================================
--- trunk/scipy/ndimage/segmenter.py 2008-02-11 23:27:05 UTC (rev 3915)
+++ trunk/scipy/ndimage/segmenter.py 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,6 +1,6 @@
import math
import numpy as N
-import scipy.ndimage.segment as S
+import scipy.ndimage._segment as S
# make sure this is local to use as default
inputname = 'slice112.raw'
Modified: trunk/scipy/ndimage/setup.py
===================================================================
--- trunk/scipy/ndimage/setup.py 2008-02-11 23:27:05 UTC (rev 3915)
+++ trunk/scipy/ndimage/setup.py 2008-02-11 23:32:11 UTC (rev 3916)
@@ -14,7 +14,13 @@
include_dirs=['src']+[get_include()],
)
- config.add_subpackage('segment')
+ config.add_extension('_segment',
+ sources=['src/segment/Segmenter_EXT.c',
+ 'src/segment/Segmenter_IMPL.c'],
+ depends = ['src/segment/ndImage_Segmenter_structs.h']
+ )
+
+ #config.add_subpackage('segment')
config.add_data_dir('tests')
config.add_subpackage('register')
Copied: trunk/scipy/ndimage/src/segment (from rev 3913, trunk/scipy/ndimage/segment)
Deleted: trunk/scipy/ndimage/src/segment/Segmenter_EXT.c
===================================================================
--- trunk/scipy/ndimage/segment/Segmenter_EXT.c 2008-02-11 13:46:25 UTC (rev 3913)
+++ trunk/scipy/ndimage/src/segment/Segmenter_EXT.c 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,460 +0,0 @@
-#include "ndImage_Segmenter_structs.h"
-#include "Python.h"
-#include "numpy/arrayobject.h"
-
-static PyObject *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_ParseTuple(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;
- itype = NPY_USHORT;
- 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 *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_ParseTuple(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
- itype = NPY_USHORT;
- 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 *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_ParseTuple(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
- itype = NPY_USHORT;
- 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 *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_ParseTuple(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 *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_ParseTuple(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 *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_ParseTuple(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 *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_ParseTuple(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 *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_ParseTuple(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 *Segmenter_RegionGrow(PyObject *self, PyObject *args)
-{
-
- int lowThreshold;
- int highThreshold;
- int closeWindow;
- int openWindow;
- int num;
- int nd;
- int type;
- int itype;
- int groups;
- npy_intp *dims;
- double *fP1;
- unsigned short *fP2;
- PyObject *iArray = NULL;
- PyObject *eArray = NULL;
-
- //
- // pass in 2D LPF coefficients
- if(!PyArg_ParseTuple(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
- itype = NPY_USHORT;
- 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 SegmenterMethods[] =
-{
- { "canny_edges", Segmenter_CannyEdges, METH_VARARGS, NULL },
- { "shen_castan_edges", Segmenter_ShenCastanEdges, METH_VARARGS, NULL },
- { "sobel_edges", Segmenter_SobelEdges, METH_VARARGS, NULL },
- { "get_object_stats", Segmenter_GetObjectStats, METH_VARARGS, NULL },
- { "morpho_thin_filt", Segmenter_MorphoThinFilt, METH_VARARGS, NULL },
- { "build_boundary", Segmenter_BuildBoundary, METH_VARARGS, NULL },
- { "voxel_measures", Segmenter_VoxelMeasures, METH_VARARGS, NULL },
- { "texture_measures", Segmenter_TextureMeasures, METH_VARARGS, NULL },
- { "region_grow", Segmenter_RegionGrow, METH_VARARGS, NULL },
- { NULL, NULL, 0, NULL},
-};
-
-void init_segmenter(void)
-{
- Py_InitModule("_segmenter", SegmenterMethods);
- import_array();
-}
-
Copied: trunk/scipy/ndimage/src/segment/Segmenter_EXT.c (from rev 3915, trunk/scipy/ndimage/segment/Segmenter_EXT.c)
===================================================================
--- trunk/scipy/ndimage/segment/Segmenter_EXT.c 2008-02-11 23:27:05 UTC (rev 3915)
+++ trunk/scipy/ndimage/src/segment/Segmenter_EXT.c 2008-02-11 23:32:11 UTC (rev 3916)
@@ -0,0 +1,460 @@
+#include "ndImage_Segmenter_structs.h"
+#include "Python.h"
+#include "numpy/arrayobject.h"
+
+static PyObject *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_ParseTuple(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;
+ itype = NPY_USHORT;
+ 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 *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_ParseTuple(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
+ itype = NPY_USHORT;
+ 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 *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_ParseTuple(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
+ itype = NPY_USHORT;
+ 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 *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_ParseTuple(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 *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_ParseTuple(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 *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_ParseTuple(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 *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_ParseTuple(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 *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_ParseTuple(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 *Segmenter_RegionGrow(PyObject *self, PyObject *args)
+{
+
+ int lowThreshold;
+ int highThreshold;
+ int closeWindow;
+ int openWindow;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int groups;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_ParseTuple(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
+ itype = NPY_USHORT;
+ 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 SegmenterMethods[] =
+{
+ { "canny_edges", Segmenter_CannyEdges, METH_VARARGS, NULL },
+ { "shen_castan_edges", Segmenter_ShenCastanEdges, METH_VARARGS, NULL },
+ { "sobel_edges", Segmenter_SobelEdges, METH_VARARGS, NULL },
+ { "get_object_stats", Segmenter_GetObjectStats, METH_VARARGS, NULL },
+ { "morpho_thin_filt", Segmenter_MorphoThinFilt, METH_VARARGS, NULL },
+ { "build_boundary", Segmenter_BuildBoundary, METH_VARARGS, NULL },
+ { "voxel_measures", Segmenter_VoxelMeasures, METH_VARARGS, NULL },
+ { "texture_measures", Segmenter_TextureMeasures, METH_VARARGS, NULL },
+ { "region_grow", Segmenter_RegionGrow, METH_VARARGS, NULL },
+ { NULL, NULL, 0, NULL},
+};
+
+void init_segment(void)
+{
+ Py_InitModule("_segment", SegmenterMethods);
+ import_array();
+}
+
Deleted: trunk/scipy/ndimage/src/segment/Segmenter_IMPL.c
===================================================================
--- trunk/scipy/ndimage/segment/Segmenter_IMPL.c 2008-02-11 13:46:25 UTC (rev 3913)
+++ trunk/scipy/ndimage/src/segment/Segmenter_IMPL.c 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,3020 +0,0 @@
-#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 */
- }
- }
- }
- }
- 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 */
- doPreProcess(samples, rows, cols, rawImage, BPHigh, apearture, lowThreshold, highThreshold);
-
- /*
- // 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;
- unsigned short *omask;
- int olapValuesC[4];
- int olapValuesO[4];
- int CloseMaskSize = 1;
- int OpenMaskSize = 1;
- int LowValue1, HighValue1;
- int LowValue2, HighValue2;
- int spadSize;
- int maskSize = 11;
- unsigned char *ImageE;
- unsigned char *ImageC;
-
- spadSize = MAX(rows, cols);
-
- ImageE = calloc(spadSize*spadSize, sizeof(unsigned char));
- ImageC = calloc(spadSize*spadSize, sizeof(unsigned char));
-
- cmask = calloc(11*11, sizeof(unsigned short));
- omask = calloc(11*11, sizeof(unsigned short));
-
- //
- // 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*maskSize+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*maskSize+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);
-
- free(cmask);
- free(omask);
-
- 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, int *K_mask){
-
- int i, j;
- int Column;
- int maskCols = 3;
-
- for(i = 0; i < 3; ++i){
- for(j = 0; j < 30; ++j){
- J_mask[i+j*maskCols] = 0;
- K_mask[i+j*maskCols] = 0;
- }
- }
-
- Column = 0;
- J_mask[0+maskCols*(Column+0)] = 1;
- J_mask[0+maskCols*(Column+1)] = 1;
- J_mask[0+maskCols*(Column+2)] = 1;
- J_mask[1+maskCols*(Column+1)] = 1;
-
- Column += 3;
- J_mask[0+maskCols*(Column+1)] = 1;
- J_mask[1+maskCols*(Column+1)] = 1;
- J_mask[1+maskCols*(Column+2)] = 1;
-
- Column += 3;
- J_mask[0+maskCols*(Column+0)] = 1;
- J_mask[1+maskCols*(Column+0)] = 1;
- J_mask[2+maskCols*(Column+0)] = 1;
- J_mask[1+maskCols*(Column+1)] = 1;
-
- Column += 3;
- J_mask[0+maskCols*(Column+1)] = 1;
- J_mask[1+maskCols*(Column+0)] = 1;
- J_mask[1+maskCols*(Column+1)] = 1;
-
- Column += 3;
- J_mask[0+maskCols*(Column+2)] = 1;
- J_mask[1+maskCols*(Column+1)] = 1;
- J_mask[1+maskCols*(Column+2)] = 1;
- J_mask[2+maskCols*(Column+2)] = 1;
-
- Column += 3;
- J_mask[1+maskCols*(Column+0)] = 1;
- J_mask[1+maskCols*(Column+1)] = 1;
- J_mask[2+maskCols*(Column+1)] = 1;
-
- Column += 3;
- J_mask[1+maskCols*(Column+1)] = 1;
- J_mask[2+maskCols*(Column+0)] = 1;
- J_mask[2+maskCols*(Column+1)] = 1;
- J_mask[2+maskCols*(Column+2)] = 1;
-
- Column += 3;
- J_mask[1+maskCols*(Column+1)] = 1;
- J_mask[1+maskCols*(Column+2)] = 1;
- J_mask[2+maskCols*(Column+1)] = 1;
-
- Column = 0;
- K_mask[2+maskCols*(Column+0)] = 1;
- K_mask[2+maskCols*(Column+1)] = 1;
- K_mask[2+maskCols*(Column+2)] = 1;
-
- Column += 3;
- K_mask[1+maskCols*(Column+0)] = 1;
- K_mask[2+maskCols*(Column+0)] = 1;
- K_mask[2+maskCols*(Column+1)] = 1;
-
- Column += 3;
- K_mask[0+maskCols*(Column+2)] = 1;
- K_mask[1+maskCols*(Column+2)] = 1;
- K_mask[2+maskCols*(Column+2)] = 1;
-
- Column += 3;
- K_mask[1+maskCols*(Column+2)] = 1;
- K_mask[2+maskCols*(Column+1)] = 1;
- K_mask[2+maskCols*(Column+2)] = 1;
-
- Column += 3;
- K_mask[0+maskCols*(Column+0)] = 1;
- K_mask[1+maskCols*(Column+0)] = 1;
- K_mask[2+maskCols*(Column+0)] = 1;
-
- Column += 3;
- K_mask[0+maskCols*(Column+1)] = 1;
- K_mask[0+maskCols*(Column+2)] = 1;
- K_mask[1+maskCols*(Column+2)] = 1;
-
- Column += 3;
- K_mask[0+maskCols*(Column+0)] = 1;
- K_mask[0+maskCols*(Column+1)] = 1;
- K_mask[0+maskCols*(Column+2)] = 1;
-
- Column += 3;
- K_mask[0+maskCols*(Column+0)] = 1;
- K_mask[0+maskCols*(Column+1)] = 1;
- K_mask[1+maskCols*(Column+0)] = 1;
-
- return;
-
-}
-
-void ThinningFilter(int regRows, int regColumns, int spadSize, int *J_mask, int *K_mask,
- 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 maskCols = 3;
- 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+maskCols*(Column+j)];
- }
- }
- for(i = 0; i < 3; ++i){
- for(j = 0; j < 3; ++j){
- k_mask[i][j] = K_mask[i+maskCols*(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 no change */
- 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 is data dependent. */
- ++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;
- int *K_mask;
-
- 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));
- J_mask = calloc(3*30, sizeof(int));
- K_mask = calloc(3*30, sizeof(int));
-
- 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);
- free(J_mask);
- free(K_mask);
-
- 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;
- 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){
-
-
- /*
- // 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;
- int morphoMaskSize = 11;
-
- 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[morphoMaskSize*(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[morphoMaskSize*(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, unsigned short *omask,
- 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;
- pBoundary[memOffset].xy.x = value.xy.x;
- pBoundary[memOffset].xy.y = value.xy.y;
- ++memOffset;
- }
-
- 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 will be user-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;
- unsigned short *omask;
- int maskSize = 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;
-
- 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));
- cmask = calloc(11*11, sizeof(unsigned short));
- omask = calloc(11*11, sizeof(unsigned short));
-
- /*
- // Close filter
- */
- CloseMaskSize = (CloseSize-1)/2;
- for(i = 0; i < 2*CloseMaskSize+1; ++i){
- for(j = 0; j < 2*CloseMaskSize+1; ++j){
- cmask[i*maskSize+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*maskSize+j] = 1;
- }
- }
- LowValue1 = 1;
- HighValue1 = 0;
- LowValue2 = 0;
- HighValue2 = 1;
- olapValuesO[0] = LowValue1;
- olapValuesO[1] = HighValue1;
- olapValuesO[2] = LowValue2;
- olapValuesO[3] = HighValue2;
-
- 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){
- ++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;
- 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
- */
-
- offset = 0;
- numBoundaries = pBoundaryIndex[0].numberPoints;
- 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;
- }
- }
-
- 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(0){
- 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);
- free(cmask);
- free(omask);
-
- 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];
- }
-
- /* L 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;
- }
- /* -- later will need to return a view of the texture images
- 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=1.0;
- 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){
- /*
- // 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){
- // debug statement
- 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;
- objectMetrics[i].voxelVar = std;
- }
- else{
- objectMetrics[i].voxelMean = 0.0;
- objectMetrics[i].voxelVar = 0.0;
- }
- }
-
- 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
- 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;
-
-}
-
-
-
Copied: trunk/scipy/ndimage/src/segment/Segmenter_IMPL.c (from rev 3915, trunk/scipy/ndimage/segment/Segmenter_IMPL.c)
Deleted: trunk/scipy/ndimage/src/segment/__init__.py
===================================================================
--- trunk/scipy/ndimage/segment/__init__.py 2008-02-11 13:46:25 UTC (rev 3913)
+++ trunk/scipy/ndimage/src/segment/__init__.py 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,5 +0,0 @@
-# Segmentation package
-# Author: Tom Waite, 2007
-
-from _segmenter import *
-from objectdata import *
Deleted: trunk/scipy/ndimage/src/segment/ndImage_Segmenter_structs.h
===================================================================
--- trunk/scipy/ndimage/segment/ndImage_Segmenter_structs.h 2008-02-11 13:46:25 UTC (rev 3913)
+++ trunk/scipy/ndimage/src/segment/ndImage_Segmenter_structs.h 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,155 +0,0 @@
-#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 *, unsigned short *,
- 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 *);
-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 *, int *, unsigned char *,
- unsigned char *, unsigned char *, unsigned char *, unsigned char *, unsigned char *);
-void initThinFilter(int *, int *);
-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
Copied: trunk/scipy/ndimage/src/segment/ndImage_Segmenter_structs.h (from rev 3915, trunk/scipy/ndimage/segment/ndImage_Segmenter_structs.h)
Deleted: trunk/scipy/ndimage/src/segment/objectdata.py
===================================================================
--- trunk/scipy/ndimage/segment/objectdata.py 2008-02-11 13:46:25 UTC (rev 3913)
+++ trunk/scipy/ndimage/src/segment/objectdata.py 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,25 +0,0 @@
-
-import numpy as N
-
-objstruct =N.dtype([('L', 'i'),
- ('R', 'i'),
- ('T', 'i'),
- ('B', 'i'),
- ('Label', 'i'),
- ('Area', 'i'),
- ('cX', 'f'),
- ('cY', 'f'),
- ('curveClose', 'i'),
- ('cXB', 'f'),
- ('cYB', 'f'),
- ('bLength', 'f'),
- ('minRadius', 'f'),
- ('maxRadius', 'f'),
- ('aveRadius', 'f'),
- ('ratio', 'f'),
- ('compactness', 'f'),
- ('voxelMean', 'f'),
- ('voxelVar', 'f'),
- ('TEM', 'f', 20)]
- )
-
Copied: trunk/scipy/ndimage/src/segment/objectdata.py (from rev 3915, trunk/scipy/ndimage/segment/objectdata.py)
Deleted: trunk/scipy/ndimage/src/segment/setup.py
===================================================================
--- trunk/scipy/ndimage/segment/setup.py 2008-02-11 13:46:25 UTC (rev 3913)
+++ trunk/scipy/ndimage/src/segment/setup.py 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,23 +0,0 @@
-
-#!/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']
- )
-
- config.add_data_dir('tests')
-
- return config
-
-if __name__ == '__main__':
- from numpy.distutils.core import setup
- setup(**configuration(top_path='').todict())
-
-
Copied: trunk/scipy/ndimage/src/segment/tests (from rev 3915, trunk/scipy/ndimage/segment/tests)
Modified: trunk/scipy/ndimage/src/segment/tests/test_segment.py
===================================================================
--- trunk/scipy/ndimage/segment/tests/test_segment.py 2008-02-11 23:27:05 UTC (rev 3915)
+++ trunk/scipy/ndimage/src/segment/tests/test_segment.py 2008-02-11 23:32:11 UTC (rev 3916)
@@ -1,313 +1,313 @@
-
-import numpy as N
+
+import numpy as N
from scipy.testing import *
-import scipy.ndimage.segment as S
-
-inputname = 'slice112.raw'
-
-import os
-filename = os.path.join(os.path.split(__file__)[0],inputname)
-
-
-def shen_castan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048,
- highThreshold=600+2048, dust=16):
- """
- labeledEdges, ROIList = shen_castan(image, [default])
-
- implements Shen-Castan edge finding
-
- Inputs - image, IIR filter, shen_castan_low, window, low_threshold, high_threshold, dust
- - image is the numarray 2D image
- - IIR filter is filter parameter for exponential filter
- - shen_castan_low is edge threshold is range (0.0, 1.0]
- - window is search window for edge detection
- - low_ and high_ threshold are density values
- - dust is blob filter. blob area (length x width of bounding box) under this
- size threshold are filtered (referred to as dust and blown away)
-
- Outputs - labeledEdges, ROIList[>dust]
- - labeledEdges is boundary (edges) of segmented 'blobs',
- numerically labeled by blob number
- - ROIList[>dust] is a blob feature list. Only values
- with bounding box area greater than dust threshold are returned
-
- """
- 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=S.objstruct)
- # 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):
- """
- labeledEdges, ROIList = sobel(image, [default])
-
- implements sobel magnitude edge finding
-
- Inputs - image, sobel_low, tMode, low_threshold, high_threshold,
- high_filter_cutoff, filter_aperature, dust
- - image is the numarray 2D image
- - sobel_low is edge threshold is range (0.0, 1.0]
- - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
- - low_ and high_ threshold are density values
- - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
- - aperature is odd filter kernel length
- - dust is blob filter. blob area (length x width of bounding box) under this
- size threshold are filtered (referred to as dust and blown away)
-
- Outputs - labeledEdges, ROIList[>dust]
- - labeledEdges is boundary (edges) of segmented 'blobs',
- numerically labeled by blob number
- - ROIList[>dust] is a blob feature list. Only values
- with bounding box area greater than dust threshold are returned
-
- """
- # 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=S.objstruct)
- # 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):
- """
- labeledEdges, ROIList = canny(image, [default])
-
- implements canny edge finding
-
- Inputs - image, DG_sigma, canny_low, canny_high, tMode, low_threshold,
- high_threshold, high_filter_cutoff, filter_aperature, dust
- - image is the numarray 2D image
- - DG_sigma is Gaussain sigma for the derivative-of-gaussian filter
- - clow is low edge threshold is range (0.0, 1.0]
- - chigh is high edge threshold is range (0.0, 1.0]
- - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
- - low_ and high_ threshold are density values
- - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
- - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
- - aperature is odd filter kernel length
- - dust is blob filter. blob area (length x width of bounding box) under this
- size threshold are filtered (referred to as dust and blown away)
-
- Outputs - labeledEdges, ROIList[>dust]
- - labeledEdges is boundary (edges) of segmented 'blobs',
- numerically labeled by blob number
- - ROIList[>dust] is a blob feature list. Only values
- with bounding box area greater than dust threshold are returned
-
- """
- # 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=S.objstruct)
- # 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):
- """
- get_shape_mask(labeledEdges, ROIList)
-
- takes labeled edge image plus ROIList (blob descriptors) and generates
- boundary shape features and builds labeled blob masks. 'labeledEdges'
- is over-written by 'labeledMask'. Adds features to ROIList structure
-
- Inputs - labeledEdges, ROIList
- - labeledEdges is boundary (edges) of segmented 'blobs',
- numerically labeled by blob number
- - ROIList is a blob feature list.
-
- Output - no return. edge image input is over-written with mask image.
- ROIList added to.
-
- """
-
- # 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):
- """
- get_voxel_measures(rawImage, labeledEdges, ROIList)
-
- takes raw 2D image, labeled blob mask and ROIList. computes voxel features
- (moments, histogram) for each blob. Adds features to ROIList structure.
-
- Inputs - rawImage, labeledEdges, ROIList
- - rawImage is the original source 2D image
- - labeledEdges is boundary (edges) of segmented 'blobs',
- numerically labeled by blob number
- - ROIList is a blob feature list.
-
- Output - no return. ROIList added to.
-
- """
- #
- # 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):
- """
- get_texture_measures(rawImage, labeledEdges, ROIList)
-
- takes raw 2D image, labeled blob mask and ROIList. computes 2D
- texture features using 7x7 Law's texture filters applied
- to segmented blobs. TEM (texture energy metric) is computed
- for each Law's filter image and stored in TEM part of ROIList.
-
- Inputs - rawImage, labeledEdges, ROIList
- - rawImage is the original source 2D image
- - labeledEdges is boundary (edges) of segmented 'blobs',
- numerically labeled by blob number
- - ROIList is a blob feature list.
-
- Output - no return. ROIList added to.
- """
- #
- # pass raw image, labeled mask and the partially filled ROIList
- # VoxelMeasures will fill the texture (Law's, sub-edges, co-occurence, Gabor) features in the list
- #
- S.texture_measures(rawImage, labeledEdges, ROIList)
- return
-
-def segment_regions():
- """
- sourceImage, labeledMask, ROIList = segment_regions()
-
- Inputs - No Input
-
- Outputs - sourceImage, labeledMask, ROIList
- - sourceImage is raw 2D image (default cardiac CT slice for demo
- - labeledMask is mask of segmented 'blobs',
- numerically labeled by blob number
- - ROIList is numerical Python structure of intensity, shape and
- texture features for each blob
-
- High level script calls Python functions:
- get_slice() - a cardiac CT slice demo file
- sobel() - sobel magnitude edge finder,
- returns connected edges
- get_shape_mask() - gets segmented blob boundary and mask
- and shape features
- get_voxel_measures() - uses masks get object voxel moment
- and histogram features
- get_texture_measures() - uses masks get object 2D texture features
- """
- # get slice from the CT volume
- image = get_slice(filename)
- # 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():
- """
- regionMask, numberRegions = region_grow()
- Inputs - No Input
- Outputs - regionMask, numberRegions
- - regionMask is the labeled segment masks from 2D image
- - numberRegions is the number of segmented blobs
-
- High level script calls Python functions:
- get_slice() - a cardiac CT slice demo file
- region_grow() - "grows" connected blobs. default threshold
- and morphological filter structuring element
- """
- # get slice from the CT volume
- image = get_slice(filename)
- regionMask, numberRegions = region_grow(image)
- return regionMask, numberRegions
-
-
-def region_grow(image, lowThreshold=220+2048, highThreshold=600+2048, open=7, close=7):
- """
- regionMask, numberRegions = region_grow(image, [defaults])
-
- Inputs - image, low_threshold, high_threshold, open, close
- - image is the numarray 2D image
- - low_ and high_ threshold are density values
- - open is open morphology structuring element
- odd size. 0 to turn off. max is 11
- - close is close morphology structuring element
- odd size. 0 to turn off. max is 11
-
- Outputs - regionMask, numberRegions
- - regionMask is the labeled segment masks from 2D image
- - numberRegions is the number of segmented blobs
- """
- # 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='slice112.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='slice112.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)
-
-
+import scipy.ndimage._segment as S
+
+inputname = 'slice112.raw'
+
+import os
+filename = os.path.join(os.path.split(__file__)[0],inputname)
+
+
+def shen_castan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048,
+ highThreshold=600+2048, dust=16):
+ """
+ labeledEdges, ROIList = shen_castan(image, [default])
+
+ implements Shen-Castan edge finding
+
+ Inputs - image, IIR filter, shen_castan_low, window, low_threshold, high_threshold, dust
+ - image is the numarray 2D image
+ - IIR filter is filter parameter for exponential filter
+ - shen_castan_low is edge threshold is range (0.0, 1.0]
+ - window is search window for edge detection
+ - low_ and high_ threshold are density values
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ 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=S.objstruct)
+ # 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):
+ """
+ labeledEdges, ROIList = sobel(image, [default])
+
+ implements sobel magnitude edge finding
+
+ Inputs - image, sobel_low, tMode, low_threshold, high_threshold,
+ high_filter_cutoff, filter_aperature, dust
+ - image is the numarray 2D image
+ - sobel_low is edge threshold is range (0.0, 1.0]
+ - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
+ - low_ and high_ threshold are density values
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - aperature is odd filter kernel length
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ # 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=S.objstruct)
+ # 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):
+ """
+ labeledEdges, ROIList = canny(image, [default])
+
+ implements canny edge finding
+
+ Inputs - image, DG_sigma, canny_low, canny_high, tMode, low_threshold,
+ high_threshold, high_filter_cutoff, filter_aperature, dust
+ - image is the numarray 2D image
+ - DG_sigma is Gaussain sigma for the derivative-of-gaussian filter
+ - clow is low edge threshold is range (0.0, 1.0]
+ - chigh is high edge threshold is range (0.0, 1.0]
+ - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
+ - low_ and high_ threshold are density values
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - aperature is odd filter kernel length
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ # 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=S.objstruct)
+ # 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):
+ """
+ get_shape_mask(labeledEdges, ROIList)
+
+ takes labeled edge image plus ROIList (blob descriptors) and generates
+ boundary shape features and builds labeled blob masks. 'labeledEdges'
+ is over-written by 'labeledMask'. Adds features to ROIList structure
+
+ Inputs - labeledEdges, ROIList
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. edge image input is over-written with mask image.
+ ROIList added to.
+
+ """
+
+ # 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):
+ """
+ get_voxel_measures(rawImage, labeledEdges, ROIList)
+
+ takes raw 2D image, labeled blob mask and ROIList. computes voxel features
+ (moments, histogram) for each blob. Adds features to ROIList structure.
+
+ Inputs - rawImage, labeledEdges, ROIList
+ - rawImage is the original source 2D image
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. ROIList added to.
+
+ """
+ #
+ # 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):
+ """
+ get_texture_measures(rawImage, labeledEdges, ROIList)
+
+ takes raw 2D image, labeled blob mask and ROIList. computes 2D
+ texture features using 7x7 Law's texture filters applied
+ to segmented blobs. TEM (texture energy metric) is computed
+ for each Law's filter image and stored in TEM part of ROIList.
+
+ Inputs - rawImage, labeledEdges, ROIList
+ - rawImage is the original source 2D image
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. ROIList added to.
+ """
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the texture (Law's, sub-edges, co-occurence, Gabor) features in the list
+ #
+ S.texture_measures(rawImage, labeledEdges, ROIList)
+ return
+
+def segment_regions():
+ """
+ sourceImage, labeledMask, ROIList = segment_regions()
+
+ Inputs - No Input
+
+ Outputs - sourceImage, labeledMask, ROIList
+ - sourceImage is raw 2D image (default cardiac CT slice for demo
+ - labeledMask is mask of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is numerical Python structure of intensity, shape and
+ texture features for each blob
+
+ High level script calls Python functions:
+ get_slice() - a cardiac CT slice demo file
+ sobel() - sobel magnitude edge finder,
+ returns connected edges
+ get_shape_mask() - gets segmented blob boundary and mask
+ and shape features
+ get_voxel_measures() - uses masks get object voxel moment
+ and histogram features
+ get_texture_measures() - uses masks get object 2D texture features
+ """
+ # get slice from the CT volume
+ image = get_slice(filename)
+ # 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():
+ """
+ regionMask, numberRegions = region_grow()
+ Inputs - No Input
+ Outputs - regionMask, numberRegions
+ - regionMask is the labeled segment masks from 2D image
+ - numberRegions is the number of segmented blobs
+
+ High level script calls Python functions:
+ get_slice() - a cardiac CT slice demo file
+ region_grow() - "grows" connected blobs. default threshold
+ and morphological filter structuring element
+ """
+ # get slice from the CT volume
+ image = get_slice(filename)
+ regionMask, numberRegions = region_grow(image)
+ return regionMask, numberRegions
+
+
+def region_grow(image, lowThreshold=220+2048, highThreshold=600+2048, open=7, close=7):
+ """
+ regionMask, numberRegions = region_grow(image, [defaults])
+
+ Inputs - image, low_threshold, high_threshold, open, close
+ - image is the numarray 2D image
+ - low_ and high_ threshold are density values
+ - open is open morphology structuring element
+ odd size. 0 to turn off. max is 11
+ - close is close morphology structuring element
+ odd size. 0 to turn off. max is 11
+
+ Outputs - regionMask, numberRegions
+ - regionMask is the labeled segment masks from 2D image
+ - numberRegions is the number of segmented blobs
+ """
+ # 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='slice112.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='slice112.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)
+
+
class TestSegment(TestCase):
- def test1(self):
- image = get_slice(filename)
- sourceImage = image.copy()
- edges, objects = sobel(image)
- get_shape_mask(edges, objects)
- get_voxel_measures(sourceImage, edges, objects)
- get_texture_measures(sourceImage, edges, objects)
-
- def test2(self):
- sourceImage, labeledMask, ROIList = segment_regions()
-
- def test3(self):
- regionMask, numberRegions = grow_regions()
- regionMask.max()
- #save_slice(regionMask, 'regionMask.raw')
-
-
-if __name__ == "__main__":
+ def test1(self):
+ image = get_slice(filename)
+ sourceImage = image.copy()
+ edges, objects = sobel(image)
+ get_shape_mask(edges, objects)
+ get_voxel_measures(sourceImage, edges, objects)
+ get_texture_measures(sourceImage, edges, objects)
+
+ def test2(self):
+ sourceImage, labeledMask, ROIList = segment_regions()
+
+ def test3(self):
+ regionMask, numberRegions = grow_regions()
+ regionMask.max()
+ #save_slice(regionMask, 'regionMask.raw')
+
+
+if __name__ == "__main__":
inittest.main()
More information about the Scipy-svn
mailing list