1 /*
2 Author: Michael Droettboom
3 mdroe@stsci.edu
4 */
5
6 #include "astropy_wcs/sip.h"
7
8 #include <assert.h>
9 #include <stdlib.h>
10 #include <string.h>
11
12 #include <wcserr.h>
13
14 #define SIP_ERRMSG(status) WCSERR_SET(status)
15
16 void
sip_clear(sip_t * sip)17 sip_clear(
18 sip_t* sip) {
19
20 assert(sip != NULL);
21
22 sip->a_order = 0;
23 sip->a = NULL;
24 sip->b_order = 0;
25 sip->b = NULL;
26 sip->ap_order = 0;
27 sip->ap = NULL;
28 sip->bp_order = 0;
29 sip->bp = NULL;
30 sip->crpix[0] = 0.0;
31 sip->crpix[1] = 0.0;
32 sip->scratch = NULL;
33 sip->err = NULL;
34 }
35
36 int
sip_init(sip_t * sip,const unsigned int a_order,const double * a,const unsigned int b_order,const double * b,const unsigned int ap_order,const double * ap,const unsigned int bp_order,const double * bp,const double * crpix)37 sip_init(
38 sip_t* sip,
39 const unsigned int a_order, const double* a,
40 const unsigned int b_order, const double* b,
41 const unsigned int ap_order, const double* ap,
42 const unsigned int bp_order, const double* bp,
43 const double* crpix /* [2] */) {
44
45 unsigned int a_size = 0;
46 unsigned int b_size = 0;
47 unsigned int ap_size = 0;
48 unsigned int bp_size = 0;
49 unsigned int scratch_size = 0;
50 int status = 0;
51 struct wcserr** err = NULL;
52 static const char *function = "sip_init";
53
54 assert(sip != NULL);
55 sip_clear(sip);
56 err = &(sip->err);
57
58 /* We we have one of A/B or AP/BP, we must have both. */
59 if ((a == NULL) ^ (b == NULL)) {
60 return wcserr_set(
61 SIP_ERRMSG(WCSERR_BAD_COORD_TRANS),
62 "Both A and B SIP transform must be defined");
63 }
64
65 if ((ap == NULL) ^ (bp == NULL)) {
66 return wcserr_set(
67 SIP_ERRMSG(WCSERR_BAD_COORD_TRANS),
68 "Both AP and BP SIP transform must be defined");
69 }
70
71
72 if (a != NULL) {
73 sip->a_order = a_order;
74 a_size = (a_order + 1) * (a_order + 1) * sizeof(double);
75 sip->a = malloc(a_size);
76 if (sip->a == NULL) {
77 sip_free(sip);
78 status = wcserr_set(
79 SIP_ERRMSG(WCSERR_MEMORY), "Memory allocation failed");
80 goto exit;
81 }
82 memcpy(sip->a, a, a_size);
83 if (a_order > scratch_size) {
84 scratch_size = a_order;
85 }
86
87 sip->b_order = b_order;
88 b_size = (b_order + 1) * (b_order + 1) * sizeof(double);
89 sip->b = malloc(b_size);
90 if (sip->b == NULL) {
91 sip_free(sip);
92 status = wcserr_set(
93 SIP_ERRMSG(WCSERR_MEMORY), "Memory allocation failed");
94 goto exit;
95 }
96 memcpy(sip->b, b, b_size);
97 if (b_order > scratch_size) {
98 scratch_size = b_order;
99 }
100 }
101
102 if (ap != NULL) {
103 sip->ap_order = ap_order;
104 ap_size = (ap_order + 1) * (ap_order + 1) * sizeof(double);
105 sip->ap = malloc(ap_size);
106 if (sip->ap == NULL) {
107 sip_free(sip);
108 status = wcserr_set(
109 SIP_ERRMSG(WCSERR_MEMORY), "Memory allocation failed");
110 goto exit;
111 }
112 memcpy(sip->ap, ap, ap_size);
113 if (ap_order > scratch_size) {
114 scratch_size = ap_order;
115 }
116
117 sip->bp_order = bp_order;
118 bp_size = (bp_order + 1) * (bp_order + 1) * sizeof(double);
119 sip->bp = malloc(bp_size);
120 if (sip->bp == NULL) {
121 sip_free(sip);
122 status = wcserr_set(
123 SIP_ERRMSG(WCSERR_MEMORY), "Memory allocation failed");
124 goto exit;
125 }
126 memcpy(sip->bp, bp, bp_size);
127 if (bp_order > scratch_size) {
128 scratch_size = bp_order;
129 }
130 }
131
132 scratch_size = (scratch_size + 1) * sizeof(double);
133 sip->scratch = malloc(scratch_size);
134 if (sip->scratch == NULL) {
135 sip_free(sip);
136 status = wcserr_set(
137 SIP_ERRMSG(WCSERR_MEMORY), "Memory allocation failed");
138 goto exit;
139 }
140
141 sip->crpix[0] = crpix[0];
142 sip->crpix[1] = crpix[1];
143
144 exit:
145
146 return status;
147 }
148
149 void
sip_free(sip_t * sip)150 sip_free(sip_t* sip) {
151 free(sip->a);
152 sip->a = NULL;
153 free(sip->b);
154 sip->b = NULL;
155 free(sip->ap);
156 sip->ap = NULL;
157 free(sip->bp);
158 sip->bp = NULL;
159 free(sip->scratch);
160 sip->scratch = NULL;
161 free(sip->err);
162 sip->err = NULL;
163 }
164
165 static INLINE double
lu(const unsigned int order,const double * const matrix,const int x,const int y)166 lu(
167 const unsigned int order,
168 const double* const matrix,
169 const int x,
170 const int y) {
171
172 int index;
173 assert(x >= 0 && x <= (int)order);
174 assert(y >= 0 && y <= (int)order);
175
176 index = x * ((int)order + 1) + y;
177 assert(index >= 0 && index < ((int)order + 1) * ((int)order + 1));
178
179 return matrix[index];
180 }
181
182 static int
sip_compute(const unsigned int naxes,const unsigned int nelem,const unsigned int m,const double * a,const unsigned int n,const double * b,const double * crpix,double * tmp,const double * input,double * output)183 sip_compute(
184 /*@unused@*/ const unsigned int naxes,
185 const unsigned int nelem,
186 const unsigned int m,
187 /*@null@*/ const double* a,
188 const unsigned int n,
189 /*@null@*/ const double* b,
190 const double* crpix /* [2] */,
191 /*@null@*/ double* tmp,
192 /*@null@*/ const double* input /* [NAXES][nelem] */,
193 /*@null@*/ double* output /* [NAXES][nelem] */) {
194
195 unsigned int i;
196 int j, k;
197 double x, y;
198 double sum;
199 const double* input_ptr;
200 double* output_ptr;
201
202 assert(a != NULL);
203 assert(b != NULL);
204 assert(crpix != NULL);
205 assert(tmp != NULL);
206 assert(input != NULL);
207 assert(output != NULL);
208
209 /* Avoid segfaults */
210 if (input == NULL || output == NULL || tmp == NULL || crpix == NULL) {
211 return 1;
212 }
213
214 /* If we have one, we must have both... */
215 if ((a == NULL) ^ (b == NULL)) {
216 return 6;
217 }
218
219 /* If no distortion, just return values */
220 if (a == NULL /* && b == NULL ... implied */) {
221 return 0;
222 }
223
224 input_ptr = input;
225 output_ptr = output;
226 for (i = 0; i < nelem; ++i) {
227 x = *input_ptr++ - crpix[0];
228 y = *input_ptr++ - crpix[1];
229
230 for (j = 0; j <= (int)m; ++j) {
231 tmp[j] = lu(m, a, (int)m-j, j);
232 for (k = j-1; k >= 0; --k) {
233 tmp[j] = (y * tmp[j]) + lu(m, a, (int)m-j, k);
234 }
235 }
236
237 sum = tmp[0];
238 for (j = (int)m; j > 0; --j) {
239 sum = x * sum + tmp[(int)m - j + 1];
240 }
241 *output_ptr++ += sum;
242
243 for (j = 0; j <= (int)n; ++j) {
244 tmp[j] = lu(n, b, (int)n-j, j);
245 for (k = j-1; k >= 0; --k) {
246 tmp[j] = (y * tmp[j]) + lu(n, b, (int)n-j, k);
247 }
248 }
249
250 sum = tmp[0];
251 for (j = (int)n; j > 0; --j) {
252 sum = x * sum + tmp[n - j + 1];
253 }
254 *output_ptr++ += sum;
255 }
256
257 return 0;
258 }
259
260 int
sip_pix2deltas(const sip_t * sip,const unsigned int naxes,const unsigned int nelem,const double * pix,double * deltas)261 sip_pix2deltas(
262 const sip_t* sip,
263 const unsigned int naxes,
264 const unsigned int nelem,
265 const double* pix /* [NAXES][nelem] */,
266 double* deltas /* [NAXES][nelem] */) {
267
268 if (sip == NULL) {
269 return 1;
270 }
271
272 return sip_compute(naxes, nelem,
273 sip->a_order, sip->a,
274 sip->b_order, sip->b,
275 sip->crpix,
276 (double *)sip->scratch,
277 pix, deltas);
278 }
279
280 int
sip_foc2deltas(const sip_t * sip,const unsigned int naxes,const unsigned int nelem,const double * foc,double * deltas)281 sip_foc2deltas(
282 const sip_t* sip,
283 const unsigned int naxes,
284 const unsigned int nelem,
285 const double* foc /* [NAXES][nelem] */,
286 double* deltas /* [NAXES][nelem] */) {
287
288 if (sip == NULL) {
289 return 1;
290 }
291
292 return sip_compute(naxes, nelem,
293 sip->ap_order, sip->ap,
294 sip->bp_order, sip->bp,
295 sip->crpix,
296 (double *)sip->scratch,
297 foc, deltas);
298 }
299
300 int
sip_pix2foc(const sip_t * sip,const unsigned int naxes,const unsigned int nelem,const double * pix,double * foc)301 sip_pix2foc(
302 const sip_t* sip,
303 const unsigned int naxes,
304 const unsigned int nelem,
305 const double* pix /* [NAXES][nelem] */,
306 double* foc /* [NAXES][nelem] */) {
307 assert(pix);
308 assert(foc);
309
310 if (pix != foc) {
311 memcpy(foc, pix, sizeof(double) * naxes * nelem);
312 }
313
314 return sip_pix2deltas(sip, naxes, nelem, pix, foc);
315 }
316
317 int
sip_foc2pix(const sip_t * sip,const unsigned int naxes,const unsigned int nelem,const double * foc,double * pix)318 sip_foc2pix(
319 const sip_t* sip,
320 const unsigned int naxes,
321 const unsigned int nelem,
322 const double* foc /* [NAXES][nelem] */,
323 double* pix /* [NAXES][nelem] */) {
324 assert(pix);
325 assert(foc);
326
327 if (pix != foc) {
328 memcpy(pix, foc, sizeof(double) * naxes * nelem);
329 }
330
331 return sip_foc2deltas(sip, naxes, nelem, foc, pix);
332 }
333