[Scipy-svn] r3511 - in trunk/scipy/ndimage/segment: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Nov 9 14:53:58 CST 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