1 /*
2 makeit.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN, IAP & Leiden observatory
9 *
10 * Contents: main program.
11 *
12 * Last modify: 14/07/2006
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <time.h>
26
27 #include "define.h"
28 #include "globals.h"
29 #include "prefs.h"
30 #include "fits/fitscat.h"
31 #include "assoc.h"
32 #include "back.h"
33 #include "check.h"
34 #include "field.h"
35 #include "filter.h"
36 #include "growth.h"
37 #include "interpolate.h"
38 #include "psf.h"
39 #include "som.h"
40 #include "weight.h"
41 #include "xml.h"
42
43 time_t thetimet, thetimet2;
44
45 /******************************** makeit *************************************/
46 /*
47 Manage the whole stuff.
48 */
makeit()49 void makeit()
50
51 {
52 checkstruct *check;
53 picstruct *dfield, *field,*pffield[MAXFLAG], *wfield,*dwfield;
54 catstruct *imacat;
55 tabstruct *imatab;
56 static time_t thetime1, thetime2;
57 struct tm *tm;
58 int i, nok, ntab, next;
59
60 /* Install error logging */
61 error_installfunc(write_error);
62
63 /* Processing start date and time */
64 thetimet = time(NULL);
65 tm = localtime(&thetimet);
66 sprintf(prefs.sdate_start,"%04d-%02d-%02d",
67 tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
68 sprintf(prefs.stime_start,"%02d:%02d:%02d",
69 tm->tm_hour, tm->tm_min, tm->tm_sec);
70
71 NFPRINTF(OUTPUT, "");
72 QPRINTF(OUTPUT, "----- %s %s started on %s at %s with %d thread%s\n\n",
73 BANNER,
74 MYVERSION,
75 prefs.sdate_start,
76 prefs.stime_start,
77 prefs.nthreads,
78 prefs.nthreads>1? "s":"");
79
80 /* Initialize globals variables */
81 initglob();
82
83 NFPRINTF(OUTPUT, "Setting catalog parameters");
84 readcatparams(prefs.param_name);
85 useprefs(); /* update things accor. to prefs parameters */
86
87 if (prefs.psf_flag)
88 {
89 NFPRINTF(OUTPUT, "Reading PSF information");
90 thepsf = psf_load(prefs.psf_name[0]);
91 if (prefs.dpsf_flag)
92 ppsf = psf_load(prefs.psf_name[1]);
93 /*-- Need to check things up because of PSF context parameters */
94 updateparamflags();
95 useprefs();
96 }
97
98 if (prefs.filter_flag)
99 {
100 NFPRINTF(OUTPUT, "Reading detection filter");
101 getfilter(prefs.filter_name); /* get the detection filter */
102 }
103
104 if (FLAG(obj2.sprob))
105 {
106 NFPRINTF(OUTPUT, "Initializing Neural Network");
107 neurinit();
108 NFPRINTF(OUTPUT, "Reading Neural Network Weights");
109 getnnw();
110 }
111
112 if (prefs.somfit_flag)
113 {
114 int margin;
115
116 thesom = som_load(prefs.som_name);
117 if ((margin=(thesom->inputsize[1]+1)/2) > prefs.cleanmargin)
118 prefs.cleanmargin = margin;
119 if (prefs.somfit_vectorsize>thesom->neurdim)
120 {
121 prefs.somfit_vectorsize = thesom->neurdim;
122 sprintf(gstr,"%d", prefs.somfit_vectorsize);
123 warning("Dimensionality of the SOM-fit vector limited to ", gstr);
124 }
125 }
126
127 /* Prepare growth-curve buffer */
128 if (prefs.growth_flag)
129 initgrowth();
130
131 /* Compute the number of valid input extensions */
132 if (!(imacat = read_cat(prefs.image_name[0])))
133 error(EXIT_FAILURE, "*Error*: cannot open ", prefs.image_name[0]);
134 close_cat(imacat);
135 imatab = imacat->tab;
136 next = 0;
137 for (ntab = 0 ; ntab<imacat->ntab; ntab++, imatab = imatab->nexttab)
138 {
139 /*-- Check for the next valid image extension */
140 if ((imatab->naxis < 2)
141 || !strncmp(imatab->xtension, "BINTABLE", 8)
142 || !strncmp(imatab->xtension, "ASCTABLE", 8))
143 continue;
144 next++;
145 }
146 thecat.next = next;
147
148 /*-- Init the CHECK-images */
149 if (prefs.check_flag)
150 {
151 checkenum c;
152
153 NFPRINTF(OUTPUT, "Initializing check-image(s)");
154 for (i=0; i<prefs.ncheck_type; i++)
155 if ((c=prefs.check_type[i]) != CHECK_NONE)
156 {
157 if (prefs.check[c])
158 error(EXIT_FAILURE,"*Error*: 2 CHECK_IMAGEs cannot have the same ",
159 " CHECK_IMAGE_TYPE");
160 prefs.check[c] = initcheck(prefs.check_name[i], prefs.check_type[i],
161 next);
162 free(prefs.check_name[i]);
163 }
164 }
165
166 NFPRINTF(OUTPUT, "Initializing catalog");
167 initcat();
168
169 /* Initialize XML data */
170 if (prefs.xml_flag || prefs.cat_type==ASCII_VO)
171 init_xml(next);
172
173 /* Go through all images */
174 nok = -1;
175 for (ntab = 0 ; ntab<imacat->ntab; ntab++, imatab = imatab->nexttab)
176 {
177 /*-- Check for the next valid image extension */
178 if ((imatab->naxis < 2)
179 || !strncmp(imatab->xtension, "BINTABLE", 8)
180 || !strncmp(imatab->xtension, "ASCTABLE", 8))
181 continue;
182 nok++;
183
184 /*-- Initial time measurement*/
185 time(&thetime1);
186 thecat.currext = nok+1;
187
188 dfield = field = wfield = dwfield = NULL;
189
190 if (prefs.dimage_flag)
191 {
192 /*---- Init the Detection and Measurement-images */
193 dfield = newfield(prefs.image_name[0], DETECT_FIELD, nok);
194 field = newfield(prefs.image_name[1], MEASURE_FIELD, nok);
195 if ((field->width!=dfield->width) || (field->height!=dfield->height))
196 error(EXIT_FAILURE, "*Error*: Frames have different sizes","");
197 /*---- Prepare interpolation */
198 if (prefs.dweight_flag && prefs.interp_type[0] == INTERP_ALL)
199 init_interpolate(dfield, -1, -1);
200 if (prefs.interp_type[1] == INTERP_ALL)
201 init_interpolate(field, -1, -1);
202 }
203 else
204 {
205 field = newfield(prefs.image_name[0], DETECT_FIELD | MEASURE_FIELD, nok);
206 /*-- Prepare interpolation */
207 if ((prefs.dweight_flag || prefs.weight_flag)
208 && prefs.interp_type[0] == INTERP_ALL)
209 init_interpolate(field, -1, -1); /* 0.0 or anything else */
210 }
211
212 /*-- Init the WEIGHT-images */
213 if (prefs.dweight_flag || prefs.weight_flag)
214 {
215 weightenum wtype;
216 PIXTYPE interpthresh;
217
218 if (prefs.nweight_type>1)
219 {
220 /*------ Double-weight-map mode */
221 if (prefs.weight_type[1] != WEIGHT_NONE)
222 {
223 /*-------- First: the "measurement" weights */
224 wfield = newweight(prefs.wimage_name[1],field,prefs.weight_type[1],
225 nok);
226 wtype = prefs.weight_type[1];
227 interpthresh = prefs.weight_thresh[1];
228 /*-------- Convert the interpolation threshold to variance units */
229 weight_to_var(wfield, &interpthresh, 1);
230 wfield->weight_thresh = interpthresh;
231 if (prefs.interp_type[1] != INTERP_NONE)
232 init_interpolate(wfield,
233 prefs.interp_xtimeout[1], prefs.interp_ytimeout[1]);
234 }
235 /*------ The "detection" weights */
236 if (prefs.weight_type[0] != WEIGHT_NONE)
237 {
238 interpthresh = prefs.weight_thresh[0];
239 if (prefs.weight_type[0] == WEIGHT_FROMINTERP)
240 {
241 dwfield=newweight(prefs.wimage_name[0],wfield,prefs.weight_type[0],
242 nok);
243 weight_to_var(wfield, &interpthresh, 1);
244 }
245 else
246 {
247 dwfield = newweight(prefs.wimage_name[0], dfield?dfield:field,
248 prefs.weight_type[0], nok);
249 weight_to_var(dwfield, &interpthresh, 1);
250 }
251 dwfield->weight_thresh = interpthresh;
252 if (prefs.interp_type[0] != INTERP_NONE)
253 init_interpolate(dwfield,
254 prefs.interp_xtimeout[0], prefs.interp_ytimeout[0]);
255 }
256 }
257 else
258 {
259 /*------ Single-weight-map mode */
260 wfield = newweight(prefs.wimage_name[0], dfield?dfield:field,
261 prefs.weight_type[0], nok);
262 wtype = prefs.weight_type[0];
263 interpthresh = prefs.weight_thresh[0];
264 /*------ Convert the interpolation threshold to variance units */
265 weight_to_var(wfield, &interpthresh, 1);
266 wfield->weight_thresh = interpthresh;
267 if (prefs.interp_type[0] != INTERP_NONE)
268 init_interpolate(wfield,
269 prefs.interp_xtimeout[0], prefs.interp_ytimeout[0]);
270 }
271 }
272
273 /*-- Init the FLAG-images */
274 for (i=0; i<prefs.nimaflag; i++)
275 {
276 pffield[i] = newfield(prefs.fimage_name[i], FLAG_FIELD, nok);
277 if ((pffield[i]->width!=field->width)
278 || (pffield[i]->height!=field->height))
279 error(EXIT_FAILURE,
280 "*Error*: Incompatible FLAG-map size in ", prefs.fimage_name[i]);
281 }
282
283 /*-- Compute background maps for `standard' fields */
284 QPRINTF(OUTPUT, dfield? "Measurement image:"
285 : "Detection+Measurement image: ");
286 makeback(field, wfield);
287 QPRINTF(OUTPUT, (dfield || (dwfield&&dwfield->flags^INTERP_FIELD))? "(M) "
288 "Background: %-10g RMS: %-10g / Threshold: %-10g \n"
289 : "(M+D) "
290 "Background: %-10g RMS: %-10g / Threshold: %-10g \n",
291 field->backmean, field->backsig, (field->flags & DETECT_FIELD)?
292 field->dthresh: field->thresh);
293 if (dfield)
294 {
295 QPRINTF(OUTPUT, "Detection image: ");
296 makeback(dfield, dwfield? dwfield
297 : (prefs.weight_type[0] == WEIGHT_NONE?NULL:wfield));
298 QPRINTF(OUTPUT, "(D) "
299 "Background: %-10g RMS: %-10g / Threshold: %-10g \n",
300 dfield->backmean, dfield->backsig, dfield->dthresh);
301 }
302 else if (dwfield && dwfield->flags^INTERP_FIELD)
303 {
304 makeback(field, dwfield);
305 QPRINTF(OUTPUT, "(D) "
306 "Background: %-10g RMS: %-10g / Threshold: %-10g \n",
307 field->backmean, field->backsig, field->dthresh);
308 }
309
310 /*-- For interpolated weight-maps, copy the background structure */
311 if (dwfield && dwfield->flags&(INTERP_FIELD|BACKRMS_FIELD))
312 copyback(dwfield->reffield, dwfield);
313 if (wfield && wfield->flags&(INTERP_FIELD|BACKRMS_FIELD))
314 copyback(wfield->reffield, wfield);
315
316 /*-- Prepare learn and/or associations */
317 if (prefs.assoc_flag)
318 init_assoc(field); /* initialize assoc tasks */
319
320 /*-- Update the CHECK-images */
321 if (prefs.check_flag)
322 for (i=0; i<MAXCHECK; i++)
323 if ((check=prefs.check[i]))
324 reinitcheck(field, check);
325
326 /*-- Initialize PSF contexts and workspace */
327 if (prefs.psf_flag)
328 {
329 psf_readcontext(thepsf, field);
330 psf_init(thepsf);
331 if (prefs.dpsf_flag)
332 {
333 psf_readcontext(thepsf, dfield);
334 psf_init(thepsf); /*?*/
335 }
336 }
337
338 /*-- Copy field structures to static ones (for catalog info) */
339 if (dfield)
340 {
341 thefield1 = *field;
342 thefield2 = *dfield;
343 }
344 else
345 thefield1 = thefield2 = *field;
346
347 if (wfield)
348 {
349 thewfield1 = *wfield;
350 thewfield2 = dwfield? *dwfield: *wfield;
351 }
352 else if (dwfield)
353 thewfield2 = *dwfield;
354
355 reinitcat(field);
356
357 /*-- Start the extraction pipeline */
358 NFPRINTF(OUTPUT, "Scanning image");
359 scanimage(field, dfield, pffield, prefs.nimaflag, wfield, dwfield);
360
361 /*-- Finish the current CHECK-image processing */
362 if (prefs.check_flag)
363 for (i=0; i<MAXCHECK; i++)
364 if ((check=prefs.check[i]))
365 reendcheck(field, check);
366
367 /*-- Final time measurements*/
368 if (time(&thetime2)!=-1)
369 {
370 if (!strftime(thecat.ext_date, 12, "%d/%m/%Y", localtime(&thetime2)))
371 error(EXIT_FAILURE, "*Internal Error*: Date string too long ","");
372 if (!strftime(thecat.ext_time, 10, "%H:%M:%S", localtime(&thetime2)))
373 error(EXIT_FAILURE, "*Internal Error*: Time/date string too long ","");
374 thecat.ext_elapsed = difftime(thetime2, thetime1);
375 }
376
377 reendcat();
378
379 /* Update XML data */
380 if (prefs.xml_flag || prefs.cat_type==ASCII_VO)
381 update_xml(&thecat, dfield? dfield:field, field,
382 dwfield? dwfield:wfield, wfield);
383
384
385 /*-- Close ASSOC routines */
386 end_assoc(field);
387
388 for (i=0; i<prefs.nimaflag; i++)
389 endfield(pffield[i]);
390 endfield(field);
391 if (dfield)
392 endfield(dfield);
393 if (wfield)
394 endfield(wfield);
395 if (dwfield)
396 endfield(dwfield);
397
398 QPRINTF(OUTPUT, "Objects: detected %-8d / sextracted %-8d \n",
399 thecat.ndetect, thecat.ntotal);
400 }
401
402 if (nok<0)
403 error(EXIT_FAILURE, "Not enough valid FITS image extensions in ",
404 prefs.image_name[0]);
405 free_cat(&imacat, 1);
406
407 NFPRINTF(OUTPUT, "Closing files");
408
409 /* End CHECK-image processing */
410 if (prefs.check_flag)
411 for (i=0; i<MAXCHECK; i++)
412 {
413 if ((check=prefs.check[i]))
414 endcheck(check);
415 prefs.check[i] = NULL;
416 }
417
418 if (prefs.filter_flag)
419 endfilter();
420
421 if (prefs.somfit_flag)
422 som_end(thesom);
423
424 if (prefs.growth_flag)
425 endgrowth();
426
427 if (prefs.psf_flag)
428 psf_end(thepsf,thepsfit); /*?*/
429
430 if (prefs.dpsf_flag)
431 psf_end(ppsf,ppsfit);
432
433 if (FLAG(obj2.sprob))
434 neurclose();
435
436 /* Processing end date and time */
437 thetimet2 = time(NULL);
438 tm = localtime(&thetimet2);
439 sprintf(prefs.sdate_end,"%04d-%02d-%02d",
440 tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
441 sprintf(prefs.stime_end,"%02d:%02d:%02d",
442 tm->tm_hour, tm->tm_min, tm->tm_sec);
443 prefs.time_diff = difftime(thetimet2, thetimet);
444
445 /* Write XML */
446 if (prefs.xml_flag)
447 write_xml(prefs.xml_name);
448
449 endcat((char *)NULL);
450
451 if (prefs.xml_flag || prefs.cat_type==ASCII_VO)
452 end_xml();
453
454 return;
455 }
456
457
458 /******************************** initglob ***********************************/
459 /*
460 Initialize a few global variables
461 */
initglob()462 void initglob()
463 {
464 int i;
465
466 for (i=0; i<37; i++)
467 {
468 ctg[i] = cos(i*PI/18);
469 stg[i] = sin(i*PI/18);
470 }
471
472
473 return;
474 }
475
476
477 /****** write_error ********************************************************
478 PROTO int write_error(char *msg1, char *msg2)
479 PURPOSE Manage files in case of a catched error
480 INPUT a character string,
481 another character string
482 OUTPUT RETURN_OK if everything went fine, RETURN_ERROR otherwise.
483 NOTES -.
484 AUTHOR E. Bertin (IAP)
485 VERSION 14/07/2006
486 ***/
write_error(char * msg1,char * msg2)487 void write_error(char *msg1, char *msg2)
488 {
489 char error[MAXCHAR];
490
491 sprintf(error, "%s%s", msg1,msg2);
492 if (prefs.xml_flag)
493 write_xmlerror(prefs.xml_name, error);
494
495 /* Also close existing catalog */
496 endcat(error);
497
498 end_xml();
499
500 return;
501 }
502
503