Skip to content

Commit

Permalink
MAINT: Added cython version of zoom_shift and geometric transform in …
Browse files Browse the repository at this point in the history 
…ni_interpolatio.pyx file
  • Loading branch information
bewithaman committed Jul 14, 2015
1 parent c7d2e10 commit a4f0973
Showing 1 changed file with 0 additions and 311 deletions.
311 changes: 0 additions & 311 deletions scipy/ndimage/src/ni_interpolation.c
View file
Expand Up @@ -666,314 +666,3 @@ NI_GeometricTransform(PyArrayObject *input, int (*map)(npy_intp*, double*,
return PyErr_Occurred() ? 0 : 1;
}

int NI_ZoomShift(PyArrayObject *input, PyArrayObject* zoom_ar,
PyArrayObject* shift_ar, PyArrayObject *output,
int order, int mode, double cval)
{
char *po, *pi;
npy_intp **zeros = NULL, **offsets = NULL, ***edge_offsets = NULL;
npy_intp ftmp[MAXDIM], *fcoordinates = NULL, *foffsets = NULL;
npy_intp jj, hh, kk, filter_size, odimensions[MAXDIM];
npy_intp idimensions[MAXDIM], istrides[MAXDIM];
npy_intp size;
double ***splvals = NULL;
NI_Iterator io;
Float64 *zooms = zoom_ar ? (Float64*)PyArray_DATA(zoom_ar) : NULL;
Float64 *shifts = shift_ar ? (Float64*)PyArray_DATA(shift_ar) : NULL;
int rank = 0, qq;
NPY_BEGIN_THREADS_DEF;

NPY_BEGIN_THREADS;

for(kk = 0; kk < input->nd; kk++) {
idimensions[kk] = input->dimensions[kk];
istrides[kk] = input->strides[kk];
odimensions[kk] = output->dimensions[kk];
}
rank = input->nd;

/* if the mode is 'constant' we need some temps later: */
if (mode == NI_EXTEND_CONSTANT) {
zeros = malloc(rank * sizeof(npy_intp*));
if (NI_UNLIKELY(!zeros)) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}
for(jj = 0; jj < rank; jj++)
zeros[jj] = NULL;
for(jj = 0; jj < rank; jj++) {
zeros[jj] = malloc(odimensions[jj] * sizeof(npy_intp));
if (NI_UNLIKELY(!zeros[jj])) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}
}
}

/* store offsets, along each axis: */
offsets = malloc(rank * sizeof(npy_intp*));
/* store spline coefficients, along each axis: */
splvals = malloc(rank * sizeof(double**));
/* store offsets at all edges: */
edge_offsets = malloc(rank * sizeof(npy_intp**));
if (NI_UNLIKELY(!offsets || !splvals || !edge_offsets)) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}
for(jj = 0; jj < rank; jj++) {
offsets[jj] = NULL;
splvals[jj] = NULL;
edge_offsets[jj] = NULL;
}
for(jj = 0; jj < rank; jj++) {
offsets[jj] = malloc(odimensions[jj] * sizeof(npy_intp));
splvals[jj] = malloc(odimensions[jj] * sizeof(double*));
edge_offsets[jj] = malloc(odimensions[jj] * sizeof(npy_intp*));
if (NI_UNLIKELY(!offsets[jj] || !splvals[jj] || !edge_offsets[jj])) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}
for(hh = 0; hh < odimensions[jj]; hh++) {
splvals[jj][hh] = NULL;
edge_offsets[jj][hh] = NULL;
}
}

/* precalculate offsets, and offsets at the edge: */
for(jj = 0; jj < rank; jj++) {
double shift = 0.0, zoom = 0.0;
if (shifts)
shift = shifts[jj];
if (zooms)
zoom = zooms[jj];
for(kk = 0; kk < odimensions[jj]; kk++) {
double cc = (double)kk;
if (shifts)
cc += shift;
if (zooms)
cc *= zoom;
cc = map_coordinate(cc, idimensions[jj], mode);
if (cc > -1.0) {
npy_intp start;
if (zeros && zeros[jj])
zeros[jj][kk] = 0;
if (order & 1) {
start = (npy_intp)floor(cc) - order / 2;
} else {
start = (npy_intp)floor(cc + 0.5) - order / 2;
}
offsets[jj][kk] = istrides[jj] * start;
if (start < 0 || start + order >= idimensions[jj]) {
edge_offsets[jj][kk] = malloc((order + 1) * sizeof(npy_intp));
if (NI_UNLIKELY(!edge_offsets[jj][kk])) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}
for(hh = 0; hh <= order; hh++) {
npy_intp idx = start + hh;
npy_intp len = idimensions[jj];
if (len <= 1) {
idx = 0;
} else {
npy_intp s2 = 2 * len - 2;
if (idx < 0) {
idx = s2 * (npy_intp)(-idx / s2) + idx;
idx = idx <= 1 - len ? idx + s2 : -idx;
} else if (idx >= len) {
idx -= s2 * (npy_intp)(idx / s2);
if (idx >= len)
idx = s2 - idx;
}
}
edge_offsets[jj][kk][hh] = istrides[jj] * (idx - start);
}
}
if (order > 0) {
splvals[jj][kk] = malloc((order + 1) * sizeof(double));
if (NI_UNLIKELY(!splvals[jj][kk])) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}
spline_coefficients(cc, order, splvals[jj][kk]);
}
} else {
zeros[jj][kk] = 1;
}
}
}

filter_size = 1;
for(jj = 0; jj < rank; jj++)
filter_size *= order + 1;

if (!NI_InitPointIterator(output, &io))
goto exit;

pi = (void *)PyArray_DATA(input);
po = (void *)PyArray_DATA(output);

/* store all coordinates and offsets with filter: */
fcoordinates = malloc(rank * filter_size * sizeof(npy_intp));
foffsets = malloc(filter_size * sizeof(npy_intp));
if (NI_UNLIKELY(!fcoordinates || !foffsets)) {
NPY_END_THREADS;
PyErr_NoMemory();
goto exit;
}

for(jj = 0; jj < rank; jj++)
ftmp[jj] = 0;
kk = 0;
for(hh = 0; hh < filter_size; hh++) {
for(jj = 0; jj < rank; jj++)
fcoordinates[jj + hh * rank] = ftmp[jj];
foffsets[hh] = kk;
for(jj = rank - 1; jj >= 0; jj--) {
if (ftmp[jj] < order) {
ftmp[jj]++;
kk += istrides[jj];
break;
} else {
ftmp[jj] = 0;
kk -= istrides[jj] * order;
}
}
}
size = 1;
for(qq = 0; qq < output->nd; qq++)
size *= output->dimensions[qq];
for(kk = 0; kk < size; kk++) {
double t = 0.0;
npy_intp edge = 0, oo = 0, zero = 0;

for(hh = 0; hh < rank; hh++) {
if (zeros && zeros[hh][io.coordinates[hh]]) {
/* we use constant border condition */
zero = 1;
break;
}
oo += offsets[hh][io.coordinates[hh]];
if (edge_offsets[hh][io.coordinates[hh]])
edge = 1;
}

if (!zero) {
npy_intp *ff = fcoordinates;
const int type_num = NI_NormalizeType(input->descr->type_num);
t = 0.0;
for(hh = 0; hh < filter_size; hh++) {
npy_intp idx = 0;
double coeff = 0.0;

if (NI_UNLIKELY(edge)) {
/* use precalculated edge offsets: */
for(jj = 0; jj < rank; jj++) {
if (edge_offsets[jj][io.coordinates[jj]])
idx += edge_offsets[jj][io.coordinates[jj]][ff[jj]];
else
idx += ff[jj] * istrides[jj];
}
idx += oo;
} else {
/* use normal offsets: */
idx += oo + foffsets[hh];
}
switch (type_num) {
CASE_INTERP_COEFF(coeff, pi, idx, Bool);
CASE_INTERP_COEFF(coeff, pi, idx, UInt8);
CASE_INTERP_COEFF(coeff, pi, idx, UInt16);
CASE_INTERP_COEFF(coeff, pi, idx, UInt32);
#if HAS_UINT64
CASE_INTERP_COEFF(coeff, pi, idx, UInt64);
#endif
CASE_INTERP_COEFF(coeff, pi, idx, Int8);
CASE_INTERP_COEFF(coeff, pi, idx, Int16);
CASE_INTERP_COEFF(coeff, pi, idx, Int32);
CASE_INTERP_COEFF(coeff, pi, idx, Int64);
CASE_INTERP_COEFF(coeff, pi, idx, Float32);
CASE_INTERP_COEFF(coeff, pi, idx, Float64);
default:
NPY_END_THREADS;
PyErr_SetString(PyExc_RuntimeError,
"data type not supported");
goto exit;
}
/* calculate interpolated value: */
for(jj = 0; jj < rank; jj++)
if (order > 0)
coeff *= splvals[jj][io.coordinates[jj]][ff[jj]];
t += coeff;
ff += rank;
}
} else {
t = cval;
}
/* store output: */
switch (NI_NormalizeType(output->descr->type_num)) {
CASE_INTERP_OUT(po, t, Bool);
CASE_INTERP_OUT_UINT(po, t, UInt8, 0, MAX_UINT8);
CASE_INTERP_OUT_UINT(po, t, UInt16, 0, MAX_UINT16);
CASE_INTERP_OUT_UINT(po, t, UInt32, 0, MAX_UINT32);
#if HAS_UINT64
/* There was a bug in numpy as of (at least) <= 1.6.1 such that
* MAX_UINT64 was incorrectly defined, leading to a compiler error.
* NPY_MAX_UINT64 is correctly defined
*/
CASE_INTERP_OUT_UINT(po, t, UInt64, 0, NPY_MAX_UINT64);
#endif
CASE_INTERP_OUT_INT(po, t, Int8, MIN_INT8, MAX_INT8);
CASE_INTERP_OUT_INT(po, t, Int16, MIN_INT16, MAX_INT16);
CASE_INTERP_OUT_INT(po, t, Int32, MIN_INT32, MAX_INT32);
CASE_INTERP_OUT_INT(po, t, Int64, MIN_INT64, MAX_INT64);
CASE_INTERP_OUT(po, t, Float32);
CASE_INTERP_OUT(po, t, Float64);
default:
NPY_END_THREADS;
PyErr_SetString(PyExc_RuntimeError, "data type not supported");
goto exit;
}
NI_ITERATOR_NEXT(io, po);
}

exit:
NPY_END_THREADS;
if (zeros) {
for(jj = 0; jj < rank; jj++)
free(zeros[jj]);
free(zeros);
}
if (offsets) {
for(jj = 0; jj < rank; jj++)
free(offsets[jj]);
free(offsets);
}
if (splvals) {
for(jj = 0; jj < rank; jj++) {
if (splvals[jj]) {
for(hh = 0; hh < odimensions[jj]; hh++)
free(splvals[jj][hh]);
free(splvals[jj]);
}
}
free(splvals);
}
if (edge_offsets) {
for(jj = 0; jj < rank; jj++) {
if (edge_offsets[jj]) {
for(hh = 0; hh < odimensions[jj]; hh++)
free(edge_offsets[jj][hh]);
free(edge_offsets[jj]);
}
}
free(edge_offsets);
}
free(foffsets);
free(fcoordinates);
return PyErr_Occurred() ? 0 : 1;
}

0 comments on commit a4f0973

Please sign in to comment.