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

scipy-svn@scip... scipy-svn@scip...
Fri Feb 29 18:44:18 CST 2008


Author: tom.waite
Date: 2008-02-29 18:44:16 -0600 (Fri, 29 Feb 2008)
New Revision: 3965

Modified:
   trunk/scipy/ndimage/src/register/Register_IMPL.c
Log:
Bug fix and enhancements

Modified: trunk/scipy/ndimage/src/register/Register_IMPL.c
===================================================================
--- trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-03-01 00:43:51 UTC (rev 3964)
+++ trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-03-01 00:44:16 UTC (rev 3965)
@@ -322,6 +322,7 @@
 	int V101;
 	int V110;
 	int V111;
+	int g[64], f[64];
 	float valueXYZ;
 
 	//
@@ -509,6 +510,180 @@
 
 
 
+int NI_VolumeResample(int layersS, int rowsS, int colsS, int layersD, int rowsD, int colsD,
+	              int scale, int mode, unsigned char *imageD, unsigned char *imageS, double *Z)
+{
+
+	int i;
+	int x, y, z;
+	int sliceSizeSrc;
+	int sliceSizeDst;
+	int status;
+	int ivf;
+	int xf, xg, yg, zg;
+	int g_slice, f_slice;
+	int g_row, f_row;
+	int g_slicesize, f_slicesize;
+	int itemp, sOffset, dOffset;
+	int XInt, YInt, ZInt;
+	float ps1, ps2, ps3;
+	float Y[4], tpoint, reSampler;
+	float XPrime, YPrime, ZPrime;
+	float C, R, L;
+	float *RLUT;
+	float *samples;
+
+	if(mode ==1){
+	    /* 
+	     * integer subsample
+	     */
+	    g_slicesize = rowsD * colsD;
+	    f_slicesize = rowsS * colsS;
+	    for(zg = 0; zg < layersD; ++zg){
+	        g_slice = zg * g_slicesize;
+	        f_slice = zg * scale * f_slicesize;
+	        for(yg = 0; yg < rowsD; ++yg){
+		    g_row = yg * colsD;
+		    f_row = yg * scale * colsS;
+	            for(xg = 0; xg < colsD; ++xg){
+		        xf = xg * scale;
+		        ivf = imageS[f_slice+f_row+xf];
+		        imageD[g_slice+g_row+xg] = ivf;
+	            }
+	        }
+	    }
+	}
+	else if(mode ==2){
+	    /*
+	     * fractional cubic convolution resample
+	     */
+
+	    /* first resample each column in all rows and all layers */
+
+	    sliceSizeSrc = colsS * rowsS;
+	    sliceSizeDst = colsD * rowsD;
+
+	    RLUT    = calloc(colsD, sizeof(float));
+	    samples = calloc(colsS+4, sizeof(float));
+	    reSampler = (float)1.0/Z[0];
+	    tpoint = (float)0.0;
+	    for(i = 0; i < colsD; ++i){
+	        RLUT[i] = tpoint;
+	        tpoint += reSampler;
+	    }
+
+	    for(z = 0; z < layersS; ++z){
+                sOffset = z * sliceSizeSrc;
+                dOffset = z * sliceSizeDst;
+	        for(y = 0; y < rowsS; ++y){
+	            for(x = 0; x < colsS; ++x){
+                        samples[x] = (float)imageS[sOffset+x];
+	            }
+	            for(x = 1; x < colsD; ++x){
+	                XPrime = RLUT[x];
+		        XInt   = (int)XPrime;
+		        C      = XPrime - (float)XInt;
+                        Y[0]   = samples[XInt-1];
+                        Y[1]   = samples[XInt];
+                        Y[2]   = samples[XInt+1];
+                        Y[3]   = samples[XInt+2];
+		        ps1    = Y[2] - Y[0];
+		        ps2    = (float)2.0*(Y[0] - Y[1]) + Y[2] - Y[3];
+		        ps3    = -Y[0] + Y[1] - Y[2] + Y[3];
+		        itemp  = (int)(Y[1]+C*(ps1+C*(ps2+C*ps3)));
+		        if(itemp < 0)   itemp = 0;
+		        if(itemp > 255) itemp = 255;
+                        imageD[dOffset+x] = itemp;
+	            }
+                    sOffset += colsS;
+                    dOffset += colsD;
+	        }
+	    }
+	    free(RLUT);
+	    free(samples);
+
+	    /* second resample each row in all columns and all layers */
+	    RLUT    = calloc(rowsD, sizeof(float));
+	    samples = calloc(rowsS+4, sizeof(float));
+	    reSampler = (float)1.0/Z[1];
+	    tpoint = (float)0.0;
+	    for(i = 0; i < rowsD; ++i){
+	        RLUT[i] = tpoint;
+	        tpoint += reSampler;
+	    }
+
+	    for(z = 0; z < layersS; ++z){
+                dOffset = z * sliceSizeDst;
+	        for(x = 0; x < colsD; ++x){
+	            for(y = 0; y < rowsS; ++y){
+                        samples[y] = (float)imageD[dOffset+x+y*colsD];
+	            }
+	            for(y = 1; y < rowsD; ++y){
+	                YPrime = RLUT[y];
+		        YInt   = (int)YPrime;
+		        R      = YPrime - (float)YInt;
+                        Y[0]   = samples[YInt-1];
+                        Y[1]   = samples[YInt];
+                        Y[2]   = samples[YInt+1];
+                        Y[3]   = samples[YInt+2];
+		        ps1    = Y[2] - Y[0];
+		        ps2    = (float)2.0*(Y[0] - Y[1]) + Y[2] - Y[3];
+		        ps3    = -Y[0] + Y[1] - Y[2] + Y[3];
+		        itemp  = (int)(Y[1]+R*(ps1+R*(ps2+R*ps3)));
+		        if(itemp < 0)   itemp = 0;
+		        if(itemp > 255) itemp = 255;
+                        imageD[dOffset+x+y*colsD] = itemp;
+	            }
+	        }
+	    }
+	    free(RLUT);
+	    free(samples);
+
+	    /* third resample each layers in all columns and all rows */
+	    RLUT    = calloc(layersD, sizeof(float));
+	    samples = calloc(layersS+4, sizeof(float));
+	    reSampler = (float)1.0/Z[2];
+	    tpoint = (float)0.0;
+	    for(i = 0; i < layersD; ++i){
+	        RLUT[i] = tpoint;
+	        tpoint += reSampler;
+	    }
+
+	    for(y = 0; y < rowsD; ++y){
+                dOffset = y * colsD;
+	        for(x = 0; x < colsD; ++x){
+	    	    for(z = 0; z < layersS; ++z){
+                        samples[z] = (float)imageD[dOffset+x+z*sliceSizeDst];
+		    }
+	    	    for(z = 1; z < layersD; ++z){
+	                ZPrime = RLUT[z];
+		        ZInt   = (int)ZPrime;
+		        L      = ZPrime - (float)ZInt;
+                        Y[0]   = samples[ZInt-1];
+                        Y[1]   = samples[ZInt];
+                        Y[2]   = samples[ZInt+1];
+                        Y[3]   = samples[ZInt+2];
+		        ps1    = Y[2] - Y[0];
+		        ps2    = (float)2.0*(Y[0] - Y[1]) + Y[2] - Y[3];
+		        ps3    = -Y[0] + Y[1] - Y[2] + Y[3];
+		        itemp  = (int)(Y[1]+R*(ps1+R*(ps2+R*ps3)));
+		        if(itemp < 0)   itemp = 0;
+		        if(itemp > 255) itemp = 255;
+                        imageD[dOffset+x+z*sliceSizeDst] = itemp;
+		    }
+		}
+	    }
+	    free(RLUT);
+	    free(samples);
+     	}
+
+	status = 1;
+
+	return status;
+
+}
+
+
 int NI_CubicResample(int layersF, int rowsF, int colsF, int layersG, int rowsG, int colsG,
 	             int *dimSteps, double *M, unsigned char *imageG, unsigned char *imageF)
 {
@@ -541,6 +716,7 @@
 				          	zp, colsG, rowsG, layersG, sliceSizeG);
 			/* clip at hard edges */
 			if(vf < 0.0) vf = 0.0;
+			if(vf > 255.0) vf = 255.0;
 			imageG[sliceG+rowG+(int)x] = (int)vf;
 		    }
 	        }
@@ -553,8 +729,6 @@
 
 }
 
-
-
 int NI_ImageThreshold(int layers, int rows, int cols, unsigned short *image, double *H,
 	               double *IH, int histogram_elements, double threshold, int *index)
 {



More information about the Scipy-svn mailing list