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 <string.h>
7 #include <stdio.h>
8 #include <unistd.h>
9 #include <assert.h>
10 #include <sys/types.h>
11 #include <unistd.h>
12 #include <math.h>
13
14 #include "os-features.h"
15 #include "image2xy-files.h"
16 #include "image2xy.h"
17 #include "fitsio.h"
18 #include "ioutils.h"
19 #include "simplexy.h"
20 #include "errors.h"
21 #include "log.h"
22 #include "cfitsutils.h"
23
image2xy_files(const char * infn,const char * outfn,anbool do_u8,int downsample,int downsample_as_required,int extension,int plane,simplexy_t * params)24 int image2xy_files(const char* infn, const char* outfn,
25 anbool do_u8, int downsample, int downsample_as_required,
26 int extension, int plane,
27 simplexy_t* params) {
28 fitsfile *fptr = NULL;
29 fitsfile *ofptr = NULL;
30 int status = 0; // FIXME should have ostatus too
31 int naxis;
32 long naxisn[3];
33 int kk;
34 int nhdus, hdutype, nimgs;
35 char* str;
36 simplexy_t myparams;
37
38 if (params == NULL) {
39 memset(&myparams, 0, sizeof(simplexy_t));
40 params = &myparams;
41 }
42
43 // QFITS to CFITSIO extension convention switch
44 extension++;
45
46 fits_open_file(&fptr, infn, READONLY, &status);
47 CFITS_CHECK("Failed to open FITS input file %s", infn);
48
49 // Are there multiple HDU's?
50 fits_get_num_hdus(fptr, &nhdus, &status);
51 CFITS_CHECK("Failed to read number of HDUs for input file %s", infn);
52 logverb("nhdus=%d\n", nhdus);
53
54 if (extension > nhdus) {
55 logerr("Requested extension %i is greater than number of extensions (%i) in file %s\n",
56 extension, nhdus, infn);
57 return -1;
58 }
59
60 // Create output file
61 fits_create_file(&ofptr, outfn, &status);
62 CFITS_CHECK("Failed to open FITS output file %s", outfn);
63
64 fits_create_img(ofptr, 8, 0, NULL, &status);
65 CFITS_CHECK("Failed to create output image");
66
67 fits_write_key(ofptr, TSTRING, "SRCFN", (char*)infn, "Source image", &status);
68 if (extension)
69 fits_write_key(ofptr, TINT, "SRCEXT", &extension, "Source image extension (1=primary)", &status);
70
71 /* Parameters for simplexy; save for debugging */
72 fits_write_comment(ofptr, "Parameters used for source extraction", &status);
73
74 fits_write_history(ofptr, "Created by Astrometry.net's image2xy program.", &status);
75 CFITS_CHECK("Failed to write HISTORY headers");
76
77 asprintf_safe(&str, "GIT URL: %s", AN_GIT_URL);
78 fits_write_history(ofptr, str, &status);
79 CFITS_CHECK("Failed to write GIT HISTORY headers");
80 free(str);
81 asprintf_safe(&str, "GIT Rev: %s", AN_GIT_REVISION);
82 fits_write_history(ofptr, str, &status);
83 CFITS_CHECK("Failed to write GIT HISTORY headers");
84 free(str);
85 asprintf_safe(&str, "GIT Date: %s", AN_GIT_DATE);
86 fits_write_history(ofptr, str, &status);
87 CFITS_CHECK("Failed to write GIT HISTORY headers");
88 free(str);
89 fits_write_history(ofptr, "Visit us on the web at http://astrometry.net/", &status);
90 CFITS_CHECK("Failed to write GIT HISTORY headers");
91
92 nimgs = 0;
93
94 // Run simplexy on each HDU
95 for (kk=1; kk <= nhdus; kk++) {
96 char* ttype[] = {"X","Y","FLUX","BACKGROUND", "LFLUX", "LBG"};
97 char* tform[] = {"E","E","E","E","E","E"};
98 char* tunit[] = {"pix","pix","unknown","unknown","unknown","unknown"};
99 long* fpixel;
100 int a;
101 int w, h;
102 int bitpix;
103 int ncols;
104
105 if (extension && kk != extension)
106 continue;
107
108 fits_movabs_hdu(fptr, kk, &hdutype, &status);
109 fits_get_hdu_type(fptr, &hdutype, &status);
110
111 if (hdutype != IMAGE_HDU) {
112 if (extension)
113 logerr("Requested extension %i in file %s is not an image.\n", extension, infn);
114 continue;
115 }
116
117 fits_get_img_dim(fptr, &naxis, &status);
118 CFITS_CHECK("Failed to find image dimensions for HDU %i", kk);
119
120 fits_get_img_size(fptr, 2, naxisn, &status);
121 CFITS_CHECK("Failed to find image dimensions for HDU %i", kk);
122
123 nimgs++;
124
125 logverb("Got naxis=%d, na1=%lu, na2=%lu\n", naxis, naxisn[0], naxisn[1]);
126
127 fits_get_img_type(fptr, &bitpix, &status);
128 CFITS_CHECK("Failed to get FITS image type");
129
130 fpixel = malloc(naxis * sizeof(long));
131 for (a=0; a<naxis; a++)
132 fpixel[a] = 1;
133
134 if (plane && naxis == 3) {
135 if (plane <= naxisn[2]) {
136 logmsg("Grabbing image plane %i\n", plane);
137 fpixel[2] = plane;
138 } else
139 logerr("Requested plane %i but only %i are available.\n", plane, (int)naxisn[2]);
140 } else if (plane)
141 logmsg("Plane %i requested but this image has NAXIS = %i (not 3).\n", plane, naxis);
142 else if (naxis > 2)
143 logmsg("This looks like a multi-color image: processing the first image plane only. (NAXIS=%i)\n", naxis);
144
145 if (bitpix == 8 && do_u8 && !downsample) {
146 simplexy_fill_in_defaults_u8(params);
147
148 // u8 image.
149 params->image_u8 = malloc(naxisn[0] * naxisn[1]);
150 if (!params->image_u8) {
151 SYSERROR("Failed to allocate u8 image array");
152 goto bailout;
153 }
154 fits_read_pix(fptr, TBYTE, fpixel, naxisn[0]*naxisn[1], NULL,
155 params->image_u8, NULL, &status);
156
157 } else {
158 simplexy_fill_in_defaults(params);
159
160 params->image = malloc(naxisn[0] * naxisn[1] * sizeof(float));
161 if (!params->image) {
162 SYSERROR("Failed to allocate image array");
163 goto bailout;
164 }
165 fits_read_pix(fptr, TFLOAT, fpixel, naxisn[0]*naxisn[1], NULL,
166 params->image, NULL, &status);
167 }
168 free(fpixel);
169 CFITS_CHECK("Failed to read image pixels");
170
171 params->nx = naxisn[0];
172 params->ny = naxisn[1];
173
174 image2xy_run(params, downsample, downsample_as_required);
175
176 if (params->Lorder)
177 ncols = 6;
178 else
179 ncols = 4;
180 fits_create_tbl(ofptr, BINARY_TBL, params->npeaks, ncols, ttype, tform,
181 tunit, "SOURCES", &status);
182 CFITS_CHECK("Failed to create output table");
183
184 fits_write_col(ofptr, TFLOAT, 1, 1, 1, params->npeaks, params->x, &status);
185 CFITS_CHECK("Failed to write X column");
186
187 fits_write_col(ofptr, TFLOAT, 2, 1, 1, params->npeaks, params->y, &status);
188 CFITS_CHECK("Failed to write Y column");
189
190 fits_write_col(ofptr, TFLOAT, 3, 1, 1, params->npeaks, params->flux, &status);
191 CFITS_CHECK("Failed to write FLUX column");
192
193 fits_write_col(ofptr, TFLOAT, 4, 1, 1, params->npeaks, params->background, &status);
194 CFITS_CHECK("Failed to write BACKGROUND column");
195
196 if (params->Lorder) {
197 fits_write_col(ofptr, TFLOAT, 5, 1, 1, params->npeaks, params->fluxL, &status);
198 CFITS_CHECK("Failed to write LFLUX column");
199
200 fits_write_col(ofptr, TFLOAT, 6, 1, 1, params->npeaks, params->backgroundL, &status);
201 CFITS_CHECK("Failed to write LBG column");
202 }
203
204 fits_modify_comment(ofptr, "TTYPE1", "X coordinate", &status);
205 CFITS_CHECK("Failed to set X TTYPE");
206
207 fits_modify_comment(ofptr, "TTYPE2", "Y coordinate", &status);
208 CFITS_CHECK("Failed to set Y TTYPE");
209
210 fits_modify_comment(ofptr, "TTYPE3", "Flux of source", &status);
211 CFITS_CHECK("Failed to set FLUX TTYPE");
212
213 fits_modify_comment(ofptr, "TTYPE4", "Sky background of source", &status);
214 CFITS_CHECK("Failed to set BACKGROUND TTYPE");
215
216 fits_write_key(ofptr, TINT, "SRCEXT", &kk,
217 "Extension number in src image", &status);
218 CFITS_CHECK("Failed to write SRCEXT");
219
220 w = naxisn[0];
221 h = naxisn[1];
222 fits_write_key(ofptr, TINT, "IMAGEW", &w, "Input image width", &status);
223 CFITS_CHECK("Failed to write IMAGEW");
224
225 fits_write_key(ofptr, TINT, "IMAGEH", &h, "Input image height", &status);
226 CFITS_CHECK("Failed to write IMAGEH");
227
228 fits_write_key(ofptr, TFLOAT, "ESTSIGMA", &(params->sigma),
229 "Estimated source image variance", &status);
230 CFITS_CHECK("Failed to write ESTSIGMA");
231
232 fits_write_key(ofptr, TFLOAT, "DPSF", &(params->dpsf), "image2xy Assumed gaussian psf width", &status);
233 fits_write_key(ofptr, TFLOAT, "PLIM", &(params->plim), "image2xy Significance to keep", &status);
234 fits_write_key(ofptr, TFLOAT, "DLIM", &(params->dlim), "image2xy Closest two peaks can be", &status);
235 fits_write_key(ofptr, TFLOAT, "SADDLE", &(params->saddle), "image2xy Saddle difference (in sig)", &status);
236 fits_write_key(ofptr, TINT, "MAXPER", &(params->maxper), "image2xy Max num of peaks per object", &status);
237 fits_write_key(ofptr, TINT, "MAXPEAKS", &(params->maxnpeaks), "image2xy Max num of peaks total", &status);
238 fits_write_key(ofptr, TINT, "MAXSIZE", &(params->maxsize), "image2xy Max size for extended objects", &status);
239 fits_write_key(ofptr, TINT, "HALFBOX", &(params->halfbox), "image2xy Half-size for sliding sky window", &status);
240
241
242 fits_write_comment(ofptr,
243 "The X and Y points are specified assuming 1,1 is "
244 "the center of the leftmost bottom pixel of the "
245 "image in accordance with the FITS standard.", &status);
246 CFITS_CHECK("Failed to write comments");
247
248 simplexy_free_contents(params);
249 }
250
251 // Put in the optional NEXTEND keywoard
252 fits_movabs_hdu(ofptr, 1, &hdutype, &status);
253 assert(hdutype == IMAGE_HDU);
254 fits_write_key(ofptr, TINT, "NEXTEND", &nimgs, "Number of extensions", &status);
255 if (status == END_OF_FILE)
256 status = 0; /* Reset after normal error */
257 CFITS_CHECK("Failed to write NEXTEND");
258
259 fits_close_file(fptr, &status);
260 CFITS_CHECK("Failed to close FITS input file");
261 fptr = NULL;
262
263 fits_close_file(ofptr, &status);
264 CFITS_CHECK("Failed to close FITS output file");
265
266 // for valgrind
267 simplexy_clean_cache();
268
269 return 0;
270
271 bailout:
272 if (fptr)
273 fits_close_file(fptr, &status);
274 if (ofptr)
275 fits_close_file(ofptr, &status);
276 return -1;
277 }
278