[Scipy-svn] r4169 - trunk/scipy/ndimage/src/register
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Apr 23 20:32:00 EDT 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