← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

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
 
 */