[Scipy-svn] r4169 - trunk/scipy/ndimage/src/register

scipy-svn@scip... scipy-svn@scip...
Wed Apr 23 19:32:00 CDT 2008


Author: tom.waite
Date: 2008-04-23 19:31:58 -0500 (Wed, 23 Apr 2008)
New Revision: 4169

Modified:
   trunk/scipy/ndimage/src/register/Register_IMPL.c
Log:
added resample_with_gradient for affine xform part of elastic registration.

Modified: trunk/scipy/ndimage/src/register/Register_IMPL.c
===================================================================
--- trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-04-24 00:31:38 UTC (rev 4168)
+++ trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-04-24 00:31:58 UTC (rev 4169)
@@ -784,3 +784,171 @@
 }
 
 
+int NI_ResampleWithGradient(int layersS, int rowsS, int colsS, int layersD, int rowsD,
+		            int colsD, int *dimSteps, double *M, unsigned char *imageD,
+			    unsigned char *imageS, double *scale, int *offset, double *gradientX,
+			    double *gradientY, double *gradientZ)
+{
+
+	int i;
+	int seed;
+	int status;
+	int sliceD;
+	int rowD;
+	int sliceSizeD;
+	int dimsS[3];
+	int dimsD[3];
+	int dims[2];
+	float vs;
+	float x, y, z;
+	float xp, yp, zp;
+	float dx1, dy1, dz1;
+	float dx2, dy2, dz2;
+	float ran_x, ran_y, ran_z;
+	float dx, dy, dz;
+	double gradX, gradY, gradZ;
+
+	int ptr_x0;
+	int ptr_y0;
+	int ptr_z0;
+	int ptr_x1;
+	int ptr_y1;
+	int ptr_z1;
+	//
+	// Vxyz for [0,1] values of x, y, z
+	//
+	int V000;
+	int V100;
+	int V010;
+	int V001;
+	int V011;
+	int V101;
+	int V110;
+	int V111;
+	float valueXYZ;
+
+	sliceSizeD = rowsD * colsD;
+	dimsD[0] = colsD;
+	dimsD[1] = rowsD;
+	dimsD[2] = layersD;
+	dimsS[0] = colsS;
+	dimsS[1] = rowsS;
+	dimsS[2] = layersS;
+
+	dims[0] = dimsS[0];
+	dims[1] = dimsS[0]*dimsS[1];
+
+	seed = 1000;
+	srand(seed);
+
+	for(z = 0.0; z < layersD-dimSteps[2]-1; z += dimSteps[2]){
+	    sliceD = (int)z * sliceSizeD;
+	    for(y = 0.0; y < rowsD-dimSteps[1]-1; y += dimSteps[1]){
+		rowD = (int)y * colsD;
+	        for(x = 0.0; x < colsD-dimSteps[0]-1; x += dimSteps[0]){
+
+		    /* jitter the coordinates to prevent aliasing */
+		    ran_x = 1.0 * rand()/((float)RAND_MAX);
+		    ran_y = 1.0 * rand()/((float)RAND_MAX);
+		    ran_z = 1.0 * rand()/((float)RAND_MAX);
+
+		    dx = x + ran_x;
+		    dy = y + ran_y;
+		    dz = z + ran_z;
+
+		    // get the 'from' coordinates 
+		    xp = M[0]*dx + M[1]*dy + M[2]*dz  + M[3];
+		    yp = M[4]*dx + M[5]*dy + M[6]*dz  + M[7];
+		    zp = M[8]*dx + M[9]*dy + M[10]*dz + M[11];
+		    // clip the resample window 
+		    if((zp >= 0.0 && zp < layersS-dimSteps[2]) && 
+		       (yp >= 0.0 && yp < rowsS-dimSteps[1]) && 
+		       (xp >= 0.0 && xp < colsS-dimSteps[0])){
+
+			// corners of the 3D unit volume cube
+	    		ptr_z0 = (int)zp * dims[1];
+	    		ptr_z1 = ptr_z0 + dims[1];
+			ptr_y0 = (int)yp * dims[0];
+			ptr_y1 = ptr_y0 + dims[0];
+		    	ptr_x0 = (int)xp;
+		    	ptr_x1 = ptr_x0 + 1;
+
+			dx1 = xp - (int)xp; 
+			dy1 = yp - (int)yp; 
+			dz1 = zp - (int)zp; 
+			dx2 = 1.0 - dx1; 
+			dy2 = 1.0 - dy1; 
+			dz2 = 1.0 - dz1; 
+
+			V000 = imageS[ptr_x0+ptr_y0+ptr_z0];
+			V100 = imageS[ptr_x1+ptr_y0+ptr_z0];
+			V010 = imageS[ptr_x0+ptr_y1+ptr_z0];
+			V001 = imageS[ptr_x0+ptr_y0+ptr_z1];
+			V011 = imageS[ptr_x0+ptr_y1+ptr_z1];
+			V101 = imageS[ptr_x1+ptr_y0+ptr_z1];
+			V110 = imageS[ptr_x1+ptr_y1+ptr_z0];
+			V111 = imageS[ptr_x1+ptr_y1+ptr_z1];
+			
+		        vs = V000 * (dx2) * (dy2) * (dz2) +
+			     V100 * (dx1) * (dy2) * (dz2) +
+			     V010 * (dx2) * (dy1) * (dz2) +
+			     V001 * (dx2) * (dy2) * (dz1) +
+			     V101 * (dx1) * (dy2) * (dz1) +
+			     V011 * (dx2) * (dy1) * (dz1) +
+			     V110 * (dx1) * (dy1) * (dz2) +
+			     V111 * (dx1) * (dy1) * (dz1);
+
+			/* resampled voxel */
+			imageD[sliceD+rowD+(int)x] = (int)(vs*scale[(int)zp]) + offset[(int)zp];
+
+			/*
+			 * x gradient voxel. for no resample dz1, dy1 = 0.0 and
+			 * dy2, dz2 = 1.0 so gradX = V100 - V000
+			*/
+
+			/* d/d(dx1) = 1.0, d/d(dx2) = -1.0 */
+		        gradX = V000 * (-1.0) * (dy2) * (dz2) +
+			        V100 * (1.0)  * (dy2) * (dz2) +
+			        V010 * (-1.0) * (dy1) * (dz2) +
+			        V001 * (-1.0) * (dy2) * (dz1) +
+			        V101 * (1.0)  * (dy2) * (dz1) +
+			        V011 * (-1.0) * (dy1) * (dz1) +
+			        V110 * (1.0)  * (dy1) * (dz2) +
+			        V111 * (1.0)  * (dy1) * (dz1);
+
+			/* d/d(dy1) = 1.0, d/d(dy2) = -1.0 */
+		        gradY = V000 * (dx2) * (-1.0) * (dz2) +
+			        V100 * (dx1) * (-1.0) * (dz2) +
+			        V010 * (dx2) * (1.0)  * (dz2) +
+			        V001 * (dx2) * (-1.0) * (dz1) +
+			        V101 * (dx1) * (-1.0) * (dz1) +
+			        V011 * (dx2) * (1.0)  * (dz1) +
+			        V110 * (dx1) * (1.0)  * (dz2) +
+			        V111 * (dx1) * (1.0)  * (dz1);
+
+			/* d/d(dz1) = 1.0, d/d(dz2) = -1.0 */
+		        gradZ = V000 * (dx2) * (dy2) * (-1.0) +
+			        V100 * (dx1) * (dy2) * (-1.0) +
+			        V010 * (dx2) * (dy1) * (-1.0) +
+			        V001 * (dx2) * (dy2) * (1.0)  +
+			        V101 * (dx1) * (dy2) * (1.0)  +
+			        V011 * (dx2) * (dy1) * (1.0)  +
+			        V110 * (dx1) * (dy1) * (-1.0) +
+			        V111 * (dx1) * (dy1) * (1.0);
+
+			gradientX[sliceD+rowD+(int)x] = (int)(gradX*scale[(int)zp]);
+			gradientY[sliceD+rowD+(int)x] = (int)(gradY*scale[(int)zp]);
+			gradientZ[sliceD+rowD+(int)x] = (int)(gradZ*scale[(int)zp]);
+
+		    }
+	        }
+	    }
+	}
+
+	status = 1;
+
+	return status;
+
+}
+
+



More information about the Scipy-svn mailing list