1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9 #include <math.h>
10 #include <assert.h>
11
12 #include "os-features.h"
13 #include "wcs-resample.h"
14 #include "sip_qfits.h"
15 #include "an-bool.h"
16 #include "anqfits.h"
17 #include "starutil.h"
18 #include "bl.h"
19 #include "boilerplate.h"
20 #include "log.h"
21 #include "errors.h"
22 #include "fitsioutils.h"
23 #include "anwcs.h"
24 #include "resample.h"
25 /** /# Modified by Robert Lancaster for the StellarSolver Internal Library
26 int resample_wcs_files(const char* infitsfn, int infitsext,
27 const char* inwcsfn, int inwcsext,
28 const char* outwcsfn, int outwcsext,
29 const char* outfitsfn, int lorder,
30 int zero_inf) {
31
32 anwcs_t* inwcs;
33 anwcs_t* outwcs;
34 anqfits_t* anqin;
35 qfitsdumper qoutimg;
36 float* inimg;
37 float* outimg;
38 qfits_header* hdr;
39
40 int outW, outH;
41 int inW, inH;
42
43 double outpixmin, outpixmax;
44
45 // read input WCS.
46 inwcs = anwcs_open(inwcsfn, inwcsext);
47 if (!inwcs) {
48 ERROR("Failed to parse WCS header from %s extension %i", inwcsfn, inwcsext);
49 return -1;
50 }
51
52 logmsg("Read input WCS from file \"%s\" ext %i\n", inwcsfn, inwcsext);
53 anwcs_print(inwcs, stdout);
54
55 // read output WCS.
56 outwcs = anwcs_open(outwcsfn, outwcsext);
57 if (!outwcs) {
58 ERROR("Failed to parse WCS header from %s extension %i", outwcsfn, outwcsext);
59 return -1;
60 }
61
62 logmsg("Read output (target) WCS from file \"%s\" ext %i\n", outwcsfn, outwcsext);
63 anwcs_print(outwcs, stdout);
64
65 outW = anwcs_imagew(outwcs);
66 outH = anwcs_imageh(outwcs);
67
68 // read input image.
69 anqin = anqfits_open(infitsfn);
70 if (!anqin) {
71 ERROR("Failed to open \"%s\"", infitsfn);
72 return -1;
73 }
74 inimg = (float*)anqfits_readpix(anqin, infitsext, 0, 0, 0, 0, 0,
75 PTYPE_FLOAT, NULL, &inW, &inH);
76 anqfits_close(anqin);
77 anqin = NULL;
78 if (!inimg) {
79 ERROR("Failed to read pixels from \"%s\"", infitsfn);
80 return -1;
81 }
82
83 if (zero_inf) {
84 int i;
85 for (i=0; i<(inW*inH); i++) {
86 if (!isfinite(inimg[i])) {
87 inimg[i] = 0.0;
88 }
89 }
90 }
91
92 logmsg("Input image is %i x %i pixels.\n", inW, inH);
93 logmsg("Output image is %i x %i pixels.\n", outW, outH);
94
95 outimg = calloc(outW * outH, sizeof(float));
96
97 if (resample_wcs(inwcs, inimg, inW, inH,
98 outwcs, outimg, outW, outH, 1, lorder)) {
99 ERROR("Failed to resample");
100 return -1;
101 }
102
103 {
104 double pmin, pmax;
105 int i;
106 /// /# Modified by Robert Lancaster for the StellarSolver Internal Library, this commented out part was interfering with commenting out the whole block
107 pmin = HUGE_VAL;
108 pmax = -HUGE_VAL;
109 for (i=0; i<(inW*inH); i++) {
110 pmin = MIN(pmin, inimg[i]);
111 pmax = MAX(pmax, inimg[i]);
112 }
113 logmsg("Input image bounds: %g to %g\n", pmin, pmax);
114 ///# Modified by Robert Lancaster for the StellarSolver Internal Library
115 pmin = HUGE_VAL;
116 pmax = -HUGE_VAL;
117 for (i=0; i<(outW*outH); i++) {
118 pmin = MIN(pmin, outimg[i]);
119 pmax = MAX(pmax, outimg[i]);
120 }
121 logmsg("Output image bounds: %g to %g\n", pmin, pmax);
122 outpixmin = pmin;
123 outpixmax = pmax;
124 }
125
126 // prepare output image.
127 memset(&qoutimg, 0, sizeof(qoutimg));
128 qoutimg.filename = outfitsfn;
129 qoutimg.npix = outW * outH;
130 qoutimg.ptype = PTYPE_FLOAT;
131 qoutimg.fbuf = outimg;
132 qoutimg.out_ptype = BPP_IEEE_FLOAT;
133
134 hdr = fits_get_header_for_image(&qoutimg, outW, NULL);
135 anwcs_add_to_header(outwcs, hdr);
136 fits_header_add_double(hdr, "DATAMIN", outpixmin, "min pixel value");
137 fits_header_add_double(hdr, "DATAMAX", outpixmax, "max pixel value");
138
139 if (fits_write_header_and_image(hdr, &qoutimg, 0)) {
140 ERROR("Failed to write image to file \"%s\"", outfitsfn);
141 return -1;
142 }
143 free(outimg);
144 qfits_header_destroy(hdr);
145
146 anwcs_free(inwcs);
147 anwcs_free(outwcs);
148
149 return 0;
150 }
151 **/
152 // Check whether output pixels overlap with input pixels,
153 // on a grid of output pixel positions.
find_overlap_grid(int B,int outW,int outH,const anwcs_t * outwcs,const anwcs_t * inwcs,int * pBW,int * pBH)154 static anbool* find_overlap_grid(int B, int outW, int outH,
155 const anwcs_t* outwcs, const anwcs_t* inwcs,
156 int* pBW, int* pBH) {
157 int BW, BH;
158 anbool* bib = NULL;
159 anbool* bib2 = NULL;
160 int i,j;
161
162 BW = (int)ceil(outW / (float)B);
163 BH = (int)ceil(outH / (float)B);
164 bib = calloc(BW*BH, sizeof(anbool));
165 for (i=0; i<BH; i++) {
166 for (j=0; j<BW; j++) {
167 int x,y;
168 double ra,dec;
169 y = MIN(outH-1, B*i);
170 x = MIN(outW-1, B*j);
171 if (anwcs_pixelxy2radec(outwcs, x+1, y+1, &ra, &dec))
172 continue;
173 bib[i*BW+j] = anwcs_radec_is_inside_image(inwcs, ra, dec);
174 }
175 }
176 if (log_get_level() >= LOG_VERB) {
177 logverb("Input image overlaps output image:\n");
178 for (i=0; i<BH; i++) {
179 for (j=0; j<BW; j++)
180 logverb((bib[i*BW+j]) ? "*" : ".");
181 logverb("\n");
182 }
183 }
184 // Grow the in-bounds area:
185 bib2 = calloc(BW*BH, sizeof(anbool));
186 for (i=0; i<BH; i++)
187 for (j=0; j<BW; j++) {
188 int di,dj;
189 if (!bib[i*BW+j])
190 continue;
191 for (di=-1; di<=1; di++)
192 for (dj=-1; dj<=1; dj++)
193 bib2[(MIN(MAX(i+di, 0), BH-1))*BW + (MIN(MAX(j+dj, 0), BW-1))] = TRUE;
194 }
195 // swap!
196 free(bib);
197 bib = bib2;
198 bib2 = NULL;
199
200 if (log_get_level() >= LOG_VERB) {
201 logverb("After growing:\n");
202 for (i=0; i<BH; i++) {
203 for (j=0; j<BW; j++)
204 logverb((bib[i*BW+j]) ? "*" : ".");
205 logverb("\n");
206 }
207 }
208
209 *pBW = BW;
210 *pBH = BH;
211 return bib;
212 }
213
214
215
resample_wcs(const anwcs_t * inwcs,const float * inimg,int inW,int inH,const anwcs_t * outwcs,float * outimg,int outW,int outH,int weighted,int lorder)216 int resample_wcs(const anwcs_t* inwcs, const float* inimg, int inW, int inH,
217 const anwcs_t* outwcs, float* outimg, int outW, int outH,
218 int weighted, int lorder) {
219 int i,j;
220 int jlo,jhi,ilo,ihi;
221 lanczos_args_t largs;
222 double xyz[3];
223 memset(&largs, 0, sizeof(largs));
224 largs.order = lorder;
225 largs.weighted = weighted;
226
227 jlo = ilo = 0;
228 ihi = outW;
229 jhi = outH;
230
231 // find the center of "outwcs", and describe the boundary as a
232 // polygon in ~IWC coordinates. Then find the extent of "inwcs".
233
234 {
235 int ok = 1;
236 double xmin, xmax, ymin, ymax;
237 int x, y, W, H;
238 double xx,yy;
239 xmin = ymin = HUGE_VAL;
240 xmax = ymax = -HUGE_VAL;
241 // HACK -- just look at the corners. Could anwcs_walk_boundary.
242 W = anwcs_imagew(inwcs);
243 H = anwcs_imageh(inwcs);
244 for (x=0; x<2; x++) {
245 for (y=0; y<2; y++) {
246 if (anwcs_pixelxy2xyz(inwcs, 1+x * (W-1), 1+y * (H-1), xyz) ||
247 anwcs_xyz2pixelxy(outwcs, xyz, &xx, &yy)) {
248 ok = 0;
249 break;
250 }
251 xmin = MIN(xmin, xx);
252 xmax = MAX(xmax, xx);
253 ymin = MIN(ymin, yy);
254 ymax = MAX(ymax, yy);
255 }
256 if (!ok)
257 break;
258 }
259 if (ok) {
260 W = anwcs_imagew(outwcs);
261 H = anwcs_imageh(outwcs);
262 if ((xmin >= W) || (xmax < 0) ||
263 (ymin >= H) || (ymax < 0)) {
264 logverb("No overlap between input and output images\n");
265 return 0;
266 }
267 ilo = MAX(0, xmin);
268 ihi = MIN(W, xmax);
269 jlo = MAX(0, ymin);
270 jhi = MIN(H, ymax);
271 //logverb("Clipped resampling output box to [%i,%i], [%i,%i]\n",
272 //ilo,ihi,jlo,jhi);
273 }
274 }
275
276 for (j=jlo; j<jhi; j++) {
277 for (i=ilo; i<ihi; i++) {
278 double inx, iny;
279 float pix;
280 // +1 for FITS pixel coordinates.
281 if (anwcs_pixelxy2xyz(outwcs, i+1, j+1, xyz) ||
282 anwcs_xyz2pixelxy(inwcs, xyz, &inx, &iny))
283 continue;
284
285 inx -= 1.0;
286 iny -= 1.0;
287
288 if (lorder == 0) {
289 int x,y;
290 // Nearest-neighbour resampling
291 x = round(inx);
292 y = round(iny);
293 if (x < 0 || x >= inW || y < 0 || y >= inH)
294 continue;
295 pix = inimg[y * inW + x];
296 } else {
297 if (inx < (-lorder) || inx >= (inW+lorder) ||
298 iny < (-lorder) || iny >= (inH+lorder))
299 continue;
300 pix = lanczos_resample_unw_sep_f(inx, iny, inimg, inW, inH, &largs);
301 }
302 outimg[j * outW + i] = pix;
303 }
304 }
305 return 0;
306 }
307
308
resample_wcs_rgba(const anwcs_t * inwcs,const unsigned char * inimg,int inW,int inH,const anwcs_t * outwcs,unsigned char * outimg,int outW,int outH)309 int resample_wcs_rgba(const anwcs_t* inwcs, const unsigned char* inimg,
310 int inW, int inH,
311 const anwcs_t* outwcs, unsigned char* outimg,
312 int outW, int outH) {
313 int i,j;
314 int B = 20;
315 int BW, BH;
316 anbool* bib;
317 int bi,bj;
318
319 bib = find_overlap_grid(B, outW, outH, outwcs, inwcs, &BW, &BH);
320
321 // We've expanded the in-bounds boxes by 1 in each direction,
322 // so this (using the lower-left corner) should be ok.
323 for (bj=0; bj<BH; bj++) {
324 for (bi=0; bi<BW; bi++) {
325 int jlo,jhi,ilo,ihi;
326 if (!bib[bj*BW + bi])
327 continue;
328 jlo = MIN(outH, bj *B);
329 jhi = MIN(outH, (bj+1)*B);
330 ilo = MIN(outW, bi *B);
331 ihi = MIN(outW, (bi+1)*B);
332 for (j=jlo; j<jhi; j++) {
333 for (i=ilo; i<ihi; i++) {
334 double xyz[3];
335 double inx, iny;
336 int x,y;
337 // +1 for FITS pixel coordinates.
338 if (anwcs_pixelxy2xyz(outwcs, i+1, j+1, xyz) ||
339 anwcs_xyz2pixelxy(inwcs, xyz, &inx, &iny))
340 continue;
341 // FIXME - Nearest-neighbour resampling!!
342 // -1 for FITS pixel coordinates.
343 x = round(inx - 1.0);
344 y = round(iny - 1.0);
345 if (x < 0 || x >= inW || y < 0 || y >= inH)
346 continue;
347 // HACK -- straight copy
348 outimg[4 * (j * outW + i) + 0] = inimg[4 * (y * inW + x) + 0];
349 outimg[4 * (j * outW + i) + 1] = inimg[4 * (y * inW + x) + 1];
350 outimg[4 * (j * outW + i) + 2] = inimg[4 * (y * inW + x) + 2];
351 outimg[4 * (j * outW + i) + 3] = inimg[4 * (y * inW + x) + 3];
352 }
353 }
354 }
355 }
356
357 free(bib);
358
359 return 0;
360
361 }
362
363