1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5 %module(package="astrometry.util") util
6
7 %include <typemaps.i>
8 %include <cstring.i>
9 %include <exception.i>
10
11 %{
12 // numpy.
13 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
14 #include <numpy/arrayobject.h>
15
16 #include <stdint.h>
17 #include <stdlib.h>
18 #include <math.h>
19
20 #include "os-features.h"
21 #include "log.h"
22 #include "healpix.h"
23 #include "healpix-utils.h"
24 #include "anwcs.h"
25 #include "sip.h"
26 #include "fitsioutils.h"
27 #include "sip-utils.h"
28 #include "sip_qfits.h"
29 #include "index.h"
30 #include "quadfile.h"
31 #include "codekd.h"
32 #include "starkd.h"
33 #include "starutil.h"
34 #include "an-bool.h"
35 #include "ioutils.h"
36
37 #include "coadd.h"
38 #include "wcs-resample.h"
39 #include "resample.h"
40 #include "keywords.h"
41
42 #include "dimage.h"
43
44 #include "fit-wcs.h"
45
46 #include "qfits_header.h"
47 #include "qfits_rw.h"
48 #include "wcs-pv2sip.h"
49
50 #define true 1
51 #define false 0
52
53 // For sip.h
checkorder(int i,int j)54 static void checkorder(int i, int j) {
55 assert(i >= 0);
56 assert(i < SIP_MAXORDER);
57 assert(j >= 0);
58 assert(j < SIP_MAXORDER);
59 }
60
61 // From index.i:
62 /**
63 For returning single codes and quads as python lists, do something like this:
64
65 %typemap(out) float [ANY] {
66 int i;
67 $result = PyList_New($1_dim0);
68 for (i = 0; i < $1_dim0; i++) {
69 PyObject *o = PyFloat_FromDouble((double) $1[i]);
70 PyList_SetItem($result,i,o);
71 }
72 }
73 **/
74
code_alloc(int DC)75 double* code_alloc(int DC) {
76 return malloc(DC * sizeof(double));
77 }
code_free(double * code)78 void code_free(double* code) {
79 free(code);
80 }
code_get(double * code,int i)81 double code_get(double* code, int i) {
82 return code[i];
83 }
84
codekd_addr(index_t * ind)85 long codekd_addr(index_t* ind) {
86 return (long)ind->codekd;
87 }
starkd_addr(index_t * ind)88 long starkd_addr(index_t* ind) {
89 return (long)ind->starkd;
90 }
91
quadfile_addr(index_t * ind)92 long quadfile_addr(index_t* ind) {
93 return (long)ind->quads;
94 }
95 /*
96 long qidxfile_addr(qidxfile* qf) {
97 return (long)qf;
98 }
99 */
100
101 %}
102
103 %init %{
104 // numpy
105 import_array();
106 %}
107
108 // Things in keywords.h (used by healpix.h)
109 #define Const
110 #define WarnUnusedResult
111 #define InlineDeclare
112 #define Flatten
113 #define ASTROMETRY_KEYWORDS_H
114 #define ATTRIB_FORMAT(x,y,z)
115
116 void log_init(int level);
117 int log_get_level();
118 void log_set_level(int lvl);
119
120 %include "coadd.h"
121 %include "resample.h"
122 %include "an-bool.h"
123 %include "fit-wcs.h"
124
125 %inline %{
126 #define ERR(x, ...) \
127 printf(x, ## __VA_ARGS__)
128
print_array(PyObject * arr)129 static void print_array(PyObject* arr) {
130 PyArrayObject *obj;
131 int i, N;
132 PyArray_Descr *desc;
133 printf("Array: %p\n", arr);
134 if (!arr) return;
135 if (!PyArray_Check(arr)) {
136 printf(" Not a Numpy Array\n");
137 if (arr == Py_None)
138 printf(" is None\n");
139 return;
140 }
141 obj = (PyArrayObject*)arr;
142
143 printf(" Contiguous: %s\n",
144 PyArray_ISCONTIGUOUS(obj) ? "yes" : "no");
145 printf(" Writeable: %s\n",
146 PyArray_ISWRITEABLE(obj) ? "yes" : "no");
147 printf(" Aligned: %s\n",
148 PyArray_ISALIGNED(obj) ? "yes" : "no");
149 printf(" C array: %s\n",
150 PyArray_ISCARRAY(obj) ? "yes" : "no");
151
152 //printf(" typeobj: %p (float is %p)\n", arr->typeobj,
153 //&PyFloat_Type);
154
155 printf(" data: %p\n", PyArray_DATA(obj));
156 printf(" N dims: %i\n", PyArray_NDIM(obj));
157 N = PyArray_NDIM(obj);
158 for (i=0; i<N; i++)
159 printf(" dim %i: %i\n", i, (int)PyArray_DIM(obj, i));
160 for (i=0; i<N; i++)
161 printf(" stride %i: %i\n", i, (int)PyArray_STRIDE(obj, i));
162 desc = PyArray_DESCR(obj);
163 printf(" descr kind: '%c'\n", desc->kind);
164 printf(" descr type: '%c'\n", desc->type);
165 printf(" descr byteorder: '%c'\n", desc->byteorder);
166 printf(" descr elsize: %i\n", desc->elsize);
167 }
168
169
an_hist2d(PyObject * py_arrx,PyObject * py_arry,PyObject * py_hist,double xlo,double xhi,double ylo,double yhi)170 static PyObject* an_hist2d(PyObject* py_arrx, PyObject* py_arry,
171 PyObject* py_hist,
172 double xlo, double xhi,
173 double ylo, double yhi) {
174 PyArray_Descr* dtype = NULL;
175 PyArray_Descr* itype = NULL;
176 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED |
177 NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ELEMENTSTRIDES;
178 int reqout = req | NPY_ARRAY_WRITEABLE | NPY_ARRAY_UPDATEIFCOPY;
179 PyArrayObject* np_arrx;
180 PyArrayObject* np_arry;
181 PyArrayObject* np_hist;
182 double *arrx;
183 double *arry;
184 int32_t *hist;
185 int nx, ny;
186 double dx, dy, idx, idy;
187 int i, N;
188
189 dtype = PyArray_DescrFromType(NPY_DOUBLE);
190 itype = PyArray_DescrFromType(NPY_INT32);
191
192 Py_INCREF(dtype);
193 np_arrx = (PyArrayObject*)PyArray_FromAny(py_arrx, dtype, 1, 1, req, NULL);
194 if (!np_arrx) {
195 PyErr_SetString(PyExc_ValueError,"Expected x array to be double");
196 Py_DECREF(dtype);
197 return NULL;
198 }
199 N = PyArray_SIZE(np_arrx);
200
201 Py_INCREF(dtype);
202 np_arry = (PyArrayObject*)PyArray_FromAny(py_arry, dtype, 1, 1, req, NULL);
203 if (!np_arry) {
204 PyErr_SetString(PyExc_ValueError,"Expected y array to be double");
205 Py_DECREF(dtype);
206 Py_DECREF(np_arrx);
207 return NULL;
208 }
209 if (PyArray_SIZE(np_arry) != N) {
210 PyErr_SetString(PyExc_ValueError,"Expected x and y arrays to be the same length");
211 Py_DECREF(dtype);
212 Py_DECREF(np_arrx);
213 return NULL;
214 }
215 Py_CLEAR(dtype);
216 Py_INCREF(itype);
217 np_hist = (PyArrayObject*)PyArray_FromAny(py_hist, itype, 2, 2, reqout, NULL);
218 if (!np_hist) {
219 PyErr_SetString(PyExc_ValueError,"Expected hist array to be int32");
220 Py_DECREF(np_arrx);
221 Py_DECREF(np_arry);
222 Py_DECREF(itype);
223 return NULL;
224 }
225 Py_CLEAR(itype);
226
227 ny = PyArray_DIM(np_hist, 0);
228 nx = PyArray_DIM(np_hist, 1);
229 dx = (xhi - xlo) / (double)nx;
230 dy = (yhi - ylo) / (double)ny;
231 idx = 1./dx;
232 idy = 1./dy;
233
234 hist = PyArray_DATA(np_hist);
235 arrx = PyArray_DATA(np_arrx);
236 arry = PyArray_DATA(np_arry);
237
238 for (i=0; i<N; i++,
239 arrx++, arry++) {
240 double x, y;
241 int binx, biny;
242 x = *arrx;
243 y = *arry;
244 if ((x < xlo) || (x > xhi) || (y < ylo) || (y > yhi))
245 continue;
246 binx = (int)((x - xlo) * idx);
247 biny = (int)((y - ylo) * idy);
248 // == upper limit
249 if (unlikely(binx == nx)) {
250 binx--;
251 }
252 if (unlikely(biny == ny)) {
253 biny--;
254 }
255 hist[biny * nx + binx]++;
256 }
257
258 Py_DECREF(np_arrx);
259 Py_DECREF(np_arry);
260 Py_DECREF(np_hist);
261
262 Py_RETURN_NONE;
263 }
264
265
flat_percentile_f(PyObject * np_arr,double pct)266 static double flat_percentile_f(PyObject* np_arr, double pct) {
267 PyArray_Descr* dtype;
268 npy_intp N;
269 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED |
270 NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ELEMENTSTRIDES;
271 float* x;
272 float med = 0;
273 int L, R;
274 int mid;
275
276 dtype = PyArray_DescrFromType(NPY_FLOAT);
277 np_arr = PyArray_CheckFromAny(np_arr, dtype, 0, 0, req, NULL);
278 if (!np_arr) {
279 ERR("flat_median_f: Failed to convert array to float\n");
280 return 0;
281 }
282 dtype = NULL;
283 N = PyArray_Size(np_arr);
284 x = (float*)malloc(sizeof(float) * N);
285 memcpy(x, PyArray_DATA((PyArrayObject*)np_arr), sizeof(float)*N);
286 Py_DECREF(np_arr);
287
288 {
289 int i;
290 for (i=0; i<N; i++) {
291 if (!isfinite(x[i])) {
292 ERR("flat_median_f cannot handle NaN values (element %i)\n", i);
293 return x[i];
294 }
295 }
296 }
297
298 // Pseudocode from wikipedia's 'Selection algorithm' page
299 L = 0;
300 R = (int)(N-1);
301 mid = (int)(pct * 0.01 * N);
302 if (mid < 0) {
303 mid = 0;
304 }
305 if (mid >= R) {
306 mid = R;
307 }
308 while (L < R) {
309 int ipivot;
310 int i,j;
311 int k;
312 float pivot;
313 //printf("L=%i, R=%i (N=%i), mid=%i\n", L, R, 1+R-L, mid);
314 ipivot = random() % (1+R-L) + L;
315 pivot = x[ipivot];
316 // partition array...
317 i = L;
318 j = R;
319 do {
320 // scan for elements out of place
321 // scan from the left:
322 while (x[i] < pivot)
323 i++;
324 // scan from the right:
325 while (x[j] >= pivot && j>i)
326 j--;
327 // now x[i] >= pivot
328 // and (x[j] < pivot) OR j == i
329 assert(x[i] >= pivot);
330 assert((x[j] < pivot) || (j == i));
331 assert(j >= i);
332 if (i < j) {
333 // swap
334 float tmp = x[i];
335 x[i] = x[j];
336 x[j] = tmp;
337 }
338 } while (i < j);
339 {
340 for (k=L; k<i; k++) {
341 assert(x[k] < pivot);
342 }
343 for (k=i; k<=R; k++) {
344 assert(x[k] >= pivot);
345 }
346 }
347 // partition the right partition into == and >
348 j = i;
349 k = R;
350 do {
351 // scan for elements out of place
352 // scan from the right:
353 while (x[k] > pivot)
354 k--;
355 // scan from the left:
356 while (x[j] == pivot && j<k)
357 j++;
358
359 assert(x[k] == pivot);
360 assert((x[j] > pivot) || (j == k));
361 assert(k >= j);
362 if (j < k) {
363 // swap
364 float tmp = x[j];
365 x[j] = x[k];
366 x[k] = tmp;
367 }
368 } while (j < k);
369
370 j = k+1;
371
372 {
373 //printf("L=%i, i=%i, j=%i, k=%i, R=%i\n", L, i, j, k, R);
374 for (k=L; k<i; k++) {
375 assert(x[k] < pivot);
376 }
377 for (k=i; k<j; k++) {
378 assert(x[k] == pivot);
379 }
380 for (k=j; k<=R; k++) {
381 assert(x[k] > pivot);
382 }
383 }
384
385
386
387 // there must be at least one element in the right partitions
388 assert(i <= R);
389
390 // there must be at least one element in the middle partition
391 assert(j-i >= 1);
392
393 if (mid < i)
394 // the median is in the left partition (< pivot)
395 R = i-1;
396 else if (mid >= j)
397 // the median is in the right partition (> pivot)
398 L = j;
399 else {
400 // the median is in the middle partition (== pivot)
401 L = R = i;
402 break;
403 }
404 assert(L <= mid);
405 assert(R >= mid);
406 }
407 med = x[mid];
408 free(x);
409 return med;
410 }
411
flat_median_f(PyObject * np_arr)412 static double flat_median_f(PyObject* np_arr) {
413 return flat_percentile_f(np_arr, 50.0);
414 }
415
median_smooth(PyObject * py_image,PyObject * py_mask,int halfbox,PyObject * py_smooth)416 static int median_smooth(PyObject* py_image,
417 PyObject* py_mask,
418 int halfbox,
419 PyObject* py_smooth) {
420 /*
421
422 image: np.float32
423 mask: bool or uint8; 1 to IGNORE.
424 smooth: np.float32; output array.
425
426 */
427 PyArrayObject* np_image = (PyArrayObject*)py_image;
428 PyArrayObject* np_mask = (PyArrayObject*)py_mask;
429 PyArrayObject* np_smooth = (PyArrayObject*)py_smooth;
430
431 if (!PyArray_Check(np_image) ||
432 !PyArray_Check(np_smooth) ||
433 !PyArray_ISNOTSWAPPED(np_image) ||
434 !PyArray_ISNOTSWAPPED(np_smooth ) ||
435 !PyArray_ISFLOAT(np_image) ||
436 !PyArray_ISFLOAT(np_smooth ) ||
437 (PyArray_ITEMSIZE(np_image) != sizeof(float)) ||
438 (PyArray_ITEMSIZE(np_smooth ) != sizeof(float)) ||
439 !(PyArray_NDIM(np_image) == 2) ||
440 !(PyArray_NDIM(np_smooth ) == 2) ||
441 !PyArray_ISCONTIGUOUS(np_image) ||
442 !PyArray_ISCONTIGUOUS(np_smooth ) ||
443 !PyArray_ISWRITEABLE(np_smooth)) {
444 ERR("median_smooth: array type checks failed for image/smooth\n");
445 ERR("check %i %i notswapped %i %i isfloat %i %i size %i %i ndim %i %i contig %i %i writable %i\n",
446 PyArray_Check(np_image), PyArray_Check(np_smooth),
447 PyArray_ISNOTSWAPPED(np_image), PyArray_ISNOTSWAPPED(np_smooth ),
448 PyArray_ISFLOAT(np_image), PyArray_ISFLOAT(np_smooth),
449 (PyArray_ITEMSIZE(np_image) == sizeof(float)),
450 (PyArray_ITEMSIZE(np_smooth) == sizeof(float)),
451 (PyArray_NDIM(np_image) == 2),
452 (PyArray_NDIM(np_smooth ) == 2),
453 PyArray_ISCONTIGUOUS(np_image),
454 PyArray_ISCONTIGUOUS(np_smooth),
455 PyArray_ISWRITEABLE(np_smooth));
456 return -1;
457 }
458 if ((PyObject*)np_mask != Py_None) {
459 if (!PyArray_Check(np_mask) ||
460 !PyArray_ISNOTSWAPPED(np_mask) ||
461 !PyArray_ISBOOL(np_mask) ||
462 (PyArray_ITEMSIZE(np_mask) != sizeof(uint8_t)) ||
463 !(PyArray_NDIM(np_mask) == 2) ||
464 !PyArray_ISCONTIGUOUS(np_mask)) {
465 ERR("median_smooth: array type checks failed for mask\n");
466 return -1;
467 }
468 }
469 npy_intp NX, NY;
470 const float* image;
471 float* smooth;
472 const uint8_t* maskimg = NULL;
473
474 NY = PyArray_DIM(np_image, 0);
475 NX = PyArray_DIM(np_image, 1);
476 if ((PyArray_DIM(np_smooth, 0) != NY) ||
477 (PyArray_DIM(np_smooth, 1) != NX)) {
478 ERR("median_smooth: 'smooth' array is wrong shape\n");
479 return -1;
480 }
481 image = PyArray_DATA(np_image);
482 smooth = PyArray_DATA(np_smooth);
483
484 if ((PyObject*)np_mask != Py_None) {
485 if ((PyArray_DIM(np_mask, 0) != NY) ||
486 (PyArray_DIM(np_mask, 1) != NX)) {
487 ERR("median_smooth: 'mask' array is wrong shape\n");
488 return -1;
489 }
490 maskimg = PyArray_DATA(np_mask);
491 }
492
493 dmedsmooth(image, maskimg, (int)NX, (int)NY, halfbox, smooth);
494
495 return 0;
496 }
497
498 #define LANCZOS_INTERP_FUNC lanczos5_interpolate
499 #define L 5
500 static int LANCZOS_INTERP_FUNC(PyObject* np_ixi, PyObject* np_iyi,
501 PyObject* np_dx, PyObject* np_dy,
502 PyObject* loutputs, PyObject* linputs);
503 #include "lanczos.i"
504 #undef LANCZOS_INTERP_FUNC
505 #undef L
506
507 #define LANCZOS_INTERP_FUNC lanczos3_interpolate
508 #define L 3
509 static int LANCZOS_INTERP_FUNC(PyObject* np_ixi, PyObject* np_iyi,
510 PyObject* np_dx, PyObject* np_dy,
511 PyObject* loutputs, PyObject* linputs);
512 #include "lanczos.i"
513 #undef LANCZOS_INTERP_FUNC
514 #undef L
515
lanczos5_filter(PyObject * py_dx,PyObject * py_f)516 static int lanczos5_filter(PyObject* py_dx, PyObject* py_f) {
517 npy_intp N;
518 npy_intp i;
519 float* dx;
520 float* f;
521
522 PyArrayObject *np_dx = (PyArrayObject*)py_dx;
523 PyArrayObject *np_f = (PyArrayObject*)py_f;
524
525 if (!PyArray_Check(np_dx) ||
526 !PyArray_Check(np_f ) ||
527 !PyArray_ISNOTSWAPPED(np_dx) ||
528 !PyArray_ISNOTSWAPPED(np_f ) ||
529 !PyArray_ISFLOAT(np_dx) ||
530 !PyArray_ISFLOAT(np_f ) ||
531 (PyArray_ITEMSIZE(np_dx) != sizeof(float)) ||
532 (PyArray_ITEMSIZE(np_f ) != sizeof(float)) ||
533 !(PyArray_NDIM(np_dx) == 1) ||
534 !(PyArray_NDIM(np_f ) == 1) ||
535 !PyArray_ISCONTIGUOUS(np_dx) ||
536 !PyArray_ISCONTIGUOUS(np_f ) ||
537 !PyArray_ISWRITEABLE(np_f)
538 ) {
539 ERR("Arrays aren't right type\n");
540 return -1;
541 }
542 N = PyArray_DIM(np_dx, 0);
543 if (PyArray_DIM(np_f, 0) != N) {
544 ERR("Input and output must have same dimensions\n");
545 return -1;
546 }
547 dx = PyArray_DATA(np_dx);
548 f = PyArray_DATA(np_f);
549 const double fifthpi = M_PI / 5.0;
550 const double pisq = M_PI * M_PI;
551 const double fiveopisq = 5. / pisq;
552 for (i=N; i>0; i--, dx++, f++) {
553 double x = *dx;
554 if (x < -5.0 || x > 5.0) {
555 *f = 0.0;
556 } else if (x == 0) {
557 *f = 1.0;
558 } else {
559 *f = fiveopisq * sin(M_PI * x) * sin(fifthpi * x) / (x * x);
560 }
561 }
562 return 0;
563 }
564
lanczos3_filter(PyObject * py_dx,PyObject * py_f)565 static int lanczos3_filter(PyObject* py_dx, PyObject* py_f) {
566 npy_intp N;
567 npy_intp i;
568 float* dx;
569 float* f;
570
571 PyArrayObject *np_dx = (PyArrayObject*)py_dx;
572 PyArrayObject *np_f = (PyArrayObject*)py_f;
573
574 if (!PyArray_Check(np_dx) ||
575 !PyArray_Check(np_f ) ||
576 !PyArray_ISNOTSWAPPED(np_dx) ||
577 !PyArray_ISNOTSWAPPED(np_f ) ||
578 !PyArray_ISFLOAT(np_dx) ||
579 !PyArray_ISFLOAT(np_f ) ||
580 (PyArray_ITEMSIZE(np_dx) != sizeof(float)) ||
581 (PyArray_ITEMSIZE(np_f ) != sizeof(float)) ||
582 !(PyArray_NDIM(np_dx) == 1) ||
583 !(PyArray_NDIM(np_f ) == 1) ||
584 !PyArray_ISCONTIGUOUS(np_dx) ||
585 !PyArray_ISCONTIGUOUS(np_f ) ||
586 !PyArray_ISWRITEABLE(np_f)
587 ) {
588 ERR("Arrays aren't right type\n");
589 return -1;
590 }
591 N = PyArray_DIM(np_dx, 0);
592 if (PyArray_DIM(np_f, 0) != N) {
593 ERR("Input and output must have same dimensions\n");
594 return -1;
595 }
596 dx = PyArray_DATA(np_dx);
597 f = PyArray_DATA(np_f);
598 const double thirdpi = M_PI / 3.0;
599 const double pisq = M_PI * M_PI;
600 const double threeopisq = 3. / pisq;
601 for (i=N; i>0; i--, dx++, f++) {
602 double x = *dx;
603 if (x < -3.0 || x > 3.0) {
604 *f = 0.0;
605 } else if (x == 0) {
606 *f = 1.0;
607 } else {
608 *f = threeopisq * sin(M_PI * x) * sin(thirdpi * x) / (x * x);
609 }
610 }
611 return 0;
612 }
613
lanczos3_filter_table(PyObject * py_dx,PyObject * py_f,int rangecheck)614 static int lanczos3_filter_table(PyObject* py_dx, PyObject* py_f, int rangecheck) {
615 npy_intp N;
616 npy_intp i;
617 float* dx;
618 float* f;
619
620 PyArrayObject *np_dx = (PyArrayObject*)py_dx;
621 PyArrayObject *np_f = (PyArrayObject*)py_f;
622
623 // Nlutunit is number of bins per unit x
624 static const int Nlutunit = 1024;
625 static const float lut0 = -4.;
626 static const int Nlut = 8192; //8 * Nlutunit;
627 // We want bins to go from -4 to 4 (Lanczos-3 range of -3 to 3, plus some buffer)
628 // [Nlut]
629 static float lut[8192];
630 static int initialized = 0;
631
632 if (!initialized) {
633 for (i=0; i<(Nlut); i++) {
634 float x,f;
635 x = (lut0 + (i / (float)Nlutunit));
636 if (x <= -3.0 || x >= 3.0) {
637 f = 0.0;
638 } else if (x == 0) {
639 f = 1.0;
640 } else {
641 f = 3. * sin(M_PI * x) * sin(M_PI / 3.0 * x) / (M_PI * M_PI * x * x);
642 }
643 lut[i] = f;
644 }
645 initialized = 1;
646 }
647
648 if (!PyArray_Check(np_dx) ||
649 !PyArray_Check(np_f )) {
650 ERR("Array check\n");
651 }
652 if (!PyArray_ISNOTSWAPPED(np_dx) ||
653 !PyArray_ISNOTSWAPPED(np_f )) {
654 ERR("Swapped\n");
655 }
656 if (!PyArray_ISFLOAT(np_dx) ||
657 !PyArray_ISFLOAT(np_f )) {
658 ERR("Float\n");
659 }
660 if ((PyArray_ITEMSIZE(np_dx) != sizeof(float)) ||
661 (PyArray_ITEMSIZE(np_f ) != sizeof(float))) {
662 ERR("sizeof float\n");
663 }
664 if ((PyArray_ITEMSIZE(np_dx) != sizeof(float))) {
665 ERR("sizeof dx %i\n", (int)PyArray_ITEMSIZE(np_dx));
666 }
667 if ((PyArray_ITEMSIZE(np_f ) != sizeof(float))) {
668 ERR("sizeof f %i\n", (int)PyArray_ITEMSIZE(np_f));
669 }
670 if (!(PyArray_NDIM(np_dx) == 1) ||
671 !(PyArray_NDIM(np_f ) == 1)) {
672 ERR("one-d\n");
673 }
674 if (!PyArray_ISCONTIGUOUS(np_dx) ||
675 !PyArray_ISCONTIGUOUS(np_f )) {
676 ERR("contig\n");
677 }
678 if (!PyArray_ISWRITEABLE(np_f)) {
679 ERR("writable\n");
680 }
681
682
683 if (!PyArray_Check(np_dx) ||
684 !PyArray_Check(np_f ) ||
685 !PyArray_ISNOTSWAPPED(np_dx) ||
686 !PyArray_ISNOTSWAPPED(np_f ) ||
687 !PyArray_ISFLOAT(np_dx) ||
688 !PyArray_ISFLOAT(np_f ) ||
689 (PyArray_ITEMSIZE(np_dx) != sizeof(float)) ||
690 (PyArray_ITEMSIZE(np_f ) != sizeof(float)) ||
691 !(PyArray_NDIM(np_dx) == 1) ||
692 !(PyArray_NDIM(np_f ) == 1) ||
693 !PyArray_ISCONTIGUOUS(np_dx) ||
694 !PyArray_ISCONTIGUOUS(np_f ) ||
695 !PyArray_ISWRITEABLE(np_f)
696 ) {
697 ERR("Arrays aren't right type\n");
698 return -1;
699 }
700 N = PyArray_DIM(np_dx, 0);
701 if (PyArray_DIM(np_f, 0) != N) {
702 ERR("Input and output must have same dimensions\n");
703 return -1;
704 }
705 dx = PyArray_DATA(np_dx);
706 f = PyArray_DATA(np_f);
707 if (rangecheck) {
708 for (i=N; i>0; i--, dx++, f++) {
709 float x = *dx;
710 int li = (int)((x - lut0) * Nlutunit);
711 if ((li < 0) || (li >= Nlut)) {
712 *f = 0.0;
713 } else {
714 *f = lut[li];
715 }
716 }
717 } else {
718 for (i=N; i>0; i--, dx++, f++) {
719 float x = *dx;
720 int li = (int)((x - lut0) * Nlutunit);
721 *f = lut[li];
722 }
723 }
724 return 0;
725 }
726
lanczos_shift_image_c(PyObject * py_img,PyObject * py_weight,PyObject * py_outimg,PyObject * py_outweight,int order,double dx,double dy)727 static int lanczos_shift_image_c(PyObject* py_img,
728 PyObject* py_weight,
729 PyObject* py_outimg,
730 PyObject* py_outweight,
731 int order, double dx, double dy) {
732 int W,H;
733 int i,j;
734
735 lanczos_args_t lanczos;
736
737 PyArray_Descr* dtype;
738 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED |
739 NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ELEMENTSTRIDES;
740 int reqout = req | NPY_ARRAY_WRITEABLE | NPY_ARRAY_UPDATEIFCOPY;
741 double *img, *weight, *outimg, *outweight;
742
743 PyArrayObject *np_img=NULL, *np_weight=NULL, *np_outimg=NULL, *np_outweight=NULL;
744
745 weight = NULL;
746 outweight = NULL;
747 lanczos.order = order;
748
749 /*
750 printf("np_img:\n");
751 print_array(np_img);
752 printf("np_weight:\n");
753 print_array(np_weight);
754 printf("np_outimg:\n");
755 print_array(np_outimg);
756 printf("np_outweight:\n");
757 print_array(np_outweight);
758 */
759
760 dtype = PyArray_DescrFromType(NPY_DOUBLE);
761 Py_INCREF(dtype);
762 np_img = (PyArrayObject*)PyArray_CheckFromAny(py_img, dtype, 2, 2, req, NULL);
763 if (py_weight != Py_None) {
764 Py_INCREF(dtype);
765 np_weight = (PyArrayObject*)PyArray_CheckFromAny(py_weight, dtype, 2, 2, req, NULL);
766 if (!np_weight) {
767 ERR("Failed to run PyArray_FromAny on np_weight\n");
768 return -1;
769 }
770 }
771 Py_INCREF(dtype);
772 np_outimg = (PyArrayObject*)PyArray_CheckFromAny(py_outimg, dtype, 2, 2, reqout, NULL);
773 if (py_outweight != Py_None) {
774 Py_INCREF(dtype);
775 np_outweight = (PyArrayObject*)PyArray_CheckFromAny(py_outweight, dtype, 2, 2, reqout, NULL);
776 }
777 Py_DECREF(dtype);
778 dtype = NULL;
779
780 if (!np_img || !np_outimg ||
781 ((py_outweight != Py_None) && !np_outweight)) {
782 ERR("Failed to PyArray_FromAny the images (np_img=%p, np_outimg=%p, np_outweight=%p)\n",
783 np_img, np_outimg, np_outweight);
784 return -1;
785 }
786
787 H = (int)PyArray_DIM(np_img, 0);
788 W = (int)PyArray_DIM(np_img, 1);
789
790 if ((PyArray_DIM(np_outimg, 0) != H) ||
791 (PyArray_DIM(np_outimg, 1) != W)) {
792 ERR("All images must have the same dimensions.\n");
793 return -1;
794 }
795 if (np_weight) {
796 if ((PyArray_DIM(np_weight, 0) != H) ||
797 (PyArray_DIM(np_weight, 1) != W)) {
798 ERR("All images must have the same dimensions.\n");
799 return -1;
800 }
801 weight = PyArray_DATA(np_weight);
802 }
803 if (np_outweight) {
804 if ((PyArray_DIM(np_outweight, 0) != H) ||
805 (PyArray_DIM(np_outweight, 1) != W)) {
806 ERR("All images must have the same dimensions.\n");
807 return -1;
808 }
809 outweight = PyArray_DATA(np_outweight);
810 }
811
812 /*
813 printf("np_img:\n");
814 print_array(np_img);
815 printf("np_weight:\n");
816 print_array(np_weight);
817 printf("np_outimg:\n");
818 print_array(np_outimg);
819 printf("np_outweight:\n");
820 print_array(np_outweight);
821 printf("weight = %p, outweight = %p\n", weight, outweight);
822 */
823
824 img = PyArray_DATA(np_img);
825 outimg = PyArray_DATA(np_outimg);
826
827 for (i=0; i<H; i++) {
828 for (j=0; j<W; j++) {
829 double wt, val;
830 double px, py;
831 px = j - dx;
832 py = i - dy;
833 val = lanczos_resample_d(px, py, img, weight, W, H, &wt,
834 &lanczos);
835 //printf("pixel %i,%i: wt %g\n", j, i, wt);
836 if (outweight) {
837 outimg[i*W + j] = val;
838 outweight[i*W + j] = wt;
839 } else {
840 outimg[i*W + j] = val / wt;
841 }
842 }
843 }
844
845 /*
846 if (np_img != Py_None) {
847 Py_XDECREF(np_img);
848 }
849 if (np_weight != Py_None) {
850 Py_XDECREF(np_weight);
851 }
852 if (np_outweight != Py_None) {
853 Py_XDECREF(np_outweight);
854 }
855 if (np_outimg != Py_None) {
856 Py_XDECREF(np_outimg);
857 }
858 */
859 return 0;
860 }
861 %}
862
863 %pythoncode %{
864
865 def lanczos_shift_image(img, dx, dy, order=3, weight=None,
866 outimg=None, outweight=None):
867 img = img.astype(float)
868 if weight is not None:
869 weight = weight.astype(float)
870 assert(img.shape == weight.shape)
871 if outimg is None:
872 outimg = np.zeros_like(img)
873 if outweight is not None:
874 assert(outweight.shape == img.shape)
875
876 # print 'outweight:', outweight
877
878 lanczos_shift_image_c(img, weight, outimg, outweight, order, dx, dy)
879 if outweight is None:
880 return outimg
881 return outimg,outweight
882 %}
883
884 // for quadfile_get_stars(quadfile* qf, int quadid, unsigned int* stars)
885 // --> list of stars
886 // swap the int* neighbours arg for tempneigh
887 %typemap(in, numinputs=0) unsigned int *stars (unsigned int tempstars[DQMAX]) {
888 $1 = tempstars;
889 }
890 // in the argout typemap we don't know about the swap (but that's ok)
891 %typemap(argout) (const quadfile_t* qf, unsigned int quadid, unsigned int *stars) {
892 int i;
893 int D;
894 if (result == -1) {
895 goto fail;
896 }
897 D = $1->dimquads;
898 $result = PyList_New(D);
899 for (i = 0; i < D; i++) {
900 PyObject *o = PyInt_FromLong($3[i]);
901 PyList_SetItem($result, i, o);
902 }
903 }
904
905
906 /**
907 double* startree_get_data_column(startree_t* s, const char* colname, const int* indices, int N);
908 -> list of doubles.
909 -> ASSUME indices = None
910 */
911 %typemap(argout) (startree_t* s, const char* colname, const int* indices, int N) {
912 int i;
913 int N;
914 if (!result) {
915 goto fail;
916 }
917 N = $4;
918 $result = PyList_New(N);
919 for (i = 0; i < N; i++) {
920 PyObject *o = PyFloat_FromDouble(result[i]);
921 PyList_SetItem($result, i, o);
922 }
923 free(result);
924 }
925
926
927 %include "index.h"
928 %include "quadfile.h"
929 %include "codekd.h"
930 %include "starkd.h"
931 //%include "qidxfile.h"
932
933 double* code_alloc(int DC);
934 void code_free(double* code);
935 double code_get(double* code, int i);
936 long codekd_addr(index_t* ind);
937 long starkd_addr(index_t* ind);
938 long quadfile_addr(index_t* ind);
939 //long qidxfile_addr(qidxfile* qf);
940
941 %apply double *OUTPUT { double *dx, double *dy };
942 %apply double *OUTPUT { double *ra, double *dec };
943
944 // healpix_to_xyz
945 %apply double *OUTPUT { double *p_x, double *p_y, double *p_z };
946
947 // for int healpix_get_neighbours(int hp, int* neigh, int nside)
948 // --> list of neigh
949 // swap the int* neighbours arg for tempneigh
950 %typemap(in, numinputs=0) int *neighbours (int tempneigh[8]) {
951 $1 = tempneigh;
952 }
953 // in the argout typemap we don't know about the swap (but that's ok)
954 %typemap(argout) int *neighbours {
955 int i;
956 int nn;
957 // convert $result to nn
958 //nn = (int)PyInt_AsLong($result);
959 nn = result;
960 $result = PyList_New(nn);
961 for (i = 0; i < nn; i++) {
962 PyObject *o = PyInt_FromLong($1[i]);
963 PyList_SetItem($result, i, o);
964 }
965 }
966
967
968 // for il* healpix_rangesearch_radec(ra, dec, double, int nside, il* hps);
969 // --> list
970 // swallow the int* hps arg
971 %typemap(in, numinputs=0) il* hps {
972 $1 = NULL;
973 }
974 %typemap(out) il* {
975 int i;
976 int N;
977 N = il_size($1);
978 $result = PyList_New(N);
979 for (i = 0; i < N; i++) {
980 PyObject *o = PyInt_FromLong(il_get($1, i));
981 PyList_SetItem($result, i, o);
982 }
983 }
984
985 // healpix_radec_bounds
986 %apply double *OUTPUT { double *ralo, double *rahi, double *declo, double *dechi };
987
988 // xyztohealpixf
989 %apply double *OUTPUT { double *p_dx, double *p_dy };
990
991 %include "healpix.h"
992 %include "healpix-utils.h"
993
994
995 // anwcs_get_radec_center_and_radius
996 %apply double *OUTPUT { double *p_ra, double *p_dec, double *p_radius };
997
998 // anwcs_get_radec_bounds
999 %apply double *OUTPUT { double* pramin, double* pramax, double* pdecmin, double* pdecmax };
1000
1001 %apply double *OUTPUT { double *p_x, double *p_y, double *p_z };
1002 %apply double *OUTPUT { double *p_ra, double *p_dec };
1003 //%apply double *OUTPUT { double *xyz };
1004
1005 // eg anwcs_radec2pixelxy
1006 %apply double *OUTPUT { double *p_x, double *p_y };
1007
1008 // anwcs_pixelxy2xyz
1009 %typemap(in, numinputs=0) double* p_xyz (double tempxyz[3]) {
1010 $1 = tempxyz;
1011 }
1012 // in the argout typemap we don't know about the swap (but that's ok)
1013 %typemap(argout) double* p_xyz {
1014 $result = Py_BuildValue("(ddd)", $1[0], $1[1], $1[2]);
1015 }
1016
1017 // anwcs_get_cd_matrix
1018 %typemap(in, numinputs=0) double* p_cd (double tempcd[4]) {
1019 $1 = tempcd;
1020 }
1021 %typemap(argout) double* p_cd {
1022 $result = Py_BuildValue("(dddd)", $1[0], $1[1], $1[2], $1[3]);
1023 }
1024
1025
1026 %typemap(in, numinputs=0) char **stringparam (char* tempstr) {
1027 $1 = &tempstr;
1028 }
1029 %typemap(in, numinputs=0) int *stringsizeparam (int slen) {
1030 $1 = &slen;
1031 }
1032 char* anwcs_wcslib_to_string(const anwcs_t* wcs,
1033 char **stringparam, int *stringsizeparam);
1034
1035 %ignore anwcs_wcslib_to_string;
1036
1037 %include "anwcs.h"
1038
1039 %extend anwcs_t {
1040 anwcs_t(char* fn, int ext=0, int slen=0) {
1041 if ((ext == -1) ||
1042 (starts_with(fn, "SIMPLE =") && !file_exists(fn))) {
1043 // assume header string
1044 if (slen == 0) {
1045 slen = (int)strlen(fn);
1046 }
1047 return anwcs_wcslib_from_string(fn, slen);
1048 }
1049 anwcs_t* w = anwcs_open(fn, ext);
1050 return w;
1051 }
1052 ~anwcs_t() { free($self); }
1053
1054 double pixel_scale() { return anwcs_pixel_scale($self); }
1055
1056 // FIXME -- this should be more like linearizeAtPoint(x,y)
1057 //void get_cd() { return anwcs_get_cd_matrix($self); }
1058
1059 void get_center(double *p_ra, double *p_dec) {
1060 anwcs_get_radec_center_and_radius($self, p_ra, p_dec, NULL);
1061 }
1062 void get_radius(double *p_radius) {
1063 anwcs_get_radec_center_and_radius($self, NULL, NULL, p_radius);
1064 }
1065
1066 anbool is_inside(double ra, double dec) {
1067 return anwcs_radec_is_inside_image($self, ra, dec);
1068 }
1069 double get_width() {
1070 return anwcs_imagew($self);
1071 }
1072 double get_height() {
1073 return anwcs_imageh($self);
1074 }
1075 void set_width(int W) {
1076 int H = anwcs_imageh($self);
1077 anwcs_set_size($self, W, H);
1078 }
1079 void set_height(int H) {
1080 int W = anwcs_imagew($self);
1081 anwcs_set_size($self, W, H);
1082 }
1083 void pixelxy2radec(double x, double y, double *p_ra, double *p_dec) {
1084 anwcs_pixelxy2radec($self, x, y, p_ra, p_dec);
1085 }
1086 int radec2pixelxy(double ra, double dec, double *p_x, double *p_y) {
1087 return anwcs_radec2pixelxy($self, ra, dec, p_x, p_y);
1088 }
1089
1090 int write_to(const char* filename) {
1091 return anwcs_write($self, filename);
1092 }
1093
1094 }
1095 %pythoncode %{
1096 anwcs = anwcs_t
1097 anwcs.imagew = property(anwcs.get_width, anwcs.set_width, None, 'image width')
1098 anwcs.imageh = property(anwcs.get_height, anwcs.set_height, None, 'image height')
1099 anwcs.writeto = anwcs.write_to
1100
1101 def anwcs_t_get_shape(self):
1102 return int(self.get_height()), int(self.get_width())
1103 anwcs_t.get_shape = anwcs_t_get_shape
1104
1105 def anwcs_t_set_shape(self, S):
1106 H,W = S
1107 self.set_height(H)
1108 self.set_width(W)
1109 anwcs_t.set_shape = anwcs_t_set_shape
1110 anwcs_t.shape = property(anwcs_t.get_shape, anwcs_t.set_shape, None, 'image shape')
1111
1112 # same API as tan_t
1113 anwcs.radec_center = anwcs.get_center
1114 anwcs.radius = anwcs.get_radius
1115
1116 def anwcs_from_string(s):
1117 return anwcs_t(s, -1, len(s))
1118
1119 def anwcs_get_header_string(self):
1120 s = anwcs_wcslib_to_string(self)
1121 return (s +
1122 'NAXIS = 2' + ' '*69 +
1123 'NAXIS1 = % 20i' % self.imagew + ' '*50 +
1124 'NAXIS2 = % 20i' % self.imageh + ' '*50 +
1125 'END'+' '*77)
1126 anwcs.getHeaderString = anwcs_get_header_string
1127
1128 def anwcs_radec_bounds(self, stepsize=1000):
1129 r0,r1,d0,d1 = anwcs_get_radec_bounds(self, stepsize)
1130 return r0,r1,d0,d1
1131 anwcs.radec_bounds = anwcs_radec_bounds
1132
1133 def anwcs_get_cd(self):
1134 return anwcs_get_cd_matrix(self)
1135 anwcs.get_cd = anwcs_get_cd
1136
1137 %}
1138
1139
1140
1141 %include "starutil.h"
1142
1143 %apply (char *STRING, int LENGTH) { (const unsigned char *, int) };
1144
1145 %include "qfits_header.h"
1146 %include "qfits_rw.h"
1147
1148 %pythondynamic qfits_header;
1149
1150 %pythoncode %{
1151 def fitsio_to_qfits_header(hdr):
1152 hdrstr = ''
1153 for rec in hdr.records():
1154 cardstr = rec.get('card', None)
1155 if cardstr is None:
1156 cardstr = rec.get('card_string', None)
1157 if cardstr is None:
1158 cardstr = hdr._record2card(rec)
1159 # pad
1160 cardstr = cardstr + ' '*(80 - len(cardstr))
1161 hdrstr += cardstr
1162 hdrstr += 'END' + ' '*77
1163 qhdr = qfits_header_read_hdr_string(hdrstr)
1164 return qhdr
1165 %}
1166
1167 %include "wcs-pv2sip.h"
1168
1169 %pythoncode %{
1170 def wcs_pv2sip_hdr(hdr, order=5, xlo=0, xhi=0, ylo=0, yhi=0,
1171 stepsize=0, W=0, H=0):
1172 qhdr = fitsio_to_qfits_header(hdr)
1173 forcetan = False
1174 doshift = 1
1175 scamp = False
1176
1177 sip = wcs_pv2sip_header(qhdr, None, 0, stepsize, xlo, xhi, ylo, yhi, W, H,
1178 order, forcetan, doshift)
1179 return sip
1180 %}
1181
1182
1183 %typemap(in) double [ANY] (double temp[$1_dim0]) {
1184 int i;
1185 if (!PySequence_Check($input)) {
1186 PyErr_SetString(PyExc_ValueError,"Expected a sequence");
1187 return NULL;
1188 }
1189 if (PySequence_Length($input) != $1_dim0) {
1190 PyErr_SetString(PyExc_ValueError,"Size mismatch. Expected $1_dim0 elements");
1191 return NULL;
1192 }
1193 for (i = 0; i < $1_dim0; i++) {
1194 PyObject *o = PySequence_GetItem($input,i);
1195 if (PyNumber_Check(o)) {
1196 temp[i] = PyFloat_AsDouble(o);
1197 } else {
1198 PyErr_SetString(PyExc_ValueError,"Sequence elements must be numbers");
1199 return NULL;
1200 }
1201 }
1202 $1 = temp;
1203 }
1204 %typemap(out) double [ANY] {
1205 int i;
1206 $result = PyList_New($1_dim0);
1207 for (i = 0; i < $1_dim0; i++) {
1208 PyObject *o = PyFloat_FromDouble($1[i]);
1209 PyList_SetItem($result,i,o);
1210 }
1211 }
1212
1213 %typemap(in) double flatmatrix[ANY][ANY] (double temp[$1_dim0][$1_dim1]) {
1214 int i;
1215 if (!PySequence_Check($input)) {
1216 PyErr_SetString(PyExc_ValueError,"Expected a sequence");
1217 return NULL;
1218 }
1219 if (PySequence_Length($input) != ($1_dim0 * $1_dim1)) {
1220 PyErr_SetString(PyExc_ValueError,"Size mismatch. Expected $1_dim0*$1_dim1 elements");
1221 return NULL;
1222 }
1223 for (i = 0; i < ($1_dim0*$1_dim1); i++) {
1224 PyObject *o = PySequence_GetItem($input,i);
1225 if (PyNumber_Check(o)) {
1226 // FIXME -- is it dim0 or dim1?
1227 temp[i / $1_dim0][i % $1_dim0] = PyFloat_AsDouble(o);
1228 } else {
1229 PyErr_SetString(PyExc_ValueError,"Sequence elements must be numbers");
1230 return NULL;
1231 }
1232 }
1233 $1 = temp;
1234 }
1235 %typemap(out) double flatmatrix[ANY][ANY] {
1236 int i;
1237 $result = PyList_New($1_dim0 * $1_dim1);
1238 for (i = 0; i < ($1_dim0)*($1_dim1); i++) {
1239 // FIXME -- dim0 or dim1?
1240 PyObject *o = PyFloat_FromDouble($1[i / $1_dim0][i % $1_dim0]);
1241 PyList_SetItem($result,i,o);
1242 }
1243 }
1244
1245
1246 %apply double [ANY] { double crval[2] };
1247 %apply double [ANY] { double crpix[2] };
1248 %apply double flatmatrix[ANY][ANY] { double cd[2][2] };
1249
1250 // SIP coefficients; array size must match SIP_MAXORDER.
1251 %apply double flatmatrix[ANY][ANY] { double a[10][10] };
1252 %apply double flatmatrix[ANY][ANY] { double b[10][10] };
1253 %apply double flatmatrix[ANY][ANY] { double ap[10][10] };
1254 %apply double flatmatrix[ANY][ANY] { double bp[10][10] };
1255
1256
1257 %include "sip.h"
1258 %include "sip_qfits.h"
1259 %include "sip-utils.h"
1260
1261 %pythondynamic sip_t;
1262
1263 %extend sip_t {
1264 sip_t(const char* fn=NULL, int ext=0) {
1265 if (fn)
1266 return sip_read_header_file_ext_only(fn, ext, NULL);
1267 sip_t* t = (sip_t*)calloc(1, sizeof(sip_t));
1268 return t;
1269 }
1270
1271 // from string -- third arg is just to distinguish this signature.
1272 sip_t(const char* s, int len, int XXX) {
1273 return sip_from_string(s, len, NULL);
1274 }
1275
1276 // copy constructor
1277 sip_t(const sip_t* other) {
1278 sip_t* t = (sip_t*)calloc(1, sizeof(sip_t));
1279 memcpy(t, other, sizeof(sip_t));
1280 return t;
1281 }
1282
1283 sip_t(const tan_t* other) {
1284 sip_t* t = (sip_t*)calloc(1, sizeof(sip_t));
1285 memcpy(&(t->wcstan), other, sizeof(tan_t));
1286 return t;
1287 }
1288
1289 sip_t(const qfits_header* hdr) {
1290 sip_t* t = sip_read_header(hdr, NULL);
1291 return t;
1292 }
1293
1294 ~sip_t() { free($self); }
1295
1296 sip_t* get_subimage(int x0, int y0, int w, int h) {
1297 sip_t* sub = malloc(sizeof(sip_t));
1298 memcpy(sub, $self, sizeof(sip_t));
1299 sub->wcstan.crpix[0] -= x0;
1300 sub->wcstan.crpix[1] -= y0;
1301 sub->wcstan.imagew = w;
1302 sub->wcstan.imageh = h;
1303 return sub;
1304 }
1305
1306 sip_t* scale(double factor) {
1307 sip_t* s = (sip_t*)calloc(1, sizeof(sip_t));
1308 sip_scale($self, s, factor);
1309 return s;
1310 }
1311
1312 double pixel_scale() { return sip_pixel_scale($self); }
1313
1314 void radec_center(double *p_ra, double *p_dec) {
1315 sip_get_radec_center($self, p_ra, p_dec);
1316 }
1317 double radius() {
1318 return sip_get_radius_deg($self);
1319 }
1320
1321 int write_to(const char* filename) {
1322 return sip_write_to_file($self, filename);
1323 }
1324
1325 int ensure_inverse_polynomials() {
1326 return sip_ensure_inverse_polynomials($self);
1327 }
1328
1329 /*
1330 double* get_cd_matrix() {
1331 return $self->wcstan.cd;
1332 }
1333 */
1334
1335 void pixelxy2xyz(double x, double y, double *p_x, double *p_y, double *p_z) {
1336 double xyz[3];
1337 sip_pixelxy2xyzarr($self, x, y, xyz);
1338 *p_x = xyz[0];
1339 *p_y = xyz[1];
1340 *p_z = xyz[2];
1341 }
1342 void pixelxy2radec(double x, double y, double *p_ra, double *p_dec) {
1343 sip_pixelxy2radec($self, x, y, p_ra, p_dec);
1344 }
1345 int radec2pixelxy(double ra, double dec, double *p_x, double *p_y) {
1346 return sip_radec2pixelxy($self, ra, dec, p_x, p_y);
1347 }
1348 void iwc2pixelxy(double u, double v, double *p_x, double *p_y) {
1349 sip_iwc2pixelxy($self, u, v, p_x, p_y);
1350 }
1351 void pixelxy2iwc(double x, double y, double *p_x, double *p_y) {
1352 sip_pixelxy2iwc($self, x, y, p_x, p_y);
1353 }
1354 void iwc2radec(double u, double v, double *p_ra, double *p_dec) {
1355 sip_iwc2radec($self, u, v, p_ra, p_dec);
1356 }
1357 int radec2iwc(double ra, double dec, double *p_x, double *p_y) {
1358 return sip_radec2iwc($self, ra, dec, p_x, p_y);
1359 }
1360 int xyz2pixelxy(double x, double y, double z, double *p_x, double *p_y) {
1361 double xyz[3];
1362 xyz[0] = x;
1363 xyz[1] = y;
1364 xyz[2] = z;
1365 return sip_xyzarr2pixelxy($self, xyz, p_x, p_y);
1366 }
1367
1368 anbool is_inside(double ra, double dec) {
1369 return sip_is_inside_image($self, ra, dec);
1370 }
1371
1372 void set_a_term(int i, int j, double val) {
1373 checkorder(i, j);
1374 $self->a[i][j] = val;
1375 }
1376 void set_b_term(int i, int j, double val) {
1377 checkorder(i, j);
1378 $self->b[i][j] = val;
1379 }
1380 void set_ap_term(int i, int j, double val) {
1381 checkorder(i, j);
1382 $self->ap[i][j] = val;
1383 }
1384 void set_bp_term(int i, int j, double val) {
1385 checkorder(i, j);
1386 $self->bp[i][j] = val;
1387 }
1388
1389 double get_a_term(int i, int j) {
1390 checkorder(i, j);
1391 return $self->a[i][j];
1392 }
1393 double get_b_term(int i, int j) {
1394 checkorder(i, j);
1395 return $self->b[i][j];
1396 }
1397 double get_ap_term(int i, int j) {
1398 checkorder(i, j);
1399 return $self->ap[i][j];
1400 }
1401 double get_bp_term(int i, int j) {
1402 checkorder(i, j);
1403 return $self->bp[i][j];
1404 }
1405
1406 void set_width(double x) {
1407 $self->wcstan.imagew = x;
1408 }
1409 void set_height(double x) {
1410 $self->wcstan.imageh = x;
1411 }
1412
1413 double get_width() {
1414 return $self->wcstan.imagew;
1415 }
1416 double get_height() {
1417 return $self->wcstan.imageh;
1418 }
1419 void get_distortion(double x, double y, double *p_x, double *p_y) {
1420 return sip_pixel_distortion($self, x, y, p_x, p_y);
1421 }
1422 void get_undistortion(double x, double y, double *p_x, double *p_y) {
1423 return sip_pixel_undistortion($self, x, y, p_x, p_y);
1424 }
1425
1426 int write_to(const char* filename) {
1427 return sip_write_to_file($self, filename);
1428 }
1429
1430
1431 }
1432 %pythoncode %{
1433
1434 def sip_t_tostring(self):
1435 tan = self.wcstan
1436 ct = 'SIN' if tan.sin else 'TAN'
1437 return (('SIP(%s): crpix (%.1f, %.1f), crval (%g, %g), cd (%g, %g, %g, %g), '
1438 + 'image %g x %g; SIP orders A=%i, B=%i, AP=%i, BP=%i') %
1439 (ct, tan.crpix[0], tan.crpix[1], tan.crval[0], tan.crval[1],
1440 tan.cd[0], tan.cd[1], tan.cd[2], tan.cd[3],
1441 tan.imagew, tan.imageh, self.a_order, self.b_order,
1442 self.ap_order, self.bp_order))
1443 sip_t.__str__ = sip_t_tostring
1444
1445 def sip_t_addtoheader(self, hdr):
1446 '''Adds this SIP WCS header to the given fitsio header'''
1447 self.wcstan.add_to_header(hdr)
1448 hdr.delete('CTYPE1')
1449 hdr.delete('CTYPE2')
1450 for k,v,c in [
1451 ('CTYPE1', 'RA---TAN-SIP', 'TANgent plane+SIP'),
1452 ('CTYPE2', 'DEC--TAN-SIP', 'TANgent plane+SIP'),
1453 ('A_ORDER', self.a_order, 'Polynomial order, axis 1'),
1454 ('B_ORDER', self.b_order, 'Polynomial order, axis 2'),
1455 ('AP_ORDER', self.ap_order, 'Inv.polynomial order, axis 1'),
1456 ('BP_ORDER', self.bp_order, 'Inv.polynomial order, axis 2'),
1457 ]:
1458 hdr.add_record(dict(name=k, value=v, comment=c))
1459 for i in range(self.a_order + 1):
1460 for j in range(self.a_order + 1):
1461 #if i + j < 1:
1462 # drop linear (CD) terms
1463 if i + j < 2:
1464 continue
1465 if i + j > self.a_order:
1466 continue
1467 hdr.add_record(dict(name='A_%i_%i' % (i,j), value=self.get_a_term(i, j),
1468 comment='SIP polynomial term'))
1469 for i in range(self.b_order + 1):
1470 for j in range(self.b_order + 1):
1471 #if i + j < 1:
1472 # drop linear (CD) terms
1473 if i + j < 2:
1474 continue
1475 if i + j > self.b_order:
1476 continue
1477 hdr.add_record(dict(name='B_%i_%i' % (i,j), value=self.get_b_term(i, j),
1478 comment='SIP polynomial term'))
1479 for i in range(self.ap_order + 1):
1480 for j in range(self.ap_order + 1):
1481 if i + j < 1:
1482 continue
1483 if i + j > self.ap_order:
1484 continue
1485 hdr.add_record(dict(name='AP_%i_%i' % (i,j), value=self.get_ap_term(i, j),
1486 comment='SIP polynomial term'))
1487 for i in range(self.bp_order + 1):
1488 for j in range(self.bp_order + 1):
1489 if i + j < 1:
1490 continue
1491 if i + j > self.bp_order:
1492 continue
1493 hdr.add_record(dict(name='BP_%i_%i' % (i,j), value=self.get_bp_term(i, j),
1494 comment='SIP polynomial term'))
1495 sip_t.add_to_header = sip_t_addtoheader
1496
1497
1498 # def sip_t_get_subimage(self, x0, y0, w, h):
1499 # wcs2 = sip_t(self)
1500 # cpx,cpy = wcs2.crpix
1501 # wcs2.set_crpix((cpx - x0, cpy - y0))
1502 # wcs2.set_width(float(w))
1503 # wcs2.set_height(float(h))
1504 # return wcs2
1505 # sip_t.get_subimage = sip_t_get_subimage
1506
1507 def sip_t_get_shape(self):
1508 return (self.wcstan.imageh, self.wcstan.imagew)
1509 sip_t.get_shape = sip_t_get_shape
1510
1511 def sip_t_set_shape(self, S):
1512 H,W = S
1513 self.set_height(H)
1514 self.set_width(W)
1515 sip_t.set_shape = sip_t_set_shape
1516
1517 sip_t.imagew = property(sip_t.get_width, sip_t.set_width, None, 'image width')
1518 sip_t.imageh = property(sip_t.get_height, sip_t.set_height, None, 'image height')
1519 sip_t.shape = property(sip_t.get_shape, sip_t.set_shape, None, 'image shape')
1520
1521 def sip_t_get_cd(self):
1522 cd = self.wcstan.cd
1523 return (cd[0], cd[1], cd[2], cd[3])
1524 def sip_t_set_cd(self, x):
1525 self.wcstan.cd = x
1526 sip_t.get_cd = sip_t_get_cd
1527 sip_t.set_cd = sip_t_set_cd
1528
1529 def sip_t_get_crval(self):
1530 return self.wcstan.crval
1531 def sip_t_set_crval(self, x):
1532 self.wcstan.crval = x
1533 sip_t.get_crval = sip_t_get_crval
1534 sip_t.set_crval = sip_t_set_crval
1535
1536 def sip_t_get_crpix(self):
1537 return self.wcstan.crpix
1538 def sip_t_set_crpix(self, x):
1539 self.wcstan.crpix = x
1540 sip_t.get_crpix = sip_t_get_crpix
1541 sip_t.set_crpix = sip_t_set_crpix
1542
1543 sip_t.crval = property(sip_t_get_crval, sip_t_set_crval, None, 'CRVAL')
1544 sip_t.crpix = property(sip_t_get_crpix, sip_t_set_crpix, None, 'CRPIX')
1545 sip_t.cd = property(sip_t_get_cd , sip_t_set_cd, None, 'CD')
1546
1547
1548 def sip_t_radec_bounds(self):
1549 # W,H = self.wcstan.imagew, self.wcstan.imageh
1550 # r,d = self.pixelxy2radec([1, W, W, 1], [1, 1, H, H])
1551 # return (r.min(), r.max(), d.min(), d.max())
1552 W,H = self.imagew, self.imageh
1553 r,d = self.pixelxy2radec([1, W/2, W, W, W, W/2, 1, 1], [1, 1, 1, H/2, H, H, H, H/2])
1554 rx = r.max()
1555 rn = r.min()
1556 # ugh, RA wrap-around. We find the largest value < 180 (ie, near zero) and smallest value > 180 (ie, near 360)
1557 # and report them with ralo > rahi so that this case can be identified
1558 if rx - rn > 180:
1559 rx = r[r < 180].max()
1560 rn = r[r > 180].min()
1561 return (rn, rx, d.min(), d.max())
1562 sip_t.radec_bounds = sip_t_radec_bounds
1563
1564 #def sip_t_fromstring(s):
1565 # sip = sip_from_string(s, len(s),
1566
1567 _real_sip_t_init = sip_t.__init__
1568 def my_sip_t_init(self, *args, **kwargs):
1569 # fitsio header: check for '.records()' function.
1570 if len(args) == 1 and hasattr(args[0], 'records'):
1571 try:
1572 hdr = args[0]
1573 qhdr = fitsio_to_qfits_header(hdr)
1574 args = [qhdr]
1575 except:
1576 pass
1577
1578 _real_sip_t_init(self, *args, **kwargs)
1579 if self.this is None:
1580 raise RuntimeError('Duck punch!')
1581 sip_t.__init__ = my_sip_t_init
1582
1583
1584 Sip = sip_t
1585 %}
1586
1587 %pythondynamic tan_t;
1588
1589 %extend tan_t {
1590 tan_t(char* fn=NULL, int ext=0, int only=0) {
1591 tan_t* t = NULL;
1592 if (fn) {
1593 if (only) {
1594 t = tan_read_header_file_ext_only(fn, ext, NULL);
1595 } else {
1596 t = tan_read_header_file_ext(fn, ext, NULL);
1597 }
1598 } else {
1599 t = (tan_t*)calloc(1, sizeof(tan_t));
1600 }
1601 // printf("tan_t: %p\n", t);
1602 if (!t) {
1603 // SWIG_exception(SWIG_RuntimeError, "Failed to read TAN WCS header");
1604 PyErr_SetString(PyExc_RuntimeError, "Failed to read TAN WCS header");
1605 return NULL;
1606 }
1607 return t;
1608 }
1609
1610 tan_t(double crval1, double crval2, double crpix1, double crpix2,
1611 double cd11, double cd12, double cd21, double cd22,
1612 double imagew, double imageh) {
1613 tan_t* t = (tan_t*)calloc(1, sizeof(tan_t));
1614 t->crval[0] = crval1;
1615 t->crval[1] = crval2;
1616 t->crpix[0] = crpix1;
1617 t->crpix[1] = crpix2;
1618 t->cd[0][0] = cd11;
1619 t->cd[0][1] = cd12;
1620 t->cd[1][0] = cd21;
1621 t->cd[1][1] = cd22;
1622 t->imagew = imagew;
1623 t->imageh = imageh;
1624 return t;
1625 }
1626 tan_t(const tan_t* other) {
1627 tan_t* t = (tan_t*)calloc(1, sizeof(tan_t));
1628 memcpy(t, other, sizeof(tan_t));
1629 return t;
1630 }
1631
1632 tan_t(const qfits_header* hdr) {
1633 tan_t* t = tan_read_header(hdr, NULL);
1634 return t;
1635 }
1636
1637 ~tan_t() { free($self); }
1638 void set(double crval1, double crval2,
1639 double crpix1, double crpix2,
1640 double cd11, double cd12, double cd21, double cd22,
1641 double imagew, double imageh) {
1642 $self->crval[0] = crval1;
1643 $self->crval[1] = crval2;
1644 $self->crpix[0] = crpix1;
1645 $self->crpix[1] = crpix2;
1646 $self->cd[0][0] = cd11;
1647 $self->cd[0][1] = cd12;
1648 $self->cd[1][0] = cd21;
1649 $self->cd[1][1] = cd22;
1650 $self->imagew = imagew;
1651 $self->imageh = imageh;
1652 }
1653
1654 anbool is_inside(double ra, double dec) {
1655 return tan_is_inside_image($self, ra, dec);
1656 }
1657
1658 tan_t* scale(double factor) {
1659 tan_t* t = (tan_t*)calloc(1, sizeof(tan_t));
1660 tan_scale($self, t, factor);
1661 return t;
1662 }
1663 tan_t* rotate(double angle_deg) {
1664 tan_t* t = (tan_t*)calloc(1, sizeof(tan_t));
1665 tan_rotate($self, t, angle_deg);
1666 return t;
1667 }
1668 double get_width() {
1669 return $self->imagew;
1670 }
1671 double get_height() {
1672 return $self->imageh;
1673 }
1674
1675 void set_width(double x) {
1676 $self->imagew = x;
1677 }
1678 void set_height(double x) {
1679 $self->imageh = x;
1680 }
1681
1682 double pixel_scale() { return tan_pixel_scale($self); }
1683 void radec_center(double *p_ra, double *p_dec) {
1684 tan_get_radec_center($self, p_ra, p_dec);
1685 }
1686 double radius() {
1687 return tan_get_radius_deg($self);
1688 }
1689 void xyzcenter(double *p_x, double *p_y, double *p_z) {
1690 double xyz[3];
1691 tan_pixelxy2xyzarr($self, 0.5+$self->imagew/2.0, 0.5+$self->imageh/2.0, xyz);
1692 *p_x = xyz[0];
1693 *p_y = xyz[1];
1694 *p_z = xyz[2];
1695 }
1696 void pixelxy2xyz(double x, double y, double *p_x, double *p_y, double *p_z) {
1697 double xyz[3];
1698 tan_pixelxy2xyzarr($self, x, y, xyz);
1699 *p_x = xyz[0];
1700 *p_y = xyz[1];
1701 *p_z = xyz[2];
1702 }
1703 void pixelxy2radec(double x, double y, double *p_ra, double *p_dec) {
1704 tan_pixelxy2radec($self, x, y, p_ra, p_dec);
1705 }
1706 int radec2pixelxy(double ra, double dec, double *p_x, double *p_y) {
1707 return tan_radec2pixelxy($self, ra, dec, p_x, p_y);
1708 }
1709 void iwc2pixelxy(double u, double v, double *p_x, double *p_y) {
1710 tan_iwc2pixelxy($self, u, v, p_x, p_y);
1711 }
1712 void pixelxy2iwc(double x, double y, double *p_x, double *p_y) {
1713 tan_pixelxy2iwc($self, x, y, p_x, p_y);
1714 }
1715 void iwc2radec(double u, double v, double *p_ra, double *p_dec) {
1716 tan_iwc2radec($self, u, v, p_ra, p_dec);
1717 }
1718 int radec2iwc(double ra, double dec, double *p_x, double *p_y) {
1719 return tan_radec2iwc($self, ra, dec, p_x, p_y);
1720 }
1721 int xyz2pixelxy(double x, double y, double z, double *p_x, double *p_y) {
1722 double xyz[3];
1723 xyz[0] = x;
1724 xyz[1] = y;
1725 xyz[2] = z;
1726 return tan_xyzarr2pixelxy($self, xyz, p_x, p_y);
1727 }
1728 int write_to(const char* filename) {
1729 return tan_write_to_file($self, filename);
1730 }
1731 void set_crval(double ra, double dec) {
1732 $self->crval[0] = ra;
1733 $self->crval[1] = dec;
1734 }
1735 void set_crpix(double x, double y) {
1736 $self->crpix[0] = x;
1737 $self->crpix[1] = y;
1738 }
1739 void set_cd(double cd11, double cd12, double cd21, double cd22) {
1740 $self->cd[0][0] = cd11;
1741 $self->cd[0][1] = cd12;
1742 $self->cd[1][0] = cd21;
1743 $self->cd[1][1] = cd22;
1744 }
1745 void set_imagesize(double w, double h) {
1746 $self->imagew = w;
1747 $self->imageh = h;
1748 }
1749
1750 /*
1751 double* get_cd_matrix() {
1752 return $self->cd;
1753 }
1754 */
1755
1756
1757 };
1758
1759
1760 %inline %{
1761
1762 // Wrapper on coadd_add_image that accepts numpy arrays.
1763
1764 static int coadd_add_numpy(coadd_t* c,
1765 PyObject* py_img, PyObject* py_weight,
1766 float fweight, const anwcs_t* wcs) {
1767 PyArray_Descr* dtype = PyArray_DescrFromType(NPY_FLOAT);
1768 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED | NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ELEMENTSTRIDES;
1769 float *img, *weight=NULL;
1770
1771 PyArrayObject *np_img=NULL, *np_weight=NULL;
1772
1773 Py_INCREF(dtype);
1774 np_img = (PyArrayObject*)PyArray_CheckFromAny(py_img, dtype, 2, 2, req, NULL);
1775 img = PyArray_DATA(np_img);
1776 if (!np_img) {
1777 ERR("Failed to PyArray_FromAny the image\n");
1778 Py_XDECREF(np_img);
1779 Py_DECREF(dtype);
1780 return -1;
1781 }
1782 if (py_weight != Py_None) {
1783 Py_INCREF(dtype);
1784 np_weight = (PyArrayObject*)PyArray_CheckFromAny(py_weight, dtype, 2, 2, req, NULL);
1785 if (!np_weight) {
1786 ERR("Failed to PyArray_FromAny the weight\n");
1787 Py_XDECREF(np_weight);
1788 Py_DECREF(dtype);
1789 return -1;
1790 }
1791 weight = PyArray_DATA(np_weight);
1792 }
1793
1794 int rtn = coadd_add_image(c, img, weight, fweight, wcs);
1795
1796 Py_DECREF(np_img);
1797 if (weight) {
1798 Py_DECREF(np_weight);
1799 }
1800 Py_DECREF(dtype);
1801 return rtn;
1802 }
1803
1804 static PyObject* coadd_get_snapshot_numpy(coadd_t* co, float badpix) {
1805 npy_intp dim[2];
1806 PyObject* npimg;
1807 dim[0] = co->H;
1808 dim[1] = co->W;
1809 npimg = PyArray_EMPTY(2, dim, NPY_FLOAT, 0);
1810 coadd_get_snapshot(co, PyArray_DATA((PyArrayObject*)npimg), badpix);
1811 return npimg;
1812 }
1813
1814 static sip_t* fit_sip_wcs_py(PyObject* py_starxyz,
1815 PyObject* py_fieldxy,
1816 PyObject* py_weights,
1817 tan_t* tanin,
1818 int sip_order,
1819 int inv_order) {
1820 PyArray_Descr* dtype = PyArray_DescrFromType(NPY_DOUBLE);
1821 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED | NPY_ARRAY_NOTSWAPPED |
1822 NPY_ARRAY_ELEMENTSTRIDES;
1823 PyArrayObject *np_starxyz=NULL, *np_fieldxy=NULL, *np_weights=NULL;
1824
1825 Py_INCREF(dtype);
1826 np_starxyz = (PyArrayObject*)PyArray_CheckFromAny(py_starxyz, dtype, 2, 2, req, NULL);
1827 Py_INCREF(dtype);
1828 np_fieldxy = (PyArrayObject*)PyArray_CheckFromAny(py_fieldxy, dtype, 2, 2, req, NULL);
1829 if (!np_starxyz || !np_fieldxy) {
1830 Py_DECREF(dtype);
1831 Py_DECREF(dtype);
1832 printf("Failed to convert starxyz or fieldxy to numpy double arrays\n");
1833 return NULL;
1834 }
1835 if (py_weights != Py_None) {
1836 Py_INCREF(dtype);
1837 np_weights = (PyArrayObject*)PyArray_CheckFromAny(py_weights, dtype, 1, 1, req,NULL);
1838 if (!np_weights) {
1839 Py_DECREF(dtype);
1840 printf("Failed to convert weights to numpy double array\n");
1841 return NULL;
1842 }
1843 }
1844 Py_DECREF(dtype);
1845 dtype = NULL;
1846
1847 int M = (int)PyArray_DIM(np_starxyz, 0);
1848 if (PyArray_DIM(np_fieldxy, 0) != M) {
1849 printf("Expected starxyz and fieldxy to have the same length\n");
1850 return NULL;
1851 }
1852 if (np_weights && (PyArray_DIM(np_weights, 0) != M)) {
1853 printf("Expected starxyz and weights to have the same length\n");
1854 return NULL;
1855 }
1856 if ((PyArray_DIM(np_starxyz, 1) != 3) ||
1857 (PyArray_DIM(np_fieldxy, 1) != 2)) {
1858 printf("Expected starxyz Mx3 and fieldxy Mx2\n");
1859 return NULL;
1860 }
1861
1862 sip_t* sipout = calloc(1, sizeof(sip_t));
1863
1864 double* weights = NULL;
1865 if (np_weights)
1866 weights = PyArray_DATA(np_weights);
1867
1868 int doshift = 1;
1869 int rtn = fit_sip_wcs(PyArray_DATA(np_starxyz),
1870 PyArray_DATA(np_fieldxy),
1871 weights, M, tanin, sip_order, inv_order,
1872 doshift, sipout);
1873 if (rtn) {
1874 free(sipout);
1875 printf("fit_sip_wcs() returned %i\n", rtn);
1876 return NULL;
1877 }
1878 return sipout;
1879 }
1880
1881
1882 %}
1883
1884
1885 %inline %{
1886
1887 typedef anbool (*f_2to2ok)(const void*, double, double, double*, double*);
1888 typedef void (*f_2to2) (const void*, double, double, double*, double*);
1889 typedef int (*f_2to2i) (const void*, double, double, double*, double*);
1890
1891 static PyObject* broadcast_2to2ok
1892 (
1893 //anbool func(const void*, double, double, double*, double*),
1894 f_2to2ok func,
1895 const void* baton,
1896 PyObject* in1, PyObject* in2);
1897
1898 static PyObject* broadcast_2to2
1899 (
1900 //void func(const void*, double, double, double*, double*),
1901 f_2to2 func,
1902 const void* baton,
1903 PyObject* in1, PyObject* in2);
1904
1905 static PyObject* broadcast_2to2i
1906 (
1907 //int func(const void*, double, double, double*, double*),
1908 f_2to2i func,
1909 const void* baton,
1910 PyObject* in1, PyObject* in2);
1911
1912 static PyObject* tan_rd2xy_wrapper(const tan_t* wcs,
1913 PyObject* in1, PyObject* in2) {
1914 return broadcast_2to2ok((f_2to2ok)tan_radec2pixelxy, wcs, in1, in2);
1915 }
1916 static PyObject* sip_rd2xy_wrapper(const sip_t* wcs,
1917 PyObject* in1, PyObject* in2) {
1918 return broadcast_2to2ok((f_2to2ok)sip_radec2pixelxy, wcs, in1, in2);
1919 }
1920 static PyObject* anwcs_rd2xy_wrapper(const anwcs_t* wcs,
1921 PyObject* in1, PyObject* in2) {
1922 return broadcast_2to2i((f_2to2i)anwcs_radec2pixelxy, wcs, in1, in2);
1923 }
1924
1925 static PyObject* tan_iwc2xy_wrapper(const tan_t* wcs,
1926 PyObject* in1, PyObject* in2) {
1927 return broadcast_2to2((f_2to2)tan_iwc2pixelxy, wcs, in1, in2);
1928 }
1929 static PyObject* sip_iwc2xy_wrapper(const sip_t* wcs,
1930 PyObject* in1, PyObject* in2) {
1931 return broadcast_2to2((f_2to2)sip_iwc2pixelxy, wcs, in1, in2);
1932 }
1933
1934 static PyObject* tan_xy2iwc_wrapper(const tan_t* wcs,
1935 PyObject* in1, PyObject* in2) {
1936 return broadcast_2to2((f_2to2)tan_pixelxy2iwc, wcs, in1, in2);
1937 }
1938 static PyObject* sip_xy2iwc_wrapper(const sip_t* wcs,
1939 PyObject* in1, PyObject* in2) {
1940 return broadcast_2to2((f_2to2)sip_pixelxy2iwc, wcs, in1, in2);
1941 }
1942
1943 static PyObject* tan_iwc2rd_wrapper(const tan_t* wcs,
1944 PyObject* in1, PyObject* in2) {
1945 return broadcast_2to2((f_2to2)tan_iwc2radec, wcs, in1, in2);
1946 }
1947 static PyObject* sip_iwc2rd_wrapper(const sip_t* wcs,
1948 PyObject* in1, PyObject* in2) {
1949 return broadcast_2to2((f_2to2)sip_iwc2radec, wcs, in1, in2);
1950 }
1951
1952 static PyObject* tan_rd2iwc_wrapper(const tan_t* wcs,
1953 PyObject* in1, PyObject* in2) {
1954 return broadcast_2to2ok((f_2to2ok)tan_radec2iwc, wcs, in1, in2);
1955 }
1956 static PyObject* sip_rd2iwc_wrapper(const sip_t* wcs,
1957 PyObject* in1, PyObject* in2) {
1958 return broadcast_2to2ok((f_2to2ok)sip_radec2iwc, wcs, in1, in2);
1959 }
1960
1961 static PyObject* tan_xy2rd_wrapper(const tan_t* wcs,
1962 PyObject* in1, PyObject* in2) {
1963 return broadcast_2to2((f_2to2)tan_pixelxy2radec, wcs, in1, in2);
1964 }
1965 static PyObject* sip_xy2rd_wrapper(const sip_t* wcs,
1966 PyObject* in1, PyObject* in2) {
1967 return broadcast_2to2((f_2to2)sip_pixelxy2radec, wcs, in1, in2);
1968 }
1969 static PyObject* anwcs_xy2rd_wrapper(const anwcs_t* wcs,
1970 PyObject* in1, PyObject* in2) {
1971 return broadcast_2to2i((f_2to2i)anwcs_pixelxy2radec, wcs, in1, in2);
1972 }
1973
1974 static PyObject* broadcast_2to2ok
1975 (
1976 //anbool func(const void*, double, double, double*, double*),
1977 f_2to2ok func,
1978 const void* baton,
1979 PyObject* in1, PyObject* in2) {
1980
1981 NpyIter *iter = NULL;
1982 NpyIter_IterNextFunc *iternext;
1983 PyArrayObject *op[5];
1984 PyObject *ret;
1985 npy_uint32 flags;
1986 npy_uint32 op_flags[5];
1987 npy_intp *innersizeptr;
1988 char **dataptrarray;
1989 npy_intp* strideptr;
1990 PyArray_Descr* dtypes[5];
1991 npy_intp i, N;
1992
1993 // we'll do the inner loop ourselves
1994 flags = NPY_ITER_EXTERNAL_LOOP;
1995 // use buffers to satisfy dtype casts
1996 flags |= NPY_ITER_BUFFERED;
1997 // grow inner loop
1998 flags |= NPY_ITER_GROWINNER;
1999
2000 op[0] = (PyArrayObject*)PyArray_FromAny(in1, NULL, 0, 0, 0, NULL);
2001 op[1] = (PyArrayObject*)PyArray_FromAny(in2, NULL, 0, 0, 0, NULL);
2002 // automatically allocate the output arrays.
2003 op[2] = NULL;
2004 op[3] = NULL;
2005 op[4] = NULL;
2006
2007 if ((PyArray_Size((PyObject*)op[0]) == 0) ||
2008 (PyArray_Size((PyObject*)op[1]) == 0)) {
2009 // empty inputs -- empty outputs
2010 npy_intp dim = 0;
2011 ret = Py_BuildValue("(NNN)",
2012 PyArray_SimpleNew(1, &dim, NPY_BOOL),
2013 PyArray_SimpleNew(1, &dim, NPY_DOUBLE),
2014 PyArray_SimpleNew(1, &dim, NPY_DOUBLE));
2015 goto cleanup;
2016 }
2017
2018 op_flags[0] = NPY_ITER_READONLY | NPY_ITER_NBO;
2019 op_flags[1] = NPY_ITER_READONLY | NPY_ITER_NBO;
2020 op_flags[2] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO;
2021 op_flags[3] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO;
2022 op_flags[4] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO;
2023
2024 dtypes[0] = PyArray_DescrFromType(NPY_DOUBLE);
2025 dtypes[1] = PyArray_DescrFromType(NPY_DOUBLE);
2026 dtypes[2] = PyArray_DescrFromType(NPY_DOUBLE);
2027 dtypes[3] = PyArray_DescrFromType(NPY_DOUBLE);
2028 dtypes[4] = PyArray_DescrFromType(NPY_BOOL);
2029
2030 iter = NpyIter_MultiNew(5, op, flags, NPY_KEEPORDER, NPY_SAFE_CASTING,
2031 op_flags, dtypes);
2032 for (i=0; i<5; i++)
2033 Py_DECREF(dtypes[i]);
2034
2035 if (!iter)
2036 return NULL;
2037
2038 iternext = NpyIter_GetIterNext(iter, NULL);
2039 strideptr = NpyIter_GetInnerStrideArray(iter);
2040 // The inner loop size and data pointers may change during the
2041 // loop, so just cache the addresses.
2042 innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);
2043 dataptrarray = NpyIter_GetDataPtrArray(iter);
2044
2045 do {
2046 // are the inputs contiguous? (Outputs will be, since we
2047 // allocated them)
2048 if ((strideptr[0] == sizeof(double)) &&
2049 (strideptr[1] == sizeof(double))) {
2050 // printf("Contiguous inputs; going fast\n");
2051 double* din1 = (double*)(dataptrarray[0]);
2052 double* din2 = (double*)(dataptrarray[1]);
2053 double* dout1 = (double*)(dataptrarray[2]);
2054 double* dout2 = (double*)(dataptrarray[3]);
2055 char* ok = dataptrarray[4];
2056 N = *innersizeptr;
2057 while (N--) {
2058 *ok = func(baton, *din1, *din2, dout1, dout2);
2059 ok++;
2060 din1++;
2061 din2++;
2062 dout1++;
2063 dout2++;
2064 }
2065 } else {
2066 // printf("Non-contiguous inputs; going slow\n");
2067 npy_intp stride1 = strideptr[0];
2068 npy_intp stride2 = strideptr[1];
2069 npy_intp size = *innersizeptr;
2070 char* src1 = dataptrarray[0];
2071 char* src2 = dataptrarray[1];
2072 double* dout1 = (double*)dataptrarray[2];
2073 double* dout2 = (double*)dataptrarray[3];
2074 char* ok = dataptrarray[4];
2075
2076 for (i=0; i<size; i++) {
2077 *ok = func(baton, *((double*)src1), *((double*)src2),
2078 dout1, dout2);
2079 ok++;
2080 src1 += stride1;
2081 src2 += stride2;
2082 dout1++;
2083 dout2++;
2084 }
2085 }
2086 } while (iternext(iter));
2087
2088 if (PyArray_IsPythonScalar(in1) && PyArray_IsPythonScalar(in2)) {
2089 PyObject* px = (PyObject*)(NpyIter_GetOperandArray(iter)[2]);
2090 PyObject* py = (PyObject*)(NpyIter_GetOperandArray(iter)[3]);
2091 PyObject* pok = (PyObject*)(NpyIter_GetOperandArray(iter)[4]);
2092 //printf("Both inputs are python scalars\n");
2093 double d;
2094 unsigned char c;
2095 d = *(double*)PyArray_DATA((PyArrayObject*)px);
2096 px = PyFloat_FromDouble(d);
2097 d = *(double*)PyArray_DATA((PyArrayObject*)py);
2098 py = PyFloat_FromDouble(d);
2099 c = *(unsigned char*)PyArray_DATA((PyArrayObject*)pok);
2100 pok = PyBool_FromLong(c);
2101 ret = Py_BuildValue("(NNN)", pok, px, py);
2102 /*
2103 // I couldn't figure this out -- ScalarAsCtype didn't work
2104 if (PyArray_CheckScalar(px)) {
2105 printf("x is scalar\n");
2106 }
2107 if (PyArray_IsScalar(px, Double)) {
2108 printf("x is PyDoubleArrType\n");
2109 }
2110 if (PyArray_IsScalar(px, CDouble)) {
2111 printf("x is PyCDoubleArrType\n");
2112 }
2113 if (PyArray_ISFLOAT(px)) {
2114 printf("x ISFLOAT\n");
2115 }
2116 //PyArray_ScalarAsCtype(px, &d);
2117 */
2118 } else {
2119 ret = Py_BuildValue("(OOO)",
2120 NpyIter_GetOperandArray(iter)[4],
2121 NpyIter_GetOperandArray(iter)[2],
2122 NpyIter_GetOperandArray(iter)[3]);
2123 }
2124
2125 cleanup:
2126 if (NpyIter_Deallocate(iter) != NPY_SUCCEED) {
2127 Py_DECREF(ret);
2128 return NULL;
2129 }
2130 Py_DECREF(op[0]);
2131 Py_DECREF(op[1]);
2132 return ret;
2133 }
2134
2135
2136 static PyObject* broadcast_2to2i
2137 (
2138 //int func(const void*, double, double, double*, double*),
2139 f_2to2i func,
2140 const void* baton,
2141 PyObject* in1, PyObject* in2) {
2142
2143 NpyIter *iter = NULL;
2144 NpyIter_IterNextFunc *iternext;
2145 PyArrayObject *op[5];
2146 PyObject *ret;
2147 npy_uint32 flags;
2148 npy_uint32 op_flags[5];
2149 npy_intp *innersizeptr;
2150 char **dataptrarray;
2151 npy_intp* strideptr;
2152 PyArray_Descr* dtypes[5];
2153 int j;
2154
2155 // we'll do the inner loop ourselves
2156 flags = NPY_ITER_EXTERNAL_LOOP;
2157 // use buffers to satisfy dtype casts
2158 flags |= NPY_ITER_BUFFERED;
2159 // grow inner loop
2160 flags |= NPY_ITER_GROWINNER;
2161
2162 op[0] = (PyArrayObject*)PyArray_FromAny(in1, NULL, 0, 0, 0, NULL);
2163 op[1] = (PyArrayObject*)PyArray_FromAny(in2, NULL, 0, 0, 0, NULL);
2164 // automatically allocate the output arrays.
2165 op[2] = NULL;
2166 op[3] = NULL;
2167 op[4] = NULL;
2168
2169 if ((PyArray_Size((PyObject*)op[0]) == 0) ||
2170 (PyArray_Size((PyObject*)op[1]) == 0)) {
2171 // empty inputs -- empty outputs
2172 npy_intp dim = 0;
2173 ret = Py_BuildValue("(NNN)",
2174 PyArray_SimpleNew(1, &dim, NPY_INT),
2175 PyArray_SimpleNew(1, &dim, NPY_DOUBLE),
2176 PyArray_SimpleNew(1, &dim, NPY_DOUBLE));
2177 goto cleanup;
2178 }
2179
2180 op_flags[0] = NPY_ITER_READONLY | NPY_ITER_NBO;
2181 op_flags[1] = NPY_ITER_READONLY | NPY_ITER_NBO;
2182 op_flags[2] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO | NPY_ITER_CONTIG | NPY_ITER_ALIGNED;
2183 op_flags[3] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO | NPY_ITER_CONTIG | NPY_ITER_ALIGNED;
2184 op_flags[4] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO | NPY_ITER_CONTIG | NPY_ITER_ALIGNED;
2185
2186 dtypes[0] = PyArray_DescrFromType(NPY_DOUBLE);
2187 dtypes[1] = PyArray_DescrFromType(NPY_DOUBLE);
2188 dtypes[2] = PyArray_DescrFromType(NPY_DOUBLE);
2189 dtypes[3] = PyArray_DescrFromType(NPY_DOUBLE);
2190 dtypes[4] = PyArray_DescrFromType(NPY_INT);
2191
2192 iter = NpyIter_MultiNew(5, op, flags, NPY_KEEPORDER, NPY_SAFE_CASTING,
2193 op_flags, dtypes);
2194 for (j=0; j<5; j++)
2195 Py_DECREF(dtypes[j]);
2196
2197 if (!iter)
2198 return NULL;
2199
2200 iternext = NpyIter_GetIterNext(iter, NULL);
2201 strideptr = NpyIter_GetInnerStrideArray(iter);
2202 // The inner loop size and data pointers may change during the
2203 // loop, so just cache the addresses.
2204 innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);
2205 dataptrarray = NpyIter_GetDataPtrArray(iter);
2206
2207 do {
2208 npy_intp i, N;
2209 char* src1 = dataptrarray[0];
2210 char* src2 = dataptrarray[1];
2211 double* dout1 = (double*)(dataptrarray[2]);
2212 double* dout2 = (double*)(dataptrarray[3]);
2213 int* ok = (int*)dataptrarray[4];
2214 N = *innersizeptr;
2215
2216 //printf("2to2i: N=%i, strides %i,%i\n", N, strideptr[0], strideptr[1]);
2217
2218 // are the inputs contiguous? (Outputs will be, since we
2219 // allocated them)
2220 if ((strideptr[0] == sizeof(double)) &&
2221 (strideptr[1] == sizeof(double))) {
2222 // printf("Contiguous inputs; going fast\n");
2223 double* din1 = (double*)src1;
2224 double* din2 = (double*)src2;
2225 while (N--) {
2226 *ok = func(baton, *din1, *din2, dout1, dout2);
2227 ok++;
2228 din1++;
2229 din2++;
2230 dout1++;
2231 dout2++;
2232 }
2233 } else {
2234 // printf("Non-contiguous inputs; going slow\n");
2235 npy_intp stride1 = strideptr[0];
2236 npy_intp stride2 = strideptr[1];
2237 for (i=0; i<N; i++) {
2238 *ok = func(baton, *((double*)src1), *((double*)src2),
2239 dout1, dout2);
2240 ok++;
2241 src1 += stride1;
2242 src2 += stride2;
2243 dout1++;
2244 dout2++;
2245 }
2246 }
2247 } while (iternext(iter));
2248
2249 if (PyArray_IsPythonScalar(in1) && PyArray_IsPythonScalar(in2)) {
2250 PyObject* px = (PyObject*)(NpyIter_GetOperandArray(iter)[2]);
2251 PyObject* py = (PyObject*)(NpyIter_GetOperandArray(iter)[3]);
2252 PyObject* pok = (PyObject*)(NpyIter_GetOperandArray(iter)[4]);
2253 //printf("Both inputs are python scalars\n");
2254 double d;
2255 int i;
2256 d = *(double*)PyArray_DATA((PyArrayObject*)px);
2257 px = PyFloat_FromDouble(d);
2258 d = *(double*)PyArray_DATA((PyArrayObject*)py);
2259 py = PyFloat_FromDouble(d);
2260 i = *(int*)PyArray_DATA((PyArrayObject*)pok);
2261 pok = PyInt_FromLong(i);
2262 ret = Py_BuildValue("(NNN)", pok, px, py);
2263 } else {
2264 // Grab the results -- note "4,2,3" order -- ok,x,y
2265 ret = Py_BuildValue("(OOO)",
2266 NpyIter_GetOperandArray(iter)[4],
2267 NpyIter_GetOperandArray(iter)[2],
2268 NpyIter_GetOperandArray(iter)[3]);
2269 }
2270 cleanup:
2271 if (NpyIter_Deallocate(iter) != NPY_SUCCEED) {
2272 Py_DECREF(ret);
2273 return NULL;
2274 }
2275 Py_DECREF(op[0]);
2276 Py_DECREF(op[1]);
2277 return ret;
2278 }
2279
2280 static PyObject* broadcast_2to2
2281 (
2282 //void func(const void*, double, double, double*, double*),
2283 f_2to2 func,
2284 const void* baton,
2285 PyObject* in1, PyObject* in2) {
2286
2287 NpyIter *iter = NULL;
2288 NpyIter_IterNextFunc *iternext;
2289 PyArrayObject *op[4];
2290 PyObject *ret;
2291 npy_uint32 flags;
2292 npy_uint32 op_flags[4];
2293 npy_intp *innersizeptr;
2294 char **dataptrarray;
2295 npy_intp* strideptr;
2296 PyArray_Descr* dtypes[4];
2297 npy_intp i;
2298
2299 // we'll do the inner loop ourselves
2300 flags = NPY_ITER_EXTERNAL_LOOP;
2301 // use buffers to satisfy dtype casts
2302 flags |= NPY_ITER_BUFFERED;
2303 // grow inner loop
2304 flags |= NPY_ITER_GROWINNER;
2305
2306 op[0] = (PyArrayObject*)PyArray_FromAny(in1, NULL, 0, 0, 0, NULL);
2307 op[1] = (PyArrayObject*)PyArray_FromAny(in2, NULL, 0, 0, 0, NULL);
2308 // automatically allocate the output arrays.
2309 op[2] = NULL;
2310 op[3] = NULL;
2311
2312 if ((PyArray_Size((PyObject*)op[0]) == 0) ||
2313 (PyArray_Size((PyObject*)op[1]) == 0)) {
2314 // empty inputs -- empty outputs
2315 npy_intp dim = 0;
2316 ret = Py_BuildValue("(NN)",
2317 PyArray_SimpleNew(1, &dim, NPY_DOUBLE),
2318 PyArray_SimpleNew(1, &dim, NPY_DOUBLE));
2319 goto cleanup;
2320 }
2321
2322 op_flags[0] = NPY_ITER_READONLY | NPY_ITER_NBO;
2323 op_flags[1] = NPY_ITER_READONLY | NPY_ITER_NBO;
2324 op_flags[2] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO;
2325 op_flags[3] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO;
2326
2327 dtypes[0] = PyArray_DescrFromType(NPY_DOUBLE);
2328 dtypes[1] = PyArray_DescrFromType(NPY_DOUBLE);
2329 dtypes[2] = PyArray_DescrFromType(NPY_DOUBLE);
2330 dtypes[3] = PyArray_DescrFromType(NPY_DOUBLE);
2331
2332 iter = NpyIter_MultiNew(4, op, flags, NPY_KEEPORDER, NPY_SAFE_CASTING,
2333 op_flags, dtypes);
2334 for (i=0; i<4; i++)
2335 Py_DECREF(dtypes[i]);
2336 if (!iter)
2337 return NULL;
2338
2339 iternext = NpyIter_GetIterNext(iter, NULL);
2340 strideptr = NpyIter_GetInnerStrideArray(iter);
2341 // The inner loop size and data pointers may change during the
2342 // loop, so just cache the addresses.
2343 innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);
2344 dataptrarray = NpyIter_GetDataPtrArray(iter);
2345
2346 do {
2347 // are the inputs contiguous? (Outputs will be, since we
2348 // allocated them)
2349 if ((strideptr[0] == sizeof(double)) &&
2350 (strideptr[1] == sizeof(double))) {
2351
2352 npy_intp N = *innersizeptr;
2353 double* din1 = (double*)(dataptrarray[0]);
2354 double* din2 = (double*)(dataptrarray[1]);
2355 double* dout1 = (double*)(dataptrarray[2]);
2356 double* dout2 = (double*)(dataptrarray[3]);
2357
2358 //printf("Contiguous inputs; going fast\n");
2359 //printf("Inner loop: %i\n", (int)N);
2360 //printf("Output strides: %i %i\n", (int)strideptr[2], (int)strideptr[3]);
2361 //printf("Strides: %i %i %i %i\n", (int)strideptr[0], (int)strideptr[1], (int)strideptr[2], (int)strideptr[3]);
2362
2363 while (N--) {
2364 //printf("Calling %i: inputs (%12g,%12g)\n", (int)N, *din1, *din2);
2365 func(baton, *din1, *din2, dout1, dout2);
2366 din1++;
2367 din2++;
2368 dout1++;
2369 dout2++;
2370 }
2371 } else {
2372 npy_intp stride1 = NpyIter_GetInnerStrideArray(iter)[0];
2373 npy_intp stride2 = NpyIter_GetInnerStrideArray(iter)[1];
2374 npy_intp size = *innersizeptr;
2375 char* src1 = dataptrarray[0];
2376 char* src2 = dataptrarray[1];
2377 double* dout1 = (double*)(dataptrarray[2]);
2378 double* dout2 = (double*)(dataptrarray[3]);
2379
2380 //printf("Non-contiguous inputs; going slow\n");
2381 //printf("%i items\n", (int)size);
2382
2383 for (i=0; i<size; i++) {
2384 //printf("Call %i: inputs (%12g,%12g)\n", (int)i, ((double*)src1)[0], ((double*)src2)[0]);
2385 func(baton, *((double*)src1), *((double*)src2),
2386 dout1, dout2);
2387 src1 += stride1;
2388 src2 += stride2;
2389 dout1++;
2390 dout2++;
2391 }
2392 }
2393 } while (iternext(iter));
2394
2395 if (PyArray_IsPythonScalar(in1) && PyArray_IsPythonScalar(in2)) {
2396 PyObject* px = (PyObject*)(NpyIter_GetOperandArray(iter)[2]);
2397 PyObject* py = (PyObject*)(NpyIter_GetOperandArray(iter)[3]);
2398 //printf("Both inputs are python scalars\n");
2399 double d;
2400 d = *(double*)PyArray_DATA((PyArrayObject*)px);
2401 px = PyFloat_FromDouble(d);
2402 d = *(double*)PyArray_DATA((PyArrayObject*)py);
2403 py = PyFloat_FromDouble(d);
2404 ret = Py_BuildValue("(NN)", px, py);
2405 } else {
2406 // Grab the results
2407 ret = Py_BuildValue("(OO)",
2408 NpyIter_GetOperandArray(iter)[2],
2409 NpyIter_GetOperandArray(iter)[3]);
2410 }
2411
2412 cleanup:
2413 if (NpyIter_Deallocate(iter) != NPY_SUCCEED) {
2414 Py_DECREF(ret);
2415 return NULL;
2416 }
2417 Py_DECREF(op[0]);
2418 Py_DECREF(op[1]);
2419 return ret;
2420 }
2421
2422 static int tan_wcs_resample(tan_t* inwcs, tan_t* outwcs,
2423 PyObject* py_inimg, PyObject* py_outimg,
2424 int weighted, int lorder) {
2425 PyArray_Descr* dtype = PyArray_DescrFromType(NPY_FLOAT);
2426 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED | NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ELEMENTSTRIDES;
2427 int reqout = req | NPY_ARRAY_WRITEABLE | NPY_ARRAY_UPDATEIFCOPY;
2428 PyArrayObject *np_inimg=NULL, *np_outimg=NULL;
2429
2430 Py_INCREF(dtype);
2431 Py_INCREF(dtype);
2432 np_inimg = (PyArrayObject*)PyArray_CheckFromAny(py_inimg, dtype, 2, 2, req, NULL);
2433 np_outimg = (PyArrayObject*)PyArray_CheckFromAny(py_outimg, dtype, 2, 2, reqout, NULL);
2434 if (!np_inimg || !np_outimg) {
2435 ERR("Failed to PyArray_FromAny the images (np_inimg=%p, np_outimg=%p)\n",
2436 np_inimg, np_outimg);
2437 Py_XDECREF(np_inimg);
2438 Py_XDECREF(np_outimg);
2439 Py_DECREF(dtype);
2440 return -1;
2441 }
2442
2443 int inW, inH, outW, outH;
2444 float *inimg, *outimg;
2445 inH = (int)PyArray_DIM(np_inimg, 0);
2446 inW = (int)PyArray_DIM(np_inimg, 1);
2447 outH = (int)PyArray_DIM(np_outimg, 0);
2448 outW = (int)PyArray_DIM(np_outimg, 1);
2449 inimg = PyArray_DATA(np_inimg);
2450 outimg = PyArray_DATA(np_outimg);
2451
2452 anwcs_t* inanwcs = anwcs_new_tan(inwcs);
2453 anwcs_t* outanwcs = anwcs_new_tan(outwcs);
2454
2455 int res = resample_wcs(inanwcs, inimg, inW, inH,
2456 outanwcs, outimg, outW, outH,
2457 weighted, lorder);
2458
2459 anwcs_free(inanwcs);
2460 anwcs_free(outanwcs);
2461
2462 Py_DECREF(dtype);
2463 Py_DECREF(np_inimg);
2464 Py_DECREF(np_outimg);
2465
2466 return res;
2467 }
2468
2469 static int tan_numpy_xyz2pixelxy(tan_t* tan, PyArrayObject* npxyz,
2470 PyArrayObject* npx, PyArrayObject* npy) {
2471 npy_intp i, N;
2472 int rtn = 0;
2473 double *x, *y;
2474
2475 if (PyArray_NDIM(npx) != 1) {
2476 PyErr_SetString(PyExc_ValueError, "arrays must be one-dimensional");
2477 return -1;
2478 }
2479 if (PyArray_TYPE(npx) != NPY_DOUBLE) {
2480 PyErr_SetString(PyExc_ValueError, "array must contain doubles");
2481 return -1;
2482 }
2483 N = PyArray_DIM(npx, 0);
2484 if ((PyArray_DIM(npy, 0) != N) ||
2485 (PyArray_DIM(npxyz, 0) != N) ||
2486 (PyArray_DIM(npxyz, 1) != 3)) {
2487 PyErr_SetString(PyExc_ValueError, "arrays must be the same size");
2488 return -1;
2489 }
2490 x = PyArray_GETPTR1(npx, 0);
2491 y = PyArray_GETPTR1(npy, 0);
2492 for (i=0; i<N; i++) {
2493 double xyz[3];
2494 anbool ok;
2495 xyz[0] = *((double*)PyArray_GETPTR2(npxyz, i, 0));
2496 xyz[1] = *((double*)PyArray_GETPTR2(npxyz, i, 1));
2497 xyz[2] = *((double*)PyArray_GETPTR2(npxyz, i, 2));
2498 ok = tan_xyzarr2pixelxy(tan, xyz, x+i, y+i);
2499 if (!ok) {
2500 x[i] = -1.0;
2501 y[i] = -1.0;
2502 rtn = -1;
2503 }
2504 }
2505 return rtn;
2506 }
2507
2508
2509 static int an_tally(PyObject* py_counts, PyObject* py_x, PyObject* py_y) {
2510 PyArray_Descr* itype;
2511 PyArrayObject *np_counts=NULL, *np_x=NULL, *np_y=NULL;
2512 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED | NPY_ARRAY_NOTSWAPPED |
2513 NPY_ARRAY_ELEMENTSTRIDES;
2514 int reqout = req | NPY_ARRAY_WRITEABLE | NPY_ARRAY_UPDATEIFCOPY;
2515 int32_t *counts, *px, *py;
2516 int W, H, i, N;
2517
2518 itype = PyArray_DescrFromType(NPY_INT32);
2519 Py_INCREF(itype);
2520 Py_INCREF(itype);
2521
2522 np_counts = (PyArrayObject*)PyArray_CheckFromAny(py_counts, itype, 2, 2, reqout, NULL);
2523 np_x = (PyArrayObject*)PyArray_CheckFromAny(py_x, itype, 1, 1, req, NULL);
2524 np_y = (PyArrayObject*)PyArray_CheckFromAny(py_y, itype, 1, 1, req, NULL);
2525
2526 if (!np_counts || !np_x || !np_y) {
2527 ERR("Failed to PyArray_FromAny the counts, x, and y arrays.\n");
2528 Py_XDECREF(np_counts);
2529 Py_XDECREF(np_x);
2530 Py_XDECREF(np_y);
2531 return -1;
2532 }
2533 N = (int)PyArray_DIM(np_x, 0);
2534 if (PyArray_DIM(np_y, 0) != N) {
2535 ERR("Expected x and y arrays to have the same lengths!\n");
2536 Py_XDECREF(np_counts);
2537 Py_XDECREF(np_x);
2538 Py_XDECREF(np_y);
2539 return -1;
2540 }
2541
2542 H = (int)PyArray_DIM(np_counts, 0);
2543 W = (int)PyArray_DIM(np_counts, 1);
2544 //printf("Counts array size %i x %i\n", W, H);
2545 counts = PyArray_DATA(np_counts);
2546 px = PyArray_DATA(np_x);
2547 py = PyArray_DATA(np_y);
2548 for (i=0; i<N; i++) {
2549 int32_t xi = (*px);
2550 int32_t yi = (*py);
2551 if (yi < 0 || yi >= H || xi < 0 || xi >= W) {
2552 printf("Warning: skipping out-of-range value: i=%i, xi,yi = %i,%i\n", i, xi, yi);
2553 } else {
2554 counts[yi*W + xi]++;
2555 }
2556 px++;
2557 py++;
2558 }
2559 Py_DECREF(np_counts);
2560 Py_DECREF(np_x);
2561 Py_DECREF(np_y);
2562 return 0;
2563 }
2564
2565 %}
2566
2567 %pythoncode %{
2568 import numpy as np
2569
2570 def tan_t_tostring(self):
2571 ct = 'SIN' if self.sin else 'TAN'
2572 return ('%s: crpix (%.1f, %.1f), crval (%g, %g), cd (%g, %g, %g, %g), image %g x %g' %
2573 (ct, self.crpix[0], self.crpix[1], self.crval[0], self.crval[1],
2574 self.cd[0], self.cd[1], self.cd[2], self.cd[3],
2575 self.imagew, self.imageh))
2576 tan_t.__str__ = tan_t_tostring
2577
2578 def tan_t_addtoheader(self, hdr):
2579 '''Adds this TAN WCS header to the given fitsio header'''
2580 hdr.add_record(dict(name='CTYPE1', value='RA---TAN', comment='TANgent plane'))
2581 hdr.add_record(dict(name='CTYPE2', value='DEC--TAN', comment='TANgent plane'))
2582 hdr.add_record(dict(name='CRVAL1', value=self.crval[0], comment='Reference RA'))
2583 hdr.add_record(dict(name='CRVAL2', value=self.crval[1], comment='Reference Dec'))
2584 hdr.add_record(dict(name='CRPIX1', value=self.crpix[0], comment='Reference x'))
2585 hdr.add_record(dict(name='CRPIX2', value=self.crpix[1], comment='Reference y'))
2586 hdr.add_record(dict(name='CD1_1', value=self.cd[0], comment='CD matrix'))
2587 hdr.add_record(dict(name='CD1_2', value=self.cd[1], comment='CD matrix'))
2588 hdr.add_record(dict(name='CD2_1', value=self.cd[2], comment='CD matrix'))
2589 hdr.add_record(dict(name='CD2_2', value=self.cd[3], comment='CD matrix'))
2590 hdr.add_record(dict(name='IMAGEW', value=self.imagew, comment='Image width'))
2591 hdr.add_record(dict(name='IMAGEH', value=self.imageh, comment='Image height'))
2592 tan_t.add_to_header = tan_t_addtoheader
2593
2594 ## picklable?
2595 def tan_t_getstate(self):
2596 return (self.crpix[0], self.crpix[1], self.crval[0], self.crval[1],
2597 self.cd[0], self.cd[1], self.cd[2], self.cd[3],
2598 self.imagew, self.imageh, self.sin)
2599 def tan_t_setstate(self, state):
2600 #print 'setstate: self', self, 'state', state
2601 #print 'state', state
2602 self.this = _util.new_tan_t()
2603 #print 'self', repr(self)
2604 p0,p1,v0,v1,cd0,cd1,cd2,cd3,w,h,sin = state
2605 self.set_crpix(p0,p1)
2606 self.set_crval(v0,v1)
2607 self.set_cd(cd0,cd1,cd2,cd3)
2608 self.set_imagesize(w,h)
2609 self.sin = sin
2610 #(self.crpix[0], self.crpix[1], self.crval[0], self.crval[1],
2611 #self.cd[0], self.cd[1], self.cd[2], self.cd[3],
2612 #self.imagew, self.imageh) = state
2613 def tan_t_getnewargs(self):
2614 return ()
2615 tan_t.__getstate__ = tan_t_getstate
2616 tan_t.__setstate__ = tan_t_setstate
2617 tan_t.__getnewargs__ = tan_t_getnewargs
2618
2619 def tan_t_getshape(self):
2620 return int(self.imageh), int(self.imagew)
2621
2622 tan_t.shape = property(tan_t_getshape)
2623
2624 def tan_t_get_cd(self):
2625 cd = self.cd
2626 return (cd[0], cd[1], cd[2], cd[3])
2627 tan_t.get_cd = tan_t_get_cd
2628
2629 def tan_t_pixelxy2radec(self, x, y):
2630 return tan_xy2rd_wrapper(self.this, x, y)
2631 tan_t.pixelxy2radec_single = tan_t.pixelxy2radec
2632 tan_t.pixelxy2radec = tan_t_pixelxy2radec
2633
2634 def tan_t_radec2pixelxy(self, r, d):
2635 return tan_rd2xy_wrapper(self.this, r, d)
2636 tan_t.radec2pixelxy_single = tan_t.radec2pixelxy
2637 tan_t.radec2pixelxy = tan_t_radec2pixelxy
2638
2639 def tan_t_iwc2pixelxy(self, r, d):
2640 return tan_iwc2xy_wrapper(self.this, r, d)
2641 tan_t.iwc2pixelxy_single = tan_t.iwc2pixelxy
2642 tan_t.iwc2pixelxy = tan_t_iwc2pixelxy
2643
2644 def tan_t_pixelxy2iwc(self, x,y):
2645 return tan_xy2iwc_wrapper(self.this, x,y)
2646 tan_t.pixelxy2iwc_single = tan_t.pixelxy2iwc
2647 tan_t.pixelxy2iwc = tan_t_pixelxy2iwc
2648
2649 def tan_t_radec2iwc(self, r, d):
2650 return tan_rd2iwc_wrapper(self.this, r, d)
2651 tan_t.radec2iwc_single = tan_t.radec2iwc
2652 tan_t.radec2iwc = tan_t_radec2iwc
2653
2654 def tan_t_iwc2radec(self, u, v):
2655 return tan_iwc2rd_wrapper(self.this, u, v)
2656 tan_t.iwc2radec_single = tan_t.iwc2radec
2657 tan_t.iwc2radec = tan_t_iwc2radec
2658
2659 def sip_t_pixelxy2radec(self, x, y):
2660 return sip_xy2rd_wrapper(self.this, x, y)
2661 sip_t.pixelxy2radec_single = sip_t.pixelxy2radec
2662 sip_t.pixelxy2radec = sip_t_pixelxy2radec
2663
2664 def sip_t_radec2pixelxy(self, r, d):
2665 return sip_rd2xy_wrapper(self.this, r, d)
2666 sip_t.radec2pixelxy_single = sip_t.radec2pixelxy
2667 sip_t.radec2pixelxy = sip_t_radec2pixelxy
2668
2669 def sip_t_iwc2pixelxy(self, r, d):
2670 return sip_iwc2xy_wrapper(self.this, r, d)
2671 sip_t.iwc2pixelxy_single = sip_t.iwc2pixelxy
2672 sip_t.iwc2pixelxy = sip_t_iwc2pixelxy
2673
2674 def sip_t_pixelxy2iwc(self, x,y):
2675 return sip_xy2iwc_wrapper(self.this, x,y)
2676 sip_t.pixelxy2iwc_single = sip_t.pixelxy2iwc
2677 sip_t.pixelxy2iwc = sip_t_pixelxy2iwc
2678
2679 def sip_t_radec2iwc(self, r, d):
2680 return sip_rd2iwc_wrapper(self.this, r, d)
2681 sip_t.radec2iwc_single = sip_t.radec2iwc
2682 sip_t.radec2iwc = sip_t_radec2iwc
2683
2684 def sip_t_iwc2radec(self, u, v):
2685 return sip_iwc2rd_wrapper(self.this, u, v)
2686 sip_t.iwc2radec_single = sip_t.iwc2radec
2687 sip_t.iwc2radec = sip_t_iwc2radec
2688
2689
2690 def anwcs_t_pixelxy2radec(self, x, y):
2691 ok,r,d = anwcs_xy2rd_wrapper(self.this, x, y)
2692 return (ok == 0),r,d
2693 anwcs_t.pixelxy2radec_single = anwcs_t.pixelxy2radec
2694 anwcs_t.pixelxy2radec = anwcs_t_pixelxy2radec
2695
2696 def anwcs_t_radec2pixelxy(self, r, d):
2697 ok,x,y = anwcs_rd2xy_wrapper(self.this, r, d)
2698 return (ok == 0),x,y
2699 anwcs_t.radec2pixelxy_single = anwcs_t.radec2pixelxy
2700 anwcs_t.radec2pixelxy = anwcs_t_radec2pixelxy
2701
2702 def tan_t_radec_bounds(self):
2703 W,H = self.imagew, self.imageh
2704 r,d = self.pixelxy2radec([1, W/2, W, W, W, W/2, 1, 1], [1, 1, 1, H/2, H, H, H, H/2])
2705 rx = r.max()
2706 rn = r.min()
2707 # ugh, RA wrap-around. We find the largest value < 180 (ie, near zero) and smallest value > 180 (ie, near 360)
2708 # and report them with ralo > rahi so that this case can be identified
2709 if rx - rn > 180:
2710 rx = r[r < 180].max()
2711 rn = r[r > 180].min()
2712 return (rn, rx, d.min(), d.max())
2713 tan_t.radec_bounds = tan_t_radec_bounds
2714
2715 _real_tan_t_init = tan_t.__init__
2716 def my_tan_t_init(self, *args, **kwargs):
2717 # fitsio header: check for '.records()' function.
2718 if len(args) == 1 and hasattr(args[0], 'records'):
2719 try:
2720 hdr = args[0]
2721 qhdr = fitsio_to_qfits_header(hdr)
2722 args = [qhdr]
2723 except:
2724 pass
2725
2726 _real_tan_t_init(self, *args, **kwargs)
2727 if self.this is None:
2728 raise RuntimeError('Duck punch!')
2729 tan_t.__init__ = my_tan_t_init
2730
2731 Tan = tan_t
2732
2733 def tan_t_get_subimage(self, x0, y0, w, h):
2734 wcs2 = tan_t(self)
2735 cpx,cpy = wcs2.crpix
2736 wcs2.set_crpix(cpx - x0, cpy - y0)
2737 wcs2.set_width(float(w))
2738 wcs2.set_height(float(h))
2739 return wcs2
2740 tan_t.get_subimage = tan_t_get_subimage
2741
2742 # Deja Vu!
2743 # def sip_t_get_subimage(self, xlo, xhi, ylo, yhi):
2744 # sipout = sip_t(self)
2745 # sip_shift(self.this, sipout.this, float(xlo), float(xhi), float(ylo), float(yhi))
2746 # return sipout
2747 # sip_t.get_subimage = sip_t_get_subimage
2748
2749 # picklable
2750 def sip_t_getstate(self):
2751 t = (self.wcstan.__getstate__(),
2752 self.a_order, self.b_order, self.a, self.b,
2753 self.ap_order, self.bp_order, self.ap, self.bp)
2754 return t
2755
2756 def sip_t_setstate(self, s):
2757 self.this = _util.new_sip_t()
2758 (t, self.a_order, self.b_order, self.a, self.b,
2759 self.ap_order, self.bp_order, self.ap, self.bp) = s
2760 #self.wcstan.__setstate__(t)
2761 # disturbingly, tan_t_setstate does not work because it resets self.this = ...
2762 p0,p1,v0,v1,cd0,cd1,cd2,cd3,w,h,sin = t
2763 self.wcstan.set_crpix(p0,p1)
2764 self.wcstan.set_crval(v0,v1)
2765 self.wcstan.set_cd(cd0,cd1,cd2,cd3)
2766 self.wcstan.set_imagesize(w,h)
2767 self.wcstan.sin = sin
2768
2769 def sip_t_getnewargs(self):
2770 return ()
2771
2772 sip_t.__getstate__ = sip_t_getstate
2773 sip_t.__setstate__ = sip_t_setstate
2774 sip_t.__getnewargs__ = sip_t_getnewargs
2775
2776 %}
2777
2778 %include "fitsioutils.h"
2779
2780 // dcen3x3
2781 %apply float *OUTPUT { float *xcen, float *ycen };
2782
2783 %include "dimage.h"
2784
2785 %inline %{
2786 int dcen3x3b(float i0, float i1, float i2,
2787 float i3, float i4, float i5,
2788 float i6, float i7, float i8,
2789 float *xcen, float *ycen) {
2790 float im[9];
2791 im[0] = i0;
2792 im[1] = i1;
2793 im[2] = i2;
2794 im[3] = i3;
2795 im[4] = i4;
2796 im[5] = i5;
2797 im[6] = i6;
2798 im[7] = i7;
2799 im[8] = i8;
2800 return dcen3x3(im, xcen, ycen);
2801 }
2802 %}
2803