[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