← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/stsci/image/src/_combinemodule.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/stsci/image/src/_combinemodule.c
 
1
#include "Python.h"
 
2
 
 
3
#include <stdio.h>
 
4
#include <math.h>
 
5
#include <signal.h>
 
6
#include <ctype.h>
 
7
 
 
8
#include "numpy/libnumarray.h"
 
9
 
 
10
#define MAX_ARRAYS 1024
 
11
 
 
12
static PyObject *_Error;
 
13
 
 
14
typedef Float64 (*combiner)(int, int, int, Float64 temp[MAX_ARRAYS]);
 
15
 
 
16
 
 
17
static int
 
18
_mask_and_sort(int ninputs, int index, Float64 **inputs, Bool **masks,
 
19
               Float64 temp[MAX_ARRAYS])
 
20
{
 
21
        int i, j, goodpix;
 
22
        if (masks) {
 
23
                for (i=j=0; i<ninputs; i++) {
 
24
                        if (masks[i][index] == 0) 
 
25
                                temp[j++] = inputs[i][index];
 
26
                }
 
27
        } else {
 
28
                for (i=j=0; i<ninputs; i++) {
 
29
                        temp[j++] = inputs[i][index];
 
30
                }
 
31
        }
 
32
        goodpix = j;
 
33
        for(i=0; i<goodpix; i++) {
 
34
                for (j=i+1; j<goodpix; j++) {
 
35
                        if (temp[j] < temp[i]) {
 
36
                                Float64 temp2 = temp[j];
 
37
                                temp[j] = temp[i];
 
38
                                temp[i] = temp2;
 
39
                        }
 
40
                }
 
41
        }
 
42
        return goodpix;
 
43
}
 
44
 
 
45
static Float64
 
46
_inner_median(int goodpix, int nlow, int nhigh, Float64 temp[MAX_ARRAYS])
 
47
{
 
48
        Float64 median;
 
49
        int midpoint, medianpix = goodpix-nhigh-nlow;
 
50
        if (medianpix <= 0) {
 
51
                median = 0;
 
52
        } else {
 
53
                midpoint = medianpix / 2;
 
54
                if (medianpix % 2) /* odd */ {
 
55
                        median = temp[ midpoint + nlow ];
 
56
                } else {
 
57
                        median = (temp[ midpoint + nlow ] + 
 
58
                                  temp[ midpoint + nlow - 1 ]) / 2.0;
 
59
                }       
 
60
        }
 
61
        return median;
 
62
}
 
63
 
 
64
static Float64
 
65
_inner_average(int goodpix, int nlow, int nhigh, Float64 temp[MAX_ARRAYS])
 
66
{
 
67
        Float64 average;
 
68
        int i, averagepix = goodpix-nhigh-nlow;
 
69
 
 
70
        if (averagepix <= 0) {
 
71
                average = 0;
 
72
        } else {
 
73
                for(i=nlow, average=0; i<averagepix+nlow;  i++)
 
74
                        average += temp[i];
 
75
                average /= averagepix;
 
76
        }
 
77
        return average;
 
78
}
 
79
 
 
80
static Float64
 
81
_inner_minimum(int goodpix, int nlow, int nhigh, Float64 temp[MAX_ARRAYS])
 
82
{
 
83
        int minimumpix = goodpix-nhigh-nlow;
 
84
        if (minimumpix <= 0) {
 
85
                return 0;
 
86
        } else {
 
87
               return temp[nlow];
 
88
        }
 
89
}
 
90
 
 
91
static int
 
92
_combine(combiner f, int dim, int maxdim, int ninputs, int nlow, int nhigh,
 
93
        PyArrayObject *inputs[], PyArrayObject *masks[], PyArrayObject *output)
 
94
{
 
95
        int i, j;
 
96
 
 
97
        if (dim == maxdim-1) {
 
98
                Float64 sorted[MAX_ARRAYS];
 
99
                Float64 *tinputs[MAX_ARRAYS];
 
100
                Bool    *tmasks[MAX_ARRAYS];
 
101
                Float64 *toutput;
 
102
                int cols = inputs[0]->dimensions[dim];
 
103
 
 
104
                /* Allocate and convert 1 temporary row at a time */
 
105
                for(i=0; i<ninputs; i++)
 
106
                        tinputs[i] = (Float64 *) inputs[i]->data;
 
107
                if (masks) {
 
108
                        for(i=0; i<ninputs; i++)
 
109
                                tmasks[i] = (Bool *) masks[i]->data;
 
110
                }
 
111
                toutput = (Float64 *) output->data;
 
112
                
 
113
                for(j=0; j<cols; j++) {
 
114
                        int goodpix = _mask_and_sort(
 
115
                                ninputs, j, tinputs, masks ? tmasks : NULL, sorted);
 
116
                        toutput[j] = f(goodpix, nlow, nhigh, sorted);
 
117
                }
 
118
        } else {
 
119
                for (i=0; i<inputs[0]->dimensions[dim]; i++) {
 
120
                        for(j=0; j<ninputs; j++) {
 
121
                                inputs[j]->data += inputs[j]->strides[dim]*i;
 
122
                                if (masks) {
 
123
                                        masks[j]->data += masks[j]->strides[dim]*i;
 
124
                                }
 
125
                        }
 
126
                        output->data += output->strides[dim]*i;
 
127
                        _combine(f, dim+1, maxdim, ninputs, nlow, nhigh, 
 
128
                                inputs, masks, output);
 
129
                        for(j=0; j<ninputs; j++) {
 
130
                                inputs[j]->data -= inputs[j]->strides[dim]*i;
 
131
                                if (masks) {
 
132
                                        masks[j]->data -= masks[j]->strides[dim]*i;
 
133
                                }
 
134
                        }
 
135
                        output->data -= output->strides[dim]*i;
 
136
                }
 
137
        }
 
138
        return 0;
 
139
}
 
140
 
 
141
typedef struct
 
142
{
 
143
        char *name;
 
144
        combiner fptr;
 
145
} fmapping;
 
146
 
 
147
static fmapping functions[] = {
 
148
        {"median", _inner_median},
 
149
        {"average", _inner_average},
 
150
        {"minimum", _inner_minimum},
 
151
};
 
152
 
 
153
 
 
154
static PyObject *
 
155
_Py_combine(PyObject *obj, PyObject *args, PyObject *kw)
 
156
{
 
157
        PyObject   *arrays, *output;
 
158
        int        nlow=0, nhigh=0, narrays;
 
159
        PyObject   *badmasks=Py_None;
 
160
        char       *keywds[] = { "arrays", "output", "nlow", "nhigh", 
 
161
                                 "badmasks", "kind", NULL };
 
162
        char *kind;
 
163
        combiner f;
 
164
        PyArrayObject  *arr[MAX_ARRAYS], *bmk[MAX_ARRAYS], *toutput;
 
165
        int i;
 
166
 
 
167
        if (!PyArg_ParseTupleAndKeywords(args, kw, "OO|iiOs:combine", keywds, 
 
168
                         &arrays, &output, &nlow, &nhigh, &badmasks, &kind))
 
169
                return NULL;
 
170
 
 
171
        narrays = PySequence_Length(arrays);
 
172
        if (narrays < 0) 
 
173
                return PyErr_Format(
 
174
                        PyExc_TypeError, "combine: arrays is not a sequence");
 
175
        if (narrays > MAX_ARRAYS)
 
176
                return PyErr_Format(
 
177
                        PyExc_TypeError, "combine: too many arrays.");
 
178
 
 
179
        for(i=0; i<narrays; i++) {
 
180
                PyObject *a = PySequence_GetItem(arrays, i);
 
181
                if (!a) return NULL;
 
182
                arr[i] = NA_InputArray(a, tFloat64, C_ARRAY);
 
183
                if (!arr[i]) return NULL;
 
184
                Py_DECREF(a);
 
185
                if (badmasks != Py_None) {
 
186
                        a =  PySequence_GetItem(badmasks, i);
 
187
                        if (!a) return NULL;    
 
188
                        bmk[i] = NA_InputArray(a, tBool, C_ARRAY);
 
189
                        if (!bmk[i]) return NULL;
 
190
                        Py_DECREF(a);
 
191
                }
 
192
        }
 
193
 
 
194
        toutput = NA_OutputArray(output, tFloat64, C_ARRAY);
 
195
        if (!toutput) return NULL;
 
196
        
 
197
        for (i=0,f=0; i<ELEM(functions); i++)
 
198
                if  (!strcmp(kind, functions[i].name)) {
 
199
                        f = functions[i].fptr;
 
200
                        break;
 
201
                }
 
202
        if (!f) return PyErr_Format(
 
203
                PyExc_ValueError, "Invalid comination function.");
 
204
 
 
205
        if (_combine( f, 0, arr[0]->nd, narrays, nlow, nhigh, 
 
206
                     arr, (badmasks != Py_None ? bmk : NULL), 
 
207
                     toutput) < 0)
 
208
                return NULL;
 
209
 
 
210
        for(i=0; i<narrays; i++) {
 
211
                Py_DECREF(arr[i]);
 
212
                if (badmasks != Py_None) {
 
213
                        Py_DECREF(bmk[i]);
 
214
                }
 
215
        }
 
216
        Py_DECREF(toutput);
 
217
 
 
218
        Py_INCREF(Py_None);
 
219
        return Py_None;
 
220
}
 
221
 
 
222
static PyMethodDef _combineMethods[] = {
 
223
    {"combine", (PyCFunction) _Py_combine, METH_VARARGS | METH_KEYWORDS}, 
 
224
    {NULL, NULL} /* Sentinel */
 
225
};
 
226
 
 
227
/* platform independent*/
 
228
#ifdef MS_WIN32
 
229
__declspec(dllexport)
 
230
#endif
 
231
 
 
232
void init_combine(void)
 
233
{
 
234
        PyObject *m, *d;
 
235
        m = Py_InitModule("_combine", _combineMethods);
 
236
        d = PyModule_GetDict(m);
 
237
        _Error = PyErr_NewException("_combine.error", NULL, NULL);
 
238
        PyDict_SetItemString(d, "error", _Error);
 
239
        import_libnumarray();
 
240
}
 
241
 
 
242
/*
 
243
 * Local Variables:
 
244
 * mode: C
 
245
 * c-file-style: "python"
 
246
 * End:
 
247
 */