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

scipy-svn@scip... scipy-svn@scip...
Thu Feb 14 20:13:20 CST 2008


Author: tom.waite
Date: 2008-02-14 20:13:15 -0600 (Thu, 14 Feb 2008)
New Revision: 3941

Modified:
   trunk/scipy/ndimage/src/register/Register_EXT.c
   trunk/scipy/ndimage/src/register/Register_IMPL.c
Log:
Added tri-cubic resampler


Modified: trunk/scipy/ndimage/src/register/Register_EXT.c
===================================================================
--- trunk/scipy/ndimage/src/register/Register_EXT.c	2008-02-14 16:46:23 UTC (rev 3940)
+++ trunk/scipy/ndimage/src/register/Register_EXT.c	2008-02-15 02:13:15 UTC (rev 3941)
@@ -140,6 +140,62 @@
 
 }
 
+
+static PyObject *Register_CubicResample(PyObject *self, PyObject *args)
+{
+
+    int num;
+    int nd;
+    int type;
+    int itype;
+    int nd_rotmatrix;
+    int nd_S;
+    npy_intp *dimsF;
+    npy_intp *dimsG;
+    npy_intp *dims_rotmatrix;
+    npy_intp *dims_S;
+    unsigned char *imageG;
+    unsigned char *imageF;
+    double        *M;
+    int           *S;
+    PyObject *imgArray1 = NULL;
+    PyObject *imgArray2 = NULL;
+    PyObject *rotArray  = NULL;
+    PyObject *SArray    = NULL;
+	
+    if(!PyArg_ParseTuple(args, "OOOO", &imgArray1, &imgArray2, &rotArray, &SArray))
+	goto exit;
+
+    /* check in the Python code that F and G are the same dims, type */
+    imageF = (unsigned char *)PyArray_DATA(imgArray1);
+    imageG = (unsigned char *)PyArray_DATA(imgArray2);
+    /* reads dims as 0 = layers, 1 = rows, 2 = cols */
+    nd     = PyArray_NDIM(imgArray1);
+    dimsF  = PyArray_DIMS(imgArray1);
+    dimsG  = PyArray_DIMS(imgArray2);
+    type   = PyArray_TYPE(imgArray1);
+    num    = PyArray_SIZE(imgArray1);
+
+    M = (double *)PyArray_DATA(rotArray);
+    nd_rotmatrix   = PyArray_NDIM(rotArray);
+    dims_rotmatrix = PyArray_DIMS(rotArray);
+
+    S = (int *)PyArray_DATA(SArray);
+    nd_S   = PyArray_NDIM(SArray);
+    dims_S = PyArray_DIMS(SArray);
+
+    if(!NI_CubicResample((int)dimsF[0], (int)dimsF[1], (int)dimsF[2], 
+                         (int)dimsG[0], (int)dimsG[1], (int)dimsG[2], 
+		          S, M, imageG, imageF))
+	    goto exit;
+
+exit:
+
+    return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue(""); 
+
+}
+
+
 static PyObject *Register_LinearResample(PyObject *self, PyObject *args)
 {
 
@@ -199,6 +255,7 @@
     { "register_histogram",       Register_Histogram,      METH_VARARGS, NULL },
     { "register_histogram_lite",  Register_HistogramLite,  METH_VARARGS, NULL },
     { "register_linear_resample", Register_LinearResample, METH_VARARGS, NULL },
+    { "register_cubic_resample",  Register_CubicResample,  METH_VARARGS, NULL },
     {  NULL, NULL, 0, NULL},
 };
 

Modified: trunk/scipy/ndimage/src/register/Register_IMPL.c
===================================================================
--- trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-02-14 16:46:23 UTC (rev 3940)
+++ trunk/scipy/ndimage/src/register/Register_IMPL.c	2008-02-15 02:13:15 UTC (rev 3941)
@@ -1,6 +1,100 @@
 #include<stdio.h>
 #include<stdlib.h>
 
+float tri_cubic_convolve(unsigned char *pVolume, int x, int y, int z, float xp, float yp,
+	                 float zp, int colsG, int rowsG, int layersG, int sliceSizeG){
+
+	int i, j, k;
+	int layerOffsets[4];
+	int rowOffsets[4];
+	float ps1, ps2, ps3;
+	float Y[4], NewRow[4], NewLayer[4];
+	float R, C, L, D, T;
+	float valueXYZ = 0.0;
+	float dataCube[4][4][4];
+	/*            [cols][rows][layers] */
+
+	rowOffsets[0]   = (y-1)*colsG;
+	rowOffsets[1]   = (y  )*colsG;
+	rowOffsets[2]   = (y+1)*colsG;
+	rowOffsets[3]   = (y+2)*colsG;
+
+	layerOffsets[0] = (z-1)*sliceSizeG;
+	layerOffsets[1] = (z  )*sliceSizeG;
+	layerOffsets[2] = (z+1)*sliceSizeG;
+	layerOffsets[3] = (z+2)*sliceSizeG;
+
+	/* get numerator for interpolation */
+	C = xp - (float)x;
+	R = yp - (float)y;
+	L = zp - (float)z;
+	D = (float)0.002;
+
+	/* get 4x4 window over all 4 layers */
+	for(i = 0; i < 4; ++i){
+	    for(j = 0; j < 4; ++j){
+		dataCube[0][j][i] = (float)pVolume[layerOffsets[i]+rowOffsets[j]+x-1];
+		dataCube[1][j][i] = (float)pVolume[layerOffsets[i]+rowOffsets[j]+x];
+		dataCube[2][j][i] = (float)pVolume[layerOffsets[i]+rowOffsets[j]+x+1];
+		dataCube[3][j][i] = (float)pVolume[layerOffsets[i]+rowOffsets[j]+x+2];
+	    }
+	}
+
+	for(i = 0; i < 4; ++i){
+	    /* interpolate 4 rows in all 4 layers */
+	    for(j = 0; j < 4; ++j){
+		if(C > D){
+		    Y[0] = dataCube[0][j][i];
+		    Y[1] = dataCube[1][j][i];
+		    Y[2] = dataCube[2][j][i];
+		    Y[3] = dataCube[3][j][i];
+		    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];
+		    NewRow[j] = Y[1]+C*(ps1+C*(ps2+C*ps3));
+		}
+		else{
+		    NewRow[j] = dataCube[1][j][i];
+		}
+	    }
+	    /* interpolate across 4 columns */
+	    if(R > D){
+		Y[0] = NewRow[0];
+		Y[1] = NewRow[1];
+		Y[2] = NewRow[2];
+		Y[3] = NewRow[3];
+		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];
+		T    = (Y[1]+R*(ps1+R*(ps2+R*ps3)));
+		NewLayer[i] = T;
+	    }
+	    else{
+		T = NewRow[1];
+		NewLayer[i] = T;
+	    } 
+	}
+	/* interpolate across 4 layers */
+	if(R > D){
+	    Y[0] = NewLayer[0];
+	    Y[1] = NewLayer[1];
+	    Y[2] = NewLayer[2];
+	    Y[3] = NewLayer[3];
+	    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];
+	    T    = (Y[1]+R*(ps1+R*(ps2+R*ps3)));
+	    valueXYZ = T;
+	}
+	else{
+	    T = NewLayer[1];
+	    valueXYZ = T;
+	} 
+
+	return(valueXYZ);
+
+}
+
 float trilinear_A(unsigned char *pVolume, int x, int y, int z, float dx, float dy, float dz, int dims[]){
 
 	// Vxyz for [0,1] values of x, y, z
@@ -414,3 +508,49 @@
 }
 
 
+
+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)
+{
+
+	int i;
+	int status;
+	int sliceG;
+	int rowG;
+	int sliceSizeG;
+	int ivf;
+	float vf;
+	float x, y, z;
+	float xp, yp, zp;
+
+	sliceSizeG = rowsG * colsG;
+	for(z = 1.0; z < layersG-dimSteps[2]-2; z += dimSteps[2]){
+	    sliceG = (int)z * sliceSizeG;
+	    for(y = 1.0; y < rowsG-dimSteps[1]-2; y += dimSteps[1]){
+		rowG = (int)y * colsG;
+	        for(x = 1.0; x < colsG-dimSteps[0]-2; x += dimSteps[0]){
+		    // 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 >= 1.0 && zp < layersF-dimSteps[2]-2) && 
+		       (yp >= 1.0 && yp < rowsF-dimSteps[1]-2) && 
+		       (xp >= 1.0 && xp < colsF-dimSteps[0]-2)){
+			vf = tri_cubic_convolve(imageF, (int)xp, (int)yp, (int)zp, xp, yp,
+				          	zp, colsG, rowsG, layersG, sliceSizeG);
+			/* clip at hard edges */
+			if(vf < 0.0) vf = 0.0;
+			imageG[sliceG+rowG+(int)x] = (int)vf;
+		    }
+	        }
+	    }
+	}
+
+	status = 1;
+
+	return status;
+
+}
+
+



More information about the Scipy-svn mailing list