[Scipy-svn] r3511 - in trunk/scipy/ndimage/segment: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Nov 9 15:53:58 EST 2007
Author: tom.waite
Date: 2007-11-09 14:53:50 -0600 (Fri, 09 Nov 2007)
New Revision: 3511
Modified:
trunk/scipy/ndimage/segment/Segmenter_IMPL.c
trunk/scipy/ndimage/segment/tests/test_segment.py
Log:
Added doc string to test_segment.py. clean up column length in Segmenter_IMPL.c
Modified: trunk/scipy/ndimage/segment/Segmenter_IMPL.c
===================================================================
--- trunk/scipy/ndimage/segment/Segmenter_IMPL.c 2007-11-08 03:42:40 UTC (rev 3510)
+++ trunk/scipy/ndimage/segment/Segmenter_IMPL.c 2007-11-09 20:53:50 UTC (rev 3511)
@@ -10,7 +10,8 @@
#define FALSE 0
#define TRUE 1
-int NI_GetObjectStats(int rows, int cols, int numberObjects, unsigned short *labeledEdges, objStruct objectMetrics[]){
+int NI_GetObjectStats(int rows, int cols, int numberObjects, unsigned short *labeledEdges,
+ objStruct objectMetrics[]){
int i, j, k, m;
int offset;
@@ -47,7 +48,7 @@
}
offset += cols;
}
- // the bounding box for the 2D blob
+ /* the bounding box for the 2D blob */
objectMetrics[k-1].L = LowX;
objectMetrics[k-1].R = HighX;
objectMetrics[k-1].B = LowY;
@@ -75,11 +76,11 @@
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;
@@ -96,7 +97,7 @@
kernel[j++] = t4;
}
- // normalize the kernel so unity gain (as is LP filter this is easy)
+ /* 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];
@@ -112,7 +113,8 @@
return;
}
-void filter2D(int HalfFilterTaps, int rows, int cols, int lowThreshold, int highThreshold, float *kernel, double *Image){
+void filter2D(int HalfFilterTaps, int rows, int cols, int lowThreshold, int highThreshold,
+ float *kernel, double *Image){
int i, j, k, n, num1;
int offset;
@@ -122,11 +124,11 @@
num1 = HalfFilterTaps + 1;
offset = 0;
for(i = 0; i < rows; ++i){
- // copy image row to local buffer
+ /* 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
+ /* constant pad the ends of the buffer */
for(j = 0; j < num1; ++j){
buffer[j] = buffer[num1];
}
@@ -134,7 +136,7 @@
buffer[j] = buffer[cols-1+num1];
}
- // Perform Symmetric Convolution in the X dimension.
+ /* 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){
@@ -147,13 +149,13 @@
offset = 0;
for(i = 0; i < cols; ++i){
- // copy image column to local buffer
+ /* 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
+ /* constant pad the ends of the buffer */
for(j = 0; j < num1; ++j){
buffer[j] = buffer[num1];
}
@@ -161,7 +163,7 @@
buffer[j] = buffer[rows-1+num1];
}
- // Perform Symmetric Convolution in the Y dimension.
+ /* Perform Symmetric Convolution in the Y dimension. */
offset = 0;
for(j = num1; j < (rows+num1); ++j){
sum = buffer[j] * kernel[num1];
@@ -173,7 +175,7 @@
}
}
- // threshold the image
+ /* threshold the image */
offset = 0;
for(i = 0; i < rows; ++i){
for(j = 0; j < cols; ++j){
@@ -189,13 +191,14 @@
}
-void doPreProcess(int samples, int rows, int cols, double *rawImage, double BPHigh, int apearture, int lowThreshold, int highThreshold){
+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;
@@ -221,9 +224,9 @@
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){
@@ -237,9 +240,9 @@
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){
@@ -267,9 +270,9 @@
}
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){
@@ -298,7 +301,7 @@
offset -= cols;
}
if(!Change) break;
- } // end while loop
+ } /* end while loop */
Classes[0] = 0;
Label = 1;
@@ -314,7 +317,7 @@
if(NewLabel){
Classes[Label++] = m;
if(Label > 4000){
- return 0; // too many labeled regions. this is a pathology in the image slice
+ return 0; /* too many labeled regions. this is a pathology */
}
}
}
@@ -322,9 +325,9 @@
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){
@@ -349,7 +352,8 @@
return (float)sqrt(X*X + Y*Y);
}
-int traceEdge(int i, int j, int rows, int cols, double cannyLow, float *magImage, float *HYSImage){
+int traceEdge(int i, int j, int rows, int cols, double cannyLow, float *magImage,
+ float *HYSImage){
int n, m;
int ptr;
@@ -357,9 +361,9 @@
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){
@@ -368,9 +372,9 @@
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;
@@ -388,7 +392,8 @@
}
-void edgeThreshold(int rows, int cols, double cannyLow, float *magImage, float *HYSImage){
+void edgeThreshold(int rows, int cols, double cannyLow, float *magImage,
+ float *HYSImage){
int i, j;
int ptr;
@@ -406,7 +411,8 @@
}
-void edgeHysteresis(int rows, int cols, double cannyLow, double cannyHigh, float *magImage, float *HYSImage){
+void edgeHysteresis(int rows, int cols, double cannyLow, double cannyHigh,
+ float *magImage, float *HYSImage){
int i, j;
int ptr;
@@ -424,8 +430,9 @@
}
-void nonMaxSupress(int rows, int cols, float aveXValue, float aveYValue, double *cannyLow, double *cannyHigh,
- int mode, float *hDGImage, float *vDGImage, float *magImage){
+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;
@@ -452,7 +459,7 @@
if((fabs(xC) < aveXValue) && (fabs(yC) < aveYValue)) continue;
G = magnitude(xC, yC);
if(fabs(yC) > fabs(xC)){
- // vertical gradient
+ /* vertical gradient */
xSlope = (float)(fabs(xC) / fabs(yC));
ySlope = (float)1.0;
G2 = magnitude(hDGImage[ptr_m1+j], vDGImage[ptr_m1+j]);
@@ -467,7 +474,7 @@
}
}
else{
- // horizontal gradient
+ /* horizontal gradient */
xSlope = (float)(fabs(yC) / fabs(xC));
ySlope = (float)1.0;
G2 = magnitude(hDGImage[ptr+j+1], vDGImage[ptr+j+1]);
@@ -481,7 +488,7 @@
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)) ){
+ 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];
@@ -517,9 +524,9 @@
++ptr;
}
}
- //
+ /*
// now get the max after skipping the low values
- //
+ */
mValue = -1;
mIndex = 0;
for(i = 10; i < 256; ++i){
@@ -530,12 +537,12 @@
}
if(mode == 1){
- // based on the mean value of edge energy
+ /* based on the mean value of edge energy */
*cannyLow = ((*cannyLow) * tAve);
*cannyHigh = ((*cannyHigh) * tAve);
}
else{
- // based on the mode value of edge energy
+ /* based on the mode value of edge energy */
*cannyLow = ((*cannyLow) * ((float)mIndex/step));
*cannyHigh = ((*cannyHigh) * ((float)mIndex/step));
}
@@ -548,9 +555,9 @@
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;
@@ -564,9 +571,9 @@
mLength = MAX(rows, cols) + 64;
tBuffer = calloc(mLength, sizeof(float));
- //
+ /*
// filter X
- //
+ */
count = 0;
for(i = 0; i < rows; ++i){
ptr = i * cols;
@@ -585,11 +592,11 @@
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
+ /* 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){
@@ -612,7 +619,7 @@
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
+ /* this is 50% of the max, hardwirred for now, and is part of the threshold */
}
free(tBuffer);
@@ -622,8 +629,10 @@
}
-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,
+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;
@@ -642,14 +651,12 @@
float *magImage = NULL;
float *tBuffer = NULL;
- // filter
- printf("do preProcess\n");
+ /* filter */
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));
@@ -657,10 +664,10 @@
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;
@@ -671,8 +678,10 @@
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);
+ 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);
}
@@ -680,17 +689,17 @@
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){
@@ -713,7 +722,8 @@
}
-void doSobel(int samples, int rows, int cols, double sobelLow, int mode, double *rawImage, unsigned short *edgeImage){
+void doSobel(int samples, int rows, int cols, double sobelLow, int mode,
+ double *rawImage, unsigned short *edgeImage){
int i, j;
int p, m, n;
@@ -745,9 +755,9 @@
offset += cols;
}
- //
+ /*
// Sobel
- //
+ */
offset = cols;
for(i = 1; i < rows-1; ++i){
offsetM1 = offset - cols;
@@ -769,7 +779,7 @@
offset += cols;
}
- // threshold based on ave
+ /* threshold based on ave */
pAve /= count;
scale = 1.0 / maxValue;
@@ -785,9 +795,9 @@
}
offset += cols;
}
- //
+ /*
// now get the max after skipping the low values
- //
+ */
maxValue = -1;
maxIndex = 0;
for(i = 10; i < 256; ++i){
@@ -798,11 +808,11 @@
}
if(mode == 1){
- // based on the mean value of edge energy
+ /* based on the mean value of edge energy */
pThreshold = (int)(sobelLow * (float)pAve);
}
else{
- // based on the mode value of edge energy
+ /* based on the mode value of edge energy */
pThreshold = (sobelLow * (minValue + ((float)maxIndex/step)));
}
@@ -827,7 +837,8 @@
}
-void estimateThreshold(float *lowThreshold, float *highThreshold, float ShenCastanLow, int rows, int cols, float *SourceImage){
+void estimateThreshold(float *lowThreshold, float *highThreshold, float ShenCastanLow,
+ int rows, int cols, float *SourceImage){
int i, j;
int offset;
@@ -862,9 +873,9 @@
offset += cols;
}
- //
+ /*
// now get the edge energy mode
- //
+ */
value = 0;
mIndex = 10;
for(i = 10; i < 256; ++i){
@@ -884,16 +895,17 @@
}
-void thresholdEdges(float *SourceImage, unsigned short *EdgeImage, double ShenCastanLow, int rows, int cols){
+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;
@@ -913,7 +925,8 @@
}
-float adaptiveGradient(float *BLImage, float *FilterImage, int nrow, int ncol, int cols, int window){
+float adaptiveGradient(float *BLImage, float *FilterImage, int nrow, int ncol,
+ int cols, int window){
int i, j;
int offset;
@@ -957,7 +970,8 @@
}
-void getZeroCrossings(float *SourceImage, float *FilterImage, float *BLImage, int rows, int cols, int window){
+void getZeroCrossings(float *SourceImage, float *FilterImage, float *BLImage,
+ int rows, int cols, int window){
int i, j;
int offset;
@@ -996,7 +1010,7 @@
}
}
if(validEdge){
- // adaptive gradeint is signed
+ /* adaptive gradeint is signed */
SourceImage[offset+j] = (float)fabs(adaptiveGradient(BLImage, FilterImage, i, j, cols, window));
}
}
@@ -1014,9 +1028,9 @@
int offset;
float t;
- //
+ /*
// like an unsharp mask
- //
+ */
offset = 0;
for(i = 0; i < rows; ++i){
for(j = 0; j < cols; ++j){
@@ -1060,7 +1074,8 @@
}
-void ISEF_Vertical(float *SourceImage, float *FilterImage, float *A, float *B, int rows, int cols, double b){
+void ISEF_Vertical(float *SourceImage, float *FilterImage, float *A, float *B,
+ int rows, int cols, double b){
int i, j;
@@ -1070,40 +1085,40 @@
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
+ /* process row 0 */
A[i] = b1 * SourceImage[i];
- // process row N-1
+ /* 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;
@@ -1114,9 +1129,9 @@
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){
@@ -1129,12 +1144,13 @@
}
-void ISEF_Horizontal(float *SourceImage, float *FilterImage, float *A, float *B, int rows, int cols, double b){
+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;
@@ -1143,9 +1159,9 @@
b1 = ((float)1.0 - b)/((float)1.0 + b);
b2 = b * b1;
- //
+ /*
// columns boundaries
- //
+ */
offset = 0;
for(i = 0; i < rows; ++i){
// col 0
@@ -1154,9 +1170,9 @@
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){
@@ -1165,9 +1181,9 @@
offset += cols;
}
- //
+ /*
// anti-causal IIR part
- //
+ */
offset = 0;
for(j = cols-2; j > 0; --j){
for(i = 0; i < rows; ++i){
@@ -1176,17 +1192,17 @@
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];
@@ -1218,7 +1234,8 @@
}
-void Shen_Castan(double b, double ShenCastanLow, int rows, int cols, int window, int lowThreshold, int highThreshold,
+void Shen_Castan(double b, double ShenCastanLow, int rows, int cols, int window,
+ int lowThreshold, int highThreshold,
double *RawImage, unsigned short *EdgeImage){
int i;
@@ -1235,10 +1252,10 @@
SourceImage[i] = RawImage[i];
}
computeISEF(SourceImage, FilterImage, rows, cols, b);
- // optional thresholding based on low, high
+ /* 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
+ /* the new source image is now the adaptive gradient */
getZeroCrossings(SourceImage, FilterImage, BinaryLaplacianImage, rows, cols, window);
thresholdEdges(SourceImage, EdgeImage, ShenCastanLow, rows, cols);
@@ -1250,8 +1267,9 @@
}
-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 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;
@@ -1281,7 +1299,8 @@
}
-void buildBinaryImage(int rows, int cols, double *rawImage, unsigned short *edgeImage, int lowThreshold, int highThreshold){
+void buildBinaryImage(int rows, int cols, double *rawImage, unsigned short *edgeImage,
+ int lowThreshold, int highThreshold){
int i, j;
int offset;
@@ -1306,7 +1325,8 @@
-void morphoFilterBinaryImage(int rows, int cols, unsigned short *edgeImage, int CloseSize, int OpenSize){
+void morphoFilterBinaryImage(int rows, int cols, unsigned short *edgeImage,
+ int CloseSize, int OpenSize){
int i, j;
@@ -1315,8 +1335,8 @@
unsigned short omask[11][11];
int olapValuesC[4];
int olapValuesO[4];
- int CloseMaskSize=1;
- int OpenMaskSize=1;
+ int CloseMaskSize = 1;
+ int OpenMaskSize = 1;
int LowValue1, HighValue1;
int LowValue2, HighValue2;
int spadSize;
@@ -1348,9 +1368,9 @@
olapValuesC[3] = HighValue2;
}
- //
+ /*
// Open filter
- //
+ */
if(OpenSize){
OpenMaskSize = (OpenSize-1)/2;
for(i = 0; i < 2*OpenMaskSize+1; ++i){
@@ -1391,11 +1411,11 @@
for(i = 0; i < rows; ++i){
for(j = 0; j < cols; ++j){
if(ImageE[offset2+j] == 1){
- // this will activate some original off-pixels
+ /* this will activate some original off-pixels */
edgeImage[offset+j] = 1;
}
else{
- // this will zero some original on-pixels
+ /* this will zero some original on-pixels */
edgeImage[offset+j] = 0;
}
}
@@ -1410,7 +1430,8 @@
}
-void doRegionGrow(int samples, int rows, int cols, double *rawImage, unsigned short *edgeImage, int lowThreshold,
+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);
@@ -1420,14 +1441,16 @@
}
-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 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);
+ doRegionGrow(samples, rows, cols, rawImage, edgeImage, lowThreshold,
+ highThreshold, closeWindow, openWindow);
*groups = ConnectedEdgePoints(rows, cols, edgeImage);
//
@@ -1448,7 +1471,8 @@
}
-int NI_SobelEdges(int samples, int rows, int cols, double sobelLow, int mode, int lowThreshold, int highThreshold, double BPHigh,
+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){
@@ -1461,9 +1485,9 @@
*groups = ConnectedEdgePoints(rows, cols, edgeImage);
- //
+ /*
// prune the isolated pixels
- //
+ */
offset = 0;
for(i = 0; i < rows; ++i){
for(j = 0; j < cols; ++j){
@@ -1611,7 +1635,7 @@
nloop = 0;
while(1){
- // erode
+ /* erode */
Column = 0;
for(n = 0; n < 8; ++n){
for(i = 0; i < 3; ++i){
@@ -1644,7 +1668,7 @@
Offset += spadSize;
}
- // dialate
+ /* dialate */
Offset = 0;
for(i = 0; i < N; ++i){
for(j = 0; j < M; ++j){
@@ -1671,7 +1695,7 @@
Offset += spadSize;
}
- // form the HMT
+ /* form the HMT */
Offset = 0;
for(i = 0; i < N; ++i){
for(j = 0; j < M; ++j){
@@ -1681,7 +1705,7 @@
Offset += spadSize;
}
- // Thin for stage n
+ /* Thin for stage n */
Offset = 0;
for(i = 0; i < N; ++i){
@@ -1701,7 +1725,7 @@
}
}
- // check for the NULL set
+ /* check for no change */
hit = 0;
Offset = 0;
for(i = 0; i < N; ++i){
@@ -1730,7 +1754,8 @@
}
-int NI_ThinFilter(int samples, int rows, int cols, int numberObjects, unsigned short *edgeImage, objStruct objectMetrics[]){
+int NI_ThinFilter(int samples, int rows, int cols, int numberObjects,
+ unsigned short *edgeImage, objStruct objectMetrics[]){
int i, j;
int loop;
@@ -1752,9 +1777,9 @@
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));
@@ -1773,9 +1798,9 @@
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){
@@ -1784,9 +1809,9 @@
srcOffset += cols;
}
- //
+ /*
// copy the ROI for MAT (medial axis transformation) filter
- //
+ */
dstOffset = inflate*rows;
for(i = bottom; i < top; ++i){
srcOffset = i*cols;
@@ -1797,11 +1822,12 @@
}
dstOffset += cols;
}
- ThinningFilter(roiRows, roiCols, cols, J_mask, K_mask, Input, CInput, ErosionStage, DialationStage, HMT, Copy);
+ 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;
@@ -1814,10 +1840,10 @@
}
}
- //
+ /*
// 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];
}
@@ -1839,11 +1865,11 @@
void generateMask(unsigned char *ImageH, bPOINT *boundary, int newSamples, int label, int cols){
- //
+ /*
// get the boundary point pairs (left, right) for each line
// if there is no pair, then the boundary is open
// then fill the image in with the current label
- //
+ */
int i, j, k, m;
int list[2048];
@@ -1873,7 +1899,7 @@
}
}
}
- // now get the closest boundary
+ /* now get the closest boundary */
if(k){
distance = maxDistance;
index = -1;
@@ -1899,9 +1925,9 @@
low = boundary[index].x;
high = boundary[i].x;
}
- //
+ /*
// do the fill
- //
+ */
offset = y * cols;
for(j = low; j <= high; ++j){
ImageH[offset+j] = label;
@@ -1909,7 +1935,7 @@
}
}
else{
- // boundary point is isolated
+ /* boundary point is isolated */
boundary[i].linkIndex = i;
}
}
@@ -1919,7 +1945,8 @@
}
-void getBoundaryMetrics(bPOINT *boundary, float *length, float *minRadius, float *maxRadius, float *aveRadius,
+void getBoundaryMetrics(bPOINT *boundary, float *length, float *minRadius,
+ float *maxRadius, float *aveRadius,
float Xcenter, float Ycenter, int newSamples){
int j;
@@ -1990,11 +2017,6 @@
p = 1;
while(p < mcount){
offset = (CurI-inflate)*spadSize;
- if(offset < 0){
- printf("offset < 0 ");
- printf("CurI [%d]. p [%d]. mcount [%d]\n", CurI, p, mcount);
- getchar();
- }
MinD = 1024;
NewI = -1;
NewJ = -1;
@@ -2002,7 +2024,7 @@
for(j = CurJ-inflate; j < CurJ+inflate; ++j){
m = Input[offset+j];
if(m == 1){
- // city block distance
+ /* city block distance */
k = abs(i-CurI) + abs(j-CurJ);
if(k < MinD){
MinD = k;
@@ -2031,10 +2053,10 @@
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;
@@ -2045,8 +2067,8 @@
LowValue2 = olapValues[2];
HighValue2 = olapValues[3];
- // close - step 1 is dialate
- // open - step 1 is erode
+ /* 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){
@@ -2065,8 +2087,8 @@
offset += spadSize;
}
- // close - step 2 is erode
- // open - step 2 is dialate
+ /* 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){
@@ -2088,7 +2110,8 @@
return;
}
-void getCompactness(unsigned char *Input, RECT roi, int label, int spadSize, float *vCompactness, float length){
+void getCompactness(unsigned char *Input, RECT roi, int label, int spadSize,
+ float *vCompactness, float length){
int i, j;
int maskOffset;
@@ -2115,8 +2138,9 @@
}
-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],
+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;
@@ -2133,9 +2157,9 @@
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;
@@ -2147,20 +2171,20 @@
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){
@@ -2178,8 +2202,10 @@
}
-void getBoundary(unsigned short *ThinEdgeImage, unsigned char *Input, blobBoundary *pBoundary, blobBoundary *lBoundary,
- boundaryIndex *pBoundaryIndex, RECT boundBox, int label, int bBox, int nextSlot, int memOffset,
+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;
@@ -2202,7 +2228,7 @@
Input[i] = 0;
}
- //copy to spad
+ /* copy to spad */
count = 0;
rows = boundBox.top-boundBox.bottom+2*inflate;
@@ -2227,7 +2253,7 @@
if(Input[srcOffset+j]){
if(first){
first = FALSE;
- // index of the seed sample
+ /* index of the seed sample */
value.xy.x = i;
value.xy.y = j;
}
@@ -2243,12 +2269,10 @@
for(i = 0; i < mcount; ++i){
value.xy.x = lBoundary[i].xy.x + boundBox.left - inflate;
value.xy.y = lBoundary[i].xy.y + boundBox.bottom - inflate + 1;
- //printf("[%d, %d]\n", value.xy.x, value.xy.y);
pBoundary[memOffset].xy.x = value.xy.x;
pBoundary[memOffset].xy.y = value.xy.y;
++memOffset;
}
- //getchar();
return;
@@ -2256,7 +2280,7 @@
void buildBoundary(objStruct objectMetrics[], int searchWindow, unsigned short *ThinEdgeImage,
- int numberObjects, int srcRows, int srcCols){
+ int numberObjects, int srcRows, int srcCols){
int i, j, k;
int count;
@@ -2267,7 +2291,7 @@
int end;
int label;
int distance;
- // these should be setup parameters
+ /* these will be user-setup parameters */
int closureDistance = 12;
int CloseSize = 5;
int OpenSize = 5;
@@ -2281,7 +2305,7 @@
float maxRadius;
float aveRadius;
float vCompactness;
- // for morphological close of mask. max structuring element is 11x11
+ /* for morphological close of mask. max structuring element is 11x11 */
unsigned short cmask[11][11];
unsigned short omask[11][11];
int olapValuesC[4];
@@ -2301,9 +2325,9 @@
unsigned char *ImageC;
unsigned char *ImageH;
- //
+ /*
// Close filter
- //
+ */
CloseMaskSize = (CloseSize-1)/2;
for(i = 0; i < 2*CloseMaskSize+1; ++i){
for(j = 0; j < 2*CloseMaskSize+1; ++j){
@@ -2319,9 +2343,9 @@
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){
@@ -2373,21 +2397,16 @@
bBox.bottom = objectMetrics[i].B;
label = objectMetrics[i].Label;
pBoundaryIndex[i+1].Label = label;
- //printf("(%d, %d, %d, %d [%d])\n", bBox.left, bBox.right, bBox.top, bBox.bottom, label);
getBoundary(ThinEdgeImage, Input, pBoundary, lBoundary, pBoundaryIndex, bBox, label,
i, pBoundaryIndex[0].numberPoints, count, spadSize, searchWindow);
}
- //
+ /*
// Input will now be used in the fill. Copy the labeled edge image
- //
+ */
- //
- // numBoundaries = numberObjects
- //
offset = 0;
numBoundaries = pBoundaryIndex[0].numberPoints;
- //printf("numBoundaries [%d]\n", numBoundaries);
for(i = 0; i < numBoundaries; ++i){
numSamples = pBoundaryIndex[i+1].numberPoints;
end = numSamples-2;
@@ -2397,9 +2416,9 @@
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){
@@ -2410,7 +2429,6 @@
break;
}
}
- //printf("[%d] newSamples [%d]\n", i, newSamples);
distance = abs(boundary[0].x-boundary[end-2].x) + abs(boundary[0].y-boundary[end-2].y);
pBoundaryIndex[i+1].curveClose = distance;
@@ -2429,7 +2447,8 @@
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);
+ (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;
@@ -2441,35 +2460,35 @@
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;
@@ -2509,9 +2528,9 @@
}
}
- //
+ /*
// 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;
@@ -2525,7 +2544,7 @@
}
// debug only
- if(1){
+ if(0){
for(i = 0; i < numBoundaries; ++i){
if(pBoundaryIndex[i+1].boundaryLength != (float)0.0){
printf("boundary %d:\n", i);
@@ -2556,9 +2575,9 @@
}
}
- //
+ /*
// need to return input which is now mask image
- //
+ */
offset = 0;
offset2 = 0;
@@ -2613,7 +2632,7 @@
lawsFilter->lawsKernel[5][i] = O7[i];
}
- // DC filter is unity gain
+ /* L filter is unity gain */
sum = (float)0.0;
for(i = 0; i < 7; ++i){
sum += lawsFilter->lawsKernel[0][i];
@@ -2633,7 +2652,7 @@
float result[7];
float sum;
- // filter rows
+ /* filter rows */
for(i = 0; i < kernelSize; ++i){
sum = (float)0.0;
offset = i * kernelSize;
@@ -2643,7 +2662,7 @@
result[i] = sum;
}
- //filter columns
+ /* filter columns */
sum = (float)0.0;
for(j = 0; j < kernelSize; ++j){
sum += (rowFilter[j]*result[j]);
@@ -2654,8 +2673,10 @@
}
-void getLawsTexture(LawsFilter7 lawsFilter, tTEM LawsFeatures[], objStruct objectMetrics[], double *sourceImage,
- unsigned short *MaskImage, int numberObjects, int srcRows, int srcCols){
+void getLawsTexture(LawsFilter7 lawsFilter, tTEM LawsFeatures[],
+ objStruct objectMetrics[], double *sourceImage,
+ unsigned short *MaskImage, int numberObjects,
+ int srcRows, int srcCols){
int i, j;
int label;
@@ -2676,15 +2697,15 @@
bBox.bottom = objectMetrics[i].B;
label = objectMetrics[i].Label;
if(objectMetrics[i].voxelMean != (float)0.0){
- //
+ /*
// valid size region
- //
+ */
computeLaws(lawsFilter, LawsFeatures, bBox, label, aperature, srcRows, srcCols, ImageH, ImageT,
MaskImage, lawsImage, sourceImage);
for(j = 1; j < lawsFilter.numberFilterLayers; ++j){
objectMetrics[i].TEM[j-1] = LawsFeatures[j].Variance;
}
- /*
+ /* -- later will need to return a view of the texture images
int index;
int offset;
int layerStep = srcRows*srcCols;
@@ -2711,12 +2732,14 @@
}
-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){
+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;
@@ -2744,7 +2767,7 @@
char dual[24];
- // zero the laws mask memory first
+ /* zero the laws mask memory first */
for(i = 0; i < srcRows*srcCols; ++i){
ImageH[i] = 0;
}
@@ -2760,9 +2783,9 @@
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;
@@ -2774,22 +2797,22 @@
}
}
if(count == fullMask){
- //
- // full mask. 100% coverage. now do the Law's texture filters
- //
+ /*
+ // 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.
+ /* 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;
}
@@ -2799,17 +2822,20 @@
}
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){
+ */
+ 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);
+ lawsImage[lawsLayer*layerStep + i*srcCols + j] =
+ (filterResult1 / lawsLL) + (filterResult2 / lawsLL);
++lawsLayer;
}
}
@@ -2837,6 +2863,7 @@
}
if(count == 0){
+ // debug statement
printf("no samples for texture\n");
return;
}
@@ -2861,9 +2888,9 @@
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;
@@ -2878,8 +2905,9 @@
}
-void getVoxelMeasures(objStruct objectMetrics[], double *sourceImage, unsigned short *MaskImage,
- int numberObjects, int srcRows, int srcCols){
+void getVoxelMeasures(objStruct objectMetrics[], double *sourceImage,
+ unsigned short *MaskImage, int numberObjects,
+ int srcRows, int srcCols){
int i, j, k;
int label;
@@ -2935,10 +2963,10 @@
}
-int NI_BuildBoundary(int samples, int rows, int cols, int numberObjects, unsigned short *edgeImage,
- objStruct objectMetrics[]){
+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 searchWindow = 5; // 5 is good value for Sobel
int status = 1;
buildBoundary(objectMetrics, searchWindow, edgeImage, numberObjects, rows, cols);
@@ -2966,7 +2994,8 @@
tTEM LawsFeatures[21];
initLaws(&lawsFilter);
- getLawsTexture(lawsFilter, LawsFeatures, objectMetrics, sourceImage, maskImage, numberObjects, rows, cols);
+ getLawsTexture(lawsFilter, LawsFeatures, objectMetrics, sourceImage,
+ maskImage, numberObjects, rows, cols);
return status;
Modified: trunk/scipy/ndimage/segment/tests/test_segment.py
===================================================================
--- trunk/scipy/ndimage/segment/tests/test_segment.py 2007-11-08 03:42:40 UTC (rev 3510)
+++ trunk/scipy/ndimage/segment/tests/test_segment.py 2007-11-09 20:53:50 UTC (rev 3511)
@@ -9,17 +9,65 @@
filename = os.path.join(os.path.split(__file__)[0],inputname)
-def shen_castan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048, highThreshold=600+2048, dust=16):
- labeledEdges, numberObjects = S.shen_castan_edges(scLow, IIRFilter, window, lowThreshold, highThreshold, image)
+def shen_castan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048,
+ highThreshold=600+2048, dust=16):
+ """
+ labeledEdges, ROIList = shen_castan(image, [default])
+
+ implements Shen-Castan edge finding
+
+ Inputs - image, IIR filter, shen_castan_low, window, low_threshold, high_threshold, dust
+ - image is the numarray 2D image
+ - IIR filter is filter parameter for exponential filter
+ - shen_castan_low is edge threshold is range (0.0, 1.0]
+ - window is search window for edge detection
+ - low_ and high_ threshold are density values
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ labeledEdges, numberObjects = S.shen_castan_edges(scLow, IIRFilter, window,
+ lowThreshold, highThreshold, image)
# allocated struct array for edge object measures. for now just the rect bounding box
ROIList = N.zeros(numberObjects, dtype=S.objstruct)
# return the bounding box for each connected edge
S.get_object_stats(labeledEdges, ROIList)
return labeledEdges, ROIList[ROIList['Area']>dust]
-def sobel(image, sLow=0.3, tMode=1, lowThreshold=220+2048, highThreshold=600+2048, BPHigh=10.0, apearture=21, dust=16):
+def sobel(image, sLow=0.3, tMode=1, lowThreshold=220+2048, highThreshold=600+2048, BPHigh=10.0,
+ apearture=21, dust=16):
+ """
+ labeledEdges, ROIList = sobel(image, [default])
+
+ implements sobel magnitude edge finding
+
+ Inputs - image, sobel_low, tMode, low_threshold, high_threshold,
+ high_filter_cutoff, filter_aperature, dust
+ - image is the numarray 2D image
+ - sobel_low is edge threshold is range (0.0, 1.0]
+ - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
+ - low_ and high_ threshold are density values
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - aperature is odd filter kernel length
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
# get sobel edge points. return edges that are labeled (1..numberObjects)
- labeledEdges, numberObjects = S.sobel_edges(sLow, tMode, lowThreshold, highThreshold, BPHigh, apearture, image)
+ labeledEdges, numberObjects = S.sobel_edges(sLow, tMode, lowThreshold,
+ highThreshold, BPHigh, apearture, image)
# allocated struct array for edge object measures. for now just the rect bounding box
ROIList = N.zeros(numberObjects, dtype=S.objstruct)
# return the bounding box for each connected edge
@@ -28,8 +76,34 @@
S.morpho_thin_filt(labeledEdges, ROIList)
return labeledEdges, ROIList[ROIList['Area']>dust]
-def canny(image, cSigma=1.0, cLow=0.5, cHigh=0.8, tMode=1, lowThreshold=220+2048, highThreshold=600+2048,
- BPHigh=10.0, apearture=21, dust=16):
+def canny(image, cSigma=1.0, cLow=0.5, cHigh=0.8, tMode=1, lowThreshold=220+2048,
+ highThreshold=600+2048, BPHigh=10.0, apearture=21, dust=16):
+ """
+ labeledEdges, ROIList = canny(image, [default])
+
+ implements canny edge finding
+
+ Inputs - image, DG_sigma, canny_low, canny_high, tMode, low_threshold,
+ high_threshold, high_filter_cutoff, filter_aperature, dust
+ - image is the numarray 2D image
+ - DG_sigma is Gaussain sigma for the derivative-of-gaussian filter
+ - clow is low edge threshold is range (0.0, 1.0]
+ - chigh is high edge threshold is range (0.0, 1.0]
+ - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
+ - low_ and high_ threshold are density values
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - aperature is odd filter kernel length
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
# get canny edge points. return edges that are labeled (1..numberObjects)
labeledEdges, numberObjects = S.canny_edges(cSigma, cLow, cHigh, tMode, lowThreshold, highThreshold,
BPHigh, apearture, image)
@@ -40,6 +114,23 @@
return labeledEdges, ROIList[ROIList['Area']>dust]
def get_shape_mask(labeledEdges, ROIList):
+ """
+ get_shape_mask(labeledEdges, ROIList)
+
+ takes labeled edge image plus ROIList (blob descriptors) and generates
+ boundary shape features and builds labeled blob masks. 'labeledEdges'
+ is over-written by 'labeledMask'. Adds features to ROIList structure
+
+ Inputs - labeledEdges, ROIList
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. edge image input is over-written with mask image.
+ ROIList added to.
+
+ """
+
# pass in Sobel morph-thinned labeled edge image (LEI) and ROIList
# GetShapeMask will augment the ROI list
# labeledEdges is the original edge image and overwritten as mask image
@@ -48,6 +139,21 @@
return
def get_voxel_measures(rawImage, labeledEdges, ROIList):
+ """
+ get_voxel_measures(rawImage, labeledEdges, ROIList)
+
+ takes raw 2D image, labeled blob mask and ROIList. computes voxel features
+ (moments, histogram) for each blob. Adds features to ROIList structure.
+
+ Inputs - rawImage, labeledEdges, ROIList
+ - rawImage is the original source 2D image
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. ROIList added to.
+
+ """
#
# pass raw image, labeled mask and the partially filled ROIList
# VoxelMeasures will fill the voxel features in the list
@@ -56,14 +162,52 @@
return
def get_texture_measures(rawImage, labeledEdges, ROIList):
+ """
+ get_texture_measures(rawImage, labeledEdges, ROIList)
+
+ takes raw 2D image, labeled blob mask and ROIList. computes 2D
+ texture features using 7x7 Law's texture filters applied
+ to segmented blobs. TEM (texture energy metric) is computed
+ for each Law's filter image and stored in TEM part of ROIList.
+
+ Inputs - rawImage, labeledEdges, ROIList
+ - rawImage is the original source 2D image
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. ROIList added to.
+ """
#
# pass raw image, labeled mask and the partially filled ROIList
- # VoxelMeasures will fill the texture (Law's, co-occurence, Gabor) features in the list
+ # VoxelMeasures will fill the texture (Law's, sub-edges, co-occurence, Gabor) features in the list
#
S.texture_measures(rawImage, labeledEdges, ROIList)
return
def segment_regions():
+ """
+ sourceImage, labeledMask, ROIList = segment_regions()
+
+ Inputs - No Input
+
+ Outputs - sourceImage, labeledMask, ROIList
+ - sourceImage is raw 2D image (default cardiac CT slice for demo
+ - labeledMask is mask of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is numerical Python structure of intensity, shape and
+ texture features for each blob
+
+ High level script calls Python functions:
+ get_slice() - a cardiac CT slice demo file
+ sobel() - sobel magnitude edge finder,
+ returns connected edges
+ get_shape_mask() - gets segmented blob boundary and mask
+ and shape features
+ get_voxel_measures() - uses masks get object voxel moment
+ and histogram features
+ get_texture_measures() - uses masks get object 2D texture features
+ """
# get slice from the CT volume
image = get_slice(filename)
# need a copy of original image as filtering will occur on the extracted slice
@@ -80,6 +224,18 @@
return sourceImage, labeledMask, ROIList
def grow_regions():
+ """
+ regionMask, numberRegions = region_grow()
+ Inputs - No Input
+ Outputs - regionMask, numberRegions
+ - regionMask is the labeled segment masks from 2D image
+ - numberRegions is the number of segmented blobs
+
+ High level script calls Python functions:
+ get_slice() - a cardiac CT slice demo file
+ region_grow() - "grows" connected blobs. default threshold
+ and morphological filter structuring element
+ """
# get slice from the CT volume
image = get_slice(filename)
regionMask, numberRegions = region_grow(image)
@@ -87,12 +243,27 @@
def region_grow(image, lowThreshold=220+2048, highThreshold=600+2048, open=7, close=7):
+ """
+ regionMask, numberRegions = region_grow(image, [defaults])
+
+ Inputs - image, low_threshold, high_threshold, open, close
+ - image is the numarray 2D image
+ - low_ and high_ threshold are density values
+ - open is open morphology structuring element
+ odd size. 0 to turn off. max is 11
+ - close is close morphology structuring element
+ odd size. 0 to turn off. max is 11
+
+ Outputs - regionMask, numberRegions
+ - regionMask is the labeled segment masks from 2D image
+ - numberRegions is the number of segmented blobs
+ """
# morphology filters need to be clipped to 11 max and be odd
regionMask, numberRegions = S.region_grow(lowThreshold, highThreshold, close, open, image)
return regionMask, numberRegions
-def get_slice(imageName='junk.raw', bytes=2, rows=512, columns=512):
+def get_slice(imageName='slice112.raw', bytes=2, rows=512, columns=512):
# get a slice alrady extracted from the CT volume
#image = open(imageName, 'rb')
#slice = image.read(rows*columns*bytes)
@@ -105,7 +276,7 @@
ImageSlice[505:512, :] = 0
return (ImageSlice).astype(float)
-def get_slice2(image_name='junk.raw', bytes=2, shape=(512,512)):
+def get_slice2(image_name='slice112.raw', bytes=2, shape=(512,512)):
import mmap
file = open(image_name, 'rb')
mm = mmap.mmap(file.fileno(), 0, access=mmap.ACCESS_READ)
More information about the Scipy-svn
mailing list