public software.sextractor

[/] [trunk/] [src/] [makeit.c] - Blame information for rev 6

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 2 bertin
/*
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 5 bertin
*       Last modify:    12/01/2006
13 2 bertin
*
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
 
42
/******************************** makeit *************************************/
43
/*
44
Manage the whole stuff.
45
*/
46
void    makeit()
47
 
48
  {
49
   checkstruct          *check;
50
   picstruct            *dfield, *field,*pffield[MAXFLAG], *wfield,*dwfield;
51
   catstruct            *imacat;
52
   tabstruct            *imatab;
53
   static time_t        thetime1, thetime2;
54
   int                  i, nok, ntab, next;
55
 
56
/* Initialize globals variables */
57
  initglob();
58
 
59
  NFPRINTF(OUTPUT, "Setting catalog parameters");
60
  readcatparams(prefs.param_name);
61
  useprefs();                   /* update things accor. to prefs parameters */
62
 
63
  if (prefs.psf_flag)
64
    {
65
    NFPRINTF(OUTPUT, "Reading PSF information");
66 5 bertin
    thepsf = psf_load(prefs.psf_name[0]);
67
    if (prefs.dpsf_flag)
68
      ppsf = psf_load(prefs.psf_name[1]);
69 2 bertin
 /*-- Need to check things up because of PSF context parameters */
70
    updateparamflags();
71
    useprefs();
72
    }
73
 
74
  if (prefs.filter_flag)
75
    {
76
    NFPRINTF(OUTPUT, "Reading detection filter");
77
    getfilter(prefs.filter_name);       /* get the detection filter */
78
    }
79
 
80
  if (FLAG(obj2.sprob))
81
    {
82
    NFPRINTF(OUTPUT, "Initializing Neural Network");
83
    neurinit();
84
    NFPRINTF(OUTPUT, "Reading Neural Network Weights");
85
    getnnw();
86
    }
87
 
88
  if (prefs.somfit_flag)
89
    {
90
     int        margin;
91
 
92
    thesom = som_load(prefs.som_name);
93
    if ((margin=(thesom->inputsize[1]+1)/2) > prefs.cleanmargin)
94
      prefs.cleanmargin = margin;
95
    if (prefs.somfit_vectorsize>thesom->neurdim)
96
      {
97
      prefs.somfit_vectorsize = thesom->neurdim;
98
      sprintf(gstr,"%d", prefs.somfit_vectorsize);
99
      warning("Dimensionality of the SOM-fit vector limited to ", gstr);
100
      }
101
    }
102
 
103
/* Prepare growth-curve buffer */
104
  if (prefs.growth_flag)
105
    initgrowth();
106
 
107
/* Compute the number of valid input extensions */
108
  if (!(imacat = read_cat(prefs.image_name[0])))
109
    error(EXIT_FAILURE, "*Error*: cannot open ", prefs.image_name[0]);
110
  close_cat(imacat);
111
  imatab = imacat->tab;
112
  next = 0;
113
  for (ntab = 0 ; ntab<imacat->ntab; ntab++, imatab = imatab->nexttab)
114
    {
115
/*--  Check for the next valid image extension */
116
    if ((imatab->naxis < 2)
117
        || !strncmp(imatab->xtension, "BINTABLE", 8)
118
        || !strncmp(imatab->xtension, "ASCTABLE", 8))
119
      continue;
120
    next++;
121
    }
122
  thecat.next = next;
123
 
124
/*-- Init the CHECK-images */
125
  if (prefs.check_flag)
126
    {
127
     checkenum  c;
128
 
129
    NFPRINTF(OUTPUT, "Initializing check-image(s)");
130
    for (i=0; i<prefs.ncheck_type; i++)
131
    if ((c=prefs.check_type[i]) != CHECK_NONE)
132
      {
133
      if (prefs.check[c])
134
         error(EXIT_FAILURE,"*Error*: 2 CHECK_IMAGEs cannot have the same ",
135
                        " CHECK_IMAGE_TYPE");
136
      prefs.check[c] = initcheck(prefs.check_name[i], prefs.check_type[i],
137
                        next);
138
      free(prefs.check_name[i]);
139
      }
140
    }
141
 
142
  NFPRINTF(OUTPUT, "Initializing catalog");
143
  initcat();
144
 
145
 
146
/* Go through all images */
147
  nok = -1;
148
  for (ntab = 0 ; ntab<imacat->ntab; ntab++, imatab = imatab->nexttab)
149
    {
150
/*--  Check for the next valid image extension */
151
    if ((imatab->naxis < 2)
152
        || !strncmp(imatab->xtension, "BINTABLE", 8)
153
        || !strncmp(imatab->xtension, "ASCTABLE", 8))
154
      continue;
155
    nok++;
156
 
157
/*-- Initial time measurement*/
158
    time(&thetime1);
159
    thecat.currext = nok+1;
160
 
161
    dfield = field = wfield = dwfield = NULL;
162
 
163
    if (prefs.dimage_flag)
164
      {
165
/*---- Init the Detection and Measurement-images */
166
      dfield = newfield(prefs.image_name[0], DETECT_FIELD, nok);
167
      field = newfield(prefs.image_name[1], MEASURE_FIELD, nok);
168
      if ((field->width!=dfield->width) || (field->height!=dfield->height))
169
        error(EXIT_FAILURE, "*Error*: Frames have different sizes","");
170
/*---- Prepare interpolation */
171
      if (prefs.dweight_flag && prefs.interp_type[0] == INTERP_ALL)
172
        init_interpolate(dfield, -1, -1);
173
      if (prefs.interp_type[1] == INTERP_ALL)
174
        init_interpolate(field, -1, -1);
175
      }
176
    else
177
      {
178
      field = newfield(prefs.image_name[0], DETECT_FIELD | MEASURE_FIELD, nok);
179
/*-- Prepare interpolation */
180
      if ((prefs.dweight_flag || prefs.weight_flag)
181
        && prefs.interp_type[0] == INTERP_ALL)
182
      init_interpolate(field, -1, -1);       /* 0.0 or anything else */
183
      }
184
 
185
/*-- Init the WEIGHT-images */
186
    if (prefs.dweight_flag || prefs.weight_flag)
187
      {
188
       weightenum       wtype;
189
       PIXTYPE  interpthresh;
190
 
191
      if (prefs.nweight_type>1)
192
        {
193
/*------ Double-weight-map mode */
194
        if (prefs.weight_type[1] != WEIGHT_NONE)
195
          {
196
/*-------- First: the "measurement" weights */
197
          wfield = newweight(prefs.wimage_name[1],field,prefs.weight_type[1],
198
                nok);
199
          wtype = prefs.weight_type[1];
200
          interpthresh = prefs.weight_thresh[1];
201
/*-------- Convert the interpolation threshold to variance units */
202
          weight_to_var(wfield, &interpthresh, 1);
203
          wfield->weight_thresh = interpthresh;
204
          if (prefs.interp_type[1] != INTERP_NONE)
205
            init_interpolate(wfield,
206
                prefs.interp_xtimeout[1], prefs.interp_ytimeout[1]);
207
          }
208
/*------ The "detection" weights */
209
        if (prefs.weight_type[0] != WEIGHT_NONE)
210
          {
211
          interpthresh = prefs.weight_thresh[0];
212
          if (prefs.weight_type[0] == WEIGHT_FROMINTERP)
213
            {
214
            dwfield=newweight(prefs.wimage_name[0],wfield,prefs.weight_type[0],
215
                nok);
216
            weight_to_var(wfield, &interpthresh, 1);
217
            }
218
          else
219
            {
220
            dwfield = newweight(prefs.wimage_name[0], dfield?dfield:field,
221
                prefs.weight_type[0], nok);
222
            weight_to_var(dwfield, &interpthresh, 1);
223
            }
224
          dwfield->weight_thresh = interpthresh;
225
          if (prefs.interp_type[0] != INTERP_NONE)
226
            init_interpolate(dwfield,
227
                prefs.interp_xtimeout[0], prefs.interp_ytimeout[0]);
228
          }
229
        }
230
      else
231
        {
232
/*------ Single-weight-map mode */
233
        wfield = newweight(prefs.wimage_name[0], dfield?dfield:field,
234
                        prefs.weight_type[0], nok);
235
        wtype = prefs.weight_type[0];
236
        interpthresh = prefs.weight_thresh[0];
237
/*------ Convert the interpolation threshold to variance units */
238
        weight_to_var(wfield, &interpthresh, 1);
239
        wfield->weight_thresh = interpthresh;
240
        if (prefs.interp_type[0] != INTERP_NONE)
241
          init_interpolate(wfield,
242
                prefs.interp_xtimeout[0], prefs.interp_ytimeout[0]);
243
        }
244
      }
245
 
246
/*-- Init the FLAG-images */
247
    for (i=0; i<prefs.nimaflag; i++)
248
      {
249
      pffield[i] = newfield(prefs.fimage_name[i], FLAG_FIELD, nok);
250
      if ((pffield[i]->width!=field->width)
251
        || (pffield[i]->height!=field->height))
252
        error(EXIT_FAILURE,
253
        "*Error*: Incompatible FLAG-map size in ", prefs.fimage_name[i]);
254
      }
255
 
256
/*-- Compute background maps for `standard' fields */
257
    QPRINTF(OUTPUT, dfield? "Measurement image:"
258
                        : "Detection+Measurement image: ");
259
    makeback(field, wfield);
260
    QPRINTF(OUTPUT, (dfield || (dwfield&&dwfield->flags^INTERP_FIELD))? "(M)   "
261
                "Background: %-10g RMS: %-10g / Threshold: %-10g \n"
262
                : "(M+D) "
263
                "Background: %-10g RMS: %-10g / Threshold: %-10g \n",
264
        field->backmean, field->backsig, (field->flags & DETECT_FIELD)?
265
        field->dthresh: field->thresh);
266
    if (dfield)
267
    {
268
      QPRINTF(OUTPUT, "Detection image: ");
269
      makeback(dfield, dwfield? dwfield
270
                        : (prefs.weight_type[0] == WEIGHT_NONE?NULL:wfield));
271
      QPRINTF(OUTPUT, "(D)   "
272
                "Background: %-10g RMS: %-10g / Threshold: %-10g \n",
273
        dfield->backmean, dfield->backsig, dfield->dthresh);
274
      }
275
    else if (dwfield && dwfield->flags^INTERP_FIELD)
276
      {
277
      makeback(field, dwfield);
278
      QPRINTF(OUTPUT, "(D)   "
279
                "Background: %-10g RMS: %-10g / Threshold: %-10g \n",
280
        field->backmean, field->backsig, field->dthresh);
281
      }
282
 
283
/*-- For interpolated weight-maps, copy the background structure */
284
    if (dwfield && dwfield->flags&(INTERP_FIELD|BACKRMS_FIELD))
285
      copyback(dwfield->reffield, dwfield);
286
    if (wfield && wfield->flags&(INTERP_FIELD|BACKRMS_FIELD))
287
      copyback(wfield->reffield, wfield);
288
 
289
/*-- Prepare learn and/or associations */
290
    if (prefs.assoc_flag)
291
      init_assoc(field);                  /* initialize assoc tasks */
292
 
293
/*-- Update the CHECK-images */
294
    if (prefs.check_flag)
295
      for (i=0; i<MAXCHECK; i++)
296
        if ((check=prefs.check[i]))
297
          reinitcheck(field, check);
298
 
299
/*-- Initialize PSF contexts and workspace */
300
    if (prefs.psf_flag)
301
      {
302
      psf_readcontext(thepsf, field);
303
      psf_init(thepsf);
304 5 bertin
      if (prefs.dpsf_flag)
305
        {
306
        psf_readcontext(thepsf, dfield);
307
        psf_init(thepsf); /*?*/
308
        }
309 2 bertin
      }
310
 
311
/*-- Copy field structures to static ones (for catalog info) */
312
    if (dfield)
313
      {
314
      thefield1 = *field;
315
      thefield2 = *dfield;
316
      }
317
    else
318
      thefield1 = thefield2 = *field;
319
 
320
    if (wfield)
321
      {
322
      thewfield1 = *wfield;
323
      thewfield2 = dwfield? *dwfield: *wfield;
324
      }
325
    else if (dwfield)
326
      thewfield2 = *dwfield;
327
 
328
    reinitcat(field);
329
 
330
/*-- Start the extraction pipeline */
331
    NFPRINTF(OUTPUT, "Scanning image");
332
    scanimage(field, dfield, pffield, prefs.nimaflag, wfield, dwfield);
333
 
334
/*-- Finish the current CHECK-image processing */
335
    if (prefs.check_flag)
336
      for (i=0; i<MAXCHECK; i++)
337
        if ((check=prefs.check[i]))
338
          reendcheck(field, check);
339
 
340
/*-- Final time measurements*/
341
    if (time(&thetime2)!=-1)
342
      {
343
      if (!strftime(thecat.ext_date, 12, "%d/%m/%Y", localtime(&thetime2)))
344
        error(EXIT_FAILURE, "*Internal Error*: Date string too long ","");
345
      if (!strftime(thecat.ext_time, 10, "%H:%M:%S", localtime(&thetime2)))
346
        error(EXIT_FAILURE, "*Internal Error*: Time/date string too long ","");
347
      thecat.ext_elapsed = difftime(thetime2, thetime1);
348
      }
349
 
350
    reendcat();
351
 
352
/*-- Close ASSOC routines */
353
    end_assoc(field);
354
 
355
    for (i=0; i<prefs.nimaflag; i++)
356
      endfield(pffield[i]);
357
    endfield(field);
358
    if (dfield)
359
      endfield(dfield);
360
    if (wfield)
361
      endfield(wfield);
362
    if (dwfield)
363
      endfield(dwfield);
364
 
365 6 bertin
    QPRINTF(OUTPUT, "Objects: detected %-8d / sextracted %-8d               \n",
366 2 bertin
        thecat.ndetect, thecat.ntotal);
367
    }
368
 
369
  if (nok<0)
370
    error(EXIT_FAILURE, "Not enough valid FITS image extensions in ",
371
        prefs.image_name[0]);
372
  free_cat(&imacat, 1);
373
 
374
  NFPRINTF(OUTPUT, "Closing files");
375
 
376
  endcat();
377
 
378
/* End CHECK-image processing */
379
  if (prefs.check_flag)
380
    for (i=0; i<MAXCHECK; i++)
381
      {
382
      if ((check=prefs.check[i]))
383
        endcheck(check);
384
      prefs.check[i] = NULL;
385
      }
386
 
387
  if (prefs.filter_flag)
388
    endfilter();
389
 
390
  if (prefs.somfit_flag)
391
    som_end(thesom);
392
 
393
  if (prefs.growth_flag)
394
    endgrowth();
395
 
396
  if (prefs.psf_flag)
397 5 bertin
    psf_end(thepsf,thepsfit); /*?*/
398 2 bertin
 
399 5 bertin
  if (prefs.dpsf_flag)
400
    psf_end(ppsf,ppsfit);
401
 
402 2 bertin
  if (FLAG(obj2.sprob))
403
    neurclose();
404
 
405
  return;
406
  }
407
 
408
 
409
/******************************** initglob ***********************************/
410
/*
411
Initialize a few global variables
412
*/
413
void    initglob()
414
  {
415
   int  i;
416
 
417
  for (i=0; i<37; i++)
418
    {
419
    ctg[i] = cos(i*PI/18);
420
    stg[i] = sin(i*PI/18);
421
    }
422
 
423
 
424
  return;
425
  }
426
 
427
/*
428
int matherr(struct exception *x)
429
{
430
printf("***MATH ERROR***: %d %s %f %f\n",
431
x->type, x->name, x->arg1, x->retval);
432
return (0);
433
}
434
 
435
*/