I have been trying to import a C code in python file. But its not working. So I've decided to convert the C program into Python so that importing the function will be easier. The C code I want to convert is given below. [I got this from github]
#include "Python.h"
#include "numpy/arrayobject.h"
#include
# define CUBE[x] [[x] * [x] * [x]]
# define SQR[x] [[x] * [x]]
static PyObject *interp3_tricubic[PyObject *self, PyObject *args];
float TriCubic [float px, float py, float pz, float *volume, int xDim, int yDim, int zDim];
// what function are exported
static PyMethodDef tricubicmethods[] = {
{"_interp3_tricubic", interp3_tricubic, METH_VARARGS},
{NULL, NULL}
};
// This function is essential for an extension for Numpy created in C
void inittricubic[] {
[void] Py_InitModule["tricubic", tricubicmethods];
import_array[];
}
// the data should be FLOAT32 and should be ensured in the wrapper
static PyObject *interp3_tricubic[PyObject *self, PyObject *args]
{
PyArrayObject *volume, *result, *C, *R, *S;
float *pr, *pc, *ps;
float *pvol, *pvc;
int xdim, ydim, zdim;
// We expect 4 arguments of the PyArray_Type
if[!PyArg_ParseTuple[args, "O!O!O!O!",
&PyArray_Type, &volume,
&PyArray_Type, &C,
&PyArray_Type, &R,
&PyArray_Type, &S]] return NULL;
if [ NULL == volume ] return NULL;
if [ NULL == C ] return NULL;
if [ NULL == R ] return NULL;
if [ NULL == S ] return NULL;
// result matrix is the same size as C and is float
result = [PyArrayObject*] PyArray_ZEROS[PyArray_NDIM[C], C->dimensions, NPY_FLOAT, 0];
// This is for reference counting [ I think ]
PyArray_FLAGS[result] |= NPY_OWNDATA;
// massive use of iterators to progress through the data
PyArrayIterObject *itr_v, *itr_r, *itr_c, *itr_s;
itr_v = [PyArrayIterObject *] PyArray_IterNew[result];
itr_r = [PyArrayIterObject *] PyArray_IterNew[R];
itr_c = [PyArrayIterObject *] PyArray_IterNew[C];
itr_s = [PyArrayIterObject *] PyArray_IterNew[S];
pvol = [float *]PyArray_DATA[volume];
xdim = PyArray_DIM[volume, 0];
ydim = PyArray_DIM[volume, 1];
zdim = PyArray_DIM[volume, 2];
while[PyArray_ITER_NOTDONE[itr_v]] {
pvc = [float *] PyArray_ITER_DATA[itr_v];
pr = [float *] PyArray_ITE R_DATA[itr_r];
pc = [float *] PyArray_ITER_DATA[itr_c];
ps = [float *] PyArray_ITER_DATA[itr_s];
*pvc = TriCubic[*pc, *pr, *ps, pvol, xdim, ydim, zdim];
PyArray_ITER_NEXT[itr_v];
PyArray_ITER_NEXT[itr_r];
PyArray_ITER_NEXT[itr_c];
PyArray_ITER_NEXT[itr_s];
}
return result;
}
/*
* TriCubic - tri-cubic interpolation at point, p.
* inputs:
* px, py, pz - the interpolation point.
* volume - a pointer to the float volume data, stored in x,
* y, then z order [x index increasing fastest].
* xDim, yDim, zDim - dimensions of the array of volume data.
* returns:
* the interpolated value at p.
* note:
* rudimentary range checking is done in this function.
*/
float TriCubic [float px, float py, float pz, float *volume, int xDim, int yDim, int zDim]
{
int x, y, z;
int i, j, k;
float dx, dy, dz;
float *pv;
float u[4], v[4], w[4];
float r[4], q[4];
float vox = 0;
int xyDim;
xyDim = xDim * yDim;
x = [int] px, y = [int] py, z = [int] pz;
// necessary evil truncating at dim-2 because tricubic needs 2 more values
// which is criminal near edges
// future work includes doing trilinear for edge cases
// range checking is extremely important here
if [x < 3 || x > xDim-3 || y < 3 || y > yDim-3 || z < 3 || z > zDim-3]
return [0];
dx = px - [float] x, dy = py - [float] y, dz = pz - [float] z;
pv = volume + [x - 1] + [y - 1] * xDim + [z - 1] * xyDim;
/* factors for Catmull-Rom interpolation */
u[0] = -0.5 * CUBE [dx] + SQR [dx] - 0.5 * dx;
u[1] = 1.5 * CUBE [dx] - 2.5 * SQR [dx] + 1;
u[2] = -1.5 * CUBE [dx] + 2 * SQR [dx] + 0.5 * dx;
u[3] = 0.5 * CUBE [dx] - 0.5 * SQR [dx];
v[0] = -0.5 * CUBE [dy] + SQR [dy] - 0.5 * dy;
v[1] = 1.5 * CUBE [dy] - 2.5 * SQR [dy] + 1;
v[2] = -1.5 * CUBE [dy] + 2 * SQR [dy] + 0.5 * dy;
v[3] = 0.5 * CUBE [dy] - 0.5 * SQR [dy];
w[0] = -0.5 * CUBE [dz] + SQR [dz] - 0.5 * dz;
w[1] = 1.5 * CUBE [dz] - 2.5 * SQR [dz] + 1;
w[2] = -1.5 * CUBE [dz] + 2 * SQR [dz] + 0.5 * dz;
w[3] = 0.5 * CUBE [dz] - 0.5 * SQR [dz];
for [k = 0; k < 4; k++]
{
q[k] = 0;
for [j = 0; j < 4; j++]
{
r[j] = 0;
for [i = 0; i < 4; i++]
{
r[j] += u[i] * *pv;
pv++;
}
q[k] += v[j] * r[j];
pv += xDim - 4;
}
vox += w[k] * q[k];
pv += xyDim - 4 * xDim;
}
return vox;
}
I have tried to convert this code to python. But the output I got is wrong. The python code I created is added below.
import numpy as N
import math
import scipy
global result
def interp3_tricubic[volume, C, R, S]:
if volume is None :
result = 0
elif C is None:
result = 0
elif R is None:
result = 0
elif S is None:
result = 0
else:
result = N.zeros[len[C], dtype=['float']]
tri_v = N.array[volume, dtype=["float"]]
tri_r = N.array[R, dtype=["float"]]
tri_c = N.array[C, dtype=["float"]]
tri_s = N.array[S, dtype=["float"]]
tri_vol = N.array[volume, dtype=["float"]]
xDim = volume.shape[0]
yDim = volume.shape[1]
zDim = volume.shape[2]
for i in range[len[C]]:
tri_v = TriCubic[tri_c[i], tri_r[i], tri_s[i], volume, xDim, yDim, zDim]
i = i + 1
# print[tri_v, "tri_v"]
return tri_v
def TriCubic [ px, py, pz, volume, xDim, yDim, zDim]:
xyDim = xDim * yDim
x = px.astype[int]
y = py.astype[int]
z = pz.astype[int]
dx = px - x
dy = py - y
dz = pz - z
pv = volume + [x - 1] + [y - 1] * xDim + [z - 1] * xyDim;
def cube[num]:
return num * num * num
def sqrt[num]:
return num * num
u = N.array[[0,0,0,0], dtype=['float']]
v = N.array[[0,0,0,0], dtype=['float']]
w = N.array[[0,0,0,0], dtype=['float']]
vox = N.zeros_like[volume, dtype=['float']]
u[0] = -0.5 * cube [dx] + sqrt [dx] - 0.5 * dx;
u[1] = 1.5 * cube [dx] - 2.5 * sqrt [dx] + 1;
u[2] = -1.5 * cube [dx] + 2 * sqrt [dx] + 0.5 * dx;
u[3] = 0.5 * cube [dx] - 0.5 * sqrt [dx];
v[0] = -0.5 * cube [dy] + sqrt [dy] - 0.5 * dy;
v[1] = 1.5 * cube [dy] - 2.5 * sqrt [dy] + 1;
v[2] = -1.5 * cube [dy] + 2 * sqrt [dy] + 0.5 * dy;
v[3] = 0.5 * cube [dy] - 0.5 * sqrt [dy];
w[0] = -0.5 * cube [dz] + sqrt [dz] - 0.5 * dz;
w[1] = 1.5 * cube [dz] - 2.5 * sqrt [dz] + 1;
w[2] = -1.5 * cube [dz] + 2 * sqrt [dz] + 0.5 * dz;
w[3] = 0.5 * cube [dz] - 0.5 * sqrt [dz];
k = 0
j = 0
i = 0
q = [0,0,0,0]
r = [0,0,0,0]
for k in range[4]:
for j in range[4]:
for i in range[4]:
r[j] += u[i] * pv[i]
i = i+1
q[k] += v[j] * r[j]
pv += xDim - 4
j = j+1
vox += w[k] * q[k]
pv += xyDim - 4 * xDim
k = k+1
return vox
I am confused on the meaning of some lines. Like these lines...
static PyObject *interp3_tricubic[PyObject *self, PyObject *args];
itr_v = [PyArrayIterObject *] PyArray_IterNew[result];
r[j] += u[i] * *pv;
Please help me correct the code. I am stuck!