[Scipy-svn] r3492 - in branches/ndimage_segmenter/scipy: . NDISegmenter
scipy-svn@scip...
scipy-svn@scip...
Fri Nov 2 15:09:21 CDT 2007
Author: tom.waite
Date: 2007-11-02 15:09:04 -0500 (Fri, 02 Nov 2007)
New Revision: 3492
Added:
branches/ndimage_segmenter/scipy/NDISegmenter/
branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter.py
branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_EXT.c
branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_IMPL.c
branches/ndimage_segmenter/scipy/NDISegmenter/__multiarray_api.h
branches/ndimage_segmenter/scipy/NDISegmenter/__ufunc_api.h
branches/ndimage_segmenter/scipy/NDISegmenter/arrayobject.h
branches/ndimage_segmenter/scipy/NDISegmenter/arrayscalars.h
branches/ndimage_segmenter/scipy/NDISegmenter/config.h
branches/ndimage_segmenter/scipy/NDISegmenter/control.py
branches/ndimage_segmenter/scipy/NDISegmenter/ndImage_Segmenter_structs.h
branches/ndimage_segmenter/scipy/NDISegmenter/ndarrayobject.h
branches/ndimage_segmenter/scipy/NDISegmenter/noprefix.h
branches/ndimage_segmenter/scipy/NDISegmenter/npy_interrupt.h
branches/ndimage_segmenter/scipy/NDISegmenter/old_defines.h
branches/ndimage_segmenter/scipy/NDISegmenter/oldnumeric.h
branches/ndimage_segmenter/scipy/NDISegmenter/setup.py
branches/ndimage_segmenter/scipy/NDISegmenter/slice112.raw
branches/ndimage_segmenter/scipy/NDISegmenter/testShenCastan.py
branches/ndimage_segmenter/scipy/NDISegmenter/ufuncobject.h
branches/ndimage_segmenter/scipy/NDISegmenter/volumeInput.py
Log:
Add NDISegmenter code.
Added: branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter.py
===================================================================
--- branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter.py 2007-11-02 19:51:44 UTC (rev 3491)
+++ branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter.py 2007-11-02 20:09:04 UTC (rev 3492)
@@ -0,0 +1,129 @@
+import numpy as N
+import NDI_Segmenter as S
+import volumeInput as V
+import struct
+
+def ShenCastan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048, highThreshold=600+2048, dust=16):
+ datatype = [('L', 'i', 1), ('R', 'i', 1), ('T', 'i', 1), ('B', 'i', 1),
+ ('Label', 'i', 1), ('Area', 'i', 1), ('cX', 'f', 1), ('cY', 'f', 1),
+ ('curveClose', 'i', 1), ('cXB', 'f', 1), ('cYB', 'f', 1), ('bLength', 'f', 1),
+ ('minRadius', 'f', 1), ('maxRadius', 'f', 1), ('aveRadius', 'f', 1), ('ratio', 'f', 1),
+ ('compactness', 'f', 1), ('voxelMean', 'f', 1), ('voxelVar', 'f', 1), ('TEM', 'f', 20)]
+ labeledEdges, numberObjects = S.ShenCastanEdges(scLow, IIRFilter, window, lowThreshold, highThreshold, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=datatype)
+ # return the bounding box for each connected edge
+ S.SetObjectStats(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def Sobel(image, sLow=0.3, tMode=1, lowThreshold=220+2048, highThreshold=600+2048, BPHigh=10.0, apearture=21, dust=16):
+ datatype = [('L', 'i', 1), ('R', 'i', 1), ('T', 'i', 1), ('B', 'i', 1),
+ ('Label', 'i', 1), ('Area', 'i', 1), ('cX', 'f', 1), ('cY', 'f', 1),
+ ('curveClose', 'i', 1), ('cXB', 'f', 1), ('cYB', 'f', 1), ('bLength', 'f', 1),
+ ('minRadius', 'f', 1), ('maxRadius', 'f', 1), ('aveRadius', 'f', 1), ('ratio', 'f', 1),
+ ('compactness', 'f', 1), ('voxelMean', 'f', 1), ('voxelVar', 'f', 1), ('TEM', 'f', 20)]
+ # get sobel edge points. return edges that are labeled (1..numberObjects)
+ labeledEdges, numberObjects = S.SobelEdges(sLow, tMode, lowThreshold, highThreshold, BPHigh, apearture, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=datatype)
+ # return the bounding box for each connected edge
+ S.GetObjectStats(labeledEdges, ROIList)
+ # thin (medial axis transform) of the sobel edges as the sobel produces a 'band edge'
+ S.MorphoThinFilt(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def Canny(image, cSigma=1.0, cLow=0.5, cHigh=0.8, tMode=1, lowThreshold=220+2048, highThreshold=600+2048,
+ BPHigh=10.0, apearture=21, dust=16):
+ datatype = [('L', 'i', 1), ('R', 'i', 1), ('T', 'i', 1), ('B', 'i', 1),
+ ('Label', 'i', 1), ('Area', 'i', 1), ('cX', 'f', 1), ('cY', 'f', 1),
+ ('curveClose', 'i', 1), ('cXB', 'f', 1), ('cYB', 'f', 1), ('bLength', 'f', 1),
+ ('minRadius', 'f', 1), ('maxRadius', 'f', 1), ('aveRadius', 'f', 1), ('ratio', 'f', 1),
+ ('compactness', 'f', 1), ('voxelMean', 'f', 1), ('voxelVar', 'f', 1), ('TEM', 'f', 20)]
+ # get canny edge points. return edges that are labeled (1..numberObjects)
+ labeledEdges, numberObjects = S.CannyEdges(cSigma, cLow, cHigh, tMode, lowThreshold, highThreshold,
+ BPHigh, apearture, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=datatype)
+ # return the bounding box for each connected edge
+ S.GetObjectStats(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def GetShapeMask(labeledEdges, ROIList):
+ # pass in Sobel morph-thinned labeled edge image (LEI) and ROIList
+ # GetShapeMask will augment the ROI list
+ # labeledEdges is the original edge image and overwritten as mask image
+ # maskImage is the mask that is used for blob texture / pixel features
+ S.BuildBoundary(labeledEdges, ROIList)
+ return
+
+def GetVoxelMeasures(rawImage, labeledEdges, ROIList):
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the voxel features in the list
+ #
+ S.VoxelMeasures(rawImage, labeledEdges, ROIList)
+ return
+
+def GetTextureMeasures(rawImage, labeledEdges, ROIList):
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the texture (Law's, co-occurence, Gabor) features in the list
+ #
+ S.TextureMeasures(rawImage, labeledEdges, ROIList)
+ return
+
+def SegmentRegions(volSlice=112):
+ # get slice from the CT volume
+ image = GetSliceFromVolume(volSlice)
+ # 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
+ GetShapeMask(labeledMask, ROIList)
+ # Use the labeled mask and source image (raw) to get voxel features
+ GetVoxelMeasures(sourceImage, labeledMask, ROIList)
+ # Use the labeled mask and source image (raw) to get texture features
+ GetTextureMeasures(sourceImage, labeledMask, ROIList)
+ return sourceImage, labeledMask, ROIList
+
+def GrowRegions(volSlice=112):
+ # get slice from the CT volume
+ image = GetSliceFromVolume(volSlice)
+ regionMask, numberRegions = RegionGrow(image)
+ return regionMask, numberRegions
+
+
+def RegionGrow(image, lowThreshold=220+2048, highThreshold=600+2048, open=7, close=7):
+ # morphology filters need to be clipped to 11 max and be odd
+ regionMask, numberRegions = S.RegionGrow(lowThreshold, highThreshold, close, open, image)
+ return regionMask, numberRegions
+
+
+def GetSlice(imageName='junk.raw', bytes=2, rows=512, columns=512):
+ # get a slice alrady extracted from the CT volume
+ image = open(imageName, 'rb')
+ slice = image.read(rows*columns*bytes)
+ values = struct.unpack('h'*rows*columns, slice)
+ ImageSlice = N.array(values, dtype=float).reshape(rows, columns)
+ return (ImageSlice)
+
+def GetSliceFromVolume(mySlice, rows=512, columns=512, bytes=2):
+ # extract a slice the CT volume. Hardwirred
+ image = open('C:\PythonStuff\CardiacCT.vol', 'rb')
+ image.seek(mySlice*rows*columns*bytes)
+ slice = image.read(rows*columns*bytes)
+ values = struct.unpack('h'*rows*columns, slice)
+ ImageSlice = N.array(values, dtype=float).reshape(rows, columns)
+ return (ImageSlice+2048)
+
+def SaveSlice(mySlice, filename='junk.raw', rows=512, columns=512, bytes=2):
+ # just save the slice to a fixed file
+ slice = mySlice.astype(int)
+ image = open(filename, 'wb')
+ image.write(slice)
+ image.close()
+ return
+
+
Property changes on: branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter.py
___________________________________________________________________
Name: svn:executable
+ *
Added: branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_EXT.c
===================================================================
--- branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_EXT.c 2007-11-02 19:51:44 UTC (rev 3491)
+++ branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_EXT.c 2007-11-02 20:09:04 UTC (rev 3492)
@@ -0,0 +1,445 @@
+#include "ndImage_Segmenter_structs.h"
+#include "Python.h"
+#include "noprefix.h"
+
+static PyObject *NDI_Segmenter_CannyEdges(PyObject *self, PyObject *args)
+{
+
+ double sigma;
+ double cannyLow;
+ double cannyHigh;
+ double BPHigh;
+ int lowThreshold;
+ int highThreshold;
+ int apearture;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int mode;
+ int groups;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_Parse(args, "(dddiiidiO)", &sigma, &cannyLow, &cannyHigh, &mode, &lowThreshold, &highThreshold,
+ &BPHigh, &apearture, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ itype = 4;
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+ if(!NI_CannyEdges(num, (int)dims[0], (int)dims[1], sigma, cannyLow, cannyHigh, mode, lowThreshold,
+ highThreshold, BPHigh, apearture, fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups);
+
+}
+
+static PyObject *NDI_Segmenter_SobelEdges(PyObject *self, PyObject *args)
+{
+
+ double sobelLow;
+ double BPHigh;
+ int lowThreshold;
+ int highThreshold;
+ int apearture;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int groups;
+ int mode;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_Parse(args, "(diiidiO)", &sobelLow, &mode, &lowThreshold, &highThreshold, &BPHigh, &apearture, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // this is int type and hard-wirred. pass this in from Python code
+ itype = 4; // unsigned short
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+
+ if(!NI_SobelEdges(num, (int)dims[0], (int)dims[1], sobelLow, mode, lowThreshold, highThreshold, BPHigh, apearture,
+ fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups-1);
+
+}
+
+
+
+static PyObject *NDI_Segmenter_ShenCastanEdges(PyObject *self, PyObject *args)
+{
+ int window;
+ int lowThreshold;
+ int highThreshold;
+ double ShenCastanLow;
+ double b;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ int groups;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ if(!PyArg_Parse(args, "(ddiiiO)", &ShenCastanLow, &b, &window, &lowThreshold, &highThreshold, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // this is int type and hard-wirred. pass this in from Python code
+ itype = 4; // unsigned short
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+ if(!NI_ShenCastanEdges(num, (int)dims[0], (int)dims[1], b, ShenCastanLow, window, lowThreshold, highThreshold,
+ fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups-1);
+
+}
+
+static PyObject *NDI_Segmenter_GetObjectStats(PyObject *self, PyObject *args)
+{
+
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ unsigned short *fP1;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ objStruct *myData;
+
+ if(!PyArg_Parse(args, "(OO)", &iArray, &nArray))
+ goto exit;
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(nArray))
+ goto exit;
+
+ //
+ // PyArray_ContiguousFromObject or PyArray_ContiguousFromAny to be explored
+ // for non-contiguous
+ //
+
+
+ // pointer to the edge-labeled image
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+ fP1 = (unsigned short *)PyArray_DATA(iArray);
+
+ // the object descriptor array that was allocated from numpy
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ myData = (objStruct*)PyArray_DATA(nArray);
+
+ if(!NI_GetObjectStats((int)dims[0], (int)dims[1], (int)objNumber[0], fP1, myData))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_MorphoThinFilt(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ unsigned short *fP1;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OO)", &iArray, &nArray))
+ goto exit;
+
+ fP1 = (unsigned short *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ if(!NI_ThinFilter(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_BuildBoundary(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ unsigned short *fP1;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OO)", &iArray, &nArray))
+ goto exit;
+
+ fP1 = (unsigned short *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+ //
+ // this is int type and hard-wirred. pass this in from Python code
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ //
+ // pass in ROI list and labeled edges
+ // return an augmented ROI list
+ // replace the edgeImage with maskImage
+ //
+ if(!NI_BuildBoundary(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+
+static PyObject *NDI_Segmenter_VoxelMeasures(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ PyObject *eArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OOO)", &iArray, &eArray, &nArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // eArray and iArray are same dims
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ //
+ // pass in ROI list and labeled edges
+ // return an augmented ROI list
+ // replace the edgeImage with maskImage
+ //
+
+ if(!NI_VoxelMeasures(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, fP2, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_TextureMeasures(PyObject *self, PyObject *args)
+{
+
+ int num;
+ int nd;
+ int type;
+ npy_intp *dims;
+ npy_intp *objNumber;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *nArray = NULL;
+ PyObject *eArray = NULL;
+ objStruct *ROIList;
+
+ if(!PyArg_Parse(args, "(OOO)", &iArray, &eArray, &nArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // eArray and iArray are same dims
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ objNumber = PyArray_DIMS(nArray); // this is the number of labels in the edge image
+ ROIList = (objStruct*)PyArray_DATA(nArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray))
+ goto exit;
+
+ //
+ // pass in ROI list and labeled edges
+ // return an augmented ROI list
+ // replace the edgeImage with maskImage
+ //
+
+ if(!NI_TextureMeasures(num, (int)dims[0], (int)dims[1], (int)objNumber[0], fP1, fP2, ROIList))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("");
+
+}
+
+static PyObject *NDI_Segmenter_RegionGrow(PyObject *self, PyObject *args)
+{
+
+ int lowThreshold;
+ int highThreshold;
+ int closeWindow;
+ int openWindow;
+ int num;
+ int nd;
+ int type;
+ int itype;
+ int groups;
+ int mode;
+ npy_intp *dims;
+ double *fP1;
+ unsigned short *fP2;
+ PyObject *iArray = NULL;
+ PyObject *eArray = NULL;
+
+ //
+ // pass in 2D LPF coefficients
+ if(!PyArg_Parse(args, "(iiiiO)", &lowThreshold, &highThreshold, &closeWindow, &openWindow, &iArray))
+ goto exit;
+
+ fP1 = (double *)PyArray_DATA(iArray);
+ nd = PyArray_NDIM(iArray);
+ dims = PyArray_DIMS(iArray);
+ type = PyArray_TYPE(iArray);
+ num = PyArray_SIZE(iArray);
+
+ // this is int type and hard-wirred. pass this in from Python code
+ itype = 4; // unsigned short
+ eArray = (PyObject*)PyArray_SimpleNew(nd, dims, itype);
+ fP2 = (unsigned short *)PyArray_DATA(eArray);
+
+ if(!PyArray_ISCONTIGUOUS(iArray) || !PyArray_ISCONTIGUOUS(eArray))
+ goto exit;
+
+
+ if(!NI_RegionGrow(num, (int)dims[0], (int)dims[1], lowThreshold, highThreshold, closeWindow, openWindow,
+ fP1, fP2, &groups))
+ goto exit;
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("Oi", eArray, groups-1);
+
+}
+
+static PyMethodDef NDI_SegmenterMethods[] =
+{
+ { "CannyEdges", NDI_Segmenter_CannyEdges, METH_VARARGS },
+ { "ShenCastanEdges", NDI_Segmenter_ShenCastanEdges, METH_VARARGS },
+ { "SobelEdges", NDI_Segmenter_SobelEdges, METH_VARARGS },
+ { "GetObjectStats", NDI_Segmenter_GetObjectStats, METH_VARARGS },
+ { "MorphoThinFilt", NDI_Segmenter_MorphoThinFilt, METH_VARARGS },
+ { "BuildBoundary", NDI_Segmenter_BuildBoundary, METH_VARARGS },
+ { "VoxelMeasures", NDI_Segmenter_VoxelMeasures, METH_VARARGS },
+ { "TextureMeasures", NDI_Segmenter_TextureMeasures, METH_VARARGS },
+ { "RegionGrow", NDI_Segmenter_RegionGrow, METH_VARARGS },
+ { NULL, NULL },
+};
+
+void initNDI_Segmenter()
+{
+ Py_InitModule("NDI_Segmenter", NDI_SegmenterMethods);
+ import_array();
+}
+
Property changes on: branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_EXT.c
___________________________________________________________________
Name: svn:executable
+ *
Added: branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_IMPL.c
===================================================================
--- branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_IMPL.c 2007-11-02 19:51:44 UTC (rev 3491)
+++ branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_IMPL.c 2007-11-02 20:09:04 UTC (rev 3492)
@@ -0,0 +1,2965 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "ndImage_Segmenter_structs.h"
+
+// these are for this standalone and come out with the full build
+//
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+#define FALSE 0
+#define TRUE 1
+
+int NI_GetObjectStats(int rows, int cols, int numberObjects, unsigned short *labeledEdges, objStruct objectMetrics[]){
+
+ int i, j, k, m;
+ int offset;
+ int count;
+ int LowX;
+ int LowY;
+ int HighX;
+ int HighY;
+ int status;
+ float centerX;
+ float centerY;
+
+ for(k = 1; k < numberObjects; ++k){
+ offset = cols;
+ LowX = 32767;
+ LowY = 32767;
+ HighX = 0;
+ HighY = 0;
+ count = 0;
+ centerX = (float)0.0;
+ centerY = (float)0.0;
+ for(i = 1; i < (rows-1); ++i){
+ for(j = 1; j < (cols-1); ++j){
+ m = labeledEdges[offset+j];
+ if(k == m){
+ if(i < LowY) LowY = i;
+ if(j < LowX) LowX = j;
+ if(i > HighY) HighY = i;
+ if(j > HighX) HighX = j;
+ centerX += (float)j;
+ centerY += (float)i;
+ ++count;
+ }
+ }
+ offset += cols;
+ }
+ // the bounding box for the 2D blob
+ objectMetrics[k-1].L = LowX;
+ objectMetrics[k-1].R = HighX;
+ objectMetrics[k-1].B = LowY;
+ objectMetrics[k-1].T = HighY;
+ objectMetrics[k-1].Area = count;
+ objectMetrics[k-1].cX = centerX/(float)count;
+ objectMetrics[k-1].cY = centerY/(float)count;
+ objectMetrics[k-1].Label = k;
+ }
+
+ status = numberObjects;
+ return status;
+
+}
+
+
+void buildKernel(double BPHigh, int HalfFilterTaps, int apearture, float *kernel){
+
+ int i, j;
+ float r, t1, t2, t3, t4;
+ float LC, HC, tLOW, tHIGH;
+ float pi = (float)3.14159, rad = (float)0.01745;
+
+ LC = (float)0.0;
+ HC = BPHigh * rad;
+ t2 = (float)2.0*pi;
+ t1 = (float)2.0*HalfFilterTaps + (float)1.0;
+ //
+ // build the Filter Kernel
+ // the kernel starts at 1 only because it is linked to the internal filter2D routine
+ // the code is not a Fortran code
+ //
+ j = 1;
+ for(i = -HalfFilterTaps; i <= HalfFilterTaps; ++i){
+ r = (float)i;
+ if(r == (float)0.0){
+ tLOW = LC;
+ tHIGH = HC;
+ }
+ else{
+ tLOW = (float)(sin(r*LC))/r;
+ tHIGH = (float)(sin(r*HC))/r;
+ }
+ t3 = (float)0.54 + (float)0.46*((float)cos(r*t2/t1));
+ t4 = t3*(tHIGH-tLOW);
+ kernel[j++] = t4;
+ }
+
+ // normalize the kernel so unity gain (as is LP filter this is easy)
+ t1 = (float)0.0;
+ for(j = 1; j <= apearture; ++j){
+ t1 += kernel[j];
+ }
+ for(j = 1; j <= apearture; ++j){
+ kernel[j] /= t1;
+ }
+
+ t1 = (float)0.0;
+ for(j = 1; j <= apearture; ++j){
+ t1 += kernel[j];
+ }
+ return;
+}
+
+void filter2D(int HalfFilterTaps, int rows, int cols, int lowThreshold, int highThreshold, float *kernel, double *Image){
+
+ int i, j, k, n, num1;
+ int offset;
+ float sum, value;
+ float buffer[1024];
+
+ num1 = HalfFilterTaps + 1;
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ // copy image row to local buffer
+ for(j = 0; j < cols; ++j){
+ buffer[num1+j] = Image[offset+j];
+ }
+ // constant pad the ends of the buffer
+ for(j = 0; j < num1; ++j){
+ buffer[j] = buffer[num1];
+ }
+ for(j = cols+num1; j < cols+2*num1; ++j){
+ buffer[j] = buffer[cols-1+num1];
+ }
+
+ // Perform Symmetric Convolution in the X dimension.
+ for(n = 0, j = num1; j < (cols+num1); ++j, ++n){
+ sum = buffer[j] * kernel[num1];
+ for(k = 1; k < num1; ++k){
+ sum += kernel[num1-k] * (buffer[j+k] + buffer[j-k]);
+ }
+ Image[offset+n] = sum;
+ }
+ offset += cols;
+ }
+
+ offset = 0;
+ for(i = 0; i < cols; ++i){
+ // copy image column to local buffer
+ offset = 0;
+ for(j = 0; j < rows; ++j){
+ buffer[num1+j] = Image[offset+i];
+ offset += cols;
+ }
+ // constant pad the ends of the buffer
+ for(j = 0; j < num1; ++j){
+ buffer[j] = buffer[num1];
+ }
+ for(j = rows+num1; j < rows+2*num1; ++j){
+ buffer[j] = buffer[rows-1+num1];
+ }
+
+ // Perform Symmetric Convolution in the Y dimension.
+ offset = 0;
+ for(j = num1; j < (rows+num1); ++j){
+ sum = buffer[j] * kernel[num1];
+ for(k = 1; k < num1; ++k){
+ sum += kernel[num1-k] * (buffer[j+k] + buffer[j-k]);
+ }
+ Image[offset+i] = sum;
+ offset += cols;
+ }
+ }
+
+ // threshold the image
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = Image[offset+j];
+ if(value < (float)lowThreshold) value = (float)0.0;
+ if(value > (float)highThreshold) value = (float)0.0;
+ Image[offset+j] = value;
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+void doPreProcess(int samples, int rows, int cols, double *rawImage, double BPHigh, int apearture, int lowThreshold, int highThreshold){
+
+ //
+ // 2D low pass filter using bisinc and threshold
+ // this specific example is on cardiac CT and focuses on segmenting the
+ // aorta and blood-filled chambers. for MRI the threshold will be different
+ //
+
+ float *kernel;
+ int HalfFilterTaps = (apearture-1)/2;
+ kernel = calloc(apearture+16, sizeof(float));
+
+ buildKernel(BPHigh, HalfFilterTaps, apearture, kernel);
+ filter2D(HalfFilterTaps, rows, cols, lowThreshold, highThreshold, kernel, rawImage);
+
+ free(kernel);
+
+ return;
+
+}
+
+
+int ConnectedEdgePoints(int rows, int cols, unsigned short *connectedEdges){
+
+ int i, j, k, l, m;
+ int offset;
+ int Label;
+ int Classes[4096];
+ bool NewLabel;
+ bool Change;
+ unsigned short T[12];
+
+ //
+ // connected components labeling. pixels touch within 3x3 mask for edge connectedness.
+ //
+ Label = 1;
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(connectedEdges[offset+j] == 1){
+ connectedEdges[offset+j] = Label++;
+ }
+ }
+ offset += cols;
+ }
+
+ while(1){
+ Change = FALSE;
+ //
+ // TOP-DOWN Pass for labeling
+ //
+ offset = cols;
+ for(i = 1; i < rows-1; ++i){
+ for(j = 1; j < cols-1; ++j){
+ if(connectedEdges[offset+j] != 0){
+ T[0] = connectedEdges[offset+j];
+ T[1] = connectedEdges[offset+j+1];
+ T[2] = connectedEdges[offset-cols+j+1];
+ T[3] = connectedEdges[offset-cols+j];
+ T[4] = connectedEdges[offset-cols+j-1];
+ T[5] = connectedEdges[offset+j-1];
+ T[6] = connectedEdges[offset+cols+j-1];
+ T[7] = connectedEdges[offset+cols+j];
+ T[8] = connectedEdges[offset+cols+j+1];
+ m = T[0];
+ for(l = 1; l < 9; ++l){
+ if(T[l] != 0){
+ if(T[l] < m) m = T[l];
+ }
+ }
+ if(m != connectedEdges[offset+j]){
+ Change = TRUE;
+ connectedEdges[offset+j] = m;
+ }
+ }
+ }
+ offset += cols;
+ }
+ //
+ // BOTTOM-UP Pass for labeling
+ //
+ offset = (rows-1)*cols;
+ for(i = (rows-1); i > 1; --i){
+ for(j = (cols-1); j > 1; --j){
+ if(connectedEdges[offset+j] != 0){
+ T[0] = connectedEdges[offset+j];
+ T[1] = connectedEdges[offset+j+1];
+ T[2] = connectedEdges[offset-cols+j+1];
+ T[3] = connectedEdges[offset-cols+j];
+ T[4] = connectedEdges[offset-cols+j-1];
+ T[5] = connectedEdges[offset+j-1];
+ T[6] = connectedEdges[offset+cols+j-1];
+ T[7] = connectedEdges[offset+cols+j];
+ T[8] = connectedEdges[offset+cols+j+1];
+ m = T[0];
+ for(l = 1; l < 9; ++l){
+ if(T[l] != 0){
+ if(T[l] < m) m = T[l];
+ }
+ }
+ if(m != connectedEdges[offset+j]){
+ Change = TRUE;
+ connectedEdges[offset+j] = m;
+ }
+ }
+ }
+ offset -= cols;
+ }
+ if(!Change) break;
+ } // end while loop
+
+ Classes[0] = 0;
+ Label = 1;
+ offset = cols;
+ for(i = 1; i < (rows-1); ++i){
+ for(j = 1; j < (cols-1); ++j){
+ m = connectedEdges[offset+j];
+ if(m > 0){
+ NewLabel = TRUE;
+ for(k = 1; k < Label; ++k){
+ if(Classes[k] == m) NewLabel = FALSE;
+ }
+ if(NewLabel){
+ Classes[Label++] = m;
+ if(Label > 4000){
+ return 0; // too many labeled regions. this is a pathology in the image slice
+ }
+ }
+ }
+ }
+ offset += cols;
+ }
+
+ //
+ // re-label the connected blobs in continuous label order
+ //
+ offset = cols;
+ for(i = 1; i < (rows-1); ++i){
+ for(j = 1; j < (cols-1); ++j){
+ m = connectedEdges[offset+j];
+ if(m > 0){
+ for(k = 1; k < Label; ++k){
+ if(Classes[k] == m){
+ connectedEdges[offset+j] = (unsigned short)k;
+ break;
+ }
+ }
+ }
+ }
+ offset += cols;
+ }
+
+ return Label;
+}
+
+float magnitude(float X, float Y){
+
+ return (float)sqrt(X*X + Y*Y);
+}
+
+int traceEdge(int i, int j, int rows, int cols, double cannyLow, float *magImage, float *HYSImage){
+
+ int n, m;
+ int ptr;
+ int flag;
+
+ ptr = i * cols;
+ if(HYSImage[ptr+j] == (float)0.0){
+ //
+ // this point is above high threshold
+ //
+ HYSImage[ptr+j] = (float)1.0;
+ flag = 0;
+ for(n = -1; n <= 1; ++n){
+ for(m = -1; m <= 1; ++m){
+ if(n == 0 && m == 0) continue;
+ if(((i+n) > 0) && ((j+m) > 0) && ((i+n) < rows) && ((j+m) < cols)){
+ ptr = (i+n) * cols;
+ if(magImage[ptr+j+m] > cannyLow){
+ //
+ // this point is above low threshold
+ //
+ if(traceEdge(i+n, j+m, rows, cols, cannyLow, magImage, HYSImage)){
+ flag = 1;
+ break;
+ }
+ }
+ }
+ }
+ if(flag) break;
+ }
+ return(1);
+ }
+
+ return(0);
+
+}
+
+
+void edgeThreshold(int rows, int cols, double cannyLow, float *magImage, float *HYSImage){
+
+ int i, j;
+ int ptr;
+
+ for(i = 0; i < rows; ++i){
+ ptr = i * cols;
+ for(j = 0; j < cols; ++j){
+ if(magImage[ptr+j] > cannyLow){
+ HYSImage[ptr+j] = (float)1.0;
+ }
+ }
+ }
+
+ return;
+
+}
+
+void edgeHysteresis(int rows, int cols, double cannyLow, double cannyHigh, float *magImage, float *HYSImage){
+
+ int i, j;
+ int ptr;
+
+ for(i = 0; i < rows; ++i){
+ ptr = i * cols;
+ for(j = 0; j < cols; ++j){
+ if(magImage[ptr+j] > cannyHigh){
+ traceEdge(i, j, rows, cols, cannyLow, magImage, HYSImage);
+ }
+ }
+ }
+
+ return;
+
+}
+
+void nonMaxSupress(int rows, int cols, float aveXValue, float aveYValue, double *cannyLow, double *cannyHigh,
+ int mode, float *hDGImage, float *vDGImage, float *magImage){
+
+ int i, j;
+ int ptr, ptr_m1, ptr_p1;
+ float xSlope, ySlope, G1, G2, G3, G4, G, xC, yC;
+ float scale;
+ float maxValue = (float)0.0;
+ float minValue = (float)-1.0;
+ int histogram[256];
+ int value;
+ int mValue;
+ int mIndex;
+ int count;
+ double step;
+ double tAve;
+
+ for(i = 1; i < rows-1; ++i){
+ ptr = i * cols;
+ ptr_m1 = ptr - cols;
+ ptr_p1 = ptr + cols;
+ for(j = 1; j < cols; ++j){
+ magImage[ptr+j] = (float)0.0;
+ xC = hDGImage[ptr+j];
+ yC = vDGImage[ptr+j];
+ if((fabs(xC) < aveXValue) && (fabs(yC) < aveYValue)) continue;
+ G = magnitude(xC, yC);
+ if(fabs(yC) > fabs(xC)){
+ // vertical gradient
+ xSlope = (float)(fabs(xC) / fabs(yC));
+ ySlope = (float)1.0;
+ G2 = magnitude(hDGImage[ptr_m1+j], vDGImage[ptr_m1+j]);
+ G4 = magnitude(hDGImage[ptr_p1+j], vDGImage[ptr_p1+j]);
+ if((xC*yC) > (float)0.0){
+ G1 = magnitude(hDGImage[ptr_m1+j-1], vDGImage[ptr_m1+j-1]);
+ G3 = magnitude(hDGImage[ptr_p1+j+1], vDGImage[ptr_p1+j+1]);
+ }
+ else{
+ G1 = magnitude(hDGImage[ptr_m1+j+1], vDGImage[ptr_m1+j+1]);
+ G3 = magnitude(hDGImage[ptr_p1+j-1], vDGImage[ptr_p1+j-1]);
+ }
+ }
+ else{
+ // horizontal gradient
+ xSlope = (float)(fabs(yC) / fabs(xC));
+ ySlope = (float)1.0;
+ G2 = magnitude(hDGImage[ptr+j+1], vDGImage[ptr+j+1]);
+ G4 = magnitude(hDGImage[ptr+j-1], vDGImage[ptr+j-1]);
+ if((xC*yC) > (float)0.0){
+ G1 = magnitude(hDGImage[ptr_p1+j+1], vDGImage[ptr_p1+j+1]);
+ G3 = magnitude(hDGImage[ptr_m1+j-1], vDGImage[ptr_m1+j-1]);
+ }
+ else{
+ G1 = magnitude(hDGImage[ptr_m1+j+1], vDGImage[ptr_m1+j+1]);
+ G3 = magnitude(hDGImage[ptr_p1+j-1], vDGImage[ptr_p1+j-1]);
+ }
+ }
+ if( (G > (xSlope*G1 + (ySlope-xSlope)*G2)) && (G > (xSlope*G3 + (ySlope-xSlope)*G4)) ){
+ magImage[ptr+j] = G;
+ }
+ if(magImage[ptr+j] > maxValue) maxValue = magImage[ptr+j];
+ if(magImage[ptr+j] < minValue) minValue = magImage[ptr+j];
+ }
+ }
+
+ scale = (float)1.0 / (maxValue-minValue);
+ ptr = 0;
+ count = 0;
+ tAve = 0.0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ magImage[ptr] = scale * (magImage[ptr]-minValue);
+ if(magImage[ptr] > 0.0){
+ tAve += magImage[ptr];
+ ++count;
+ }
+ ++ptr;
+ }
+ }
+ tAve /= (float)count;
+
+ step = 255.0;
+ for(i = 0; i < 256; ++i){
+ histogram[i] = 0;
+ }
+ ptr = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = (int)(step*(magImage[ptr]));
+ ++histogram[value];
+ ++ptr;
+ }
+ }
+ //
+ // now get the max after skipping the low values
+ //
+ mValue = -1;
+ mIndex = 0;
+ for(i = 10; i < 256; ++i){
+ if(histogram[i] > mValue){
+ mValue = histogram[i];
+ mIndex = i;
+ }
+ }
+
+ if(mode == 1){
+ // based on the mean value of edge energy
+ *cannyLow = ((*cannyLow) * tAve);
+ *cannyHigh = ((*cannyHigh) * tAve);
+ }
+ else{
+ // based on the mode value of edge energy
+ *cannyLow = ((*cannyLow) * ((float)mIndex/step));
+ *cannyHigh = ((*cannyHigh) * ((float)mIndex/step));
+ }
+
+ return;
+
+}
+
+void DGFilters(int samples, int rows, int cols, double cannySigma, int gWidth,
+ float *aveXValue, float *aveYValue, double *rawImage,
+ double *dgKernel, float *hDGImage, float *vDGImage){
+
+ //
+ // implements the derivative of Gaussian filter. kernel set by CannyEdges
+ //
+ int i, j, k;
+ int ptr;
+ int mLength;
+ int count;
+ float *tBuffer = NULL;
+ double sum;
+
+ *aveXValue = (float)0.0;
+ *aveYValue = (float)0.0;
+
+ mLength = MAX(rows, cols) + 64;
+ tBuffer = calloc(mLength, sizeof(float));
+
+ //
+ // filter X
+ //
+ count = 0;
+ for(i = 0; i < rows; ++i){
+ ptr = i * cols;
+ for(j = gWidth; j < cols-gWidth; ++j){
+ sum = dgKernel[0] * rawImage[ptr+j];
+ for(k = 1; k < gWidth; ++k){
+ sum += dgKernel[k] * (-rawImage[ptr+j+k] + rawImage[ptr+j-k]);
+ }
+ hDGImage[ptr+j] = (float)sum;
+ if(sum != (float)0.0){
+ ++count;
+ *aveXValue += (float)fabs(sum);
+ }
+ }
+ }
+ if(count){
+ *aveXValue /= (float)count;
+ *aveXValue = (float)0.5 * (*aveXValue);
+ // this is 50% of the max, hardwirred for now, and is part of the threshold
+ }
+ //
+ // filter Y
+ //
+ count = 0;
+ for(i = 0; i < cols; ++i){
+ for(j = 0; j < rows; ++j){
+ ptr = j * cols;
+ tBuffer[j] = rawImage[ptr+i];
+ }
+ for(j = gWidth; j < rows-gWidth; ++j){
+ ptr = j * cols;
+ sum = dgKernel[0] * tBuffer[j];
+ for(k = 1; k < gWidth; ++k){
+ sum += dgKernel[k] * (-tBuffer[j+k] + tBuffer[j-k]);
+ }
+ vDGImage[ptr+i] = sum;
+ if(sum != (float)0.0){
+ ++count;
+ *aveYValue += (float)fabs(sum);
+ }
+ }
+ }
+ if(count){
+ *aveYValue /= (float)count;
+ *aveYValue = (float)0.5 * (*aveYValue);
+ // this is 50% of the max, hardwirred for now, and is part of the threshold
+ }
+
+ free(tBuffer);
+
+ return;
+
+}
+
+
+int NI_CannyEdges(int samples, int rows, int cols, double cannySigma, double cannyLow, double cannyHigh, int mode,
+ int lowThreshold, int highThreshold, double BPHigh, int apearture, double *rawImage,
+ unsigned short *edgeImage, int *groups){
+
+ int i, j;
+ int offset;
+ int doHysteresis = 0;
+ int gWidth;
+ int mLength;
+ int status;
+ float aveXValue;
+ float aveYValue;
+ double t;
+ double dgKernel[20];
+ float *HYSImage = NULL;
+ float *hDGImage = NULL;
+ float *vDGImage = NULL;
+ float *magImage = NULL;
+ float *tBuffer = NULL;
+
+ // filter
+ printf("do preProcess\n");
+ doPreProcess(samples, rows, cols, rawImage, BPHigh, apearture, lowThreshold, highThreshold);
+ printf("do Canny\n");
+
+ //
+ // memory for magnitude, horizontal and vertical derivative of Gaussian filter
+ //
+ mLength = MAX(rows, cols) + 64;
+ HYSImage = calloc(samples, sizeof(float));
+ hDGImage = calloc(samples, sizeof(float));
+ vDGImage = calloc(samples, sizeof(float));
+ magImage = calloc(samples, sizeof(float));
+ tBuffer = calloc(mLength, sizeof(float));
+
+ //
+ // build derivative of Gaussian filter kernel
+ // kernel is anti-symmetric so convolution is k[j]*(v[i+j] - v[i-j])
+ //
+ gWidth = 20;
+ for(i = 0; i < gWidth; ++i){
+ t = (float)i;
+ dgKernel[i] = (float)exp((double)((-t*t)/((float)2.0 * cannySigma * cannySigma)));
+ dgKernel[i] *= -(t / (cannySigma * cannySigma));
+ }
+ for(i = 0; i < samples; ++i){
+ HYSImage[i] = (float)0.0;
+ }
+
+ DGFilters(samples, rows, cols, cannySigma, gWidth, &aveXValue, &aveYValue, rawImage, dgKernel, hDGImage, vDGImage);
+ nonMaxSupress(rows, cols, aveXValue, aveYValue, &cannyLow, &cannyHigh, mode, hDGImage, vDGImage, magImage);
+ if(doHysteresis){
+ edgeHysteresis(rows, cols, cannyLow, cannyHigh, magImage, HYSImage);
+ }
+ else{
+ edgeThreshold(rows, cols, cannyLow, magImage, HYSImage);
+ }
+
+ //
+ // edge image
+ //
+ for(i = 0; i < samples; ++i){
+ edgeImage[i] = (unsigned short)HYSImage[i];
+ }
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+
+ free(tBuffer);
+ free(hDGImage);
+ free(vDGImage);
+ free(magImage);
+ free(HYSImage);
+
+ status = *groups;
+ return status;
+
+}
+
+void doSobel(int samples, int rows, int cols, double sobelLow, int mode, double *rawImage, unsigned short *edgeImage){
+
+ int i, j;
+ int p, m, n;
+ int offset;
+ int offsetM1;
+ int offsetP1;
+ int minValue, maxValue;
+ int pAve = 0;
+ int count = 0;
+ int histogram[256];
+ int value;
+ int maxIndex;
+ float pThreshold;
+ double scale;
+ double step;
+ float *filteredImage = NULL;
+
+ filteredImage = calloc(samples, sizeof(float));
+
+ minValue = 10000;
+ maxValue = -10000;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ filteredImage[offset+j] = 0;
+ edgeImage[offset+j] = 0;
+ }
+ offset += cols;
+ }
+
+ //
+ // Sobel
+ //
+ offset = cols;
+ for(i = 1; i < rows-1; ++i){
+ offsetM1 = offset - cols;
+ offsetP1 = offset + cols;
+ for(j = 1; j < cols-1; ++j){
+ n = 2*rawImage[offsetM1+j] + rawImage[offsetM1+j-1] + rawImage[offsetM1+j+1] -
+ 2*rawImage[offsetP1+j] - rawImage[offsetP1+j-1] - rawImage[offsetP1+j+1];
+ m = 2*rawImage[offset+j-1] + rawImage[offsetM1+j-1] + rawImage[offsetP1+j-1] -
+ 2*rawImage[offset+j+1] - rawImage[offsetM1+j+1] - rawImage[offsetP1+j+1];
+ p = (int)sqrt((float)(m*m) + (float)(n*n));
+ if(p > 0){
+ pAve += p;
+ ++count;
+ if(p > maxValue) maxValue = p;
+ if(p < minValue) minValue = p;
+ }
+ filteredImage[offset+j] = p;
+ }
+ offset += cols;
+ }
+
+ // threshold based on ave
+ pAve /= count;
+ scale = 1.0 / maxValue;
+
+ step = 255.0/(maxValue-minValue);
+ for(i = 0; i < 256; ++i){
+ histogram[i] = 0;
+ }
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = (int)(step*(filteredImage[offset+j]-minValue));
+ ++histogram[value];
+ }
+ offset += cols;
+ }
+ //
+ // now get the max after skipping the low values
+ //
+ maxValue = -1;
+ maxIndex = 0;
+ for(i = 10; i < 256; ++i){
+ if(histogram[i] > maxValue){
+ maxValue = histogram[i];
+ maxIndex = i;
+ }
+ }
+
+ if(mode == 1){
+ // based on the mean value of edge energy
+ pThreshold = (int)(sobelLow * (float)pAve);
+ }
+ else{
+ // based on the mode value of edge energy
+ pThreshold = (sobelLow * (minValue + ((float)maxIndex/step)));
+ }
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(filteredImage[offset+j] > pThreshold){
+ edgeImage[offset+j] = 1;
+ }
+ else{
+ edgeImage[offset+j] = 0;
+ }
+ filteredImage[offset+j] *= scale;
+ }
+ offset += cols;
+ }
+
+ free(filteredImage);
+
+ return;
+
+
+}
+
+void estimateThreshold(float *lowThreshold, float *highThreshold, float ShenCastanLow, int rows, int cols, float *SourceImage){
+
+ int i, j;
+ int offset;
+ int value;
+ int mIndex;
+ int histogram[256];
+ float low, high;
+ float scale;
+
+ low = (float)1000.0;
+ high = (float)-1000.0;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(fabs(SourceImage[offset+j]) > high) high = fabs(SourceImage[offset+j]);
+ if(fabs(SourceImage[offset+j]) < low) low = fabs(SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ scale = (float)255.0 / (high-low);
+ for(i = 0; i < 256; ++i){
+ histogram[i] = 0;
+ }
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = (int)(scale*(fabs(SourceImage[offset+j]) - low));
+ ++histogram[value];
+ }
+ offset += cols;
+ }
+
+ //
+ // now get the edge energy mode
+ //
+ value = 0;
+ mIndex = 10;
+ for(i = 10; i < 256; ++i){
+ if(histogram[i] > value){
+ value = histogram[i];
+ mIndex = i;
+ }
+ }
+
+ *highThreshold = ((float)mIndex / scale) + low;
+ *lowThreshold = ((float)mIndex / scale) + low;
+
+ *highThreshold *= ShenCastanLow;
+ *lowThreshold *= ShenCastanLow;
+
+ return;
+
+}
+
+void thresholdEdges(float *SourceImage, unsigned short *EdgeImage, double ShenCastanLow, int rows, int cols){
+
+ int i, j;
+ int offset;
+ float tLow, tHigh;
+
+ //
+ // SourceImage contains the adaptive gradient
+ // get threshold from the mode of the edge energy
+ //
+ estimateThreshold(&tLow, &tHigh, ShenCastanLow, rows, cols, SourceImage);
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(SourceImage[offset+j] > tLow){
+ EdgeImage[offset+j] = 1;
+ }
+ else{
+ EdgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+float adaptiveGradient(float *BLImage, float *FilterImage, int nrow, int ncol, int cols, int window){
+
+ int i, j;
+ int offset;
+ int numOn, numOff;
+ int hWindow = window/2;
+ float sumOn, sumOff;
+ float aveOn, aveOff;
+
+ numOn = 0;
+ numOff = 0;
+
+ sumOn = (float)0.0;
+ sumOff = (float)0.0;
+
+ aveOn = (float)0.0;
+ aveOff = (float)0.0;
+
+ offset = nrow * cols;
+ for(i = -hWindow; i < hWindow; ++i){
+ for(j = -hWindow; j < hWindow; ++j){
+ if(BLImage[offset+(i*cols)+(j+ncol)] == 1){
+ sumOn += FilterImage[offset+(i*cols)+(j+ncol)];
+ ++numOn;
+ }
+ else{
+ sumOff += FilterImage[offset+(i*cols)+(j+ncol)];
+ ++numOff;
+ }
+ }
+ }
+
+ if(numOn){
+ aveOn = sumOn / numOn;
+ }
+
+ if(numOff){
+ aveOff = sumOff / numOff;
+ }
+
+ return (aveOff-aveOn);
+
+}
+
+void getZeroCrossings(float *SourceImage, float *FilterImage, float *BLImage, int rows, int cols, int window){
+
+ int i, j;
+ int offset;
+ bool validEdge;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ SourceImage[offset+j] = 0.0;
+ }
+ offset += cols;
+ }
+
+ offset = window*cols;
+ for(i = window; i < rows-window; ++i){
+ for(j = window; j < cols-window; ++j){
+ validEdge = FALSE;
+ if((BLImage[offset+j] == 1) && (BLImage[offset+cols+j] == 0)){
+ if((FilterImage[offset+cols+j] - FilterImage[offset-cols+j]) > 0.0){
+ validEdge = TRUE;
+ }
+ }
+ else if((BLImage[offset+j] == 1) && (BLImage[offset+j+1] == 0)){
+ if((FilterImage[offset+j+1] - FilterImage[offset+j-1]) > 0.0){
+ validEdge = TRUE;
+ }
+ }
+ else if((BLImage[offset+j] == 1) && (BLImage[offset-cols+j] == 0)){
+ if((FilterImage[offset+cols+j] - FilterImage[offset-cols+j]) < 0.0){
+ validEdge = TRUE;
+ }
+ }
+ else if((BLImage[offset+j] == 1) && (BLImage[offset+j-1] == 0)){
+ if((FilterImage[offset+j+1] - FilterImage[offset+j-1]) < 0.0){
+ validEdge = TRUE;
+ }
+ }
+ if(validEdge){
+ // adaptive gradeint is signed
+ SourceImage[offset+j] = (float)fabs(adaptiveGradient(BLImage, FilterImage, i, j, cols, window));
+ }
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+
+void computeBandedLaplacian(float *image1, float *image2, float *BLImage, int rows, int cols){
+
+ int i, j;
+ int offset;
+ float t;
+
+ //
+ // like an unsharp mask
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ t = image1[offset+j] - image2[offset+j];
+ if(t < (float)0.0){
+ t = (float)0.0;
+ }
+ else{
+ t = (float)1.0;
+ }
+ BLImage[offset+j] = t;
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+void thresholdImage(float *Raw, float *Filtered, int rows, int cols, int tLow, int tHigh){
+
+ int i, j;
+ int ptr;
+
+ ptr = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(Raw[ptr] > tHigh){
+ Raw[ptr] = 0.0;
+ Filtered[ptr] = 0.0;
+ }
+ if(Raw[ptr] < tLow){
+ Raw[ptr] = 0.0;
+ Filtered[ptr] = 0.0;
+ }
+ ++ptr;
+ }
+ }
+
+ return;
+
+}
+
+void ISEF_Vertical(float *SourceImage, float *FilterImage, float *A, float *B, int rows, int cols, double b){
+
+
+ int i, j;
+ int offset;
+ float b1, b2;
+
+ b1 = ((float)1.0 - b)/((float)1.0 + b);
+ b2 = b * b1;
+
+ //
+ // set the boundaries
+ //
+ offset = (rows-1)*cols;
+ for(i = 0; i < cols; ++i){
+ // process row 0
+ A[i] = b1 * SourceImage[i];
+ // process row N-1
+ B[offset+i] = b2 * SourceImage[offset+i];
+ }
+
+ //
+ // causal component of IIR filter
+ //
+ offset = cols;
+ for(i = 1; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ //
+ // IIR ISEF filter applied across rows
+ //
+ A[offset+j] = (b * A[offset-cols+j]) + (b1 * SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ //
+ // anti-causal component of IIR filter
+ //
+ offset = (rows-2)*cols;
+ for(i = rows-2; i >= 0; --i){
+ for(j = 0; j < cols; ++j){
+ //
+ // IIR ISEF filter applied across rows
+ //
+ B[offset+j] = (b * B[offset+cols+j]) + (b2 * SourceImage[offset+j]);
+ }
+ offset -= cols;
+ }
+
+ offset = (rows-1)*cols;
+ for(j = 0; j < cols-1; ++j){
+ FilterImage[offset+j] = A[offset+j];
+ }
+
+ //
+ // add causal and anti-causal IIR parts
+ //
+ offset = 0;
+ for(i = 1; i < rows-2; ++i){
+ for(j = 0; j < cols-1; ++j){
+ FilterImage[offset+j] = A[offset+j] + B[offset+cols+j];
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+void ISEF_Horizontal(float *SourceImage, float *FilterImage, float *A, float *B, int rows, int cols, double b){
+
+
+ //
+ // source and smooth are the same in this pass of the 2D IIR
+ //
+
+ int i, j;
+ int offset;
+ float b1, b2;
+
+ b1 = ((float)1.0 - b)/((float)1.0 + b);
+ b2 = b * b1;
+
+ //
+ // columns boundaries
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ // col 0
+ A[offset] = b1 * SourceImage[offset];
+ // col N-1
+ B[offset+cols-1] = b2 * SourceImage[offset+cols-1];
+ }
+
+ //
+ // causal IIR part
+ //
+ offset = 0;
+ for(j = 1; j < cols; ++j){
+ for(i = 0; i < rows; ++i){
+ A[offset+j] = (b * A[offset+j-1]) + (b1 * SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ //
+ // anti-causal IIR part
+ //
+ offset = 0;
+ for(j = cols-2; j > 0; --j){
+ for(i = 0; i < rows; ++i){
+ B[offset+j] = (b * B[offset+j+1]) + (b2 * SourceImage[offset+j]);
+ }
+ offset += cols;
+ }
+
+ //
+ // filtered output. this is 2-pass IIR and pass 1 is vertical
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ FilterImage[offset+cols-1] = A[offset+cols-1];
+ }
+
+ //
+ // add causal and anti-causal IIR parts
+ //
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols-1; ++j){
+ FilterImage[offset+j] = A[offset+j] + B[offset+j+1];
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+
+void computeISEF(float *SourceImage, float *FilterImage, int rows, int cols, double b){
+
+ int imageSize = rows*cols;
+ float *A;
+ float *B;
+
+ A = calloc(imageSize, sizeof(float));
+ B = calloc(imageSize, sizeof(float));
+
+ ISEF_Vertical(SourceImage, FilterImage, A, B, rows, cols, b);
+ ISEF_Horizontal(FilterImage, FilterImage, A, B, rows, cols, b);
+
+ free(A);
+ free(B);
+
+ return;
+
+}
+
+void Shen_Castan(double b, double ShenCastanLow, int rows, int cols, int window, int lowThreshold, int highThreshold,
+ double *RawImage, unsigned short *EdgeImage){
+
+ int i;
+ int imageSize = rows*cols;
+ float *FilterImage;
+ float *BinaryLaplacianImage;
+ float *SourceImage;
+
+ FilterImage = calloc(imageSize, sizeof(float));
+ BinaryLaplacianImage = calloc(imageSize, sizeof(float));
+ SourceImage = calloc(imageSize, sizeof(float));
+
+ for(i = 0; i < imageSize; ++i){
+ SourceImage[i] = RawImage[i];
+ }
+ computeISEF(SourceImage, FilterImage, rows, cols, b);
+ // optional thresholding based on low, high
+ thresholdImage(SourceImage, FilterImage, rows, cols, lowThreshold, highThreshold);
+ computeBandedLaplacian(FilterImage, SourceImage, BinaryLaplacianImage, rows, cols);
+ // the new source image is now the adaptive gradient
+ getZeroCrossings(SourceImage, FilterImage, BinaryLaplacianImage, rows, cols, window);
+ thresholdEdges(SourceImage, EdgeImage, ShenCastanLow, rows, cols);
+
+ free(FilterImage);
+ free(BinaryLaplacianImage);
+ free(SourceImage);
+
+ return;
+
+}
+
+int NI_ShenCastanEdges(int samples, int rows, int cols, double b, double ShenCastanLow, int window, int lowThreshold,
+ int highThreshold, double *rawImage, unsigned short *edgeImage, int *groups){
+
+
+ int i, j;
+ int offset;
+ int status = 0;
+
+ Shen_Castan(b, ShenCastanLow, rows, cols, window, lowThreshold, highThreshold, rawImage, edgeImage);
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ status = *groups;
+
+ return status;
+
+}
+
+void buildBinaryImage(int rows, int cols, double *rawImage, unsigned short *edgeImage, int lowThreshold, int highThreshold){
+
+ int i, j;
+ int offset;
+ double value;
+ int maskValue;
+
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ value = rawImage[offset+j];
+ maskValue = 1;
+ if(value < (double)lowThreshold) maskValue = 0;
+ if(value > (double)highThreshold) maskValue = 0;
+ edgeImage[offset+j] = maskValue;
+ }
+ offset += cols;
+ }
+
+ return;
+
+}
+
+
+
+void morphoFilterBinaryImage(int rows, int cols, unsigned short *edgeImage, int CloseSize, int OpenSize){
+
+
+ int i, j;
+ int offset, offset2;
+ unsigned short cmask[11][11];
+ unsigned short omask[11][11];
+ int olapValuesC[4];
+ int olapValuesO[4];
+ int CloseMaskSize;
+ int OpenMaskSize;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+ int spadSize;
+ unsigned char *ImageE;
+ unsigned char *ImageC;
+
+ spadSize = MAX(rows, cols);
+
+ ImageE = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageC = calloc(spadSize*spadSize, sizeof(unsigned char));
+
+ //
+ // Close filter
+ //
+ if(CloseSize){
+ CloseMaskSize = (CloseSize-1)/2;
+ for(i = 0; i < 2*CloseMaskSize+1; ++i){
+ for(j = 0; j < 2*CloseMaskSize+1; ++j){
+ cmask[i][j] = 1;
+ }
+ }
+ LowValue1 = 0;
+ HighValue1 = 1;
+ LowValue2 = 1;
+ HighValue2 = 0;
+ olapValuesC[0] = LowValue1;
+ olapValuesC[1] = HighValue1;
+ olapValuesC[2] = LowValue2;
+ olapValuesC[3] = HighValue2;
+ }
+
+ //
+ // Open filter
+ //
+ if(OpenSize){
+ OpenMaskSize = (OpenSize-1)/2;
+ for(i = 0; i < 2*OpenMaskSize+1; ++i){
+ for(j = 0; j < 2*OpenMaskSize+1; ++j){
+ omask[i][j] = 1;
+ }
+ }
+ LowValue1 = 1;
+ HighValue1 = 0;
+ LowValue2 = 0;
+ HighValue2 = 1;
+ olapValuesO[0] = LowValue1;
+ olapValuesO[1] = HighValue1;
+ olapValuesO[2] = LowValue2;
+ olapValuesO[3] = HighValue2;
+ }
+
+ offset = 0;
+ offset2 = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ ImageE[offset2+j] = (unsigned char)edgeImage[offset+j];
+ }
+ offset2 += spadSize;
+ offset += cols;
+ }
+
+ if(OpenSize){
+ OpenCloseFilter(olapValuesO, OpenMaskSize, rows, cols, spadSize, ImageE, ImageC, omask);
+ }
+
+ if(CloseSize){
+ OpenCloseFilter(olapValuesC, CloseMaskSize, rows, cols, spadSize, ImageE, ImageC, cmask);
+ }
+
+ offset = 0;
+ offset2 = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(ImageE[offset2+j] == 1){
+ // this will activate some original off-pixels
+ edgeImage[offset+j] = 1;
+ }
+ else{
+ // this will zero some original on-pixels
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset2 += spadSize;
+ offset += cols;
+ }
+
+ free(ImageE);
+ free(ImageC);
+
+ return;
+
+}
+
+void doRegionGrow(int samples, int rows, int cols, double *rawImage, unsigned short *edgeImage, int lowThreshold,
+ int highThreshold, int closeWindow, int openWindow){
+
+ buildBinaryImage(rows, cols, rawImage, edgeImage, lowThreshold, highThreshold);
+ morphoFilterBinaryImage(rows, cols, edgeImage, closeWindow, openWindow);
+
+ return;
+
+}
+
+int NI_RegionGrow(int samples, int rows, int cols, int lowThreshold, int highThreshold, int closeWindow,
+ int openWindow, double *rawImage, unsigned short *edgeImage, int *groups){
+
+ int i, j;
+ int offset;
+ int status;
+
+ doRegionGrow(samples, rows, cols, rawImage, edgeImage, lowThreshold, highThreshold, closeWindow, openWindow);
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ status = *groups;
+ return status;
+
+}
+
+int NI_SobelEdges(int samples, int rows, int cols, double sobelLow, int mode, int lowThreshold, int highThreshold, double BPHigh,
+ int apearture, double *rawImage, unsigned short *edgeImage, int *groups){
+
+
+ int i, j;
+ int offset;
+ int status;
+
+ doPreProcess(samples, rows, cols, rawImage, BPHigh, apearture, lowThreshold, highThreshold);
+ doSobel(samples, rows, cols, sobelLow, mode, rawImage, edgeImage);
+ *groups = ConnectedEdgePoints(rows, cols, edgeImage);
+
+
+ //
+ // prune the isolated pixels
+ //
+ offset = 0;
+ for(i = 0; i < rows; ++i){
+ for(j = 0; j < cols; ++j){
+ if(edgeImage[offset+j] > (*groups)){
+ edgeImage[offset+j] = 0;
+ }
+ }
+ offset += cols;
+ }
+
+ status = *groups;
+ return status;
+
+}
+
+void initThinFilter(int J_mask[3][30], int K_mask[3][30]){
+
+ int i, j;
+ int Column;
+
+ for(i = 0; i < 3; ++i){
+ for(j = 0; j < 30; ++j){
+ J_mask[i][j] = 0;
+ K_mask[i][j] = 0;
+ }
+ }
+
+ Column = 0;
+ J_mask[0][Column+0] = 1;
+ J_mask[0][Column+1] = 1;
+ J_mask[0][Column+2] = 1;
+ J_mask[1][Column+1] = 1;
+
+ Column += 3;
+ J_mask[0][Column+1] = 1;
+ J_mask[1][Column+1] = 1;
+ J_mask[1][Column+2] = 1;
+
+ Column += 3;
+ J_mask[0][Column+0] = 1;
+ J_mask[1][Column+0] = 1;
+ J_mask[2][Column+0] = 1;
+ J_mask[1][Column+1] = 1;
+
+ Column += 3;
+ J_mask[0][Column+1] = 1;
+ J_mask[1][Column+0] = 1;
+ J_mask[1][Column+1] = 1;
+
+ Column += 3;
+ J_mask[0][Column+2] = 1;
+ J_mask[1][Column+1] = 1;
+ J_mask[1][Column+2] = 1;
+ J_mask[2][Column+2] = 1;
+
+ Column += 3;
+ J_mask[1][Column+0] = 1;
+ J_mask[1][Column+1] = 1;
+ J_mask[2][Column+1] = 1;
+
+ Column += 3;
+ J_mask[1][Column+1] = 1;
+ J_mask[2][Column+0] = 1;
+ J_mask[2][Column+1] = 1;
+ J_mask[2][Column+2] = 1;
+
+ Column += 3;
+ J_mask[1][Column+1] = 1;
+ J_mask[1][Column+2] = 1;
+ J_mask[2][Column+1] = 1;
+
+ Column = 0;
+ K_mask[2][Column+0] = 1;
+ K_mask[2][Column+1] = 1;
+ K_mask[2][Column+2] = 1;
+
+ Column += 3;
+ K_mask[1][Column+0] = 1;
+ K_mask[2][Column+0] = 1;
+ K_mask[2][Column+1] = 1;
+
+ Column += 3;
+ K_mask[0][Column+2] = 1;
+ K_mask[1][Column+2] = 1;
+ K_mask[2][Column+2] = 1;
+
+ Column += 3;
+ K_mask[1][Column+2] = 1;
+ K_mask[2][Column+1] = 1;
+ K_mask[2][Column+2] = 1;
+
+ Column += 3;
+ K_mask[0][Column+0] = 1;
+ K_mask[1][Column+0] = 1;
+ K_mask[2][Column+0] = 1;
+
+ Column += 3;
+ K_mask[0][Column+1] = 1;
+ K_mask[0][Column+2] = 1;
+ K_mask[1][Column+2] = 1;
+
+ Column += 3;
+ K_mask[0][Column+0] = 1;
+ K_mask[0][Column+1] = 1;
+ K_mask[0][Column+2] = 1;
+
+ Column += 3;
+ K_mask[0][Column+0] = 1;
+ K_mask[0][Column+1] = 1;
+ K_mask[1][Column+0] = 1;
+
+ return;
+
+}
+
+void ThinningFilter(int regRows, int regColumns, int spadSize, int J_mask[3][30], int K_mask[3][30],
+ unsigned char *Input, unsigned char *CInput, unsigned char *ErosionStage,
+ unsigned char *DialationStage, unsigned char *HMT, unsigned char *Copy){
+
+ int i, j, k, l, m, n, overlap, hit;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+ int Column, T, nloop;
+ int Offset;
+ int N, M;
+ int j_mask[3][3];
+ int k_mask[3][3];
+
+ N = regRows;
+ M = regColumns;
+
+ LowValue1 = 1;
+ HighValue1 = 0;
+
+ LowValue2 = 0;
+ HighValue2 = 1;
+
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ Copy[Offset+j] = Input[Offset+j];
+ }
+ Offset += spadSize;
+ }
+
+ nloop = 0;
+ while(1){
+ // erode
+ Column = 0;
+ for(n = 0; n < 8; ++n){
+ for(i = 0; i < 3; ++i){
+ for(j = 0; j < 3; ++j){
+ j_mask[i][j] = J_mask[i][Column+j];
+ }
+ }
+ for(i = 0; i < 3; ++i){
+ for(j = 0; j < 3; ++j){
+ k_mask[i][j] = K_mask[i][Column+j];
+ }
+ }
+ Column += 3;
+
+ Offset = spadSize;
+ for(i = 1; i < N-1; ++i){
+ for(j = 1; j < M-1; ++j){
+ hit = LowValue1;
+ for(k = -1; k < 2; ++k){
+ for(l = -1; l < 2; ++l){
+ T = j_mask[k+1][l+1];
+ if(T == 1){
+ overlap = T*Input[Offset+(k*spadSize)+j+l];
+ if(overlap == HighValue1) hit = HighValue1;
+ }
+ }
+ }
+ ErosionStage[Offset+j] = hit;
+ }
+ Offset += spadSize;
+ }
+
+ // dialate
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ CInput[Offset+j] = (~Input[Offset+j]) & 0x1;
+ }
+ Offset += spadSize;
+ }
+
+ Offset = spadSize;
+ for(i = 1; i < N-1; ++i){
+ for(j = 1; j < M-1; ++j){
+ hit = LowValue1;
+ for(k = -1; k < 2; ++k){
+ for(l = -1; l < 2; ++l){
+ T = k_mask[k+1][l+1];
+ if(T == 1){
+ overlap = T*CInput[Offset+(k*spadSize)+j+l];
+ if(overlap == HighValue1) hit = HighValue1;
+ }
+ }
+ }
+ DialationStage[Offset+j] = hit;
+ }
+ Offset += spadSize;
+ }
+
+ // form the HMT
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ m = (ErosionStage[Offset+j]*DialationStage[Offset+j]);
+ HMT[Offset+j] = m;
+ }
+ Offset += spadSize;
+ }
+
+ // Thin for stage n
+
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ HMT[Offset+j] = (~HMT[Offset+j]) & 0x1;
+ }
+ Offset += spadSize;
+ }
+
+ Offset = 0;
+ for (i = 0; i < N; ++i){
+ for (j = 0; j < M; ++j){
+ m = (Input[Offset+j]*HMT[Offset+j]);
+ Input[Offset+j] = m;
+ }
+ Offset += spadSize;
+ }
+ }
+
+ // check for the NULL set
+ hit = 0;
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ hit += abs(Copy[Offset+j]-Input[Offset+j]);
+ }
+ Offset += spadSize;
+ }
+ if(!hit) break;
+
+ hit = 0;
+ Offset = 0;
+ for(i = 0; i < N; ++i){
+ for(j = 0; j < M; ++j){
+ Copy[Offset+j] = Input[Offset+j];
+ if(Input[Offset+j]) ++hit;
+ }
+ Offset += spadSize;
+ }
+ ++nloop;
+ }
+
+
+ return;
+
+}
+
+
+int NI_ThinFilter(int samples, int rows, int cols, int numberObjects, unsigned short *edgeImage, objStruct objectMetrics[]){
+
+ int i, j;
+ int loop;
+ int label;
+ int left, right, top, bottom;
+ int roiRows, roiCols;
+ int srcOffset;
+ int dstOffset;
+ int status;
+ int inflate = 1;
+ int J_mask[3][30];
+ int K_mask[3][30];
+
+ unsigned char *Input;
+ unsigned char *CInput;
+ unsigned char *ErosionStage;
+ unsigned char *DialationStage;
+ unsigned char *HMT;
+ unsigned char *Copy;
+ unsigned short *thinEdgeImage;
+
+ //
+ // scratch pad (spad) memory
+ //
+ Input = calloc(samples, sizeof(unsigned char));
+ CInput = calloc(samples, sizeof(unsigned char));
+ ErosionStage = calloc(samples, sizeof(unsigned char));
+ DialationStage = calloc(samples, sizeof(unsigned char));
+ HMT = calloc(samples, sizeof(unsigned char));
+ Copy = calloc(samples, sizeof(unsigned char));
+ thinEdgeImage = calloc(samples, sizeof(unsigned short));
+
+ initThinFilter(J_mask, K_mask);
+ for(loop = 0; loop < numberObjects; ++loop){
+ label = objectMetrics[loop].Label;
+ left = objectMetrics[loop].L;
+ right = objectMetrics[loop].R;
+ top = objectMetrics[loop].T;
+ bottom = objectMetrics[loop].B;
+ roiRows = top-bottom+2*inflate;
+ roiCols = right-left+2*inflate;
+
+ //
+ // clear the scratch pad
+ //
+ srcOffset = 0;
+ for(i = 0; i < roiRows; ++i){
+ for(j = 0; j < roiCols; ++j){
+ Input[srcOffset+j] = 0;
+ }
+ srcOffset += cols;
+ }
+
+ //
+ // copy the ROI for MAT (medial axis transformation) filter
+ //
+ dstOffset = inflate*rows;
+ for(i = bottom; i < top; ++i){
+ srcOffset = i*cols;
+ for(j = left; j < right; ++j){
+ if(edgeImage[srcOffset+j] == label){
+ Input[dstOffset+j-left+inflate] = 1;
+ }
+ }
+ dstOffset += cols;
+ }
+ ThinningFilter(roiRows, roiCols, cols, J_mask, K_mask, Input, CInput, ErosionStage, DialationStage, HMT, Copy);
+
+ //
+ // copy the MAT roi to the new edgeImage (clip the inflate border)
+ //
+ dstOffset = inflate*rows;
+ for(i = bottom; i < top; ++i){
+ srcOffset = i*cols;
+ for(j = left; j < right; ++j){
+ if(Input[dstOffset+j-left+inflate]){
+ thinEdgeImage[srcOffset+j] = label;
+ }
+ }
+ dstOffset += cols;
+ }
+ }
+
+ //
+ // copy the MAT edges and return the thinned edges
+ // this will prune the isolated edge points from the edgeImage source
+ //
+ for(i = 0; i < rows*cols; ++i){
+ edgeImage[i] = thinEdgeImage[i];
+ }
+
+ free(Input);
+ free(CInput);
+ free(ErosionStage);
+ free(DialationStage);
+ free(HMT);
+ free(Copy);
+ free(thinEdgeImage);
+
+ status = 1;
+
+ return status;
+
+}
+
+
+void generateMask(unsigned char *ImageH, bPOINT *boundary, int newSamples, int label, int cols){
+
+ //
+ // get the boundary point pairs (left, right) for each line
+ // if there is no pair, then the boundary is open
+ // then fill the image in with the current label
+ //
+
+ int i, j, k, m;
+ int list[512];
+ 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 = 0;
+ NewJ = 0;
+ 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;
+ }
+ CurI = NewI;
+ CurJ = NewJ;
+ offset = CurI * spadSize;
+ Input[offset+CurJ] = 0;
+ lBoundary[p].xy.x = CurJ;
+ lBoundary[p].xy.y = CurI;
+ ++p;
+ }
+
+ return;
+
+}
+
+
+void OpenCloseFilter(int olapValues[], int maskSize, int rows, int columns, int spadSize,
+ unsigned char *input, unsigned char *output, unsigned short mask[11][11]){
+
+
+ //
+ // do morphological open/close image filtering. the olapValues array determines
+ // if the filter is Open or Close.
+ //
+ int i, j, k, l, m, overlap, hit;
+ int offset;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+
+ LowValue1 = olapValues[0];
+ HighValue1 = olapValues[1];
+ LowValue2 = olapValues[2];
+ HighValue2 = olapValues[3];
+
+ // close - step 1 is dialate
+ // open - step 1 is erode
+ offset = maskSize*spadSize;
+ for(i = maskSize; i < rows-maskSize; ++i){
+ for(j = maskSize; j < columns-maskSize; ++j){
+ hit = LowValue1;
+ for(k = -maskSize; k < maskSize; ++k){
+ m = k*spadSize;
+ for(l = -maskSize; l < maskSize; ++l){
+ overlap = mask[k+maskSize][l+maskSize]*input[offset+m+j+l];
+ if(overlap == HighValue1){
+ hit = HighValue1;
+ }
+ }
+ }
+ output[offset+j] = hit;
+ }
+ offset += spadSize;
+ }
+
+ // close - step 2 is erode
+ // open - step 2 is dialate
+ offset = maskSize*spadSize;
+ for(i = maskSize; i < rows-maskSize; ++i){
+ for(j = maskSize; j < columns-maskSize; ++j){
+ hit = LowValue2;
+ for(k = -maskSize; k < maskSize; ++k){
+ m = k*spadSize;
+ for(l = -maskSize; l < maskSize; ++l){
+ overlap = mask[k+maskSize][l+maskSize]*output[offset+m+j+l];
+ if(overlap == HighValue2){
+ hit = HighValue2;
+ }
+ }
+ }
+ input[offset+j] = hit;
+ }
+ offset += spadSize;
+ }
+
+ return;
+}
+
+void getCompactness(unsigned char *Input, RECT roi, int label, int spadSize, float *vCompactness, float length){
+
+ int i, j;
+ int maskOffset;
+ int area;
+ static float fpi = (float)(4.0 * 3.14159);
+
+ area = 0;
+ for(i = roi.bottom; i < roi.top; ++i){
+ maskOffset = i*spadSize;
+ for(j = roi.left; j < roi.right; ++j){
+ if(Input[maskOffset+j] == label){
+ ++area;
+ }
+ }
+ }
+ if(area && (length != (float)0.0)){
+ *vCompactness = (fpi * (float)area) / (length*length);
+ }
+ else{
+ *vCompactness = (float)0.0;
+ }
+
+ return;
+}
+
+
+void doMorphology(unsigned char *Input, unsigned char *ImageE, unsigned char *ImageC, unsigned char *ImageH,
+ int olapValuesC[],int olapValuesO[], unsigned short cmask[11][11], unsigned short omask[11][11],
+ RECT roi, int label, int CloseMaskSize, int OpenMaskSize, int spadSize){
+
+ int i, j;
+ int rows, cols;
+ int srcOffset;
+ int dstOffset;
+ int maskSize;
+
+ cols = roi.right - roi.left;
+ rows = roi.top - roi.bottom;
+
+ for(i = 0; i < spadSize*spadSize; ++i){
+ ImageE[i] = 0;
+ ImageC[i] = 0;
+ }
+
+ //
+ // put the ROI in the ImageE array centered in ULC
+ //
+ dstOffset = 0;
+ for(i = roi.bottom; i < roi.top; ++i){
+ srcOffset = i*spadSize;
+ for(j = roi.left; j < roi.right; ++j){
+ if(ImageH[srcOffset+j] == label){
+ ImageE[dstOffset+j-roi.left] = 1;
+ }
+ }
+ dstOffset += spadSize;
+ }
+
+ //
+ // open
+ //
+ maskSize = OpenMaskSize;
+ OpenCloseFilter(olapValuesO, maskSize, rows, cols, spadSize, ImageE, ImageC, omask);
+ //
+ // close
+ //
+ maskSize = CloseMaskSize;
+ OpenCloseFilter(olapValuesC, maskSize, rows, cols, spadSize, ImageE, ImageC, cmask);
+
+ //
+ // put the closed ROI (in ImageE) back in its roi space
+ //
+
+ srcOffset = 0;
+ for(i = roi.bottom; i < roi.top+2*maskSize+1; ++i){
+ dstOffset = (i-(2*maskSize+1))*spadSize;
+ for(j = roi.left-maskSize-1; j < roi.right+maskSize+1; ++j){
+ if(ImageE[srcOffset+j-roi.left] == 1){
+ Input[dstOffset+j-maskSize+1] = label;
+ }
+ }
+ srcOffset += spadSize;
+ }
+
+ return;
+
+}
+
+
+void getBoundary(unsigned short *ThinEdgeImage, unsigned char *Input, blobBoundary *pBoundary, blobBoundary *lBoundary,
+ boundaryIndex *pBoundaryIndex, RECT boundBox, int label, int bBox, int nextSlot, int memOffset,
+ int spadSize, int searchWindow){
+
+ int i, j;
+ int dstOffset;
+ int srcOffset;
+ int mcount;
+ int rows;
+ int columns;
+ bool first;
+ blobBoundary value;
+ int inflate = 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 should be setup parameters
+ int closureDistance = 12;
+ int CloseSize = 5;
+ int OpenSize = 5;
+ int threshold = 3;
+ int newSamples;
+ int spadSize;
+ POINT rectPoint[4];
+ int in[4];
+ float length;
+ float minRadius;
+ float maxRadius;
+ float aveRadius;
+ float vCompactness;
+ // for morphological close of mask. max structuring element is 11x11
+ unsigned short cmask[11][11];
+ unsigned short omask[11][11];
+ int olapValuesC[4];
+ int olapValuesO[4];
+ int CloseMaskSize;
+ int OpenMaskSize;
+ int LowValue1, HighValue1;
+ int LowValue2, HighValue2;
+
+ //
+ // Close filter
+ //
+ CloseMaskSize = (CloseSize-1)/2;
+ for(i = 0; i < 2*CloseMaskSize+1; ++i){
+ for(j = 0; j < 2*CloseMaskSize+1; ++j){
+ cmask[i][j] = 1;
+ }
+ }
+ LowValue1 = 0;
+ HighValue1 = 1;
+ LowValue2 = 1;
+ HighValue2 = 0;
+ olapValuesC[0] = LowValue1;
+ olapValuesC[1] = HighValue1;
+ olapValuesC[2] = LowValue2;
+ olapValuesC[3] = HighValue2;
+
+ //
+ // Open filter
+ //
+ OpenMaskSize = (OpenSize-1)/2;
+ for(i = 0; i < 2*OpenMaskSize+1; ++i){
+ for(j = 0; j < 2*OpenMaskSize+1; ++j){
+ omask[i][j] = 1;
+ }
+ }
+ LowValue1 = 1;
+ HighValue1 = 0;
+ LowValue2 = 0;
+ HighValue2 = 1;
+ olapValuesO[0] = LowValue1;
+ olapValuesO[1] = HighValue1;
+ olapValuesO[2] = LowValue2;
+ olapValuesO[3] = HighValue2;
+
+ RECT bBox;
+ boundaryIndex *pBoundaryIndex;
+ blobBoundary *pBoundary;
+ blobBoundary *lBoundary;
+ bPOINT *boundary;
+ unsigned char *Input;
+ unsigned char *ImageE;
+ unsigned char *ImageC;
+ unsigned char *ImageH;
+
+ spadSize = MAX(srcRows, srcCols);
+ pBoundaryIndex = calloc(srcRows+srcCols, sizeof(boundaryIndex));
+ Input = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageE = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageC = calloc(spadSize*spadSize, sizeof(unsigned char));
+ ImageH = calloc(spadSize*spadSize, sizeof(unsigned char));
+ pBoundary = calloc(srcRows*srcCols, sizeof(blobBoundary));
+ lBoundary = calloc(32767, sizeof(blobBoundary));
+ boundary = calloc(32767, sizeof(POINT));
+
+ for(i = 0; i < (srcRows+srcCols); ++i){
+ pBoundaryIndex[i].numberPoints = 0;
+ pBoundaryIndex[i].curveClose = 0;
+ pBoundaryIndex[i].isWithin = FALSE;
+ pBoundaryIndex[i].criticalSize = FALSE;
+ pBoundaryIndex[i].closedCurve = FALSE;
+ }
+
+ for(i = 0; i < numberObjects; ++i){
+ ++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
+ //
+
+ //
+ // numBoundaries = numberObjects
+ //
+ 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);
+
+ return;
+
+}
+
+
+void initLaws(LawsFilter7 *lawsFilter){
+
+ int i;
+ float sum;
+ float L7[7] = { 1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0};
+ float E7[7] = {-1.0, -4.0, -5.0, 0.0, 5.0, 4.0, 1.0};
+ float S7[7] = {-1.0, -2.0, 1.0, 4.0, 1.0, -2.0, -1.0};
+ float W7[7] = {-1.0, 0.0, 3.0, 0.0, -3.0, 0.0, 1.0};
+ float R7[7] = { 1.0, -2.0, -1.0, 4.0, -1.0, -2.0, 1.0};
+ float O7[7] = {-1.0, 6.0, -15.0, 20.0, -15.0, 6.0, -1.0};
+
+ lawsFilter->numberKernels = 6;
+ lawsFilter->kernelLength = 7;
+ lawsFilter->numberFilterLayers = 21;
+ lawsFilter->name[0] = 'L';
+ lawsFilter->name[1] = 'E';
+ lawsFilter->name[2] = 'S';
+ lawsFilter->name[3] = 'W';
+ lawsFilter->name[4] = 'R';
+ lawsFilter->name[5] = 'O';
+ for(i = 0; i < 7; ++i){
+ lawsFilter->lawsKernel[0][i] = L7[i];
+ lawsFilter->lawsKernel[1][i] = E7[i];
+ lawsFilter->lawsKernel[2][i] = S7[i];
+ lawsFilter->lawsKernel[3][i] = W7[i];
+ lawsFilter->lawsKernel[4][i] = R7[i];
+ lawsFilter->lawsKernel[5][i] = O7[i];
+ }
+
+ // DC filter is unity gain
+ sum = (float)0.0;
+ for(i = 0; i < 7; ++i){
+ sum += lawsFilter->lawsKernel[0][i];
+ }
+ for(i = 0; i < 7; ++i){
+ lawsFilter->lawsKernel[0][i] /= sum;
+ }
+
+ return;
+
+}
+
+float lawsConvolution(float *image, float *rowFilter, float *colFilter, int kernelSize){
+
+ int i, j;
+ int offset;
+ float result[7];
+ float sum;
+
+ // filter rows
+ for(i = 0; i < kernelSize; ++i){
+ sum = (float)0.0;
+ offset = i * kernelSize;
+ for(j = 0; j < kernelSize; ++j){
+ sum += (rowFilter[j]*image[offset+j]);
+ }
+ result[i] = sum;
+ }
+
+ // filter columns
+ sum = (float)0.0;
+ for(j = 0; j < kernelSize; ++j){
+ sum += (rowFilter[j]*result[j]);
+ }
+
+ return(sum);
+
+}
+
+
+void getLawsTexture(LawsFilter7 lawsFilter, tTEM LawsFeatures[], objStruct objectMetrics[], double *sourceImage,
+ unsigned short *MaskImage, int numberObjects, int srcRows, int srcCols){
+
+ int i, j;
+ int label;
+ RECT bBox;
+ int aperature = (lawsFilter.kernelLength-1)/2;
+ unsigned char *ImageH;
+ float *ImageT;
+ float *lawsImage;
+
+ /*
+ // TO DO. return all 21 filter slices from a user selected blob
+ int debugBlob = 8;
+ int index;
+ int offset;
+ int layerStep = srcRows*srcCols;
+ */
+
+ 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;
+ }
+ /*
+ 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 with the texture filter that has the largest variance
+ offset = index * layerStep;
+ for(j = 0; j < layerStep; ++j){
+ sourceImage[j] = lawsImage[offset+j];
+ }
+ }
+ */
+ }
+ }
+
+ free(ImageH);
+ free(ImageT);
+ free(lawsImage);
+
+ return;
+
+}
+
+void computeLaws(LawsFilter7 lawsFilter, tTEM LawsFeatures[], RECT roi, int label, int aperature, int srcRows, int srcCols,
+ unsigned char *ImageH, float *ImageT, unsigned short *MaskImage, float *lawsImage, double *sourceImage){
+
+ //
+ // hard-wirred to Law's 7 kernels
+ //
+ int i, j, k;
+ int lawsLayer;
+ int column, row;
+ int offset;
+ int maskOffset[7];
+ int dataOffset[7];
+ float myImage[49];
+ int count;
+ int outerKernelNumber;
+ int innerKernelNumber;
+ int rowNumber;
+ int kernelSize = lawsFilter.kernelLength;
+ int fullMask = kernelSize*kernelSize;
+ int layerStep = srcRows*srcCols;
+ float *rowFilter;
+ float *colFilter;
+ float filterResult1;
+ float filterResult2;
+ float lawsLL;
+ float t;
+ float maxValue;
+ float scale;
+ char I, J;
+ char combo[24];
+ char dual[24];
+
+
+ // zero the laws mask memory first
+ for(i = 0; i < srcRows*srcCols; ++i){
+ ImageH[i] = 0;
+ }
+ for(j = 0; j < lawsFilter.numberFilterLayers; ++j){
+ LawsFeatures[j].Mean = (float)0.0;
+ LawsFeatures[j].Variance = (float)0.0;
+ }
+
+ for(i = roi.bottom+aperature; i < roi.top-aperature; ++i){
+ // get the row array offset for mask and data source.
+ for(row = -aperature; row <= aperature; ++row){
+ maskOffset[row+aperature] = (i+row)*srcCols;
+ dataOffset[row+aperature] = maskOffset[row+aperature];
+ }
+ for(j = roi.left+aperature; j < roi.right-aperature; ++j){
+ //
+ // get 7x7 segment and make sure have 100% mask coverage
+ //
+ count = 0;
+ for(row = -aperature; row <= aperature; ++row){
+ rowNumber = (row+aperature)*kernelSize;
+ for(column = -aperature; column <= aperature; ++column){
+ if(MaskImage[maskOffset[row+aperature]+j+column] == label){
+ myImage[rowNumber+column+aperature] = sourceImage[dataOffset[row+aperature]+j+column];
+ ++count;
+ }
+ }
+ }
+ if(count == fullMask){
+ //
+ // full mask. 100% coverage. now do the Law's texture filters
+ //
+ ImageH[i*srcCols+j] = 1;
+ lawsLayer = 0;
+ for(outerKernelNumber = 0; outerKernelNumber < lawsFilter.numberKernels; ++outerKernelNumber){
+ //
+ // outer loop pulls the i'th kernel. kernel 0 is the LP kernel
+ // the outer loop is the iso-kernel
+ //
+ I = lawsFilter.name[outerKernelNumber];
+ sprintf(dual, "%c_%c", I, I);
+ rowFilter = &lawsFilter.lawsKernel[outerKernelNumber][0];
+ colFilter = &lawsFilter.lawsKernel[outerKernelNumber][0];
+ filterResult1 = lawsConvolution(myImage, rowFilter, colFilter, kernelSize);
+ // lawsLayer 0 is the LP and needs to be used to scale.
+ if(outerKernelNumber){
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] = (float)2.0 * filterResult1 / lawsLL;
+ }
+ else{
+ lawsLL = (float)2.0 * filterResult1;
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] = (float)2.0 * filterResult1;
+ }
+ strcpy(&LawsFeatures[lawsLayer].filterName[0], dual);
+ ++lawsLayer;
+ //
+ // now do the inner loop and get the column filters for the other laws kernels
+ //
+ for(innerKernelNumber = outerKernelNumber+1; innerKernelNumber < lawsFilter.numberKernels; ++innerKernelNumber){
+ J = lawsFilter.name[innerKernelNumber];
+ sprintf(combo, "%c_%c", I, J);
+ strcpy(&LawsFeatures[lawsLayer].filterName[0], combo);
+ colFilter = &lawsFilter.lawsKernel[innerKernelNumber][0];
+ filterResult1 = lawsConvolution(myImage, rowFilter, colFilter, kernelSize);
+ filterResult2 = lawsConvolution(myImage, colFilter, rowFilter, kernelSize);
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] = (filterResult1 / lawsLL) + (filterResult2 / lawsLL);
+ ++lawsLayer;
+ }
+ }
+ }
+ }
+ }
+
+ for(i = 0; i < lawsFilter.numberFilterLayers; ++i){
+ LawsFeatures[i].Mean = (float)0.0;
+ LawsFeatures[i].Variance = (float)0.0;
+ }
+
+ count = 0;
+ for(i = roi.bottom+aperature; i < roi.top-aperature; ++i){
+ row = i * srcCols;
+ for(j = roi.left+aperature; j < roi.right-aperature; ++j){
+ if(ImageH[row+j]){
+ ++count;
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ offset = k * layerStep + row;
+ LawsFeatures[k].Mean += lawsImage[offset+j];
+ }
+ }
+ }
+ }
+
+ if(count == 0){
+ //printf("no samples for texture\n");
+ return;
+ }
+
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ LawsFeatures[k].Mean /= (float)count;
+ }
+ for(i = roi.bottom+aperature; i < roi.top-aperature; ++i){
+ row = i * srcCols;
+ for(j = roi.left+aperature; j < roi.right-aperature; ++j){
+ if(ImageH[row+j]){
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ offset = k * layerStep + row;
+ t = lawsImage[offset+j] - LawsFeatures[k].Mean;
+ LawsFeatures[k].Variance += (t * t);
+ }
+ }
+ }
+ }
+ for(k = 0; k < lawsFilter.numberFilterLayers; ++k){
+ LawsFeatures[k].Variance /= (float)count;
+ LawsFeatures[k].Variance = (float)(sqrt(LawsFeatures[k].Variance));
+ }
+
+ //
+ // now normalize the variance feature (TEM)
+ //
+ maxValue = (float)0.0;
+ for(i = 1; i < lawsFilter.numberFilterLayers; ++i){
+ if((LawsFeatures[i].Variance) > maxValue) maxValue = LawsFeatures[i].Variance;
+ }
+ scale = (float)1.0 / maxValue;
+ for(i = 1; i < lawsFilter.numberFilterLayers; ++i){
+ LawsFeatures[i].Variance = scale * LawsFeatures[i].Variance;
+ }
+
+
+ return;
+
+}
+
+void getVoxelMeasures(objStruct objectMetrics[], double *sourceImage, unsigned short *MaskImage,
+ int numberObjects, int srcRows, int srcCols){
+
+ int i, j, k;
+ int label;
+ int offset;
+ int count;
+ float mean, std, t;
+ RECT bBox;
+
+ for(i = 0; i < numberObjects; ++i){
+ bBox.left = objectMetrics[i].L;
+ bBox.right = objectMetrics[i].R;
+ bBox.top = objectMetrics[i].T;
+ bBox.bottom = objectMetrics[i].B;
+ label = objectMetrics[i].Label;
+ count = 0;
+ mean = (float)0.0;
+ for(j = bBox.bottom; j < bBox.top; ++j){
+ offset = j * srcCols;
+ for(k = bBox.left; k < bBox.right; ++k){
+ if(MaskImage[offset+k] == label){
+ mean += sourceImage[offset+k];
+ ++count;
+ }
+ }
+ }
+ if(count){
+ mean /= (float)count;
+ std = (float)0.0;
+ for(j = bBox.bottom; j < bBox.top; ++j){
+ offset = j * srcCols;
+ for(k = bBox.left; k < bBox.right; ++k){
+ if(MaskImage[offset+k] == label){
+ t = (sourceImage[offset+k]-mean);
+ std += (t * t);
+ }
+ }
+ }
+ }
+ if(count){
+ std /= (float)count;
+ std = sqrt(std);
+ objectMetrics[i].voxelMean = mean - 2048.0; // the 2048 is only for the cardiac CT volume
+ objectMetrics[i].voxelVar = std;
+ }
+ else{
+ objectMetrics[i].voxelMean = 0.0;
+ objectMetrics[i].voxelVar = 0.0;
+ }
+ if(0) printf("(%d) mean %f, std %f\n", label, objectMetrics[i].voxelMean, objectMetrics[i].voxelVar);
+ }
+
+ return;
+
+}
+
+int NI_BuildBoundary(int samples, int rows, int cols, int numberObjects, unsigned short *edgeImage,
+ objStruct objectMetrics[]){
+
+ int searchWindow = 5; // 5 is good value for Sobel - (should be 13 for Canny edges)
+ int status = 1;
+
+ buildBoundary(objectMetrics, searchWindow, edgeImage, numberObjects, rows, cols);
+
+ return status;
+
+}
+
+int NI_VoxelMeasures(int samples, int rows, int cols, int numberObjects, double *sourceImage,
+ unsigned short *maskImage, objStruct objectMetrics[]){
+
+ int status = 1;
+ getVoxelMeasures(objectMetrics, sourceImage, maskImage, numberObjects, rows, cols);
+
+ return status;
+
+}
+
+
+int NI_TextureMeasures(int samples, int rows, int cols, int numberObjects, double *sourceImage,
+ unsigned short *maskImage, objStruct objectMetrics[]){
+
+ int status = 1;
+ LawsFilter7 lawsFilter;
+ tTEM LawsFeatures[21];
+
+ initLaws(&lawsFilter);
+ getLawsTexture(lawsFilter, LawsFeatures, objectMetrics, sourceImage, maskImage, numberObjects, rows, cols);
+
+ return status;
+
+}
+
+
+
Property changes on: branches/ndimage_segmenter/scipy/NDISegmenter/Segmenter_IMPL.c
___________________________________________________________________
Name: svn:executable
+ *
Added: branches/ndimage_segmenter/scipy/NDISegmenter/__multiarray_api.h
===================================================================
--- branches/ndimage_segmenter/scipy/NDISegmenter/__multiarray_api.h 2007-11-02 19:51:44 UTC (rev 3491)
+++ branches/ndimage_segmenter/scipy/NDISegmenter/__multiarray_api.h 2007-11-02 20:09:04 UTC (rev 3492)
@@ -0,0 +1,974 @@
+
+#ifdef _MULTIARRAYMODULE
+
+typedef struct {
+ PyObject_HEAD
+ npy_bool obval;
+} PyBoolScalarObject;
+
+
+static unsigned int PyArray_GetNDArrayCVersion (void);
+static PyTypeObject PyBigArray_Type;
+static PyTypeObject PyArray_Type;
+static PyTypeObject PyArrayDescr_Type;
+static PyTypeObject PyArrayFlags_Type;
+static PyTypeObject PyArrayIter_Type;
+static PyTypeObject PyArrayMapIter_Type;
+static PyTypeObject PyArrayMultiIter_Type;
+static int NPY_NUMUSERTYPES=0;
+static PyTypeObject PyBoolArrType_Type;
+static PyBoolScalarObject _PyArrayScalar_BoolValues[2];
+
+static PyTypeObject PyGenericArrType_Type;
+static PyTypeObject PyNumberArrType_Type;
+static PyTypeObject PyIntegerArrType_Type;
+static PyTypeObject PySignedIntegerArrType_Type;
+static PyTypeObject PyUnsignedIntegerArrType_Type;
+static PyTypeObject PyInexactArrType_Type;
+static PyTypeObject PyFloatingArrType_Type;
+static PyTypeObject PyComplexFloatingArrType_Type;
+static PyTypeObject PyFlexibleArrType_Type;
+static PyTypeObject PyCharacterArrType_Type;
+static PyTypeObject PyByteArrType_Type;
+static PyTypeObject PyShortArrType_Type;
+static PyTypeObject PyIntArrType_Type;
+static PyTypeObject PyLongArrType_Type;
+static PyTypeObject PyLongLongArrType_Type;
+static PyTypeObject PyUByteArrType_Type;
+static PyTypeObject PyUShortArrType_Type;
+static PyTypeObject PyUIntArrType_Type;
+static PyTypeObject PyULongArrType_Type;
+static PyTypeObject PyULongLongArrType_Type;
+static PyTypeObject PyFloatArrType_Type;
+static PyTypeObject PyDoubleArrType_Type;
+static PyTypeObject PyLongDoubleArrType_Type;
+static PyTypeObject PyCFloatArrType_Type;
+static PyTypeObject PyCDoubleArrType_Type;
+static PyTypeObject PyCLongDoubleArrType_Type;
+static PyTypeObject PyObjectArrType_Type;
+static PyTypeObject PyStringArrType_Type;
+static PyTypeObject PyUnicodeArrType_Type;
+static PyTypeObject PyVoidArrType_Type;
+static int PyArray_SetNumericOps \
+ (PyObject *);
+static PyObject * PyArray_GetNumericOps \
+ (void);
+static int PyArray_INCREF \
+ (PyArrayObject *);
+static int PyArray_XDECREF \
+ (PyArrayObject *);
+static void PyArray_SetStringFunction \
+ (PyObject *, int);
+static PyArray_Descr * PyArray_DescrFromType \
+ (int);
+static PyObject * PyArray_TypeObjectFromType \
+ (int);
+static char * PyArray_Zero \
+ (PyArrayObject *);
+static char * PyArray_One \
+ (PyArrayObject *);
+static PyObject * PyArray_CastToType \
+ (PyArrayObject *, PyArray_Descr *, int);
+static int PyArray_CastTo \
+ (PyArrayObject *, PyArrayObject *);
+static int PyArray_CastAnyTo \
+ (PyArrayObject *, PyArrayObject *);
+static int PyArray_CanCastSafely \
+ (int, int);
+static npy_bool PyArray_CanCastTo \
+ (PyArray_Descr *, PyArray_Descr *);
+static int PyArray_ObjectType \
+ (PyObject *, int);
+static PyArray_Descr * PyArray_DescrFromObject \
+ (PyObject *, PyArray_Descr *);
+static PyArrayObject ** PyArray_ConvertToCommonType \
+ (PyObject *, int *);
+static PyArray_Descr * PyArray_DescrFromScalar \
+ (PyObject *);
+static PyArray_Descr * PyArray_DescrFromTypeObject \
+ (PyObject *);
+static npy_intp PyArray_Size \
+ (PyObject *);
+static PyObject * PyArray_Scalar \
+ (void *, PyArray_Descr *, PyObject *);
+static PyObject * PyArray_FromScalar \
+ (PyObject *, PyArray_Descr *);
+static void PyArray_ScalarAsCtype \
+ (PyObject *, void *);
+static int PyArray_CastScalarToCtype \
+ (PyObject *, void *, PyArray_Descr *);
+static int PyArray_CastScalarDirect \
+ (PyObject *, PyArray_Descr *, void *, int);
+static PyObject * PyArray_ScalarFromObject \
+ (PyObject *);
+static PyArray_VectorUnaryFunc * PyArray_GetCastFunc \
+ (PyArray_Descr *, int);
+static PyObject * PyArray_FromDims \
+ (int, int *, int);
+static PyObject * PyArray_FromDimsAndDataAndDescr \
+ (int, int *, PyArray_Descr *, char *);
+static PyObject * PyArray_FromAny \
+ (PyObject *, PyArray_Descr *, int, int, int, PyObject *);
+static PyObject * PyArray_EnsureArray \
+ (PyObject *);
+static PyObject * PyArray_EnsureAnyArray \
+ (PyObject *);
+static PyObject * PyArray_FromFile \
+ (FILE *, PyArray_Descr *, npy_intp, char *);
+static PyObject * PyArray_FromString \
+ (char *, npy_intp, PyArray_Descr *, npy_intp, char *);
+static PyObject * PyArray_FromBuffer \
+ (PyObject *, PyArray_Descr *, npy_intp, npy_intp);
+static PyObject * PyArray_FromIter \
+ (PyObject *, PyArray_Descr *, npy_intp);
+static PyObject * PyArray_Return \
+ (PyArrayObject *);
+static PyObject * PyArray_GetField \
+ (PyArrayObject *, PyArray_Descr *, int);
+static int PyArray_SetField \
+ (PyArrayObject *, PyArray_Descr *, int, PyObject *);
+static PyObject * PyArray_Byteswap \
+ (PyArrayObject *, npy_bool);
+static PyObject * PyArray_Resize \
+ (PyArrayObject *, PyArray_Dims *, int, NPY_ORDER);
+static int PyArray_MoveInto \
+ (PyArrayObject *, PyArrayObject *);
+static int PyArray_CopyInto \
+ (PyArrayObject *, PyArrayObject *);
+static int PyArray_CopyAnyInto \
+ (PyArrayObject *, PyArrayObject *);
+static int PyArray_CopyObject \
+ (PyArrayObject *, PyObject *);
+static PyObject * PyArray_NewCopy \
+ (PyArrayObject *, NPY_ORDER);
+static PyObject * PyArray_ToList \
+ (PyArrayObject *);
+static PyObject * PyArray_ToString \
+ (PyArrayObject *, NPY_ORDER);
+static int PyArray_ToFile \
+ (PyArrayObject *, FILE *, char *, char *);
+static int PyArray_Dump \
+ (PyObject *, PyObject *, int);
+static PyObject * PyArray_Dumps \
+ (PyObject *, int);
+static int PyArray_ValidType \
+ (int);
+static void PyArray_UpdateFlags \
+ (PyArrayObject *, int);
+static PyObject * PyArray_New \
+ (PyTypeObject *, int, npy_intp *, int, npy_intp *, void *, int, int, PyObject *);
+static PyObject * PyArray_NewFromDescr \
+ (PyTypeObject *, PyArray_Descr *, int, npy_intp *, npy_intp *, void *, int, PyObject *);
+static PyArray_Descr * PyArray_DescrNew \
+ (PyArray_Descr *);
+static PyArray_Descr * PyArray_DescrNewFromType \
+ (int);
+static double PyArray_GetPriority \
+ (PyObject *, double);
+static PyObject * PyArray_IterNew \
+ (PyObject *);
+static PyObject * PyArray_MultiIterNew \
+ (int, ...);
+static int PyArray_PyIntAsInt \
+ (PyObject *);
+static npy_intp PyArray_PyIntAsIntp \
+ (PyObject *);
+static int PyArray_Broadcast \
+ (PyArrayMultiIterObject *);
+static void PyArray_FillObjectArray \
+ (PyArrayObject *, PyObject *);
+static int PyArray_FillWithScalar \
+ (PyArrayObject *, PyObject *);
+static npy_bool PyArray_CheckStrides \
+ (int, int, npy_intp, npy_intp, npy_intp *, npy_intp *);
+static PyArray_Descr * PyArray_DescrNewByteorder \
+ (PyArray_Descr *, char);
+static PyObject * PyArray_IterAllButAxis \
+ (PyObject *, int *);
+static PyObject * PyArray_CheckFromAny \
+ (PyObject *, PyArray_Descr *, int, int, int, PyObject *);
+static PyObject * PyArray_FromArray \
+ (PyArrayObject *, PyArray_Descr *, int);
+static PyObject * PyArray_FromInterface \
+ (PyObject *);
+static PyObject * PyArray_FromStructInterface \
+ (PyObject *);
+static PyObject * PyArray_FromArrayAttr \
+ (PyObject *, PyArray_Descr *, PyObject *);
+static NPY_SCALARKIND PyArray_ScalarKind \
+ (int, PyArrayObject **);
+static int PyArray_CanCoerceScalar \
+ (int, int, NPY_SCALARKIND);
+static PyObject * PyArray_NewFlagsObject \
+ (PyObject *);
+static npy_bool PyArray_CanCastScalar \
+ (PyTypeObject *, PyTypeObject *);
+static int PyArray_CompareUCS4 \
+ (npy_ucs4 *, npy_ucs4 *, register size_t);
+static int PyArray_RemoveSmallest \
+ (PyArrayMultiIterObject *);
+static int PyArray_ElementStrides \
+ (PyObject *);
+static void PyArray_Item_INCREF \
+ (char *, PyArray_Descr *);
+static void PyArray_Item_XDECREF \
+ (char *, PyArray_Descr *);
+static PyObject * PyArray_FieldNames \
+ (PyObject *);
+static PyObject * PyArray_Transpose \
+ (PyArrayObject *, PyArray_Dims *);
+static PyObject * PyArray_TakeFrom \
+ (PyArrayObject *, PyObject *, int, PyArrayObject *, NPY_CLIPMODE);
+static PyObject * PyArray_PutTo \
+ (PyArrayObject *, PyObject*, PyObject *, NPY_CLIPMODE);
+static PyObject * PyArray_PutMask \
+ (PyArrayObject *, PyObject*, PyObject*);
+static PyObject * PyArray_Repeat \
+ (PyArrayObject *, PyObject *, int);
+static PyObject * PyArray_Choose \
+ (PyArrayObject *, PyObject *, PyArrayObject *, NPY_CLIPMODE);
+static int PyArray_Sort \
+ (PyArrayObject *, int, NPY_SORTKIND);
+static PyObject * PyArray_ArgSort \
+ (PyArrayObject *, int, NPY_SORTKIND);
+static PyObject * PyArray_SearchSorted \
+ (PyArrayObject *, PyObject *, NPY_SEARCHSIDE);
+static PyObject * PyArray_ArgMax \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_ArgMin \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_Reshape \
+ (PyArrayObject *, PyObject *);
+static PyObject * PyArray_Newshape \
+ (PyArrayObject *, PyArray_Dims *, NPY_ORDER);
+static PyObject * PyArray_Squeeze \
+ (PyArrayObject *);
+static PyObject * PyArray_View \
+ (PyArrayObject *, PyArray_Descr *, PyTypeObject *);
+static PyObject * PyArray_SwapAxes \
+ (PyArrayObject *, int, int);
+static PyObject * PyArray_Max \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_Min \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_Ptp \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_Mean \
+ (PyArrayObject *, int, int, PyArrayObject *);
+static PyObject * PyArray_Trace \
+ (PyArrayObject *, int, int, int, int, PyArrayObject *);
+static PyObject * PyArray_Diagonal \
+ (PyArrayObject *, int, int, int);
+static PyObject * PyArray_Clip \
+ (PyArrayObject *, PyObject *, PyObject *, PyArrayObject *);
+static PyObject * PyArray_Conjugate \
+ (PyArrayObject *, PyArrayObject *);
+static PyObject * PyArray_Nonzero \
+ (PyArrayObject *);
+static PyObject * PyArray_Std \
+ (PyArrayObject *, int, int, PyArrayObject *, int);
+static PyObject * PyArray_Sum \
+ (PyArrayObject *, int, int, PyArrayObject *);
+static PyObject * PyArray_CumSum \
+ (PyArrayObject *, int, int, PyArrayObject *);
+static PyObject * PyArray_Prod \
+ (PyArrayObject *, int, int, PyArrayObject *);
+static PyObject * PyArray_CumProd \
+ (PyArrayObject *, int, int, PyArrayObject *);
+static PyObject * PyArray_All \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_Any \
+ (PyArrayObject *, int, PyArrayObject *);
+static PyObject * PyArray_Compress \
+ (PyArrayObject *, PyObject *, int, PyArrayObject *);
+static PyObject * PyArray_Flatten \
+ (PyArrayObject *, NPY_ORDER);
+static PyObject * PyArray_Ravel \
+ (PyArrayObject *, NPY_ORDER);
+static npy_intp PyArray_MultiplyList \
+ (register npy_intp *, register int);
+static int PyArray_MultiplyIntList \
+ (register int *, register int);
+static void * PyArray_GetPtr \
+ (PyArrayObject *, register npy_intp*);
+static int PyArray_CompareLists \
+ (npy_intp *, npy_intp *, int);
+static int PyArray_AsCArray \
+ (PyObject **, void *, npy_intp *, int, PyArray_Descr*);
+static int PyArray_As1D \
+ (PyObject **, char **, int *, int);
+static int PyArray_As2D \
+ (PyObject **, char ***, int *, int *, int);
+static int PyArray_Free \
+ (PyObject *, void *);
+static int PyArray_Converter \
+ (PyObject *, PyObject **);
+static int PyArray_IntpFromSequence \
+ (PyObject *, npy_intp *, int);
+static PyObject * PyArray_Concatenate \
+ (PyObject *, int);
+static PyObject * PyArray_InnerProduct \
+ (PyObject *, PyObject *);
+static PyObject * PyArray_MatrixProduct \
+ (PyObject *, PyObject *);
+static PyObject * PyArray_CopyAndTranspose \
+ (PyObject *);
+static PyObject * PyArray_Correlate \
+ (PyObject *, PyObject *, int);
+static int PyArray_TypestrConvert \
+ (int, int);
+static int PyArray_DescrConverter \
+ (PyObject *, PyArray_Descr **);
+static int PyArray_DescrConverter2 \
+ (PyObject *, PyArray_Descr **);
+static int PyArray_IntpConverter \
+ (PyObject *, PyArray_Dims *);
+static int PyArray_BufferConverter \
+ (PyObject *, PyArray_Chunk *);
+static int PyArray_AxisConverter \
+ (PyObject *, int *);
+static int PyArray_BoolConverter \
+ (PyObject *, npy_bool *);
+static int PyArray_ByteorderConverter \
+ (PyObject *, char *);
+static int PyArray_OrderConverter \
+ (PyObject *, NPY_ORDER *);
+static unsigned char PyArray_EquivTypes \
+ (PyArray_Descr *, PyArray_Descr *);
+static PyObject * PyArray_Zeros \
+ (int, npy_intp *, PyArray_Descr *, int);
+static PyObject * PyArray_Empty \
+ (int, npy_intp *, PyArray_Descr *, int);
+static PyObject * PyArray_Where \
+ (PyObject *, PyObject *, PyObject *);
+static PyObject * PyArray_Arange \
+ (double, double, double, int);
+static PyObject * PyArray_ArangeObj \
+ (PyObject *, PyObject *, PyObject *, PyArray_Descr *);
+static int PyArray_SortkindConverter \
+ (PyObject *, NPY_SORTKIND *);
+static PyObject * PyArray_LexSort \
+ (PyObject *, int);
+static PyObject * PyArray_Round \
+ (PyArrayObject *, int, PyArrayObject *);
+static unsigned char PyArray_EquivTypenums \
+ (int, int);
+static int PyArray_RegisterDataType \
+ (PyArray_Descr *);
+static int PyArray_RegisterCastFunc \
+ (PyArray_Descr *, int, PyArray_VectorUnaryFunc *);
+static int PyArray_RegisterCanCast \
+ (PyArray_Descr *, int, NPY_SCALARKIND);
+static void PyArray_InitArrFuncs \
+ (PyArray_ArrFuncs *);
+static PyObject * PyArray_IntTupleFromIntp \
+ (int, npy_intp *);
+static int PyArray_TypeNumFromName \
+ (char *);
+static int PyArray_ClipmodeConverter \
+ (PyObject *, NPY_CLIPMODE *);
+static int PyArray_OutputConverter \
+ (PyObject *, PyArrayObject **);
+static PyObject * PyArray_BroadcastToShape \
+ (PyObject *, npy_intp *, int);
+static void _PyArray_SigintHandler \
+ (int);
+static void* _PyArray_GetSigintBuf \
+ (void);
+static int PyArray_DescrAlignConverter \
+ (PyObject *, PyArray_Descr **);
+static int PyArray_DescrAlignConverter2 \
+ (PyObject *, PyArray_Descr **);
+static int PyArray_SearchsideConverter \
+ (PyObject *, void *);
+
+#else
+
+#if defined(PY_ARRAY_UNIQUE_SYMBOL)
+#define PyArray_API PY_ARRAY_UNIQUE_SYMBOL
+#endif
+
+#if defined(NO_IMPORT) || defined(NO_IMPORT_ARRAY)
+extern void **PyArray_API;
+#else
+#if defined(PY_ARRAY_UNIQUE_SYMBOL)
+void **PyArray_API;
+#else
+static void **PyArray_API=NULL;
+#endif
+#endif
+
+#define PyArray_GetNDArrayCVersion (*(unsigned int (*)(void)) PyArray_API[0])
+#define PyBigArray_Type (*(PyTypeObject *)PyArray_API[1])
+#define PyArray_Type (*(PyTypeObject *)PyArray_API[2])
+#define PyArrayDescr_Type (*(PyTypeObject *)PyArray_API[3])
+#define PyArrayFlags_Type (*(PyTypeObject *)PyArray_API[4])
+#define PyArrayIter_Type (*(PyTypeObject *)PyArray_API[5])
+#define PyArrayMultiIter_Type (*(PyTypeObject *)PyArray_API[6])
+#define NPY_NUMUSERTYPES (*(int *)PyArray_API[7])
+#define PyBoolArrType_Type (*(PyTypeObject *)PyArray_API[8])
+#define _PyArrayScalar_BoolValues ((PyBoolScalarObject *)PyArray_API[9])
+
+#define PyGenericArrType_Type (*(PyTypeObject *)PyArray_API[10])
+#define PyNumberArrType_Type (*(PyTypeObject *)PyArray_API[11])
+#define PyIntegerArrType_Type (*(PyTypeObject *)PyArray_API[12])
+#define PySignedIntegerArrType_Type (*(PyTypeObject *)PyArray_API[13])
+#define PyUnsignedIntegerArrType_Type (*(PyTypeObject *)PyArray_API[14])
+#define PyInexactArrType_Type (*(PyTypeObject *)PyArray_API[15])
+#define PyFloatingArrType_Type (*(PyTypeObject *)PyArray_API[16])
+#define PyComplexFloatingArrType_Type (*(PyTypeObject *)PyArray_API[17])
+#define PyFlexibleArrType_Type (*(PyTypeObject *)PyArray_API[18])
+#define PyCharacterArrType_Type (*(PyTypeObject *)PyArray_API[19])
+#define PyByteArrType_Type (*(PyTypeObject *)PyArray_API[20])
+#define PyShortArrType_Type (*(PyTypeObject *)PyArray_API[21])
+#define PyIntArrType_Type (*(PyTypeObject *)PyArray_API[22])
+#define PyLongArrType_Type (*(PyTypeObject *)PyArray_API[23])
+#define PyLongLongArrType_Type (*(PyTypeObject *)PyArray_API[24])
+#define PyUByteArrType_Type (*(PyTypeObject *)PyArray_API[25])
+#define PyUShortArrType_Type (*(PyTypeObject *)PyArray_API[26])
+#define PyUIntArrType_Type (*(PyTypeObject *)PyArray_API[27])
+#define PyULongArrType_Type (*(PyTypeObject *)PyArray_API[28])
+#define PyULongLongArrType_Type (*(PyTypeObject *)PyArray_API[29])
+#define PyFloatArrType_Type (*(PyTypeObject *)PyArray_API[30])
+#define PyDoubleArrType_Type (*(PyTypeObject *)PyArray_API[31])
+#define PyLongDoubleArrType_Type (*(PyTypeObject *)PyArray_API[32])
+#define PyCFloatArrType_Type (*(PyTypeObject *)PyArray_API[33])
+#define PyCDoubleArrType_Type (*(PyTypeObject *)PyArray_API[34])
+#define PyCLongDoubleArrType_Type (*(PyTypeObject *)PyArray_API[35])
+#define PyObjectArrType_Type (*(PyTypeObject *)PyArray_API[36])
+#define PyStringArrType_Type (*(PyTypeObject *)PyArray_API[37])
+#define PyUnicodeArrType_Type (*(PyTypeObject *)PyArray_API[38])
+#define PyVoidArrType_Type (*(PyTypeObject *)PyArray_API[39])
+#define PyArray_SetNumericOps \
+ (*(int (*)(PyObject *)) \
+ PyArray_API[40])
+#define PyArray_GetNumericOps \
+ (*(PyObject * (*)(void)) \
+ PyArray_API[41])
+#define PyArray_INCREF \
+ (*(int (*)(PyArrayObject *)) \
+ PyArray_API[42])
+#define PyArray_XDECREF \
+ (*(int (*)(PyArrayObject *)) \
+ PyArray_API[43])
+#define PyArray_SetStringFunction \
+ (*(void (*)(PyObject *, int)) \
+ PyArray_API[44])
+#define PyArray_DescrFromType \
+ (*(PyArray_Descr * (*)(int)) \
+ PyArray_API[45])
+#define PyArray_TypeObjectFromType \
+ (*(PyObject * (*)(int)) \
+ PyArray_API[46])
+#define PyArray_Zero \
+ (*(char * (*)(PyArrayObject *)) \
+ PyArray_API[47])
+#define PyArray_One \
+ (*(char * (*)(PyArrayObject *)) \
+ PyArray_API[48])
+#define PyArray_CastToType \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Descr *, int)) \
+ PyArray_API[49])
+#define PyArray_CastTo \
+ (*(int (*)(PyArrayObject *, PyArrayObject *)) \
+ PyArray_API[50])
+#define PyArray_CastAnyTo \
+ (*(int (*)(PyArrayObject *, PyArrayObject *)) \
+ PyArray_API[51])
+#define PyArray_CanCastSafely \
+ (*(int (*)(int, int)) \
+ PyArray_API[52])
+#define PyArray_CanCastTo \
+ (*(npy_bool (*)(PyArray_Descr *, PyArray_Descr *)) \
+ PyArray_API[53])
+#define PyArray_ObjectType \
+ (*(int (*)(PyObject *, int)) \
+ PyArray_API[54])
+#define PyArray_DescrFromObject \
+ (*(PyArray_Descr * (*)(PyObject *, PyArray_Descr *)) \
+ PyArray_API[55])
+#define PyArray_ConvertToCommonType \
+ (*(PyArrayObject ** (*)(PyObject *, int *)) \
+ PyArray_API[56])
+#define PyArray_DescrFromScalar \
+ (*(PyArray_Descr * (*)(PyObject *)) \
+ PyArray_API[57])
+#define PyArray_DescrFromTypeObject \
+ (*(PyArray_Descr * (*)(PyObject *)) \
+ PyArray_API[58])
+#define PyArray_Size \
+ (*(npy_intp (*)(PyObject *)) \
+ PyArray_API[59])
+#define PyArray_Scalar \
+ (*(PyObject * (*)(void *, PyArray_Descr *, PyObject *)) \
+ PyArray_API[60])
+#define PyArray_FromScalar \
+ (*(PyObject * (*)(PyObject *, PyArray_Descr *)) \
+ PyArray_API[61])
+#define PyArray_ScalarAsCtype \
+ (*(void (*)(PyObject *, void *)) \
+ PyArray_API[62])
+#define PyArray_CastScalarToCtype \
+ (*(int (*)(PyObject *, void *, PyArray_Descr *)) \
+ PyArray_API[63])
+#define PyArray_CastScalarDirect \
+ (*(int (*)(PyObject *, PyArray_Descr *, void *, int)) \
+ PyArray_API[64])
+#define PyArray_ScalarFromObject \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[65])
+#define PyArray_GetCastFunc \
+ (*(PyArray_VectorUnaryFunc * (*)(PyArray_Descr *, int)) \
+ PyArray_API[66])
+#define PyArray_FromDims \
+ (*(PyObject * (*)(int, int *, int)) \
+ PyArray_API[67])
+#define PyArray_FromDimsAndDataAndDescr \
+ (*(PyObject * (*)(int, int *, PyArray_Descr *, char *)) \
+ PyArray_API[68])
+#define PyArray_FromAny \
+ (*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) \
+ PyArray_API[69])
+#define PyArray_EnsureArray \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[70])
+#define PyArray_EnsureAnyArray \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[71])
+#define PyArray_FromFile \
+ (*(PyObject * (*)(FILE *, PyArray_Descr *, npy_intp, char *)) \
+ PyArray_API[72])
+#define PyArray_FromString \
+ (*(PyObject * (*)(char *, npy_intp, PyArray_Descr *, npy_intp, char *)) \
+ PyArray_API[73])
+#define PyArray_FromBuffer \
+ (*(PyObject * (*)(PyObject *, PyArray_Descr *, npy_intp, npy_intp)) \
+ PyArray_API[74])
+#define PyArray_FromIter \
+ (*(PyObject * (*)(PyObject *, PyArray_Descr *, npy_intp)) \
+ PyArray_API[75])
+#define PyArray_Return \
+ (*(PyObject * (*)(PyArrayObject *)) \
+ PyArray_API[76])
+#define PyArray_GetField \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Descr *, int)) \
+ PyArray_API[77])
+#define PyArray_SetField \
+ (*(int (*)(PyArrayObject *, PyArray_Descr *, int, PyObject *)) \
+ PyArray_API[78])
+#define PyArray_Byteswap \
+ (*(PyObject * (*)(PyArrayObject *, npy_bool)) \
+ PyArray_API[79])
+#define PyArray_Resize \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Dims *, int, NPY_ORDER)) \
+ PyArray_API[80])
+#define PyArray_MoveInto \
+ (*(int (*)(PyArrayObject *, PyArrayObject *)) \
+ PyArray_API[81])
+#define PyArray_CopyInto \
+ (*(int (*)(PyArrayObject *, PyArrayObject *)) \
+ PyArray_API[82])
+#define PyArray_CopyAnyInto \
+ (*(int (*)(PyArrayObject *, PyArrayObject *)) \
+ PyArray_API[83])
+#define PyArray_CopyObject \
+ (*(int (*)(PyArrayObject *, PyObject *)) \
+ PyArray_API[84])
+#define PyArray_NewCopy \
+ (*(PyObject * (*)(PyArrayObject *, NPY_ORDER)) \
+ PyArray_API[85])
+#define PyArray_ToList \
+ (*(PyObject * (*)(PyArrayObject *)) \
+ PyArray_API[86])
+#define PyArray_ToString \
+ (*(PyObject * (*)(PyArrayObject *, NPY_ORDER)) \
+ PyArray_API[87])
+#define PyArray_ToFile \
+ (*(int (*)(PyArrayObject *, FILE *, char *, char *)) \
+ PyArray_API[88])
+#define PyArray_Dump \
+ (*(int (*)(PyObject *, PyObject *, int)) \
+ PyArray_API[89])
+#define PyArray_Dumps \
+ (*(PyObject * (*)(PyObject *, int)) \
+ PyArray_API[90])
+#define PyArray_ValidType \
+ (*(int (*)(int)) \
+ PyArray_API[91])
+#define PyArray_UpdateFlags \
+ (*(void (*)(PyArrayObject *, int)) \
+ PyArray_API[92])
+#define PyArray_New \
+ (*(PyObject * (*)(PyTypeObject *, int, npy_intp *, int, npy_intp *, void *, int, int, PyObject *)) \
+ PyArray_API[93])
+#define PyArray_NewFromDescr \
+ (*(PyObject * (*)(PyTypeObject *, PyArray_Descr *, int, npy_intp *, npy_intp *, void *, int, PyObject *)) \
+ PyArray_API[94])
+#define PyArray_DescrNew \
+ (*(PyArray_Descr * (*)(PyArray_Descr *)) \
+ PyArray_API[95])
+#define PyArray_DescrNewFromType \
+ (*(PyArray_Descr * (*)(int)) \
+ PyArray_API[96])
+#define PyArray_GetPriority \
+ (*(double (*)(PyObject *, double)) \
+ PyArray_API[97])
+#define PyArray_IterNew \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[98])
+#define PyArray_MultiIterNew \
+ (*(PyObject * (*)(int, ...)) \
+ PyArray_API[99])
+#define PyArray_PyIntAsInt \
+ (*(int (*)(PyObject *)) \
+ PyArray_API[100])
+#define PyArray_PyIntAsIntp \
+ (*(npy_intp (*)(PyObject *)) \
+ PyArray_API[101])
+#define PyArray_Broadcast \
+ (*(int (*)(PyArrayMultiIterObject *)) \
+ PyArray_API[102])
+#define PyArray_FillObjectArray \
+ (*(void (*)(PyArrayObject *, PyObject *)) \
+ PyArray_API[103])
+#define PyArray_FillWithScalar \
+ (*(int (*)(PyArrayObject *, PyObject *)) \
+ PyArray_API[104])
+#define PyArray_CheckStrides \
+ (*(npy_bool (*)(int, int, npy_intp, npy_intp, npy_intp *, npy_intp *)) \
+ PyArray_API[105])
+#define PyArray_DescrNewByteorder \
+ (*(PyArray_Descr * (*)(PyArray_Descr *, char)) \
+ PyArray_API[106])
+#define PyArray_IterAllButAxis \
+ (*(PyObject * (*)(PyObject *, int *)) \
+ PyArray_API[107])
+#define PyArray_CheckFromAny \
+ (*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) \
+ PyArray_API[108])
+#define PyArray_FromArray \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Descr *, int)) \
+ PyArray_API[109])
+#define PyArray_FromInterface \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[110])
+#define PyArray_FromStructInterface \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[111])
+#define PyArray_FromArrayAttr \
+ (*(PyObject * (*)(PyObject *, PyArray_Descr *, PyObject *)) \
+ PyArray_API[112])
+#define PyArray_ScalarKind \
+ (*(NPY_SCALARKIND (*)(int, PyArrayObject **)) \
+ PyArray_API[113])
+#define PyArray_CanCoerceScalar \
+ (*(int (*)(int, int, NPY_SCALARKIND)) \
+ PyArray_API[114])
+#define PyArray_NewFlagsObject \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[115])
+#define PyArray_CanCastScalar \
+ (*(npy_bool (*)(PyTypeObject *, PyTypeObject *)) \
+ PyArray_API[116])
+#define PyArray_CompareUCS4 \
+ (*(int (*)(npy_ucs4 *, npy_ucs4 *, register size_t)) \
+ PyArray_API[117])
+#define PyArray_RemoveSmallest \
+ (*(int (*)(PyArrayMultiIterObject *)) \
+ PyArray_API[118])
+#define PyArray_ElementStrides \
+ (*(int (*)(PyObject *)) \
+ PyArray_API[119])
+#define PyArray_Item_INCREF \
+ (*(void (*)(char *, PyArray_Descr *)) \
+ PyArray_API[120])
+#define PyArray_Item_XDECREF \
+ (*(void (*)(char *, PyArray_Descr *)) \
+ PyArray_API[121])
+#define PyArray_FieldNames \
+ (*(PyObject * (*)(PyObject *)) \
+ PyArray_API[122])
+#define PyArray_Transpose \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Dims *)) \
+ PyArray_API[123])
+#define PyArray_TakeFrom \
+ (*(PyObject * (*)(PyArrayObject *, PyObject *, int, PyArrayObject *, NPY_CLIPMODE)) \
+ PyArray_API[124])
+#define PyArray_PutTo \
+ (*(PyObject * (*)(PyArrayObject *, PyObject*, PyObject *, NPY_CLIPMODE)) \
+ PyArray_API[125])
+#define PyArray_PutMask \
+ (*(PyObject * (*)(PyArrayObject *, PyObject*, PyObject*)) \
+ PyArray_API[126])
+#define PyArray_Repeat \
+ (*(PyObject * (*)(PyArrayObject *, PyObject *, int)) \
+ PyArray_API[127])
+#define PyArray_Choose \
+ (*(PyObject * (*)(PyArrayObject *, PyObject *, PyArrayObject *, NPY_CLIPMODE)) \
+ PyArray_API[128])
+#define PyArray_Sort \
+ (*(int (*)(PyArrayObject *, int, NPY_SORTKIND)) \
+ PyArray_API[129])
+#define PyArray_ArgSort \
+ (*(PyObject * (*)(PyArrayObject *, int, NPY_SORTKIND)) \
+ PyArray_API[130])
+#define PyArray_SearchSorted \
+ (*(PyObject * (*)(PyArrayObject *, PyObject *, NPY_SEARCHSIDE)) \
+ PyArray_API[131])
+#define PyArray_ArgMax \
+ (*(PyObject * (*)(PyArrayObject *, int, PyArrayObject *)) \
+ PyArray_API[132])
+#define PyArray_ArgMin \
+ (*(PyObject * (*)(PyArrayObject *, int, PyArrayObject *)) \
+ PyArray_API[133])
+#define PyArray_Reshape \
+ (*(PyObject * (*)(PyArrayObject *, PyObject *)) \
+ PyArray_API[134])
+#define PyArray_Newshape \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Dims *, NPY_ORDER)) \
+ PyArray_API[135])
+#define PyArray_Squeeze \
+ (*(PyObject * (*)(PyArrayObject *)) \
+ PyArray_API[136])
+#define PyArray_View \
+ (*(PyObject * (*)(PyArrayObject *, PyArray_Descr *, PyTypeObject *)) \
+ PyArray_API[137])
+#define PyArray_SwapAxes \
+ (*(PyObject * (*)(PyArrayObject *, int, int)) \
+ PyArray_API[138])
+#define PyArray_Max \
+ (*(PyObject * (*)(PyArrayObject *, int, PyArrayObject *)) \
+ PyArray_API[139])
+#define PyArray_Min \
+ (*(PyObject * (*)(PyArrayObject *, int, PyArrayObject *)) \
+ PyArray_API[140])
+#define PyArray_Ptp \
+ (*(PyObject * (*)(PyArrayObject *, int, PyArrayObject *)) \
+ PyArray_API[141])
+#define PyArray_Mean \
+ (*(PyObject * (*)(PyArrayObject *, int, int, PyArrayObject *)) \
+ PyArray_API[142])
+#define PyArray_Trace \
+ (*(PyObject * (*)(PyArrayObject *, int, int, int, int, PyArrayObject *)) \
+ PyArray_API[143])
+#define PyArray_Diagonal \
+ (*(PyObject * (*)(PyArrayObject *, int, int, int)) \
+ PyArray_API[144])
+#define PyArray_Clip \
+ (*(PyObject * (*)(PyArrayObject *, PyObject *, PyObject *, PyArrayObject *)) \
+ PyArray_API[145])
+#define PyArray_Conjugate \
+ (*(PyObject * (*)(PyArrayObject *, PyArrayObject *)) \
+ PyArr