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

scipy-svn@scip... scipy-svn@scip...
Thu Apr 24 17:53:46 CDT 2008


Author: tom.waite
Date: 2008-04-24 17:53:44 -0500 (Thu, 24 Apr 2008)
New Revision: 4174

Modified:
   trunk/scipy/ndimage/src/register/Register_IMPL.c
Log:
resample_with_gradient to work with numpy.mgrid coordinates

Modified: trunk/scipy/ndimage/src/register/Register_IMPL.c
===================================================================
--- trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-04-24 22:53:28 UTC (rev 4173)
+++ trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-04-24 22:53:44 UTC (rev 4174)
@@ -952,3 +952,156 @@
 }
 
 
+int NI_ResampleWGradientWCoords(int size, int layersS, int rowsS, int colsS, int layersD, int rowsD,
+		                int colsD, int *dimSteps, double *X, double *Y, double *Z, double *M,
+			        unsigned char *imageD, unsigned char *imageS, double *scale, int *offset,
+			        double *gradientX, double *gradientY, double *gradientZ)
+{
+
+	int i;
+	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;
+	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];
+
+	for(i = 0; i < size; ++i){
+	    z = Z[i];
+	    y = Y[i];
+	    x = X[i];
+
+	    sliceD = (int)z * sliceSizeD;
+	    rowD   = (int)y * colsD;
+
+	    // get the 'from' coordinates 
+	    xp = M[0]*x + M[1]*y + M[2]*z  + M[3];
+	    yp = M[4]*x + M[5]*y + M[6]*z  + M[7];
+	    zp = M[8]*x + M[9]*y + M[10]*z + 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