[Scipy-svn] r4601 - in branches/Interpolate1D: . ndimage ndimage/register ndimage/segment tests
scipy-svn@scip...
scipy-svn@scip...
Tue Aug 5 14:07:54 CDT 2008
Author: fcady
Date: 2008-08-05 14:07:52 -0500 (Tue, 05 Aug 2008)
New Revision: 4601
Added:
branches/Interpolate1D/ndimage/
branches/Interpolate1D/ndimage/nd_image.c
branches/Interpolate1D/ndimage/nd_image.h
branches/Interpolate1D/ndimage/ni_filters.c
branches/Interpolate1D/ndimage/ni_filters.h
branches/Interpolate1D/ndimage/ni_fourier.c
branches/Interpolate1D/ndimage/ni_fourier.h
branches/Interpolate1D/ndimage/ni_interpolation.c
branches/Interpolate1D/ndimage/ni_interpolation.h
branches/Interpolate1D/ndimage/ni_measure.c
branches/Interpolate1D/ndimage/ni_measure.h
branches/Interpolate1D/ndimage/ni_morphology.c
branches/Interpolate1D/ndimage/ni_morphology.h
branches/Interpolate1D/ndimage/ni_support.c
branches/Interpolate1D/ndimage/ni_support.h
branches/Interpolate1D/ndimage/register/
branches/Interpolate1D/ndimage/register/Register_EXT.c
branches/Interpolate1D/ndimage/register/Register_IMPL.c
branches/Interpolate1D/ndimage/segment/
branches/Interpolate1D/ndimage/segment/Segmenter_EXT.c
branches/Interpolate1D/ndimage/segment/Segmenter_IMPL.c
branches/Interpolate1D/ndimage/segment/ndImage_Segmenter_structs.h
branches/Interpolate1D/ndimage_wrapper.py
branches/Interpolate1D/tests/test_ndimage.py
Modified:
branches/Interpolate1D/setup.py
Log:
more ND work, including unit tests. They're all failing right now, though, so I need to work on that.
Added: branches/Interpolate1D/ndimage/nd_image.c
===================================================================
--- branches/Interpolate1D/ndimage/nd_image.c 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/nd_image.c 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,1320 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+ * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#define ND_IMPORT_ARRAY
+#include "nd_image.h"
+#undef ND_IMPORT_ARRAY
+#include "ni_support.h"
+#include "ni_filters.h"
+#include "ni_fourier.h"
+#include "ni_morphology.h"
+#include "ni_interpolation.h"
+#include "ni_measure.h"
+
+typedef struct {
+ PyObject *function;
+ PyObject *extra_arguments;
+ PyObject *extra_keywords;
+} NI_PythonCallbackData;
+
+/* Convert an input array of any type, not necessarily contiguous */
+static int
+NI_ObjectToInputArray(PyObject *object, PyArrayObject **array)
+{
+ *array = NA_InputArray(object, tAny, NPY_ALIGNED|NPY_NOTSWAPPED);
+ return *array ? 1 : 0;
+}
+
+/* Convert an input array of any type, not necessarily contiguous */
+static int
+NI_ObjectToOptionalInputArray(PyObject *object, PyArrayObject **array)
+{
+ if (object == Py_None) {
+ *array = NULL;
+ return 1;
+ } else {
+ *array = NA_InputArray(object, tAny, NPY_ALIGNED|NPY_NOTSWAPPED);
+ return *array ? 1 : 0;
+ }
+}
+
+/* Convert an output array of any type, not necessarily contiguous */
+static int
+NI_ObjectToOutputArray(PyObject *object, PyArrayObject **array)
+{
+ *array = NA_OutputArray(object, tAny, NPY_ALIGNED|NPY_NOTSWAPPED);
+ return *array ? 1 : 0;
+}
+
+/* Convert an output array of any type, not necessarily contiguous */
+static int
+NI_ObjectToOptionalOutputArray(PyObject *object, PyArrayObject **array)
+{
+ if (object == Py_None) {
+ *array = NULL;
+ return 1;
+ } else {
+ *array = NA_OutputArray(object, tAny, NPY_ALIGNED|NPY_NOTSWAPPED);
+ return *array ? 1 : 0;
+ }
+}
+
+/* Convert an input/output array of any type, not necessarily contiguous */
+static int
+NI_ObjectToIoArray(PyObject *object, PyArrayObject **array)
+{
+ *array = NA_IoArray(object, tAny, NPY_ALIGNED|NPY_NOTSWAPPED);
+ return *array ? 1 : 0;
+}
+
+/* Convert an Long sequence */
+static maybelong
+NI_ObjectToLongSequenceAndLength(PyObject *object, maybelong **sequence)
+{
+ long *pa, ii;
+ PyArrayObject *array = NA_InputArray(object, PyArray_LONG, NPY_CARRAY);
+ maybelong length = PyArray_SIZE(array);
+
+ *sequence = (maybelong*)malloc(length * sizeof(maybelong));
+ if (!*sequence) {
+ PyErr_NoMemory();
+ Py_XDECREF(array);
+ return -1;
+ }
+ pa = (long*)PyArray_DATA(array);
+ for(ii = 0; ii < length; ii++)
+ (*sequence)[ii] = pa[ii];
+ Py_XDECREF(array);
+ return length;
+}
+
+static int
+NI_ObjectToLongSequence(PyObject *object, maybelong **sequence)
+{
+ return NI_ObjectToLongSequenceAndLength(object, sequence) >= 0;
+}
+
+/*********************************************************************/
+/* wrapper functions: */
+/*********************************************************************/
+
+static PyObject *Py_Correlate1D(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *weights = NULL;
+ int axis, mode;
+ long origin;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&O&iO&idl", NI_ObjectToInputArray, &input,
+ NI_ObjectToInputArray, &weights, &axis,
+ NI_ObjectToOutputArray, &output, &mode, &cval, &origin))
+ goto exit;
+ if (!NI_Correlate1D(input, weights, axis, output,
+ (NI_ExtendMode)mode, cval, origin))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(weights);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_Correlate(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *weights = NULL;
+ maybelong *origin = NULL;
+ int mode;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&idO&", NI_ObjectToInputArray, &input,
+ NI_ObjectToInputArray, &weights, NI_ObjectToOutputArray, &output,
+ &mode, &cval, NI_ObjectToLongSequence, &origin))
+ goto exit;
+ if (!NI_Correlate(input, weights, output, (NI_ExtendMode)mode, cval,
+ origin))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(weights);
+ Py_XDECREF(output);
+ if (origin)
+ free(origin);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_UniformFilter1D(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL;
+ int axis, mode;
+ long filter_size, origin;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&liO&idl", NI_ObjectToInputArray, &input,
+ &filter_size, &axis, NI_ObjectToOutputArray, &output,
+ &mode, &cval, &origin))
+ goto exit;
+ if (!NI_UniformFilter1D(input, filter_size, axis, output,
+ (NI_ExtendMode)mode, cval, origin))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_MinOrMaxFilter1D(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL;
+ int axis, mode, minimum;
+ long filter_size, origin;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&liO&idli", NI_ObjectToInputArray, &input,
+ &filter_size, &axis, NI_ObjectToOutputArray, &output,
+ &mode, &cval, &origin, &minimum))
+ goto exit;
+ if (!NI_MinOrMaxFilter1D(input, filter_size, axis, output,
+ (NI_ExtendMode)mode, cval, origin, minimum))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_MinOrMaxFilter(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *footprint = NULL;
+ PyArrayObject *structure = NULL;
+ maybelong *origin = NULL;
+ int mode, minimum;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&O&idO&i", NI_ObjectToInputArray,
+ &input, NI_ObjectToInputArray, &footprint,
+ NI_ObjectToOptionalInputArray, &structure,
+ NI_ObjectToOutputArray, &output, &mode, &cval,
+ NI_ObjectToLongSequence, &origin, &minimum))
+ goto exit;
+ if (!NI_MinOrMaxFilter(input, footprint, structure, output,
+ (NI_ExtendMode)mode, cval, origin, minimum))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(footprint);
+ Py_XDECREF(structure);
+ Py_XDECREF(output);
+ if (origin)
+ free(origin);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_RankFilter(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *footprint = NULL;
+ maybelong *origin = NULL;
+ int mode, rank;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&iO&O&idO&", NI_ObjectToInputArray,
+ &input, &rank, NI_ObjectToInputArray, &footprint,
+ NI_ObjectToOutputArray, &output, &mode, &cval,
+ NI_ObjectToLongSequence, &origin))
+ goto exit;
+ if (!NI_RankFilter(input, rank, footprint, output, (NI_ExtendMode)mode,
+ cval, origin))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(footprint);
+ Py_XDECREF(output);
+ if (origin)
+ free(origin);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static int Py_Filter1DFunc(double *iline, maybelong ilen,
+ double *oline, maybelong olen, void *data)
+{
+ PyArrayObject *py_ibuffer = NULL, *py_obuffer = NULL;
+ PyObject *rv = NULL, *args = NULL, *tmp = NULL;
+ maybelong ii;
+ double *po = NULL;
+ NI_PythonCallbackData *cbdata = (NI_PythonCallbackData*)data;
+
+ py_ibuffer = NA_NewArray(iline, PyArray_DOUBLE, 1, &ilen);
+ py_obuffer = NA_NewArray(NULL, PyArray_DOUBLE, 1, &olen);
+ if (!py_ibuffer || !py_obuffer)
+ goto exit;
+ tmp = Py_BuildValue("(OO)", py_ibuffer, py_obuffer);
+ if (!tmp)
+ goto exit;
+ args = PySequence_Concat(tmp, cbdata->extra_arguments);
+ if (!args)
+ goto exit;
+ rv = PyObject_Call(cbdata->function, args, cbdata->extra_keywords);
+ if (!rv)
+ goto exit;
+ po = (double*)PyArray_DATA(py_obuffer);
+ for(ii = 0; ii < olen; ii++)
+ oline[ii] = po[ii];
+exit:
+ Py_XDECREF(py_ibuffer);
+ Py_XDECREF(py_obuffer);
+ Py_XDECREF(rv);
+ Py_XDECREF(args);
+ Py_XDECREF(tmp);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+static PyObject *Py_GenericFilter1D(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL;
+ PyObject *fnc = NULL, *extra_arguments = NULL, *extra_keywords = NULL;
+ void *func = Py_Filter1DFunc, *data = NULL;
+ NI_PythonCallbackData cbdata;
+ int axis, mode;
+ long origin, filter_size;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&OliO&idlOO", NI_ObjectToInputArray,
+ &input, &fnc, &filter_size, &axis, NI_ObjectToOutputArray,
+ &output, &mode, &cval, &origin, &extra_arguments, &extra_keywords))
+ goto exit;
+ if (!PyTuple_Check(extra_arguments)) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "extra_arguments must be a tuple");
+ goto exit;
+ }
+ if (!PyDict_Check(extra_keywords)) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "extra_keywords must be a dictionary");
+ goto exit;
+ }
+ if (PyCObject_Check(fnc)) {
+ func = PyCObject_AsVoidPtr(fnc);
+ data = PyCObject_GetDesc(fnc);
+ } else if (PyCallable_Check(fnc)) {
+ cbdata.function = fnc;
+ cbdata.extra_arguments = extra_arguments;
+ cbdata.extra_keywords = extra_keywords;
+ data = (void*)&cbdata;
+ } else {
+ PyErr_SetString(PyExc_RuntimeError,
+ "function parameter is not callable");
+ goto exit;
+ }
+ if (!NI_GenericFilter1D(input, func, data, filter_size, axis, output,
+ (NI_ExtendMode)mode, cval, origin))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static int Py_FilterFunc(double *buffer, maybelong filter_size,
+ double *output, void *data)
+{
+ PyArrayObject *py_buffer = NULL;
+ PyObject *rv = NULL, *args = NULL, *tmp = NULL;
+ NI_PythonCallbackData *cbdata = (NI_PythonCallbackData*)data;
+
+ py_buffer = NA_NewArray(buffer, PyArray_DOUBLE, 1, &filter_size);
+ if (!py_buffer)
+ goto exit;
+ tmp = Py_BuildValue("(O)", py_buffer);
+ if (!tmp)
+ goto exit;
+ args = PySequence_Concat(tmp, cbdata->extra_arguments);
+ if (!args)
+ goto exit;
+ rv = PyObject_Call(cbdata->function, args, cbdata->extra_keywords);
+ if (!rv)
+ goto exit;
+ *output = PyFloat_AsDouble(rv);
+exit:
+ Py_XDECREF(py_buffer);
+ Py_XDECREF(rv);
+ Py_XDECREF(args);
+ Py_XDECREF(tmp);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+static PyObject *Py_GenericFilter(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *footprint = NULL;
+ PyObject *fnc = NULL, *extra_arguments = NULL, *extra_keywords = NULL;
+ void *func = Py_FilterFunc, *data = NULL;
+ NI_PythonCallbackData cbdata;
+ int mode;
+ maybelong *origin = NULL;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&OO&O&idO&OO", NI_ObjectToInputArray,
+ &input, &fnc, NI_ObjectToInputArray, &footprint,
+ NI_ObjectToOutputArray, &output, &mode, &cval,
+ NI_ObjectToLongSequence, &origin,
+ &extra_arguments, &extra_keywords))
+ goto exit;
+ if (!PyTuple_Check(extra_arguments)) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "extra_arguments must be a tuple");
+ goto exit;
+ }
+ if (!PyDict_Check(extra_keywords)) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "extra_keywords must be a dictionary");
+ goto exit;
+ }
+ if (PyCObject_Check(fnc)) {
+ func = PyCObject_AsVoidPtr(fnc);
+ data = PyCObject_GetDesc(fnc);
+ } else if (PyCallable_Check(fnc)) {
+ cbdata.function = fnc;
+ cbdata.extra_arguments = extra_arguments;
+ cbdata.extra_keywords = extra_keywords;
+ data = (void*)&cbdata;
+ } else {
+ PyErr_SetString(PyExc_RuntimeError,
+ "function parameter is not callable");
+ goto exit;
+ }
+ if (!NI_GenericFilter(input, func, data, footprint, output,
+ (NI_ExtendMode)mode, cval, origin))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(output);
+ Py_XDECREF(footprint);
+ if (origin)
+ free(origin);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_FourierFilter(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *parameters = NULL;
+ int axis, filter_type;
+ long n;
+
+ if (!PyArg_ParseTuple(args, "O&O&liO&i", NI_ObjectToInputArray, &input,
+ NI_ObjectToInputArray, ¶meters, &n, &axis,
+ NI_ObjectToOutputArray, &output, &filter_type))
+ goto exit;
+
+ if (!NI_FourierFilter(input, parameters, n, axis, output, filter_type))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(parameters);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_FourierShift(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *shifts = NULL;
+ int axis;
+ long n;
+
+ if (!PyArg_ParseTuple(args, "O&O&liO&", NI_ObjectToInputArray, &input,
+ NI_ObjectToInputArray, &shifts, &n, &axis,
+ NI_ObjectToOutputArray, &output))
+ goto exit;
+
+ if (!NI_FourierShift(input, shifts, n, axis, output))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(shifts);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_SplineFilter1D(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL;
+ int axis, order;
+
+ if (!PyArg_ParseTuple(args, "O&iiO&", NI_ObjectToInputArray, &input,
+ &order, &axis, NI_ObjectToOutputArray, &output))
+ goto exit;
+
+ if (!NI_SplineFilter1D(input, order, axis, output))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static int Py_Map(maybelong *ocoor, double* icoor, int orank, int irank,
+ void *data)
+{
+ PyObject *coors = NULL, *rets = NULL, *args = NULL, *tmp = NULL;
+ maybelong ii;
+ NI_PythonCallbackData *cbdata = (NI_PythonCallbackData*)data;
+
+ coors = PyTuple_New(orank);
+ if (!coors)
+ goto exit;
+ for(ii = 0; ii < orank; ii++) {
+ PyTuple_SetItem(coors, ii, PyInt_FromLong(ocoor[ii]));
+ if (PyErr_Occurred())
+ goto exit;
+ }
+ tmp = Py_BuildValue("(O)", coors);
+ if (!tmp)
+ goto exit;
+ args = PySequence_Concat(tmp, cbdata->extra_arguments);
+ if (!args)
+ goto exit;
+ rets = PyObject_Call(cbdata->function, args, cbdata->extra_keywords);
+ if (!rets)
+ goto exit;
+ for(ii = 0; ii < irank; ii++) {
+ icoor[ii] = PyFloat_AsDouble(PyTuple_GetItem(rets, ii));
+ if (PyErr_Occurred())
+ goto exit;
+ }
+exit:
+ Py_XDECREF(coors);
+ Py_XDECREF(tmp);
+ Py_XDECREF(rets);
+ Py_XDECREF(args);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+
+static PyObject *Py_GeometricTransform(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL;
+ PyArrayObject *coordinates = NULL, *matrix = NULL, *shift = NULL;
+ PyObject *fnc = NULL, *extra_arguments = NULL, *extra_keywords = NULL;
+ int mode, order;
+ double cval;
+ void *func = NULL, *data = NULL;
+ NI_PythonCallbackData cbdata;
+
+ if (!PyArg_ParseTuple(args, "O&OO&O&O&O&iidOO", NI_ObjectToInputArray,
+ &input, &fnc, NI_ObjectToOptionalInputArray,
+ &coordinates, NI_ObjectToOptionalInputArray,
+ &matrix, NI_ObjectToOptionalInputArray, &shift,
+ NI_ObjectToOutputArray, &output, &order, &mode,
+ &cval, &extra_arguments, &extra_keywords))
+ goto exit;
+
+ if (fnc != Py_None) {
+ if (!PyTuple_Check(extra_arguments)) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "extra_arguments must be a tuple");
+ goto exit;
+ }
+ if (!PyDict_Check(extra_keywords)) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "extra_keywords must be a dictionary");
+ goto exit;
+ }
+ if (PyCObject_Check(fnc)) {
+ func = PyCObject_AsVoidPtr(fnc);
+ data = PyCObject_GetDesc(fnc);
+ } else if (PyCallable_Check(fnc)) {
+ func = Py_Map;
+ cbdata.function = fnc;
+ cbdata.extra_arguments = extra_arguments;
+ cbdata.extra_keywords = extra_keywords;
+ data = (void*)&cbdata;
+ } else {
+ PyErr_SetString(PyExc_RuntimeError,
+ "function parameter is not callable");
+ goto exit;
+ }
+ }
+
+ if (!NI_GeometricTransform(input, func, data, matrix, shift, coordinates,
+ output, order, (NI_ExtendMode)mode, cval))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(output);
+ Py_XDECREF(coordinates);
+ Py_XDECREF(matrix);
+ Py_XDECREF(shift);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_ZoomShift(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *shift = NULL;
+ PyArrayObject *zoom = NULL;
+ int mode, order;
+ double cval;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&O&iid", NI_ObjectToInputArray,
+ &input, NI_ObjectToOptionalInputArray, &zoom,
+ NI_ObjectToOptionalInputArray, &shift, NI_ObjectToOutputArray,
+ &output, &order, &mode, &cval))
+ goto exit;
+
+ if (!NI_ZoomShift(input, zoom, shift, output, order, (NI_ExtendMode)mode,
+ cval))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(shift);
+ Py_XDECREF(zoom);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_Label(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *strct = NULL;
+ maybelong max_label;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&", NI_ObjectToInputArray, &input,
+ NI_ObjectToInputArray, &strct, NI_ObjectToOutputArray, &output))
+ goto exit;
+
+ if (!NI_Label(input, strct, &max_label, output))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(strct);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("l", (long)max_label);
+}
+
+static PyObject *Py_FindObjects(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL;
+ PyObject *result = NULL, *tuple = NULL, *start = NULL, *end = NULL;
+ PyObject *slc = NULL;
+ int jj;
+ long max_label;
+ maybelong ii, *regions = NULL;
+
+ if (!PyArg_ParseTuple(args, "O&l", NI_ObjectToInputArray, &input,
+ &max_label))
+ goto exit;
+
+ if (max_label < 0)
+ max_label = 0;
+ if (max_label > 0) {
+ if (input->nd > 0) {
+ regions = (maybelong*)malloc(2 * max_label * input->nd *
+ sizeof(maybelong));
+ } else {
+ regions = (maybelong*)malloc(max_label * sizeof(maybelong));
+ }
+ if (!regions) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+
+ if (!NI_FindObjects(input, max_label, regions))
+ goto exit;
+
+ result = PyList_New(max_label);
+ if (!result) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+
+ for(ii = 0; ii < max_label; ii++) {
+ maybelong idx = input->nd > 0 ? 2 * input->nd * ii : ii;
+ if (regions[idx] >= 0) {
+ PyObject *tuple = PyTuple_New(input->nd);
+ if (!tuple) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(jj = 0; jj < input->nd; jj++) {
+ start = PyInt_FromLong(regions[idx + jj]);
+ end = PyInt_FromLong(regions[idx + jj + input->nd]);
+ if (!start || !end) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ slc = PySlice_New(start, end, NULL);
+ if (!slc) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ Py_XDECREF(start);
+ Py_XDECREF(end);
+ start = end = NULL;
+ PyTuple_SetItem(tuple, jj, slc);
+ slc = NULL;
+ }
+ PyList_SetItem(result, ii, tuple);
+ tuple = NULL;
+ } else {
+ Py_INCREF(Py_None);
+ PyList_SetItem(result, ii, Py_None);
+ }
+ }
+
+ Py_INCREF(result);
+
+ exit:
+ Py_XDECREF(input);
+ Py_XDECREF(result);
+ Py_XDECREF(tuple);
+ Py_XDECREF(start);
+ Py_XDECREF(end);
+ Py_XDECREF(slc);
+ if (regions)
+ free(regions);
+ if (PyErr_Occurred()) {
+ Py_XDECREF(result);
+ return NULL;
+ } else {
+ return result;
+ }
+}
+
+static PyObject *Py_WatershedIFT(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *markers = NULL;
+ PyArrayObject *strct = NULL;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&O&", NI_ObjectToInputArray, &input,
+ NI_ObjectToInputArray, &markers, NI_ObjectToInputArray,
+ &strct, NI_ObjectToOutputArray, &output))
+ goto exit;
+
+ if (!NI_WatershedIFT(input, markers, strct, output))
+ goto exit;
+
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(markers);
+ Py_XDECREF(strct);
+ Py_XDECREF(output);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static int _NI_GetIndices(PyObject* indices_object,
+ maybelong** result_indices, maybelong* min_label,
+ maybelong* max_label, maybelong* n_results)
+{
+ maybelong *indices = NULL, n_indices, ii;
+
+ if (indices_object == Py_None) {
+ *min_label = -1;
+ *n_results = 1;
+ } else {
+ n_indices = NI_ObjectToLongSequenceAndLength(indices_object, &indices);
+ if (n_indices < 0)
+ goto exit;
+ if (n_indices < 1) {
+ PyErr_SetString(PyExc_RuntimeError, "no correct indices provided");
+ goto exit;
+ } else {
+ *min_label = *max_label = indices[0];
+ if (*min_label < 0) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "negative indices not allowed");
+ goto exit;
+ }
+ for(ii = 1; ii < n_indices; ii++) {
+ if (indices[ii] < 0) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "negative indices not allowed");
+ goto exit;
+ }
+ if (indices[ii] < *min_label)
+ *min_label = indices[ii];
+ if (indices[ii] > *max_label)
+ *max_label = indices[ii];
+ }
+ *result_indices = (maybelong*)malloc((*max_label - *min_label + 1) *
+ sizeof(maybelong));
+ if (!*result_indices) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(ii = 0; ii < *max_label - *min_label + 1; ii++)
+ (*result_indices)[ii] = -1;
+ *n_results = 0;
+ for(ii = 0; ii < n_indices; ii++) {
+ if ((*result_indices)[indices[ii] - *min_label] >= 0) {
+ PyErr_SetString(PyExc_RuntimeError, "duplicate index");
+ goto exit;
+ }
+ (*result_indices)[indices[ii] - *min_label] = ii;
+ ++(*n_results);
+ }
+ }
+ }
+ exit:
+ if (indices)
+ free(indices);
+ return PyErr_Occurred() == NULL;
+}
+
+
+PyObject* _NI_BuildMeasurementResultArrayObject(maybelong n_results,
+ PyArrayObject** values)
+{
+ PyObject *result = NULL;
+ if (n_results > 1) {
+ result = PyList_New(n_results);
+ if (result) {
+ maybelong ii;
+ for(ii = 0; ii < n_results; ii++) {
+ PyList_SET_ITEM(result, ii, (PyObject*)values[ii]);
+ Py_XINCREF(values[ii]);
+ }
+ }
+ } else {
+ result = (PyObject*)values[0];
+ Py_XINCREF(values[0]);
+ }
+ return result;
+}
+
+
+PyObject* _NI_BuildMeasurementResultDouble(maybelong n_results,
+ double* values)
+{
+ PyObject *result = NULL;
+ if (n_results > 1) {
+ result = PyList_New(n_results);
+ if (result) {
+ int ii;
+ for(ii = 0; ii < n_results; ii++) {
+ PyObject* val = PyFloat_FromDouble(values[ii]);
+ if (!val) {
+ Py_XDECREF(result);
+ return NULL;
+ }
+ PyList_SET_ITEM(result, ii, val);
+ }
+ }
+ } else {
+ result = Py_BuildValue("d", values[0]);
+ }
+ return result;
+}
+
+
+PyObject* _NI_BuildMeasurementResultDoubleTuple(maybelong n_results,
+ int tuple_size, double* values)
+{
+ PyObject *result = NULL;
+ maybelong ii;
+ int jj;
+
+ if (n_results > 1) {
+ result = PyList_New(n_results);
+ if (result) {
+ for(ii = 0; ii < n_results; ii++) {
+ PyObject* val = PyTuple_New(tuple_size);
+ if (!val) {
+ Py_XDECREF(result);
+ return NULL;
+ }
+ for(jj = 0; jj < tuple_size; jj++) {
+ maybelong idx = jj + ii * tuple_size;
+ PyTuple_SetItem(val, jj, PyFloat_FromDouble(values[idx]));
+ if (PyErr_Occurred()) {
+ Py_XDECREF(result);
+ return NULL;
+ }
+ }
+ PyList_SET_ITEM(result, ii, val);
+ }
+ }
+ } else {
+ result = PyTuple_New(tuple_size);
+ if (result) {
+ for(ii = 0; ii < tuple_size; ii++) {
+ PyTuple_SetItem(result, ii, PyFloat_FromDouble(values[ii]));
+ if (PyErr_Occurred()) {
+ Py_XDECREF(result);
+ return NULL;
+ }
+ }
+ }
+ }
+ return result;
+}
+
+
+PyObject* _NI_BuildMeasurementResultInt(maybelong n_results,
+ maybelong* values)
+{
+ PyObject *result = NULL;
+ if (n_results > 1) {
+ result = PyList_New(n_results);
+ if (result) {
+ maybelong ii;
+ for(ii = 0; ii < n_results; ii++) {
+ PyObject* val = PyInt_FromLong(values[ii]);
+ if (!val) {
+ Py_XDECREF(result);
+ return NULL;
+ }
+ PyList_SET_ITEM(result, ii, val);
+ }
+ }
+ } else {
+ result = Py_BuildValue("l", values[0]);
+ }
+ return result;
+}
+
+
+static PyObject *Py_Statistics(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *labels = NULL;
+ PyObject *indices_object, *result = NULL;
+ PyObject *res1 = NULL, *res2 = NULL, *res3 = NULL, *res4 = NULL;
+ double *dresult1 = NULL, *dresult2 = NULL;
+ maybelong *lresult1 = NULL, *lresult2 = NULL;
+ maybelong min_label, max_label, *result_indices = NULL, n_results, ii;
+ int type;
+
+ if (!PyArg_ParseTuple(args, "O&O&Oi", NI_ObjectToInputArray, &input,
+ NI_ObjectToOptionalInputArray, &labels, &indices_object, &type))
+ goto exit;
+
+ if (!_NI_GetIndices(indices_object, &result_indices, &min_label,
+ &max_label, &n_results))
+ goto exit;
+
+ if (type >= 0 && type <= 7) {
+ dresult1 = (double*)malloc(n_results * sizeof(double));
+ if (!dresult1) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+ if (type == 2 || type == 7) {
+ dresult2 = (double*)malloc(n_results * sizeof(double));
+ if (!dresult2) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+ if (type == 1 || type == 2 || (type >= 5 && type <= 7)) {
+ lresult1 = (maybelong*)malloc(n_results * sizeof(maybelong));
+ if (!lresult1) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+ if (type == 7) {
+ lresult2 = (maybelong*)malloc(n_results * sizeof(maybelong));
+ if (!lresult2) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+ switch(type) {
+ case 0:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, dresult1, NULL, NULL, NULL, NULL, NULL, NULL))
+ goto exit;
+ result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
+ break;
+ case 1:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, dresult1, lresult1, NULL, NULL, NULL, NULL, NULL))
+ goto exit;
+ for(ii = 0; ii < n_results; ii++)
+ dresult1[ii] = lresult1[ii] > 0 ? dresult1[ii] / lresult1[ii] : 0.0;
+
+ result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
+ break;
+ case 2:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, dresult1, lresult1, dresult2, NULL, NULL, NULL, NULL))
+ goto exit;
+ result = _NI_BuildMeasurementResultDouble(n_results, dresult2);
+ break;
+ case 3:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, NULL, NULL, NULL, dresult1, NULL, NULL, NULL))
+ goto exit;
+ result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
+ break;
+ case 4:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, NULL, NULL, NULL, NULL, dresult1, NULL, NULL))
+ goto exit;
+ result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
+ break;
+ case 5:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, NULL, NULL, NULL, dresult1, NULL, lresult1, NULL))
+ goto exit;
+ result = _NI_BuildMeasurementResultInt(n_results, lresult1);
+ break;
+ case 6:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, NULL, NULL, NULL, NULL, dresult1, NULL, lresult1))
+ goto exit;
+ result = _NI_BuildMeasurementResultInt(n_results, lresult1);
+ break;
+ case 7:
+ if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
+ n_results, NULL, NULL, NULL, dresult1, dresult2,
+ lresult1, lresult2))
+ goto exit;
+ res1 = _NI_BuildMeasurementResultDouble(n_results, dresult1);
+ res2 = _NI_BuildMeasurementResultDouble(n_results, dresult2);
+ res3 = _NI_BuildMeasurementResultInt(n_results, lresult1);
+ res4 = _NI_BuildMeasurementResultInt(n_results, lresult2);
+ if (!res1 || !res2 || !res3 || !res4)
+ goto exit;
+ result = Py_BuildValue("OOOO", res1, res2, res3, res4);
+ break;
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "operation not supported");
+ goto exit;
+ }
+
+ exit:
+ Py_XDECREF(input);
+ Py_XDECREF(labels);
+ if (result_indices)
+ free(result_indices);
+ if (dresult1)
+ free(dresult1);
+ if (dresult2)
+ free(dresult2);
+ if (lresult1)
+ free(lresult1);
+ if (lresult2)
+ free(lresult2);
+ return result;
+}
+
+
+static PyObject *Py_CenterOfMass(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *labels = NULL;
+ PyObject *indices_object, *result = NULL;
+ double *center_of_mass = NULL;
+ maybelong min_label, max_label, *result_indices = NULL, n_results;
+
+ if (!PyArg_ParseTuple(args, "O&O&O", NI_ObjectToInputArray, &input,
+ NI_ObjectToOptionalInputArray, &labels, &indices_object))
+ goto exit;
+
+ if (!_NI_GetIndices(indices_object, &result_indices, &min_label,
+ &max_label, &n_results))
+ goto exit;
+
+ center_of_mass = (double*)malloc(input->nd * n_results *
+ sizeof(double));
+ if (!center_of_mass) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+
+ if (!NI_CenterOfMass(input, labels, min_label, max_label,
+ result_indices, n_results, center_of_mass))
+ goto exit;
+
+ result = _NI_BuildMeasurementResultDoubleTuple(n_results, input->nd,
+ center_of_mass);
+
+ exit:
+ Py_XDECREF(input);
+ Py_XDECREF(labels);
+ if (result_indices)
+ free(result_indices);
+ if (center_of_mass)
+ free(center_of_mass);
+ return result;
+}
+
+static PyObject *Py_Histogram(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *labels = NULL, **histograms = NULL;
+ PyObject *indices_object, *result = NULL;
+ maybelong min_label, max_label, *result_indices = NULL, n_results;
+ maybelong jj, nbins;
+ long nbins_in;
+ double min, max;
+
+ if (!PyArg_ParseTuple(args, "O&ddlO&O", NI_ObjectToInputArray, &input,
+ &min, &max, &nbins_in, NI_ObjectToOptionalInputArray,
+ &labels, &indices_object))
+ goto exit;
+ nbins = nbins_in;
+
+ if (!_NI_GetIndices(indices_object, &result_indices, &min_label,
+ &max_label, &n_results))
+ goto exit;
+
+ /* Set all pointers to NULL, so that freeing the memory */
+ /* doesn't cause problems. */
+ histograms = (PyArrayObject**)calloc(input->nd * n_results,
+ sizeof(PyArrayObject*));
+ if (!histograms) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(jj = 0; jj < n_results; jj++) {
+ histograms[jj] = NA_NewArray(NULL, tInt32, 1, &nbins);
+ if (!histograms[jj]) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+
+ if (!NI_Histogram(input, labels, min_label, max_label, result_indices,
+ n_results, histograms, min, max, nbins))
+ goto exit;
+
+ result = _NI_BuildMeasurementResultArrayObject(n_results, histograms);
+
+ exit:
+ Py_XDECREF(input);
+ Py_XDECREF(labels);
+ if (result_indices)
+ free(result_indices);
+ if (histograms) {
+ for(jj = 0; jj < n_results; jj++) {
+ Py_XDECREF(histograms[jj]);
+ }
+ free(histograms);
+ }
+ return result;
+}
+
+static PyObject *Py_DistanceTransformBruteForce(PyObject *obj,
+ PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *features = NULL;
+ PyArrayObject *sampling = NULL;
+ int metric;
+
+ if (!PyArg_ParseTuple(args, "O&iO&O&O&", NI_ObjectToInputArray, &input,
+ &metric, NI_ObjectToOptionalInputArray, &sampling,
+ NI_ObjectToOptionalOutputArray, &output,
+ NI_ObjectToOptionalOutputArray, &features))
+ goto exit;
+ if (!NI_DistanceTransformBruteForce(input, metric, sampling,
+ output, features))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(sampling);
+ Py_XDECREF(output);
+ Py_XDECREF(features);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_DistanceTransformOnePass(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *strct = NULL, *distances = NULL, *features = NULL;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&", NI_ObjectToInputArray, &strct,
+ NI_ObjectToIoArray, &distances,
+ NI_ObjectToOptionalOutputArray, &features))
+ goto exit;
+ if (!NI_DistanceTransformOnePass(strct, distances, features))
+ goto exit;
+exit:
+ Py_XDECREF(strct);
+ Py_XDECREF(distances);
+ Py_XDECREF(features);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyObject *Py_EuclideanFeatureTransform(PyObject *obj,
+ PyObject *args)
+{
+ PyArrayObject *input = NULL, *features = NULL, *sampling = NULL;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&", NI_ObjectToInputArray, &input,
+ NI_ObjectToOptionalInputArray, &sampling,
+ NI_ObjectToOutputArray, &features))
+ goto exit;
+ if (!NI_EuclideanFeatureTransform(input, sampling, features))
+ goto exit;
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(sampling);
+ Py_XDECREF(features);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static void _FreeCoordinateList(void* ptr)
+{
+ NI_FreeCoordinateList((NI_CoordinateList*)ptr);
+}
+
+static PyObject *Py_BinaryErosion(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *input = NULL, *output = NULL, *strct = NULL;
+ PyArrayObject *mask = NULL;
+ PyObject *cobj = NULL;
+ int border_value, invert, center_is_true;
+ int changed = 0, return_coordinates;
+ NI_CoordinateList *coordinate_list = NULL;
+ maybelong *origins = NULL;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&O&iO&iii", NI_ObjectToInputArray,
+ &input, NI_ObjectToInputArray, &strct,
+ NI_ObjectToOptionalInputArray, &mask,
+ NI_ObjectToOutputArray, &output, &border_value,
+ NI_ObjectToLongSequence, &origins, &invert,
+ ¢er_is_true, &return_coordinates))
+ goto exit;
+ if (!NI_BinaryErosion(input, strct, mask, output, border_value,
+ origins, invert, center_is_true, &changed,
+ return_coordinates ? &coordinate_list : NULL))
+ goto exit;
+ if (return_coordinates) {
+ cobj = PyCObject_FromVoidPtr(coordinate_list, _FreeCoordinateList);
+ }
+exit:
+ Py_XDECREF(input);
+ Py_XDECREF(strct);
+ Py_XDECREF(mask);
+ Py_XDECREF(output);
+ if (origins)
+ free(origins);
+ if (PyErr_Occurred()) {
+ Py_XDECREF(cobj);
+ return NULL;
+ } else {
+ if (return_coordinates) {
+ return Py_BuildValue("iN", changed, cobj);
+ } else {
+ return Py_BuildValue("i", changed);
+ }
+ }
+}
+
+static PyObject *Py_BinaryErosion2(PyObject *obj, PyObject *args)
+{
+ PyArrayObject *array = NULL, *strct = NULL, *mask = NULL;
+ PyObject *cobj = NULL;
+ int invert, niter;
+ maybelong *origins = NULL;
+
+ if (!PyArg_ParseTuple(args, "O&O&O&iO&iO", NI_ObjectToIoArray, &array,
+ NI_ObjectToInputArray, &strct, NI_ObjectToOptionalInputArray,
+ &mask, &niter, NI_ObjectToLongSequence, &origins, &invert,
+ &cobj))
+ goto exit;
+
+ if (PyCObject_Check(cobj)) {
+ NI_CoordinateList *cobj_data = PyCObject_AsVoidPtr(cobj);
+ if (!NI_BinaryErosion2(array, strct, mask, niter, origins, invert,
+ &cobj_data))
+ goto exit;
+ } else {
+ PyErr_SetString(PyExc_RuntimeError, "cannot convert CObject");
+ goto exit;
+ }
+exit:
+ Py_XDECREF(array);
+ Py_XDECREF(strct);
+ Py_XDECREF(mask);
+ if (origins) free(origins);
+ return PyErr_Occurred() ? NULL : Py_BuildValue("");
+}
+
+static PyMethodDef methods[] = {
+ {"correlate1d", (PyCFunction)Py_Correlate1D,
+ METH_VARARGS, NULL},
+ {"correlate", (PyCFunction)Py_Correlate,
+ METH_VARARGS, NULL},
+ {"uniform_filter1d", (PyCFunction)Py_UniformFilter1D,
+ METH_VARARGS, NULL},
+ {"min_or_max_filter1d", (PyCFunction)Py_MinOrMaxFilter1D,
+ METH_VARARGS, NULL},
+ {"min_or_max_filter", (PyCFunction)Py_MinOrMaxFilter,
+ METH_VARARGS, NULL},
+ {"rank_filter", (PyCFunction)Py_RankFilter,
+ METH_VARARGS, NULL},
+ {"generic_filter", (PyCFunction)Py_GenericFilter,
+ METH_VARARGS, NULL},
+ {"generic_filter1d", (PyCFunction)Py_GenericFilter1D,
+ METH_VARARGS, NULL},
+ {"fourier_filter", (PyCFunction)Py_FourierFilter,
+ METH_VARARGS, NULL},
+ {"fourier_shift", (PyCFunction)Py_FourierShift,
+ METH_VARARGS, NULL},
+ {"spline_filter1d", (PyCFunction)Py_SplineFilter1D,
+ METH_VARARGS, NULL},
+ {"geometric_transform", (PyCFunction)Py_GeometricTransform,
+ METH_VARARGS, NULL},
+ {"zoom_shift", (PyCFunction)Py_ZoomShift,
+ METH_VARARGS, NULL},
+ {"label", (PyCFunction)Py_Label,
+ METH_VARARGS, NULL},
+ {"find_objects", (PyCFunction)Py_FindObjects,
+ METH_VARARGS, NULL},
+ {"watershed_ift", (PyCFunction)Py_WatershedIFT,
+ METH_VARARGS, NULL},
+ {"statistics", (PyCFunction)Py_Statistics,
+ METH_VARARGS, NULL},
+ {"center_of_mass", (PyCFunction)Py_CenterOfMass,
+ METH_VARARGS, NULL},
+ {"histogram", (PyCFunction)Py_Histogram,
+ METH_VARARGS, NULL},
+ {"distance_transform_bf", (PyCFunction)Py_DistanceTransformBruteForce,
+ METH_VARARGS, NULL},
+ {"distance_transform_op", (PyCFunction)Py_DistanceTransformOnePass,
+ METH_VARARGS, NULL},
+ {"euclidean_feature_transform",
+ (PyCFunction)Py_EuclideanFeatureTransform,
+ METH_VARARGS, NULL},
+ {"binary_erosion", (PyCFunction)Py_BinaryErosion,
+ METH_VARARGS, NULL},
+ {"binary_erosion2", (PyCFunction)Py_BinaryErosion2,
+ METH_VARARGS, NULL},
+ {NULL, NULL, 0, NULL}
+};
+
+PyMODINIT_FUNC init_nd_image(void)
+{
+ Py_InitModule("_nd_image", methods);
+ import_array();
+}
Added: branches/Interpolate1D/ndimage/nd_image.h
===================================================================
--- branches/Interpolate1D/ndimage/nd_image.h 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/nd_image.h 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,278 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+ * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef ND_IMAGE_H
+#define ND_IMAGE_H
+
+#include "Python.h"
+
+#ifndef ND_IMPORT_ARRAY
+#define NO_IMPORT_ARRAY
+#endif
+
+#include <numpy/noprefix.h>
+#undef NO_IMPORT_ARRAY
+
+/* Eventually get rid of everything below this line */
+
+typedef enum
+{
+ tAny=-1,
+ tBool=PyArray_BOOL,
+ tInt8=PyArray_INT8,
+ tUInt8=PyArray_UINT8,
+ tInt16=PyArray_INT16,
+ tUInt16=PyArray_UINT16,
+ tInt32=PyArray_INT32,
+ tUInt32=PyArray_UINT32,
+ tInt64=PyArray_INT64,
+ tUInt64=PyArray_UINT64,
+ tFloat32=PyArray_FLOAT32,
+ tFloat64=PyArray_FLOAT64,
+ tComplex64=PyArray_COMPLEX64,
+ tComplex128=PyArray_COMPLEX128,
+ tObject=PyArray_OBJECT, /* placeholder... does nothing */
+ tMaxType=PyArray_NTYPES,
+ tDefault=PyArray_FLOAT64,
+ tLong=PyArray_LONG,
+} NumarrayType;
+
+#define NI_MAXDIM NPY_MAXDIMS
+
+typedef npy_intp maybelong;
+#define MAXDIM NPY_MAXDIMS
+
+#define HAS_UINT64 1
+
+
+
+#ifdef ND_IMPORT_ARRAY
+
+/* Numarray Helper Functions */
+
+static PyArrayObject*
+NA_InputArray(PyObject *a, NumarrayType t, int requires)
+{
+ PyArray_Descr *descr;
+ if (t == tAny) descr = NULL;
+ else descr = PyArray_DescrFromType(t);
+ return (PyArrayObject *) \
+ PyArray_CheckFromAny(a, descr, 0, 0, requires, NULL);
+}
+
+/* satisfies ensures that 'a' meets a set of requirements and matches
+the specified type.
+*/
+static int
+satisfies(PyArrayObject *a, int requirements, NumarrayType t)
+{
+ int type_ok = (a->descr->type_num == t) || (t == tAny);
+
+ if (PyArray_ISCARRAY(a))
+ return type_ok;
+ if (PyArray_ISBYTESWAPPED(a) && (requirements & NPY_NOTSWAPPED))
+ return 0;
+ if (!PyArray_ISALIGNED(a) && (requirements & NPY_ALIGNED))
+ return 0;
+ if (!PyArray_ISCONTIGUOUS(a) && (requirements & NPY_CONTIGUOUS))
+ return 0;
+ if (!PyArray_ISWRITEABLE(a) && (requirements & NPY_WRITEABLE))
+ return 0;
+ if (requirements & NPY_ENSURECOPY)
+ return 0;
+ return type_ok;
+}
+
+static PyArrayObject *
+NA_OutputArray(PyObject *a, NumarrayType t, int requires)
+{
+ PyArray_Descr *dtype;
+ PyArrayObject *ret;
+
+ if (!PyArray_Check(a) || !PyArray_ISWRITEABLE(a)) {
+ PyErr_Format(PyExc_TypeError,
+ "NA_OutputArray: only writeable arrays work for output.");
+ return NULL;
+ }
+
+ if (satisfies((PyArrayObject *)a, requires, t)) {
+ Py_INCREF(a);
+ return (PyArrayObject *)a;
+ }
+ if (t == tAny) {
+ dtype = PyArray_DESCR(a);
+ Py_INCREF(dtype);
+ }
+ else {
+ dtype = PyArray_DescrFromType(t);
+ }
+ ret = (PyArrayObject *)PyArray_Empty(PyArray_NDIM(a), PyArray_DIMS(a),
+ dtype, 0);
+ ret->flags |= NPY_UPDATEIFCOPY;
+ ret->base = a;
+ PyArray_FLAGS(a) &= ~NPY_WRITEABLE;
+ Py_INCREF(a);
+ return ret;
+}
+
+/* NA_IoArray is a combination of NA_InputArray and NA_OutputArray.
+
+Unlike NA_OutputArray, if a temporary is required it is initialized to a copy
+of the input array.
+
+Unlike NA_InputArray, deallocating any resulting temporary array results in a
+copy from the temporary back to the original.
+*/
+static PyArrayObject *
+NA_IoArray(PyObject *a, NumarrayType t, int requires)
+{
+ PyArrayObject *shadow = NA_InputArray(a, t, requires | NPY_UPDATEIFCOPY );
+
+ if (!shadow) return NULL;
+
+ /* Guard against non-writable, but otherwise satisfying requires.
+ In this case, shadow == a.
+ */
+ if (!PyArray_ISWRITEABLE(shadow)) {
+ PyErr_Format(PyExc_TypeError,
+ "NA_IoArray: I/O array must be writable array");
+ PyArray_XDECREF_ERR(shadow);
+ return NULL;
+ }
+
+ return shadow;
+}
+
+#define NUM_LITTLE_ENDIAN 0
+#define NUM_BIG_ENDIAN 1
+
+static int
+NA_ByteOrder(void)
+{
+ unsigned long byteorder_test;
+ byteorder_test = 1;
+ if (*((char *) &byteorder_test))
+ return NUM_LITTLE_ENDIAN;
+ else
+ return NUM_BIG_ENDIAN;
+}
+
+/* ignores bytestride */
+static PyArrayObject *
+NA_NewAllFromBuffer(int ndim, maybelong *shape, NumarrayType type,
+ PyObject *bufferObject, maybelong byteoffset, maybelong bytestride,
+ int byteorder, int aligned, int writeable)
+{
+ PyArrayObject *self = NULL;
+ PyArray_Descr *dtype;
+
+ if (type == tAny)
+ type = tDefault;
+
+ dtype = PyArray_DescrFromType(type);
+ if (dtype == NULL) return NULL;
+
+ if (byteorder != NA_ByteOrder()) {
+ PyArray_Descr *temp;
+ temp = PyArray_DescrNewByteorder(dtype, PyArray_SWAP);
+ Py_DECREF(dtype);
+ if (temp == NULL) return NULL;
+ dtype = temp;
+ }
+
+ if (bufferObject == Py_None || bufferObject == NULL) {
+ self = (PyArrayObject *) \
+ PyArray_NewFromDescr(&PyArray_Type, dtype,
+ ndim, shape, NULL, NULL,
+ 0, NULL);
+ }
+ else {
+ npy_intp size = 1;
+ int i;
+ PyArrayObject *newself;
+ PyArray_Dims newdims;
+ for(i=0; i<ndim; i++) {
+ size *= shape[i];
+ }
+ self = (PyArrayObject *) \
+ PyArray_FromBuffer(bufferObject, dtype,
+ size, byteoffset);
+ if (self == NULL) return self;
+ newdims.len = ndim;
+ newdims.ptr = shape;
+ newself = (PyArrayObject *) \
+ PyArray_Newshape(self, &newdims, PyArray_CORDER);
+ Py_DECREF(self);
+ self = newself;
+ }
+
+ return self;
+}
+
+static PyArrayObject *
+NA_NewAll(int ndim, maybelong *shape, NumarrayType type,
+ void *buffer, maybelong byteoffset, maybelong bytestride,
+ int byteorder, int aligned, int writeable)
+{
+ PyArrayObject *result = NA_NewAllFromBuffer(
+ ndim, shape, type, Py_None,
+ byteoffset, bytestride,
+ byteorder, aligned, writeable);
+ if (result) {
+ if (!PyArray_Check((PyObject *) result)) {
+ PyErr_Format( PyExc_TypeError,
+ "NA_NewAll: non-NumArray result");
+ result = NULL;
+ } else {
+ if (buffer) {
+ memcpy(result->data, buffer, PyArray_NBYTES(result));
+ } else {
+ memset(result->data, 0, PyArray_NBYTES(result));
+ }
+ }
+ }
+ return result;
+}
+
+/* Create a new numarray which is initially a C_array, or which
+references a C_array: aligned, !byteswapped, contiguous, ...
+Call with buffer==NULL to allocate storage.
+*/
+static PyArrayObject *
+NA_NewArray(void *buffer, NumarrayType type, int ndim, maybelong *shape)
+{
+ return (PyArrayObject *) NA_NewAll(ndim, shape, type, buffer, 0, 0,
+ NA_ByteOrder(), 1, 1);
+}
+
+#endif /* ND_IMPORT_ARRAY */
+
+#endif /* ND_IMAGE_H */
Added: branches/Interpolate1D/ndimage/ni_filters.c
===================================================================
--- branches/Interpolate1D/ndimage/ni_filters.c 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/ni_filters.c 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,888 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+ * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "ni_support.h"
+#include "ni_filters.h"
+#include <stdlib.h>
+#include <math.h>
+
+#define BUFFER_SIZE 256000
+
+int NI_Correlate1D(PyArrayObject *input, PyArrayObject *weights,
+ int axis, PyArrayObject *output, NI_ExtendMode mode,
+ double cval, maybelong origin)
+{
+ int symmetric = 0, ii, jj, more;
+ maybelong ll, lines, length, size1, size2, filter_size;
+ double *ibuffer = NULL, *obuffer = NULL;
+ Float64 *fw;
+ NI_LineBuffer iline_buffer, oline_buffer;
+
+ /* test for symmetry or anti-symmetry: */
+ filter_size = weights->dimensions[0];
+ size1 = filter_size / 2;
+ size2 = filter_size - size1 - 1;
+ fw = (void *)PyArray_DATA(weights);
+ if (filter_size & 0x1) {
+ symmetric = 1;
+ for(ii = 1; ii <= filter_size / 2; ii++) {
+ if (fabs(fw[ii + size1] - fw[size1 - ii]) > DBL_EPSILON) {
+ symmetric = 0;
+ break;
+ }
+ }
+ if (symmetric == 0) {
+ symmetric = -1;
+ for(ii = 1; ii <= filter_size / 2; ii++) {
+ if (fabs(fw[size1 + ii] + fw[size1 - ii]) > DBL_EPSILON) {
+ symmetric = 0;
+ break;
+ }
+ }
+ }
+ }
+ /* allocate and initialize the line buffers: */
+ lines = -1;
+ if (!NI_AllocateLineBuffer(input, axis, size1 + origin, size2 - origin,
+ &lines, BUFFER_SIZE, &ibuffer))
+ goto exit;
+ if (!NI_AllocateLineBuffer(output, axis, 0, 0, &lines, BUFFER_SIZE,
+ &obuffer))
+ goto exit;
+ if (!NI_InitLineBuffer(input, axis, size1 + origin, size2 - origin,
+ lines, ibuffer, mode, cval, &iline_buffer))
+ goto exit;
+ if (!NI_InitLineBuffer(output, axis, 0, 0, lines, obuffer, mode, 0.0,
+ &oline_buffer))
+ goto exit;
+ length = input->nd > 0 ? input->dimensions[axis] : 1;
+ fw += size1;
+ /* iterate over all the array lines: */
+ do {
+ /* copy lines from array to buffer: */
+ if (!NI_ArrayToLineBuffer(&iline_buffer, &lines, &more))
+ goto exit;
+ /* iterate over the lines in the buffers: */
+ for(ii = 0; ii < lines; ii++) {
+ /* get lines: */
+ double *iline = NI_GET_LINE(iline_buffer, ii) + size1;
+ double *oline = NI_GET_LINE(oline_buffer, ii);
+ /* the correlation calculation: */
+ if (symmetric > 0) {
+ for(ll = 0; ll < length; ll++) {
+ oline[ll] = iline[0] * fw[0];
+ for(jj = -size1 ; jj < 0; jj++)
+ oline[ll] += (iline[jj] + iline[-jj]) * fw[jj];
+ ++iline;
+ }
+ } else if (symmetric < 0) {
+ for(ll = 0; ll < length; ll++) {
+ oline[ll] = iline[0] * fw[0];
+ for(jj = -size1 ; jj < 0; jj++)
+ oline[ll] += (iline[jj] - iline[-jj]) * fw[jj];
+ ++iline;
+ }
+ } else {
+ for(ll = 0; ll < length; ll++) {
+ oline[ll] = iline[size2] * fw[size2];
+ for(jj = -size1; jj < size2; jj++)
+ oline[ll] += iline[jj] * fw[jj];
+ ++iline;
+ }
+ }
+ }
+ /* copy lines from buffer to array: */
+ if (!NI_LineBufferToArray(&oline_buffer))
+ goto exit;
+ } while(more);
+exit:
+ if (ibuffer) free(ibuffer);
+ if (obuffer) free(obuffer);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+#define CASE_CORRELATE_POINT(_pi, _weights, _offsets, _filter_size, \
+ _cvalue, _type, _res, _mv) \
+case t ## _type: \
+{ \
+ maybelong _ii, _offset; \
+ for(_ii = 0; _ii < _filter_size; _ii++) { \
+ _offset = _offsets[_ii]; \
+ if (_offset == _mv) \
+ _res += _weights[_ii] * _cvalue; \
+ else \
+ _res += _weights[_ii] * (double)*(_type*)(_pi + _offset); \
+ } \
+} \
+break
+
+#define CASE_FILTER_OUT(_po, _tmp, _type) \
+case t ## _type: \
+ *(_type*)_po = (_type)_tmp; \
+ break
+
+int NI_Correlate(PyArrayObject* input, PyArrayObject* weights,
+ PyArrayObject* output, NI_ExtendMode mode,
+ double cvalue, maybelong *origins)
+{
+ Bool *pf = NULL;
+ maybelong fsize, jj, kk, filter_size = 0, border_flag_value;
+ maybelong *offsets = NULL, *oo, size;
+ NI_FilterIterator fi;
+ NI_Iterator ii, io;
+ char *pi, *po;
+ Float64 *pw;
+ Float64 *ww = NULL;
+ int ll;
+
+ /* get the the footprint: */
+ fsize = 1;
+ for(ll = 0; ll < weights->nd; ll++)
+ fsize *= weights->dimensions[ll];
+ pw = (Float64*)PyArray_DATA(weights);
+ pf = (Bool*)malloc(fsize * sizeof(Bool));
+ if (!pf) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(jj = 0; jj < fsize; jj++) {
+ if (fabs(pw[jj]) > DBL_EPSILON) {
+ pf[jj] = 1;
+ ++filter_size;
+ } else {
+ pf[jj] = 0;
+ }
+ }
+ /* copy the weights to contiguous memory: */
+ ww = (Float64*)malloc(filter_size * sizeof(Float64));
+ if (!ww) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ jj = 0;
+ for(kk = 0; kk < fsize; kk++) {
+ if (pf[kk]) {
+ ww[jj++] = pw[kk];
+ }
+ }
+ /* initialize filter offsets: */
+ if (!NI_InitFilterOffsets(input, pf, weights->dimensions, origins,
+ mode, &offsets, &border_flag_value, NULL))
+ goto exit;
+ /* initialize filter iterator: */
+ if (!NI_InitFilterIterator(input->nd, weights->dimensions, filter_size,
+ input->dimensions, origins, &fi))
+ goto exit;
+ /* initialize input element iterator: */
+ if (!NI_InitPointIterator(input, &ii))
+ goto exit;
+ /* initialize output element iterator: */
+ if (!NI_InitPointIterator(output, &io))
+ goto exit;
+ /* get data pointers an array size: */
+ pi = (void *)PyArray_DATA(input);
+ po = (void *)PyArray_DATA(output);
+ size = 1;
+ for(ll = 0; ll < input->nd; ll++)
+ size *= input->dimensions[ll];
+ /* iterator over the elements: */
+ oo = offsets;
+ for(jj = 0; jj < size; jj++) {
+ double tmp = 0.0;
+ switch (input->descr->type_num) {
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Bool,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, UInt8,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, UInt16,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, UInt32,
+ tmp, border_flag_value);
+#if HAS_UINT64
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, UInt64,
+ tmp, border_flag_value);
+#endif
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Int8,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Int16,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Int32,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Int64,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Float32,
+ tmp, border_flag_value);
+ CASE_CORRELATE_POINT(pi, ww, oo, filter_size, cvalue, Float64,
+ tmp, border_flag_value);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FILTER_OUT(po, tmp, Bool);
+ CASE_FILTER_OUT(po, tmp, UInt8);
+ CASE_FILTER_OUT(po, tmp, UInt16);
+ CASE_FILTER_OUT(po, tmp, UInt32);
+#if HAS_UINT64
+ CASE_FILTER_OUT(po, tmp, UInt64);
+#endif
+ CASE_FILTER_OUT(po, tmp, Int8);
+ CASE_FILTER_OUT(po, tmp, Int16);
+ CASE_FILTER_OUT(po, tmp, Int32);
+ CASE_FILTER_OUT(po, tmp, Int64);
+ CASE_FILTER_OUT(po, tmp, Float32);
+ CASE_FILTER_OUT(po, tmp, Float64);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ NI_FILTER_NEXT2(fi, ii, io, oo, pi, po);
+ }
+exit:
+ if (offsets) free(offsets);
+ if (ww) free(ww);
+ if (pf) free(pf);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+int
+NI_UniformFilter1D(PyArrayObject *input, long filter_size,
+ int axis, PyArrayObject *output, NI_ExtendMode mode,
+ double cval, long origin)
+{
+ maybelong lines, kk, ll, length, size1, size2;
+ int more;
+ double *ibuffer = NULL, *obuffer = NULL;
+ NI_LineBuffer iline_buffer, oline_buffer;
+
+ size1 = filter_size / 2;
+ size2 = filter_size - size1 - 1;
+ /* allocate and initialize the line buffers: */
+ lines = -1;
+ if (!NI_AllocateLineBuffer(input, axis, size1 + origin, size2 - origin,
+ &lines, BUFFER_SIZE, &ibuffer))
+ goto exit;
+ if (!NI_AllocateLineBuffer(output, axis, 0, 0, &lines, BUFFER_SIZE,
+ &obuffer))
+ goto exit;
+ if (!NI_InitLineBuffer(input, axis, size1 + origin, size2 - origin,
+ lines, ibuffer, mode, cval, &iline_buffer))
+ goto exit;
+ if (!NI_InitLineBuffer(output, axis, 0, 0, lines, obuffer, mode, 0.0,
+ &oline_buffer))
+ goto exit;
+ length = input->nd > 0 ? input->dimensions[axis] : 1;
+
+ /* iterate over all the array lines: */
+ do {
+ /* copy lines from array to buffer: */
+ if (!NI_ArrayToLineBuffer(&iline_buffer, &lines, &more))
+ goto exit;
+ /* iterate over the lines in the buffers: */
+ for(kk = 0; kk < lines; kk++) {
+ /* get lines: */
+ double *iline = NI_GET_LINE(iline_buffer, kk);
+ double *oline = NI_GET_LINE(oline_buffer, kk);
+ /* do the uniform filter: */
+ double tmp = 0.0;
+ double *l1 = iline;
+ double *l2 = iline + filter_size;
+ for(ll = 0; ll < filter_size; ll++)
+ tmp += iline[ll];
+ tmp /= (double)filter_size;
+ oline[0] = tmp;
+ for(ll = 1; ll < length; ll++) {
+ tmp += (*l2++ - *l1++) / (double)filter_size;
+ oline[ll] = tmp;
+ }
+ }
+ /* copy lines from buffer to array: */
+ if (!NI_LineBufferToArray(&oline_buffer))
+ goto exit;
+ } while(more);
+
+ exit:
+ if (ibuffer) free(ibuffer);
+ if (obuffer) free(obuffer);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+int
+NI_MinOrMaxFilter1D(PyArrayObject *input, long filter_size,
+ int axis, PyArrayObject *output, NI_ExtendMode mode,
+ double cval, long origin, int minimum)
+{
+ maybelong lines, kk, jj, ll, length, size1, size2;
+ int more;
+ double *ibuffer = NULL, *obuffer = NULL;
+ NI_LineBuffer iline_buffer, oline_buffer;
+
+ size1 = filter_size / 2;
+ size2 = filter_size - size1 - 1;
+ /* allocate and initialize the line buffers: */
+ lines = -1;
+ if (!NI_AllocateLineBuffer(input, axis, size1 + origin, size2 - origin,
+ &lines, BUFFER_SIZE, &ibuffer))
+ goto exit;
+ if (!NI_AllocateLineBuffer(output, axis, 0, 0, &lines, BUFFER_SIZE,
+ &obuffer))
+ goto exit;
+ if (!NI_InitLineBuffer(input, axis, size1 + origin, size2 - origin,
+ lines, ibuffer, mode, cval, &iline_buffer))
+ goto exit;
+ if (!NI_InitLineBuffer(output, axis, 0, 0, lines, obuffer, mode, 0.0,
+ &oline_buffer))
+ goto exit;
+ length = input->nd > 0 ? input->dimensions[axis] : 1;
+
+ /* iterate over all the array lines: */
+ do {
+ /* copy lines from array to buffer: */
+ if (!NI_ArrayToLineBuffer(&iline_buffer, &lines, &more))
+ goto exit;
+ /* iterate over the lines in the buffers: */
+ for(kk = 0; kk < lines; kk++) {
+ /* get lines: */
+ double *iline = NI_GET_LINE(iline_buffer, kk) + size1;
+ double *oline = NI_GET_LINE(oline_buffer, kk);
+ for(ll = 0; ll < length; ll++) {
+ /* find minimum or maximum filter: */
+ double val = iline[ll - size1];
+ for(jj = -size1 + 1; jj <= size2; jj++) {
+ double tmp = iline[ll + jj];
+ if (minimum) {
+ if (tmp < val)
+ val = tmp;
+ } else {
+ if (tmp > val)
+ val = tmp;
+ }
+ }
+ oline[ll] = val;
+ }
+ }
+ /* copy lines from buffer to array: */
+ if (!NI_LineBufferToArray(&oline_buffer))
+ goto exit;
+ } while(more);
+
+ exit:
+ if (ibuffer) free(ibuffer);
+ if (obuffer) free(obuffer);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+
+#define CASE_MIN_OR_MAX_POINT(_pi, _offsets, _filter_size, _cval, \
+ _type, _minimum, _res, _mv, _ss) \
+case t ## _type: \
+{ \
+ maybelong _ii, _oo = *_offsets; \
+ _type _cv = (_type)_cval, _tmp; \
+ _res = _oo == _mv ? _cv : *(_type*)(_pi + _oo); \
+ if (_ss) \
+ _res += *_ss; \
+ for(_ii = 1; _ii < _filter_size; _ii++) { \
+ _oo = _offsets[_ii]; \
+ _tmp = _oo == _mv ? _cv : *(_type*)(_pi + _oo); \
+ if (_ss) \
+ _tmp += _ss[_ii]; \
+ if (_minimum) { \
+ if (_tmp < _res) \
+ _res = (_type)_tmp; \
+ } else { \
+ if (_tmp > _res) \
+ _res = (_type)_tmp; \
+ } \
+ } \
+} \
+break
+
+int NI_MinOrMaxFilter(PyArrayObject* input, PyArrayObject* footprint,
+ PyArrayObject* structure, PyArrayObject* output,
+ NI_ExtendMode mode, double cvalue, maybelong *origins, int minimum)
+{
+ Bool *pf = NULL;
+ maybelong fsize, jj, kk, filter_size = 0, border_flag_value;
+ maybelong *offsets = NULL, *oo, size;
+ NI_FilterIterator fi;
+ NI_Iterator ii, io;
+ char *pi, *po;
+ int ll;
+ double *ss = NULL;
+ Float64 *ps;
+
+ /* get the the footprint: */
+ fsize = 1;
+ for(ll = 0; ll < footprint->nd; ll++)
+ fsize *= footprint->dimensions[ll];
+ pf = (Bool*)PyArray_DATA(footprint);
+ for(jj = 0; jj < fsize; jj++) {
+ if (pf[jj]) {
+ ++filter_size;
+ }
+ }
+ /* get the structure: */
+ if (structure) {
+ ss = (double*)malloc(filter_size * sizeof(double));
+ if (!ss) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ /* copy the weights to contiguous memory: */
+ ps = (Float64*)PyArray_DATA(structure);
+ jj = 0;
+ for(kk = 0; kk < fsize; kk++)
+ if (pf[kk])
+ ss[jj++] = minimum ? -ps[kk] : ps[kk];
+ }
+ /* initialize filter offsets: */
+ if (!NI_InitFilterOffsets(input, pf, footprint->dimensions, origins,
+ mode, &offsets, &border_flag_value, NULL))
+ goto exit;
+ /* initialize filter iterator: */
+ if (!NI_InitFilterIterator(input->nd, footprint->dimensions,
+ filter_size, input->dimensions, origins, &fi))
+ goto exit;
+ /* initialize input element iterator: */
+ if (!NI_InitPointIterator(input, &ii))
+ goto exit;
+ /* initialize output element iterator: */
+ if (!NI_InitPointIterator(output, &io))
+ goto exit;
+ /* get data pointers an array size: */
+ pi = (void *)PyArray_DATA(input);
+ po = (void *)PyArray_DATA(output);
+ size = 1;
+ for(ll = 0; ll < input->nd; ll++)
+ size *= input->dimensions[ll];
+ /* iterator over the elements: */
+ oo = offsets;
+ for(jj = 0; jj < size; jj++) {
+ double tmp = 0.0;
+ switch (input->descr->type_num) {
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Bool,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, UInt8,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, UInt16,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, UInt32,
+ minimum, tmp, border_flag_value, ss);
+#if HAS_UINT64
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, UInt64,
+ minimum, tmp, border_flag_value, ss);
+#endif
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Int8,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Int16,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Int32,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Int64,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Float32,
+ minimum, tmp, border_flag_value, ss);
+ CASE_MIN_OR_MAX_POINT(pi, oo, filter_size, cvalue, Float64,
+ minimum, tmp, border_flag_value, ss);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FILTER_OUT(po, tmp, Bool);
+ CASE_FILTER_OUT(po, tmp, UInt8);
+ CASE_FILTER_OUT(po, tmp, UInt16);
+ CASE_FILTER_OUT(po, tmp, UInt32);
+#if HAS_UINT64
+ CASE_FILTER_OUT(po, tmp, UInt64);
+#endif
+ CASE_FILTER_OUT(po, tmp, Int8);
+ CASE_FILTER_OUT(po, tmp, Int16);
+ CASE_FILTER_OUT(po, tmp, Int32);
+ CASE_FILTER_OUT(po, tmp, Int64);
+ CASE_FILTER_OUT(po, tmp, Float32);
+ CASE_FILTER_OUT(po, tmp, Float64);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ NI_FILTER_NEXT2(fi, ii, io, oo, pi, po);
+ }
+exit:
+ if (offsets) free(offsets);
+ if (ss) free(ss);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+static double NI_Select(double *buffer, int min, int max, int rank)
+{
+ int ii, jj;
+ double x, t;
+
+ if (min == max)
+ return buffer[min];
+
+ x = buffer[min];
+ ii = min - 1;
+ jj = max + 1;
+ for(;;) {
+ do
+ jj--;
+ while(buffer[jj] > x);
+ do
+ ii++;
+ while(buffer[ii] < x);
+ if (ii < jj) {
+ t = buffer[ii];
+ buffer[ii] = buffer[jj];
+ buffer[jj] = t;
+ } else {
+ break;
+ }
+ }
+
+ ii = jj - min + 1;
+ if (rank < ii)
+ return NI_Select(buffer, min, jj, rank);
+ else
+ return NI_Select(buffer, jj + 1, max, rank - ii);
+}
+
+#define CASE_RANK_POINT(_pi, _offsets, _filter_size, _cval, _type, \
+ _rank, _buffer, _res, _mv) \
+case t ## _type: \
+{ \
+ maybelong _ii; \
+ for(_ii = 0; _ii < _filter_size; _ii++) { \
+ maybelong _offset = _offsets[_ii]; \
+ if (_offset == _mv) \
+ _buffer[_ii] = (_type)_cval; \
+ else \
+ _buffer[_ii] = *(_type*)(_pi + _offsets[_ii]); \
+ } \
+ _res = (_type)NI_Select(_buffer, 0, _filter_size - 1, _rank); \
+} \
+break
+
+int NI_RankFilter(PyArrayObject* input, int rank,
+ PyArrayObject* footprint, PyArrayObject* output,
+ NI_ExtendMode mode, double cvalue, maybelong *origins)
+{
+ maybelong fsize, jj, filter_size = 0, border_flag_value;
+ maybelong *offsets = NULL, *oo, size;
+ NI_FilterIterator fi;
+ NI_Iterator ii, io;
+ char *pi, *po;
+ Bool *pf = NULL;
+ double *buffer = NULL;
+ int ll;
+
+ /* get the the footprint: */
+ fsize = 1;
+ for(ll = 0; ll < footprint->nd; ll++)
+ fsize *= footprint->dimensions[ll];
+ pf = (Bool*)PyArray_DATA(footprint);
+ for(jj = 0; jj < fsize; jj++) {
+ if (pf[jj]) {
+ ++filter_size;
+ }
+ }
+ /* buffer for rank calculation: */
+ buffer = (double*)malloc(filter_size * sizeof(double));
+ if (!buffer) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ /* iterator over the elements: */
+ oo = offsets;
+ /* initialize filter offsets: */
+ if (!NI_InitFilterOffsets(input, pf, footprint->dimensions, origins,
+ mode, &offsets, &border_flag_value, NULL))
+ goto exit;
+ /* initialize filter iterator: */
+ if (!NI_InitFilterIterator(input->nd, footprint->dimensions,
+ filter_size, input->dimensions, origins, &fi))
+ goto exit;
+ /* initialize input element iterator: */
+ if (!NI_InitPointIterator(input, &ii))
+ goto exit;
+ /* initialize output element iterator: */
+ if (!NI_InitPointIterator(output, &io))
+ goto exit;
+ /* get data pointers an array size: */
+ pi = (void *)PyArray_DATA(input);
+ po = (void *)PyArray_DATA(output);
+ size = 1;
+ for(ll = 0; ll < input->nd; ll++)
+ size *= input->dimensions[ll];
+ /* iterator over the elements: */
+ oo = offsets;
+ for(jj = 0; jj < size; jj++) {
+ double tmp = 0.0;
+ switch (input->descr->type_num) {
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Bool,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, UInt8,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, UInt16,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, UInt32,
+ rank, buffer, tmp, border_flag_value);
+#if HAS_UINT64
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, UInt64,
+ rank, buffer, tmp, border_flag_value);
+#endif
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Int8,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Int16,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Int32,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Int64,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Float32,
+ rank, buffer, tmp, border_flag_value);
+ CASE_RANK_POINT(pi, oo, filter_size, cvalue, Float64,
+ rank, buffer, tmp, border_flag_value);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FILTER_OUT(po, tmp, Bool);
+ CASE_FILTER_OUT(po, tmp, UInt8);
+ CASE_FILTER_OUT(po, tmp, UInt16);
+ CASE_FILTER_OUT(po, tmp, UInt32);
+#if HAS_UINT64
+ CASE_FILTER_OUT(po, tmp, UInt64);
+#endif
+ CASE_FILTER_OUT(po, tmp, Int8);
+ CASE_FILTER_OUT(po, tmp, Int16);
+ CASE_FILTER_OUT(po, tmp, Int32);
+ CASE_FILTER_OUT(po, tmp, Int64);
+ CASE_FILTER_OUT(po, tmp, Float32);
+ CASE_FILTER_OUT(po, tmp, Float64);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ NI_FILTER_NEXT2(fi, ii, io, oo, pi, po);
+ }
+exit:
+ if (offsets) free(offsets);
+ if (buffer) free(buffer);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+int NI_GenericFilter1D(PyArrayObject *input,
+ int (*function)(double*, maybelong, double*, maybelong, void*),
+ void* data, long filter_size, int axis, PyArrayObject *output,
+ NI_ExtendMode mode, double cval, long origin)
+{
+ int more;
+ maybelong ii, lines, length, size1, size2;
+ double *ibuffer = NULL, *obuffer = NULL;
+ NI_LineBuffer iline_buffer, oline_buffer;
+
+ /* allocate and initialize the line buffers: */
+ size1 = filter_size / 2;
+ size2 = filter_size - size1 - 1;
+ lines = -1;
+ if (!NI_AllocateLineBuffer(input, axis, size1 + origin, size2 - origin,
+ &lines, BUFFER_SIZE, &ibuffer))
+ goto exit;
+ if (!NI_AllocateLineBuffer(output, axis, 0, 0, &lines, BUFFER_SIZE,
+ &obuffer))
+ goto exit;
+ if (!NI_InitLineBuffer(input, axis, size1 + origin, size2 - origin,
+ lines, ibuffer, mode, cval, &iline_buffer))
+ goto exit;
+ if (!NI_InitLineBuffer(output, axis, 0, 0, lines, obuffer, mode, 0.0,
+ &oline_buffer))
+ goto exit;
+ length = input->nd > 0 ? input->dimensions[axis] : 1;
+ /* iterate over all the array lines: */
+ do {
+ /* copy lines from array to buffer: */
+ if (!NI_ArrayToLineBuffer(&iline_buffer, &lines, &more))
+ goto exit;
+ /* iterate over the lines in the buffers: */
+ for(ii = 0; ii < lines; ii++) {
+ /* get lines: */
+ double *iline = NI_GET_LINE(iline_buffer, ii);
+ double *oline = NI_GET_LINE(oline_buffer, ii);
+ if (!function(iline, length + size1 + size2, oline, length, data)) {
+ if (!PyErr_Occurred())
+ PyErr_SetString(PyExc_RuntimeError,
+ "unknown error in line processing function");
+ goto exit;
+ }
+ }
+ /* copy lines from buffer to array: */
+ if (!NI_LineBufferToArray(&oline_buffer))
+ goto exit;
+ } while(more);
+exit:
+ if (ibuffer) free(ibuffer);
+ if (obuffer) free(obuffer);
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+#define CASE_FILTER_POINT(_pi, _offsets, _filter_size, _cvalue, _type, \
+ _res, _mv, _function, _data, _buffer) \
+case t ## _type: \
+{ \
+ maybelong _ii, _offset; \
+ for(_ii = 0; _ii < _filter_size; _ii++) { \
+ _offset = _offsets[_ii]; \
+ if (_offset == _mv) \
+ _buffer[_ii] = (double)_cvalue; \
+ else \
+ _buffer[_ii] = (double)*(_type*)(_pi + _offset); \
+ } \
+ if (!_function(_buffer, _filter_size, &_res, _data)) { \
+ if (!PyErr_Occurred()) \
+ PyErr_SetString(PyExc_RuntimeError, \
+ "unknown error in filter function"); \
+ goto exit; \
+ } \
+} \
+break
+
+
+int NI_GenericFilter(PyArrayObject* input,
+ int (*function)(double*, maybelong, double*, void*), void *data,
+ PyArrayObject* footprint, PyArrayObject* output,
+ NI_ExtendMode mode, double cvalue, maybelong *origins)
+{
+ Bool *pf = NULL;
+ maybelong fsize, jj, filter_size = 0, border_flag_value;
+ maybelong *offsets = NULL, *oo, size;
+ NI_FilterIterator fi;
+ NI_Iterator ii, io;
+ char *pi, *po;
+ double *buffer = NULL;
+ int ll;
+
+ /* get the the footprint: */
+ fsize = 1;
+ for(ll = 0; ll < footprint->nd; ll++)
+ fsize *= footprint->dimensions[ll];
+ pf = (Bool*)PyArray_DATA(footprint);
+ for(jj = 0; jj < fsize; jj++) {
+ if (pf[jj])
+ ++filter_size;
+ }
+ /* initialize filter offsets: */
+ if (!NI_InitFilterOffsets(input, pf, footprint->dimensions, origins,
+ mode, &offsets, &border_flag_value, NULL))
+ goto exit;
+ /* initialize filter iterator: */
+ if (!NI_InitFilterIterator(input->nd, footprint->dimensions,
+ filter_size, input->dimensions, origins, &fi))
+ goto exit;
+ /* initialize input element iterator: */
+ if (!NI_InitPointIterator(input, &ii))
+ goto exit;
+ /* initialize output element iterator: */
+ if (!NI_InitPointIterator(output, &io))
+ goto exit;
+ /* get data pointers an array size: */
+ pi = (void *)PyArray_DATA(input);
+ po = (void *)PyArray_DATA(output);
+ size = 1;
+ for(ll = 0; ll < input->nd; ll++)
+ size *= input->dimensions[ll];
+ /* buffer for filter calculation: */
+ buffer = (double*)malloc(filter_size * sizeof(double));
+ if (!buffer) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ /* iterate over the elements: */
+ oo = offsets;
+ for(jj = 0; jj < size; jj++) {
+ double tmp = 0.0;
+ switch (input->descr->type_num) {
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Bool,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, UInt8,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, UInt16,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, UInt32,
+ tmp, border_flag_value, function, data, buffer);
+#if HAS_UINT64
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, UInt64,
+ tmp, border_flag_value, function, data, buffer);
+#endif
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Int8,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Int16,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Int32,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Int64,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Float32,
+ tmp, border_flag_value, function, data, buffer);
+ CASE_FILTER_POINT(pi, oo, filter_size, cvalue, Float64,
+ tmp, border_flag_value, function, data, buffer);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FILTER_OUT(po, tmp, Bool);
+ CASE_FILTER_OUT(po, tmp, UInt8);
+ CASE_FILTER_OUT(po, tmp, UInt16);
+ CASE_FILTER_OUT(po, tmp, UInt32);
+#if HAS_UINT64
+ CASE_FILTER_OUT(po, tmp, UInt64);
+#endif
+ CASE_FILTER_OUT(po, tmp, Int8);
+ CASE_FILTER_OUT(po, tmp, Int16);
+ CASE_FILTER_OUT(po, tmp, Int32);
+ CASE_FILTER_OUT(po, tmp, Int64);
+ CASE_FILTER_OUT(po, tmp, Float32);
+ CASE_FILTER_OUT(po, tmp, Float64);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "array type not supported");
+ goto exit;
+ }
+ NI_FILTER_NEXT2(fi, ii, io, oo, pi, po);
+ }
+exit:
+ if (offsets) free(offsets);
+ if (buffer) free(buffer);
+ return PyErr_Occurred() ? 0 : 1;
+}
Added: branches/Interpolate1D/ndimage/ni_filters.h
===================================================================
--- branches/Interpolate1D/ndimage/ni_filters.h 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/ni_filters.h 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,54 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+ * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef NI_FILTERS_H
+#define NI_FILTERS_H
+
+int NI_Correlate1D(PyArrayObject*, PyArrayObject*, int, PyArrayObject*,
+ NI_ExtendMode, double, maybelong);
+int NI_Correlate(PyArrayObject*, PyArrayObject*, PyArrayObject*,
+ NI_ExtendMode, double, maybelong*);
+int NI_UniformFilter1D(PyArrayObject*, long, int, PyArrayObject*,
+ NI_ExtendMode, double, long);
+int NI_MinOrMaxFilter1D(PyArrayObject*, long, int, PyArrayObject*,
+ NI_ExtendMode, double, long, int);
+int NI_MinOrMaxFilter(PyArrayObject*, PyArrayObject*, PyArrayObject*,
+ PyArrayObject*, NI_ExtendMode, double, maybelong*,
+ int);
+int NI_RankFilter(PyArrayObject*, int, PyArrayObject*, PyArrayObject*,
+ NI_ExtendMode, double, maybelong*);
+int NI_GenericFilter1D(PyArrayObject*, int (*)(double*, maybelong,
+ double*, maybelong, void*), void*, long, int,
+ PyArrayObject*, NI_ExtendMode, double, long);
+int NI_GenericFilter(PyArrayObject*, int (*)(double*, maybelong, double*,
+ void*), void*, PyArrayObject*, PyArrayObject*,
+ NI_ExtendMode, double, maybelong*);
+#endif
Added: branches/Interpolate1D/ndimage/ni_fourier.c
===================================================================
--- branches/Interpolate1D/ndimage/ni_fourier.c 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/ni_fourier.c 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,548 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+ * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "ni_support.h"
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>
+
+#if !defined(M_PI)
+#define M_PI 3.14159265358979323846
+#endif
+
+#define _NI_GAUSSIAN 0
+#define _NI_UNIFORM 1
+#define _NI_ELLIPSOID 2
+
+static double polevl(double x, const double coef[], int N)
+{
+ double ans;
+ const double *p = coef;
+ int i = N;
+
+ ans = *p++;
+ do
+ ans = ans * x + *p++;
+ while(--i);
+
+ return ans ;
+}
+
+double p1evl(double x, const double coef[], int N)
+{
+ double ans;
+ const double *p = coef;
+ int i = N - 1;
+
+ ans = x + *p++;
+ do
+ ans = ans * x + *p++;
+ while(--i);
+
+ return ans;
+}
+
+#define THPIO4 2.35619449019234492885
+#define SQ2OPI .79788456080286535588
+#define Z1 1.46819706421238932572E1
+#define Z2 4.92184563216946036703E1
+
+static double _bessel_j1(double x)
+{
+ double w, z, p, q, xn;
+ const double RP[4] = {
+ -8.99971225705559398224E8,
+ 4.52228297998194034323E11,
+ -7.27494245221818276015E13,
+ 3.68295732863852883286E15,
+ };
+ const double RQ[8] = {
+ 6.20836478118054335476E2,
+ 2.56987256757748830383E5,
+ 8.35146791431949253037E7,
+ 2.21511595479792499675E10,
+ 4.74914122079991414898E12,
+ 7.84369607876235854894E14,
+ 8.95222336184627338078E16,
+ 5.32278620332680085395E18,
+ };
+ const double PP[7] = {
+ 7.62125616208173112003E-4,
+ 7.31397056940917570436E-2,
+ 1.12719608129684925192E0,
+ 5.11207951146807644818E0,
+ 8.42404590141772420927E0,
+ 5.21451598682361504063E0,
+ 1.00000000000000000254E0,
+ };
+ const double PQ[7] = {
+ 5.71323128072548699714E-4,
+ 6.88455908754495404082E-2,
+ 1.10514232634061696926E0,
+ 5.07386386128601488557E0,
+ 8.39985554327604159757E0,
+ 5.20982848682361821619E0,
+ 9.99999999999999997461E-1,
+ };
+ const double QP[8] = {
+ 5.10862594750176621635E-2,
+ 4.98213872951233449420E0,
+ 7.58238284132545283818E1,
+ 3.66779609360150777800E2,
+ 7.10856304998926107277E2,
+ 5.97489612400613639965E2,
+ 2.11688757100572135698E2,
+ 2.52070205858023719784E1,
+ };
+ const double QQ[7] = {
+ 7.42373277035675149943E1,
+ 1.05644886038262816351E3,
+ 4.98641058337653607651E3,
+ 9.56231892404756170795E3,
+ 7.99704160447350683650E3,
+ 2.82619278517639096600E3,
+ 3.36093607810698293419E2,
+ };
+
+ w = x;
+ if (x < 0)
+ w = -x;
+
+ if (w <= 5.0) {
+ z = x * x;
+ w = polevl(z, RP, 3) / p1evl(z, RQ, 8);
+ w = w * x * (z - Z1) * (z - Z2);
+ return w ;
+ }
+
+ w = 5.0 / x;
+ z = w * w;
+ p = polevl(z, PP, 6) / polevl(z, PQ, 6);
+ q = polevl(z, QP, 7) / p1evl(z, QQ, 7);
+ xn = x - THPIO4;
+ p = p * cos(xn) - w * q * sin(xn);
+ return p * SQ2OPI / sqrt(x);
+}
+
+#define CASE_FOURIER_OUT_RR(_po, _tmp, _type) \
+case t ## _type: \
+ *(_type*)_po = _tmp; \
+ break
+
+#define CASE_FOURIER_OUT_RC(_po, _tmp, _type) \
+case t ## _type: \
+ (*(_type*)_po).real = tmp; \
+ (*(_type*)_po).imag = 0.0; \
+ break
+
+#define CASE_FOURIER_OUT_CC(_po, _tmp_r, _tmp_i, _type) \
+case t ## _type: \
+ (*(_type*)_po).real = _tmp_r; \
+ (*(_type*)_po).imag = _tmp_i; \
+ break
+
+#define CASE_FOURIER_FILTER_RC(_pi, _tmp, _tmp_r, _tmp_i, _type) \
+case t ## _type: \
+ _tmp_r = (*(_type*)_pi).real * _tmp; \
+ _tmp_i = (*(_type*)_pi).imag * _tmp; \
+ break;
+
+#define CASE_FOURIER_FILTER_RR(_pi, _tmp, _type) \
+case t ## _type: \
+ _tmp *= *(_type*)_pi; \
+ break;
+
+int NI_FourierFilter(PyArrayObject *input, PyArrayObject* parameter_array,
+ maybelong n, int axis, PyArrayObject* output, int filter_type)
+{
+ NI_Iterator ii, io;
+ char *pi, *po;
+ double *parameters = NULL, **params = NULL;
+ maybelong kk, hh, size;
+ Float64 *iparameters = (void *)PyArray_DATA(parameter_array);
+ int ll;
+
+ /* precalculate the parameters: */
+ parameters = (double*)malloc(input->nd * sizeof(double));
+ if (!parameters) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(kk = 0; kk < input->nd; kk++) {
+ /* along the direction of the real transform we must use the given
+ length of that dimensons, unless a complex transform is assumed
+ (n < 0): */
+ int shape = kk == axis ?
+ (n < 0 ? input->dimensions[kk] : n) : input->dimensions[kk];
+ switch (filter_type) {
+ case _NI_GAUSSIAN:
+ parameters[kk] = *iparameters++ * M_PI / (double)shape;
+ parameters[kk] = -2.0 * parameters[kk] * parameters[kk];
+ break;
+ case _NI_ELLIPSOID:
+ case _NI_UNIFORM:
+ parameters[kk] = *iparameters++;
+ break;
+ }
+ }
+ /* allocate memory for tables: */
+ params = (double**) malloc(input->nd * sizeof(double*));
+ if (!params) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(kk = 0; kk < input->nd; kk++)
+ params[kk] = NULL;
+ for(kk = 0; kk < input->nd; kk++) {
+ if (input->dimensions[kk] > 1) {
+ params[kk] = (double*)malloc(input->dimensions[kk] * sizeof(double));
+ if (!params[kk]) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+ }
+ switch (filter_type) {
+ case _NI_GAUSSIAN:
+ /* calculate the tables of exponentials: */
+ for (hh = 0; hh < input->nd; hh++) {
+ if (params[hh]) {
+ if (hh == axis && n >= 0) {
+ for(kk = 0; kk < input->dimensions[hh]; kk++) {
+ double tmp = parameters[hh] * kk * kk;
+ params[hh][kk] = fabs(tmp) > 50.0 ? 0.0 : exp(tmp);
+ }
+ } else {
+ int jj = 0;
+ for(kk = 0; kk < (input->dimensions[hh] + 1) / 2; kk++) {
+ double tmp = parameters[hh] * kk * kk;
+ params[hh][jj++] = fabs(tmp) > 50.0 ? 0.0 : exp(tmp);
+ }
+ for(kk = -(input->dimensions[hh] / 2); kk < 0; kk++) {
+ double tmp = parameters[hh] * kk * kk;
+ params[hh][jj++] = fabs(tmp) > 50.0 ? 0.0 : exp(tmp);
+ }
+ }
+ }
+ }
+ break;
+ case _NI_UNIFORM:
+ /* calculate the tables of parameters: */
+ for (hh = 0; hh < input->nd; hh++) {
+ if (params[hh]) {
+ params[hh][0] = 1.0;
+ if (hh == axis && n >= 0) {
+ double tmp = M_PI * parameters[hh] / n;
+ for(kk = 1; kk < input->dimensions[hh]; kk++)
+ params[hh][kk] = tmp > 0.0 ?
+ sin(tmp * kk) / (tmp * kk) : 0.0;
+ } else {
+ double tmp = M_PI * parameters[hh] / input->dimensions[hh];
+ int jj = 1;
+ for(kk = 1; kk < (input->dimensions[hh] + 1) / 2; kk++)
+ params[hh][jj++] = tmp > 0.0 ?
+ sin(tmp * kk) / (tmp * kk) : 0.0;
+ for(kk = -(input->dimensions[hh] / 2); kk < 0; kk++)
+ params[hh][jj++] = tmp > 0.0 ?
+ sin(tmp * kk) / (tmp * kk) : 0.0;
+ }
+ }
+ }
+ break;
+ case _NI_ELLIPSOID:
+ /* calculate the tables of parameters: */
+ for (hh = 0; hh < input->nd; hh++) {
+ if (params[hh]) {
+ params[hh][0] = 1.0;
+ if (hh == axis && n >= 0) {
+ double tmp = M_PI * parameters[hh] / n;
+ for(kk = 0; kk < input->dimensions[hh]; kk++)
+ params[hh][kk] = (double)kk * tmp;
+ } else {
+ double tmp = M_PI * parameters[hh] / input->dimensions[hh];
+ int jj = 0;
+ for(kk = 0; kk < (input->dimensions[hh] + 1) / 2; kk++)
+ params[hh][jj++] = (double)kk * tmp;
+ for(kk = -(input->dimensions[hh] / 2); kk < 0; kk++)
+ params[hh][jj++] = (double)kk * tmp;
+ }
+ } else if (input->dimensions[hh] > 0) {
+ params[hh][0] = 1.0;
+ }
+ }
+ if (input->nd > 1)
+ for(hh = 0; hh < input->nd; hh++)
+ for(kk = 0; kk < input->dimensions[hh]; kk++)
+ params[hh][kk] = params[hh][kk] * params[hh][kk];
+ break;
+ default:
+ break;
+ }
+ /* initialize input element iterator: */
+ if (!NI_InitPointIterator(input, &ii))
+ goto exit;
+ /* initialize output element iterator: */
+ if (!NI_InitPointIterator(output, &io))
+ goto exit;
+ pi = (void *)PyArray_DATA(input);
+ po = (void *)PyArray_DATA(output);
+ size = 1;
+ for(ll = 0; ll < input->nd; ll++)
+ size *= input->dimensions[ll];
+ /* iterator over the elements: */
+ for(hh = 0; hh < size; hh++) {
+ double tmp = 1.0;
+ switch (filter_type) {
+ case _NI_GAUSSIAN:
+ case _NI_UNIFORM:
+ for(kk = 0; kk < input->nd; kk++)
+ if (params[kk])
+ tmp *= params[kk][ii.coordinates[kk]];
+ break;
+ case _NI_ELLIPSOID:
+ switch (input->nd) {
+ case 1:
+ tmp = params[0][ii.coordinates[0]];
+ tmp = tmp > 0.0 ? sin(tmp) / (tmp) : 1.0;
+ break;
+ case 2:
+ tmp = 0.0;
+ for(kk = 0; kk < 2; kk++)
+ tmp += params[kk][ii.coordinates[kk]];
+ tmp = sqrt(tmp);
+ tmp = tmp > 0.0 ? 2.0 * _bessel_j1(tmp) / tmp : 1.0;
+ break;
+ case 3:
+ {
+ double r = 0.0;
+ for(kk = 0; kk < 3; kk++)
+ r += params[kk][ii.coordinates[kk]];
+ r = sqrt(r);
+ if (r > 0.0) {
+ tmp = 3.0 * (sin(r) - r * cos(r));
+ tmp /= r * r * r;
+ } else {
+ tmp = 1.0;
+ }
+ }
+ break;
+ }
+ break;
+ default:
+ break;
+ }
+ if (input->descr->type_num == tComplex64 ||
+ input->descr->type_num == tComplex128) {
+ double tmp_r = 0.0, tmp_i = 0.0;
+ switch (input->descr->type_num) {
+ CASE_FOURIER_FILTER_RC(pi, tmp, tmp_r, tmp_i, Complex64);
+ CASE_FOURIER_FILTER_RC(pi, tmp, tmp_r, tmp_i, Complex128);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "data type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FOURIER_OUT_CC(po, tmp_r, tmp_i, Complex64);
+ CASE_FOURIER_OUT_CC(po, tmp_r, tmp_i, Complex128);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "data type not supported");
+ goto exit;
+ }
+ } else {
+ switch (input->descr->type_num) {
+ CASE_FOURIER_FILTER_RR(pi, tmp, Bool)
+ CASE_FOURIER_FILTER_RR(pi, tmp, UInt8)
+ CASE_FOURIER_FILTER_RR(pi, tmp, UInt16)
+ CASE_FOURIER_FILTER_RR(pi, tmp, UInt32)
+#if HAS_UINT64
+ CASE_FOURIER_FILTER_RR(pi, tmp, UInt64)
+#endif
+ CASE_FOURIER_FILTER_RR(pi, tmp, Int8)
+ CASE_FOURIER_FILTER_RR(pi, tmp, Int16)
+ CASE_FOURIER_FILTER_RR(pi, tmp, Int32)
+ CASE_FOURIER_FILTER_RR(pi, tmp, Int64)
+ CASE_FOURIER_FILTER_RR(pi, tmp, Float32)
+ CASE_FOURIER_FILTER_RR(pi, tmp, Float64)
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "data type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FOURIER_OUT_RR(po, tmp, Float32);
+ CASE_FOURIER_OUT_RR(po, tmp, Float64);
+ CASE_FOURIER_OUT_RC(po, tmp, Complex64);
+ CASE_FOURIER_OUT_RC(po, tmp, Complex128);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "data type not supported");
+ goto exit;
+ }
+ }
+ NI_ITERATOR_NEXT2(ii, io, pi, po);
+ }
+
+ exit:
+ if (parameters) free(parameters);
+ if (params) {
+ for(kk = 0; kk < input->nd; kk++)
+ if (params[kk]) free(params[kk]);
+ free(params);
+ }
+ return PyErr_Occurred() ? 0 : 1;
+}
+
+#define CASE_FOURIER_SHIFT_R(_pi, _tmp, _r, _i, _cost, _sint, _type) \
+case t ## _type: \
+ _tmp = *(_type*)_pi; \
+ _r = _tmp * _cost; \
+ _i = _tmp * _sint; \
+ break;
+
+#define CASE_FOURIER_SHIFT_C(_pi, _r, _i, _cost, _sint, _type) \
+case t ## _type: \
+ _r = (*(_type*)_pi).real * _cost - (*(_type*)_pi).imag * _sint; \
+ _i = (*(_type*)_pi).real * _sint + (*(_type*)_pi).imag * _cost; \
+ break;
+
+int NI_FourierShift(PyArrayObject *input, PyArrayObject* shift_array,
+ maybelong n, int axis, PyArrayObject* output)
+{
+ NI_Iterator ii, io;
+ char *pi, *po;
+ double *shifts = NULL, **params = NULL;
+ maybelong kk, hh, size;
+ Float64 *ishifts = (void *)PyArray_DATA(shift_array);
+ int ll;
+
+ /* precalculate the shifts: */
+ shifts = (double*)malloc(input->nd * sizeof(double));
+ if (!shifts) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(kk = 0; kk < input->nd; kk++) {
+ /* along the direction of the real transform we must use the given
+ length of that dimensons, unless a complex transform is assumed
+ (n < 0): */
+ int shape = kk == axis ?
+ (n < 0 ? input->dimensions[kk] : n) : input->dimensions[kk];
+ shifts[kk] = -2.0 * M_PI * *ishifts++ / (double)shape;
+ }
+ /* allocate memory for tables: */
+ params = (double**) malloc(input->nd * sizeof(double*));
+ if (!params) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ for(kk = 0; kk < input->nd; kk++)
+ params[kk] = NULL;
+ for(kk = 0; kk < input->nd; kk++) {
+ if (input->dimensions[kk] > 1) {
+ params[kk] = (double*)malloc(input->dimensions[kk] * sizeof(double));
+ if (!params[kk]) {
+ PyErr_NoMemory();
+ goto exit;
+ }
+ }
+ }
+ for (hh = 0; hh < input->nd; hh++) {
+ if (params[hh]) {
+ if (hh == axis && n >= 0) {
+ for(kk = 0; kk < input->dimensions[hh]; kk++)
+ params[hh][kk] = shifts[hh] * kk;
+ } else {
+ int jj = 0;
+ for(kk = 0; kk < (input->dimensions[hh] + 1) / 2; kk++) {
+ params[hh][jj++] = shifts[hh] * kk;
+ }
+ for(kk = -(input->dimensions[hh] / 2); kk < 0; kk++) {
+ params[hh][jj++] = shifts[hh] * kk;
+ }
+ }
+ }
+ }
+ /* initialize input element iterator: */
+ if (!NI_InitPointIterator(input, &ii))
+ goto exit;
+ /* initialize output element iterator: */
+ if (!NI_InitPointIterator(output, &io))
+ goto exit;
+ pi = (void *)PyArray_DATA(input);
+ po = (void *)PyArray_DATA(output);
+ size = 1;
+ for(ll = 0; ll < input->nd; ll++)
+ size *= input->dimensions[ll];
+ /* iterator over the elements: */
+ for(hh = 0; hh < size; hh++) {
+ double tmp = 0.0, sint, cost, r = 0.0, i = 0.0;
+ for(kk = 0; kk < input->nd; kk++)
+ if (params[kk])
+ tmp += params[kk][ii.coordinates[kk]];
+ sint = sin(tmp);
+ cost = cos(tmp);
+ switch (input->descr->type_num) {
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Bool)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, UInt8)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, UInt16)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, UInt32)
+#if HAS_UINT64
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, UInt64)
+#endif
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Int8)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Int16)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Int32)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Int64)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Float32)
+ CASE_FOURIER_SHIFT_R(pi, tmp, r, i, cost, sint, Float64)
+ CASE_FOURIER_SHIFT_C(pi, r, i, cost, sint, Complex64)
+ CASE_FOURIER_SHIFT_C(pi, r, i, cost, sint, Complex128)
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "data type not supported");
+ goto exit;
+ }
+ switch (output->descr->type_num) {
+ CASE_FOURIER_OUT_CC(po, r, i, Complex64);
+ CASE_FOURIER_OUT_CC(po, r, i, Complex128);
+ default:
+ PyErr_SetString(PyExc_RuntimeError, "data type not supported");
+ goto exit;
+ }
+ NI_ITERATOR_NEXT2(ii, io, pi, po);
+ }
+
+ exit:
+ if (shifts) free(shifts);
+ if (params) {
+ for(kk = 0; kk < input->nd; kk++)
+ if (params[kk]) free(params[kk]);
+ free(params);
+ }
+ return PyErr_Occurred() ? 0 : 1;
+}
Added: branches/Interpolate1D/ndimage/ni_fourier.h
===================================================================
--- branches/Interpolate1D/ndimage/ni_fourier.h 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/ni_fourier.h 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,40 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+ * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef NI_FOURIER_H
+#define NI_FOURIER_H
+
+int NI_FourierFilter(PyArrayObject*, PyArrayObject*, maybelong, int,
+ PyArrayObject*, int);
+int NI_FourierShift(PyArrayObject*, PyArrayObject*, maybelong, int,
+ PyArrayObject*);
+
+#endif
Added: branches/Interpolate1D/ndimage/ni_interpolation.c
===================================================================
--- branches/Interpolate1D/ndimage/ni_interpolation.c 2008-08-04 20:27:24 UTC (rev 4600)
+++ branches/Interpolate1D/ndimage/ni_interpolation.c 2008-08-05 19:07:52 UTC (rev 4601)
@@ -0,0 +1,966 @@
+/* Copyright (C) 2003-2005 Peter J. Verveer
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ *
+ * 3. The name of the author may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``A