← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/stsci/convolve/src/_correlatemodule.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:
scipy/stsci/convolve/src/_correlatemodule.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
typedef enum
 
11
{
 
12
        PIX_NEAREST,
 
13
        PIX_REFLECT,
 
14
        PIX_WRAP,
 
15
        PIX_CONSTANT
 
16
} PixMode;
 
17
 
 
18
typedef struct
 
19
{
 
20
        PixMode mode;
 
21
        long    rows, cols;
 
22
        Float64 constval;
 
23
        Float64 *data;
 
24
} PixData;
 
25
 
 
26
static long 
 
27
SlowCoord(long x, long maxx, PixMode m)
 
28
{
 
29
        switch(m) {
 
30
        case PIX_NEAREST:
 
31
                if (x < 0) x = 0;
 
32
                if (x >= maxx) x = maxx-1;
 
33
                return x;
 
34
        case PIX_REFLECT:
 
35
                if (x < 0) x = -x-1;
 
36
                if (x >= maxx) x = maxx - (x - maxx) - 1;
 
37
                return x;
 
38
        case PIX_WRAP:
 
39
                if (x < 0) x += maxx;
 
40
                if (x >= maxx) x -= maxx;
 
41
                return x;
 
42
        case PIX_CONSTANT:  /* handled in SlowPix, suppress warnings */
 
43
                break;
 
44
        }
 
45
        return x;
 
46
}
 
47
 
 
48
static Float64
 
49
SlowPix(long r, long c, PixData *p)
 
50
{
 
51
        long fr, fc;
 
52
        if (p->mode == PIX_CONSTANT) {
 
53
                if ((r <  0) || (r >= p->rows) || (c < 0) || (c >= p->cols))
 
54
                        return p->constval;
 
55
                else {
 
56
                        fr = r;
 
57
                        fc = c;
 
58
                }
 
59
        } else {
 
60
                fr = SlowCoord(r, p->rows, p->mode);
 
61
                fc = SlowCoord(c, p->cols, p->mode);
 
62
        }
 
63
        return p->data[fr*p->cols + fc];
 
64
}
 
65
 
 
66
static int
 
67
_reject_complex(PyObject *a)
 
68
{
 
69
        NumarrayType t;
 
70
        if ((a == Py_None) || (a == NULL)) 
 
71
                return 0;
 
72
        t = NA_NumarrayType(a);
 
73
        if (t < 0) {
 
74
                PyErr_Clear();
 
75
                return 0;
 
76
        }
 
77
        if (t == tComplex32 || t == tComplex64) {
 
78
                PyErr_Format(PyExc_TypeError,
 
79
                             "function doesn't support complex arrays.");
 
80
                return 1;
 
81
        }
 
82
        return 0;
 
83
}
 
84
 
 
85
static void 
 
86
Correlate1d(long ksizex, Float64 *kernel, 
 
87
           long dsizex, Float64 *data, 
 
88
           Float64 *correlated)
 
89
{
 
90
        long xc;
 
91
        long halfk = ksizex/2;
 
92
 
 
93
        for(xc=0; xc<halfk; xc++)
 
94
                correlated[xc] = data[xc];
 
95
 
 
96
        for(xc=halfk; xc<dsizex-halfk; xc++) {
 
97
                long xk;
 
98
                double temp = 0;
 
99
                for (xk=0; xk<ksizex; xk++)
 
100
                        temp += kernel[xk]*data[xc-halfk+xk];
 
101
                correlated[xc] = temp;
 
102
        }
 
103
                     
 
104
        for(xc=dsizex-halfk; xc<dsizex; xc++)
 
105
                correlated[xc] = data[xc];
 
106
}
 
107
 
 
108
static PyObject *
 
109
Py_Correlate1d(PyObject *obj, PyObject *args)
 
110
{
 
111
        PyObject   *okernel, *odata, *ocorrelated=NULL;
 
112
        PyArrayObject *kernel, *data, *correlated;
 
113
 
 
114
        if (!PyArg_ParseTuple(args, "OO|O:Correlate1d", 
 
115
                              &okernel, &odata, &ocorrelated))
 
116
                return NULL;
 
117
 
 
118
        /* Align, Byteswap, Contiguous, Typeconvert */
 
119
        kernel  = NA_InputArray(okernel, tFloat64, C_ARRAY);
 
120
        data    = NA_InputArray(odata, tFloat64, C_ARRAY);
 
121
        correlated = NA_OptionalOutputArray(ocorrelated, tFloat64, C_ARRAY, 
 
122
                                            data);
 
123
 
 
124
        if (!kernel || !data || !correlated)
 
125
                goto _fail;
 
126
 
 
127
        if (_reject_complex(okernel) || _reject_complex(odata) || 
 
128
            _reject_complex(ocorrelated))
 
129
                goto _fail;
 
130
 
 
131
        if ((kernel->nd != 1) || (data->nd != 1)) {
 
132
                PyErr_Format(PyExc_ValueError,
 
133
                             "Correlate1d: numarray must have exactly 1 dimension.");
 
134
                goto _fail;
 
135
        }
 
136
 
 
137
        if (!NA_ShapeEqual(data, correlated)) {
 
138
                PyErr_Format(PyExc_ValueError,
 
139
                             "Correlate1d: data and output must have identical length.");
 
140
                goto _fail;
 
141
        }
 
142
 
 
143
        Correlate1d(kernel->dimensions[0], NA_OFFSETDATA(kernel),
 
144
                    data->dimensions[0],   NA_OFFSETDATA(data),
 
145
                    NA_OFFSETDATA(correlated));
 
146
 
 
147
        Py_DECREF(kernel);
 
148
        Py_DECREF(data);
 
149
 
 
150
        /* Align, Byteswap, Contiguous, Typeconvert */
 
151
        return NA_ReturnOutput(ocorrelated, correlated);
 
152
 
 
153
  _fail:
 
154
        Py_XDECREF(kernel);
 
155
        Py_XDECREF(data);
 
156
        Py_XDECREF(correlated);
 
157
        return NULL;
 
158
}
 
159
 
 
160
/* SlowCorrelate computes 2D correlation near the boundaries of an array.
 
161
The output array shares the same dimensions as the input array, the latter
 
162
fully described by PixData.
 
163
 
 
164
The region defined by rmin,rmax,cmin,cmax is assumed to contain only valid
 
165
coordinates.  However, access to the input array is performed using SlowPix
 
166
because pixels reachable via "kernel offsets" may be at invalid coordinates.
 
167
*/
 
168
static void 
 
169
SlowCorrelate2d(long rmin, long rmax, long cmin, long cmax, 
 
170
              long krows, long kcols, Float64 *kernel, 
 
171
              PixData *pix, Float64 *output)
 
172
{
 
173
        long kr, kc, r, c;
 
174
        long halfkrows = krows/2;
 
175
        long halfkcols = kcols/2;
 
176
 
 
177
        for(r=rmin; r<rmax; r++) { 
 
178
                for(c=cmin; c<cmax; c++) {
 
179
                        Float64 temp = 0;
 
180
                        for(kr=0; kr<krows; kr++) {
 
181
                                long pr = r + kr - halfkrows;
 
182
                                for(kc=0; kc<kcols; kc++) {
 
183
                                        long pc = c + kc - halfkcols;
 
184
                                        temp += SlowPix(pr, pc, pix) * 
 
185
                                                kernel[kr*kcols+kc];
 
186
                                }
 
187
                        }
 
188
                        output[r*pix->cols+c] = temp;
 
189
                }
 
190
        }
 
191
}
 
192
 
 
193
static void 
 
194
Correlate2d(long krows, long kcols, Float64 *kernel, 
 
195
            long drows, long dcols, Float64 *data, Float64 *correlated, 
 
196
            PixMode mode, Float64 cval)
 
197
{
 
198
        long ki, kj, di, dj;
 
199
        long halfkrows = krows/2;
 
200
        long halfkcols = kcols/2;
 
201
        
 
202
        PixData pix;
 
203
        pix.mode = mode;
 
204
        pix.data = data;
 
205
        pix.constval = cval;
 
206
        pix.rows = drows;
 
207
        pix.cols = dcols;
 
208
 
 
209
        /* Compute the boundaries using SlowPix */
 
210
 
 
211
        SlowCorrelate2d(0, halfkrows, 0, dcols, 
 
212
                      krows, kcols, kernel, &pix, correlated); /* top */
 
213
        SlowCorrelate2d(drows-halfkrows, drows, 0, dcols, 
 
214
                      krows, kcols, kernel, &pix, correlated); /* bottom */
 
215
        SlowCorrelate2d(halfkrows, drows-halfkrows, 0, halfkcols, 
 
216
                      krows, kcols, kernel, &pix, correlated); /* left */
 
217
        SlowCorrelate2d(halfkrows, drows-halfkrows, dcols-halfkcols, dcols,
 
218
                      krows, kcols, kernel, &pix, correlated); /* right */
 
219
 
 
220
        /* Correlate the center data using unchecked array access */
 
221
        for(di=halfkrows; di<drows-halfkrows; di++) {
 
222
                for(dj=halfkcols; dj<dcols-halfkcols; dj++) {
 
223
                        Float64 temp = 0;
 
224
                        for(ki=0; ki<krows; ki++) {
 
225
                                long pi = di + ki - halfkrows;
 
226
                                for(kj=0; kj<kcols; kj++) {
 
227
                                        long pj = dj + kj - halfkcols;
 
228
                                        temp += data[pi*dcols+pj] * 
 
229
                                                kernel[ki*kcols+kj];
 
230
                                }
 
231
                        }
 
232
                        correlated[di*dcols+dj] = temp;
 
233
                }
 
234
        }
 
235
}
 
236
 
 
237
static PyObject *
 
238
Py_Correlate2d(PyObject *obj, PyObject *args, PyObject *kw)
 
239
{
 
240
        PyObject      *okernel, *odata, *ocorrelated=NULL;
 
241
        PyArrayObject *kernel,  *data,  *correlated;
 
242
        Float64  cval = 0;
 
243
        int      mode = PIX_NEAREST;
 
244
        char       *keywds[] = { "kernel", "data", "output", "mode", "cval", NULL };
 
245
 
 
246
        if (!PyArg_ParseTupleAndKeywords(
 
247
                    args, kw, "OO|Oid:Correlate2d", keywds, 
 
248
                    &okernel, &odata, &ocorrelated, &mode, &cval))
 
249
                return NULL;
 
250
 
 
251
        if ((mode < PIX_NEAREST) || (mode > PIX_CONSTANT))
 
252
                return PyErr_Format(PyExc_ValueError,
 
253
                            "Correlate2d: mode value not in range(%d,%d)", 
 
254
                                    PIX_NEAREST, PIX_CONSTANT);
 
255
 
 
256
        /* Align, Byteswap, Contiguous, Typeconvert */
 
257
        kernel  = NA_InputArray(okernel, tFloat64, C_ARRAY);
 
258
        data    = NA_InputArray(odata, tFloat64, C_ARRAY);
 
259
        correlated = NA_OptionalOutputArray(ocorrelated, tFloat64, C_ARRAY, 
 
260
                                            data);
 
261
 
 
262
        if (!kernel || !data || !correlated)
 
263
                goto _fail;
 
264
 
 
265
        if ((kernel->nd != 2) || (data->nd != 2) || (correlated->nd != 2)) {
 
266
                PyErr_Format(PyExc_ValueError, "Correlate2d: inputs must have 2 dimensions.");
 
267
                goto _fail;
 
268
        }
 
269
 
 
270
        if (!NA_ShapeEqual(data, correlated)) {
 
271
                PyErr_Format(PyExc_ValueError,
 
272
                             "Correlate2d: data and output numarray need identical shapes.");
 
273
                goto _fail;
 
274
        }       
 
275
 
 
276
        if (_reject_complex(okernel) || _reject_complex(odata) || 
 
277
            _reject_complex(ocorrelated))
 
278
                goto _fail;
 
279
                
 
280
        Correlate2d(kernel->dimensions[0], kernel->dimensions[1], 
 
281
                    NA_OFFSETDATA(kernel),
 
282
                    data->dimensions[0], data->dimensions[1], 
 
283
                    NA_OFFSETDATA(data),
 
284
                    NA_OFFSETDATA(correlated), 
 
285
                    mode, cval);
 
286
 
 
287
        Py_DECREF(kernel);
 
288
        Py_DECREF(data);
 
289
 
 
290
        /* Align, Byteswap, Contiguous, Typeconvert */
 
291
        return NA_ReturnOutput(ocorrelated, correlated);
 
292
 
 
293
  _fail:
 
294
        Py_XDECREF(kernel);
 
295
        Py_XDECREF(data);
 
296
        Py_XDECREF(correlated);
 
297
        return NULL;
 
298
}
 
299
 
 
300
void Shift2d( long rows, long cols, Float64 *data, long dx, long dy, Float64 *output, int mode, Float64 cval)
 
301
{
 
302
        long r, c;
 
303
        PixData pix;
 
304
        pix.mode = mode;
 
305
        pix.constval = cval;
 
306
        pix.rows = rows;
 
307
        pix.cols = cols;
 
308
        pix.data = data;
 
309
 
 
310
        for(r=0; r<rows; r++)
 
311
                for(c=0; c<cols; c++)
 
312
                        output[ r*cols + c] = SlowPix(r+dy, c+dx, &pix);
 
313
}
 
314
 
 
315
static PyObject *
 
316
Py_Shift2d(PyObject *obj, PyObject *args, PyObject *kw)
 
317
{
 
318
        PyObject      *odata, *ooutput=NULL;
 
319
        PyArrayObject *data,  *output;
 
320
        int           dx, dy;
 
321
        Float64       cval = 0;
 
322
        int           mode = PIX_NEAREST;
 
323
        char          *keywds[] = { "data", "dx", "dy", 
 
324
                                    "output", "mode", "cval", NULL };
 
325
 
 
326
        if (!PyArg_ParseTupleAndKeywords(args, kw, "Oii|Oid:Shift2d", keywds,
 
327
                         &odata, &dx, &dy, &ooutput, &mode, &cval))
 
328
                return NULL;
 
329
 
 
330
        if ((mode < PIX_NEAREST) || (mode > PIX_CONSTANT))
 
331
                return PyErr_Format(PyExc_ValueError,
 
332
                            "Shift2d: mode value not in range(%d,%d)", 
 
333
                                    PIX_NEAREST, PIX_CONSTANT);
 
334
 
 
335
        /* Align, Byteswap, Contiguous, Typeconvert */
 
336
        data    = NA_InputArray(odata, tFloat64, C_ARRAY);
 
337
        output  = NA_OptionalOutputArray(ooutput, tFloat64, C_ARRAY, 
 
338
                                         data);
 
339
 
 
340
        if (!data || !output)
 
341
                goto _fail;
 
342
 
 
343
        if (_reject_complex(odata) || _reject_complex(ooutput))
 
344
                goto _fail;
 
345
 
 
346
        if ((data->nd != 2)) {
 
347
                PyErr_Format(PyExc_ValueError,
 
348
                             "Shift2d: numarray must have 2 dimensions.");
 
349
                goto _fail;
 
350
        }
 
351
 
 
352
        if (!NA_ShapeEqual(data, output)) {
 
353
                PyErr_Format(PyExc_ValueError,
 
354
                             "Shift2d: data and output numarray need identical shapes.");
 
355
                goto _fail;
 
356
        }
 
357
 
 
358
        /* Invert sign of deltas to match sense of 2x2 correlation. */
 
359
        Shift2d( data->dimensions[0], data->dimensions[1], NA_OFFSETDATA(data),
 
360
                 -dx, -dy, NA_OFFSETDATA(output), mode, cval);
 
361
        
 
362
        Py_XDECREF(data);
 
363
 
 
364
        /* Align, Byteswap, Contiguous, Typeconvert */
 
365
        return NA_ReturnOutput(ooutput, output);
 
366
  _fail:
 
367
        Py_XDECREF(data);
 
368
        Py_XDECREF(output);
 
369
        return NULL;
 
370
}
 
371
 
 
372
typedef struct s_BoxData BoxData;
 
373
 
 
374
typedef Float64 (*SumColFunc)(long,long,BoxData*);
 
375
typedef Float64 (*SumBoxFunc)(long,long,BoxData*);
 
376
 
 
377
struct s_BoxData {
 
378
        PixData    pix;
 
379
        long       krows, kcols;
 
380
        SumColFunc sumcol;
 
381
        SumBoxFunc sumbox;
 
382
};
 
383
 
 
384
static Float64
 
385
SlowSumCol(long r, long c, BoxData *D) 
 
386
{
 
387
        Float64 sum = 0;
 
388
        long i, krows = D->krows;
 
389
        for(i=0; i<krows; i++) {
 
390
                sum += SlowPix(r+i, c, &D->pix);
 
391
        }
 
392
        return sum;
 
393
}
 
394
 
 
395
static Float64
 
396
SlowSumBox(long r, long c, BoxData *D)
 
397
{
 
398
        long i, j;
 
399
        Float64 sum = 0;
 
400
        for(i=0; i<D->krows; i++)
 
401
                for(j=0; j<D->kcols; j++)
 
402
                        sum += SlowPix(r+i, c+j, &D->pix);
 
403
        return sum;
 
404
}
 
405
 
 
406
static Float64
 
407
FastSumCol(long r, long c, BoxData *D) 
 
408
{
 
409
        Float64 sum = 0;
 
410
        long krows = D->krows;
 
411
        long cols = D->pix.cols;
 
412
        Float64 *data = D->pix.data;
 
413
 
 
414
        data += r*cols + c;
 
415
        for(; krows--; data += cols) {
 
416
                sum += *data;
 
417
        }
 
418
        return sum;
 
419
}
 
420
 
 
421
static Float64
 
422
FastSumBox(long r, long c, BoxData *D)
 
423
{
 
424
        long i, j;
 
425
        Float64 sum = 0;
 
426
        long cols = D->pix.cols;
 
427
        Float64 *data = D->pix.data;
 
428
        data += r*cols + c;
 
429
        for(i=0; i<D->krows; i++, data += cols-D->kcols)
 
430
                for(j=0; j<D->kcols; j++, data++)
 
431
                        sum += *data;
 
432
        return sum;
 
433
}
 
434
 
 
435
static long bound(long x, long max)
 
436
{
 
437
        if (x < 0) return 0;
 
438
        else if (x > max) return max;
 
439
        else return x;
 
440
}
 
441
 
 
442
static void
 
443
BoxFunc(long rmin, long rmax, long cmin, long cmax, Float64 *output, BoxData *D)
 
444
{
 
445
        long r, c;
 
446
        long    krows2    = D->krows/2;
 
447
        long    kcols2    = D->kcols/2;
 
448
        long    kcolseven = !(D->kcols & 1);
 
449
        long    rows      = D->pix.rows;
 
450
        long    cols      = D->pix.cols;
 
451
 
 
452
        rmin = bound(rmin, rows);
 
453
        rmax = bound(rmax, rows);
 
454
        cmin = bound(cmin, cols);
 
455
        cmax = bound(cmax, cols);
 
456
 
 
457
        for(r=rmin; r<rmax; r++) {  
 
458
                Float64 sum = D->sumbox(r - krows2, cmin - kcols2, D);
 
459
                for(c=cmin; c<cmax; c++) {
 
460
                        output[r*cols + c] = sum;
 
461
                        sum -= D->sumcol(r - krows2, c - kcols2, D);
 
462
                        sum += D->sumcol(r - krows2, c + kcols2 - kcolseven + 1, D);
 
463
                }
 
464
        }
 
465
}
 
466
 
 
467
/*  BoxFuncI computes a boxcar incrementally, using a formula independent of
 
468
    the size of the boxcar.  Each incremental step is based on dropping a
 
469
    whole column of the "back" of the boxcar, and adding in a new column in
 
470
    the "front".  The sums of these columns are further optimized by realizing
 
471
    they can be computed from their counterparts one element above by adding in
 
472
    bottom corners and subtracting out top corners.
 
473
 
 
474
    incremental pixel layout:      B  C     where S is the unknown, and A, B, C are known neighbors
 
475
                                   A  S           each of these refers to the output array
 
476
 
 
477
    S = A + a1 - a0            where a0 and a1 are column vectors with *bottom* elements { a, d }
 
478
    C - B = b1 - b0            where b0 and b1 are column vectors with *top* elements { b, g }
 
479
                                    column vectors and corner elements refer to the input array
 
480
 
 
481
    offset matrix layout:          b      g                   where b is actually in b0
 
482
                                  [b0]   [b1]                       g    "           b1
 
483
                                  [a0] S [a1]                       a    "           a0
 
484
                                   a      d                         d is actually in a1
 
485
 
 
486
    a0 = b0 - b + a            column vector a0 is b0 dropping top element b and adding bottom a
 
487
    a1 = b1 - g + d            column vector a1 is b1 dropping top element g and adding bottom d
 
488
 
 
489
    S = A + (b1 - g + f) - (b0 - b + a)    by substitution
 
490
    S = A + (b1 - b0) - g + d + b - a      rearranging additions
 
491
    S = A + C - B - g + d + b - a          by substitution
 
492
*/
 
493
static void
 
494
BoxFuncI(long rmin, long rmax, long cmin, long cmax, Float64 *output, BoxData *D)
 
495
{
 
496
        long r, c;
 
497
        long krows2    = D->krows/2;
 
498
        long kcols2    = D->kcols/2;
 
499
        long krowseven = !(D->krows & 1);
 
500
        long kcolseven = !(D->kcols & 1);
 
501
        long rows      = D->pix.rows;
 
502
        long cols      = D->pix.cols;
 
503
        Float64 *input = D->pix.data;
 
504
 
 
505
        rmin = bound(rmin, rows);
 
506
        rmax = bound(rmax, rows);
 
507
        cmin = bound(cmin, cols);
 
508
        cmax = bound(cmax, cols);
 
509
 
 
510
        for(r=rmin; r<rmax; r++) {  
 
511
                long top    = r - krows2 - 1;
 
512
                long bottom = r + krows2 - krowseven;
 
513
 
 
514
                for(c=cmin; c<cmax; c++) {
 
515
                        long left   = c - kcols2 - 1;
 
516
                        long right  = c + kcols2 - kcolseven;
 
517
 
 
518
                        Float64 A = output [     r  * cols + (c-1) ];
 
519
                        Float64 B = output [  (r-1) * cols + (c-1) ];
 
520
                        Float64 C = output [  (r-1) * cols + c     ];
 
521
                        Float64 a = input  [ bottom * cols + left  ];
 
522
                        Float64 b = input  [    top * cols + left  ];
 
523
                        Float64 g = input  [    top * cols + right ];
 
524
                        Float64 d = input  [ bottom * cols + right ];
 
525
 
 
526
                        output[r*cols + c] = A + C - B - g + d + b - a;
 
527
                }
 
528
        }
 
529
}
 
530
 
 
531
static void 
 
532
Boxcar2d(long krows, long kcols, long rows, long cols, Float64 *data, 
 
533
         Float64 *output, PixMode mode, Float64 constval)
 
534
{
 
535
        long krows2 = krows/2;
 
536
        long krowseven = !(krows&1);
 
537
        long kcols2 = kcols/2;
 
538
        long kcolseven = !(kcols&1);
 
539
        long r, c;
 
540
        Float64 karea;
 
541
 
 
542
        BoxData D;
 
543
        D.pix.mode = mode;
 
544
        D.pix.constval = constval;
 
545
        D.pix.rows = rows;
 
546
        D.pix.cols = cols;
 
547
        D.pix.data = data;
 
548
        D.krows = krows;
 
549
        D.kcols = kcols;
 
550
        D.sumcol = SlowSumCol;
 
551
        D.sumbox = SlowSumBox;
 
552
        
 
553
        /* The next 4 calls compute boxcars on the boundary pixels with
 
554
           different modes detemining what data values are fetched when "out
 
555
           of bounds".  Presumably, this is a small minority of the data and
 
556
           therefore the implementation inefficiency is not *too* damaging.
 
557
           It should at least beat making 2 outsized copies of the data array. */
 
558
 
 
559
        /* top whole plus one */
 
560
        BoxFunc(0, krows2+2, 0, cols, output, &D);
 
561
        
 
562
        /* bottom whole */
 
563
        BoxFunc(rows-krows2+krowseven, rows, 0, cols, output, &D);
 
564
        
 
565
        /* left whole plus one */
 
566
        BoxFunc(0, rows, 0, kcols2+2, output, &D);
 
567
        
 
568
        /* right whole */
 
569
        BoxFunc(0, rows, cols-kcols2+kcolseven, cols, output, &D);
 
570
 
 
571
        /* Do the boxcar on the "center" data, the data where the boxcar is
 
572
           always "in bounds".  Presumably, this is the bulk of the data so
 
573
           this should be fast.  
 
574
        */
 
575
        D.sumcol = FastSumCol;
 
576
        D.sumbox = FastSumBox;
 
577
 
 
578
        BoxFuncI( krows2+2, rows-krows2+krowseven, 
 
579
                  kcols2+2, cols-kcols2+kcolseven, 
 
580
                  output, &D);
 
581
 
 
582
        karea = kcols * krows;
 
583
        for(r=0; r<rows; r++)
 
584
                for(c=0; c<cols; c++)
 
585
                        output[r*cols + c] /= karea;
 
586
}
 
587
 
 
588
static PyObject *
 
589
Py_Boxcar2d(PyObject *obj, PyObject *args, PyObject *kw)
 
590
{
 
591
        PyObject   *odata, *ooutput=NULL;
 
592
        PyArrayObject *data, *output;
 
593
        int        krows, kcols, mode=PIX_NEAREST;
 
594
        Float64    cval = 0;
 
595
        char       *keywds[] = { "data", "krows", "kcols", 
 
596
                                 "output", "mode", "cval", NULL };
 
597
 
 
598
        if (!PyArg_ParseTupleAndKeywords(args, kw, "Oii|Oid:Boxcar2d", keywds, 
 
599
                         &odata, &krows, &kcols, &ooutput, &mode, &cval))
 
600
                return NULL;
 
601
 
 
602
        /* Align, Byteswap, Contiguous, Typeconvert */
 
603
        data    = NA_InputArray(odata, tFloat64, C_ARRAY);
 
604
        output = NA_OptionalOutputArray(ooutput, tFloat64, C_ARRAY, data);
 
605
        if (!data || !output)
 
606
                goto _fail;
 
607
 
 
608
        if (_reject_complex(odata) || _reject_complex(ooutput))
 
609
                goto _fail;
 
610
 
 
611
        if ((krows < 0) || (kcols < 0)) {
 
612
                PyErr_Format(PyExc_ValueError, "krows and kcols must be > 0.");
 
613
                goto _fail;
 
614
        }
 
615
 
 
616
        if ((mode < PIX_NEAREST) || (mode > PIX_CONSTANT)) {
 
617
                PyErr_Format(PyExc_ValueError,
 
618
                             "Boxcar2d: mode value not in range(%d,%d)", 
 
619
                             PIX_NEAREST, PIX_CONSTANT);
 
620
                goto _fail;
 
621
        }
 
622
 
 
623
        if ((data->nd != 2)|| (output->nd != 2)) {
 
624
                PyErr_Format(PyExc_ValueError,
 
625
                             "Boxcar2d: numarray must have 2 dimensions.");
 
626
                goto _fail;
 
627
        }
 
628
 
 
629
        if (!NA_ShapeEqual(data, output)) {
 
630
                PyErr_Format(PyExc_ValueError,
 
631
                             "Boxcar2d: data and output numarray need identical shapes.");
 
632
                goto _fail;
 
633
        }
 
634
 
 
635
        if ((kcols <=0) || (krows <= 0)) {
 
636
                PyErr_Format(PyExc_ValueError,
 
637
                             "Boxcar2d: invalid data shape.");
 
638
                goto _fail;
 
639
        }
 
640
        if ((kcols > data->dimensions[1]) || (krows > data->dimensions[0])) {
 
641
                PyErr_Format(PyExc_ValueError, "Boxcar2d: boxcar shape incompatible with"
 
642
                             " data shape.");
 
643
                goto _fail;
 
644
        }
 
645
 
 
646
        Boxcar2d(krows, kcols, data->dimensions[0], data->dimensions[1], 
 
647
                 NA_OFFSETDATA(data), NA_OFFSETDATA(output), mode, cval);
 
648
 
 
649
        Py_XDECREF(data);
 
650
 
 
651
        /* Align, Byteswap, Contiguous, Typeconvert */
 
652
        return NA_ReturnOutput(ooutput, output);
 
653
  _fail:
 
654
        Py_XDECREF(data);
 
655
        Py_XDECREF(output);
 
656
        return NULL;
 
657
}
 
658
 
 
659
static PyMethodDef _correlateMethods[] = {
 
660
    {"Correlate1d", Py_Correlate1d, METH_VARARGS}, 
 
661
    {"Correlate2d", (PyCFunction) Py_Correlate2d, METH_VARARGS | METH_KEYWORDS},
 
662
    {"Shift2d", (PyCFunction) Py_Shift2d, METH_VARARGS | METH_KEYWORDS, 
 
663
     "Shift2d shifts and image by an integer number of pixels, and uses IRAF compatible modes for the boundary pixels."},
 
664
    {"Boxcar2d", (PyCFunction) Py_Boxcar2d, METH_VARARGS | METH_KEYWORDS,
 
665
    "Boxcar2d computes a sliding 2D boxcar average on a 2D array"},
 
666
    {NULL, NULL} /* Sentinel */
 
667
};
 
668
 
 
669
PyMODINIT_FUNC init_correlate(void)
 
670
{
 
671
        PyObject *m, *d;
 
672
        m = Py_InitModule("_correlate", _correlateMethods);
 
673
        d = PyModule_GetDict(m);
 
674
        import_libnumarray();
 
675
}
 
676
 
 
677
/*
 
678
 * Local Variables:
 
679
 * mode: C
 
680
 * c-file-style: "python"
 
681
 * End:
 
682
 */