public software.psfex

[/] [branches/] [rhl/] [src/] [sample.c] - Blame information for rev 191

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 157 bertin
/*
2
*                               sample.c
3 2 bertin
*
4 157 bertin
* Read and filter input samples from catalogues.
5 2 bertin
*
6 157 bertin
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 2 bertin
*
8 157 bertin
*       This file part of:      PSFEx
9 2 bertin
*
10 184 bertin
*       Copyright:              (C) 1997-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
11 2 bertin
*
12 157 bertin
*       License:                GNU General Public License
13
*
14
*       PSFEx is free software: you can redistribute it and/or modify
15
*       it under the terms of the GNU General Public License as published by
16
*       the Free Software Foundation, either version 3 of the License, or
17
*       (at your option) any later version.
18
*       PSFEx is distributed in the hope that it will be useful,
19
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
20
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21
*       GNU General Public License for more details.
22
*       You should have received a copy of the GNU General Public License
23
*       along with PSFEx.  If not, see <http://www.gnu.org/licenses/>.
24
*
25 184 bertin
*       Last modified:          19/07/2012
26 157 bertin
*
27
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
28 2 bertin
 
29
#ifdef HAVE_CONFIG_H
30
#include        "config.h"
31
#endif
32
 
33 189 rhl
#include <assert.h>
34 2 bertin
#include <math.h>
35
#include <stdio.h>
36
#include <stdlib.h>
37
#include <string.h>
38
 
39
#include "define.h"
40
#include "types.h"
41
#include "globals.h"
42
#include "fits/fitscat.h"
43
#include "prefs.h"
44 72 bertin
#include "context.h"
45 159 bertin
#include "misc.h"
46 2 bertin
#include "sample.h"
47
#include "vignet.h"
48
 
49 79 bertin
static float    compute_fwhmrange(float *fwhm, int nfwhm, float maxvar,
50
                float minin, float maxin, float *minout, float *maxout);
51
 
52 170 bertin
/****** load_samples *********************************************************
53
PROTO   setstruct *load_samples(char **filename, int catindex, int ncat,
54
                int ext, int next, contextstruct *context)
55
PURPOSE Examine and load point source data.
56
INPUT   Array of catalog filenames,
57
        catalog index,
58
        current extension,
59
        number of extensions,
60
        pointer to the context.
61
OUTPUT  Pointer to a set containing samples that match acceptance criteria.
62
NOTES   -.
63
AUTHOR  E. Bertin (IAP)
64 184 bertin
VERSION 19/07/2012
65 2 bertin
*/
66 119 bertin
setstruct *load_samples(char **filename, int catindex, int ncat, int ext,
67
                        int next, contextstruct *context)
68 2 bertin
  {
69
   setstruct            *set;
70
   catstruct            *cat;
71
   tabstruct            *tab;
72 170 bertin
   keystruct            *fkey, *key;
73
   char                 keynames[][32]={"ELONGATION", "FLAGS", "FLUX_RADIUS",
74
                                "SNR_WIN", ""};
75 56 bertin
   char                 str[MAXCHAR];
76 170 bertin
   char                 **pkeynames,
77
                        *head;
78 79 bertin
   float                *fwhmmin,*fwhmmax,*fwhmmode,
79 171 bertin
                        *fwhm,*fwhmt, *elong, *hl, *snr,
80 79 bertin
                        backnoise, minsn, maxelong, min,max, mode,  fval;
81 184 bertin
   unsigned int         *imaflags;
82
   unsigned short       *flags, *wflags;
83 79 bertin
   int                  *fwhmindex,
84 170 bertin
                        e,i,j,n, icat, nobj,nobjmax, ldflag, ext2, nkeys;
85 2 bertin
 
86 112 bertin
//  NFPRINTF(OUTPUT,"Loading samples...");
87 2 bertin
  minsn = (float)prefs.minsn;
88 86 bertin
  maxelong = (float)(prefs.maxellip < 1.0?
89
        (prefs.maxellip + 1.0)/(1.0 - prefs.maxellip)
90
        : 100.0);
91 79 bertin
  min = (float)prefs.fwhmrange[0];
92
  max = (float)prefs.fwhmrange[1];
93 2 bertin
  fwhm = NULL;  /* To avoid gcc -Wall warnings */
94 82 bertin
/* Allocate memory */
95
  QMALLOC(fwhmmin, float, ncat);
96
  QMALLOC(fwhmmax, float, ncat);
97
  QMALLOC(fwhmmode, float, ncat);
98 2 bertin
 
99
  if (prefs.autoselect_flag)
100
    {
101
/*-- Allocate memory */
102 82 bertin
    nobjmax = LSAMPLE_DEFSIZE;
103
    QMALLOC(fwhm, float, nobjmax);
104 79 bertin
    QMALLOC(fwhmindex, int, ncat+1);
105
    fwhmindex[0] = nobj = 0;
106 2 bertin
    fwhmt=fwhm;
107
 
108
/*-- Initialize string array */
109 170 bertin
    for (i=0; (*keynames[i]); i++);
110
    nkeys = i;
111
    QMALLOC(pkeynames, char *, nkeys);
112
    for (i=0; i<nkeys; i++)
113 2 bertin
      pkeynames[i] = keynames[i];
114
 
115
/*-- Try to estimate the most appropriate Half-light Radius range */
116
/*-- Get the Half-light radii */
117
    nobj = 0;
118
    for (i=0; i<ncat; i++)
119
      {
120
      sprintf(str,"Examining Catalog #%d", i+1);
121 112 bertin
//      NFPRINTF(OUTPUT, str);
122 2 bertin
/*---- Read input catalog */
123 119 bertin
      icat = catindex + i;
124
      if (!(cat = read_cat(filename[icat])))
125
        error(EXIT_FAILURE, "*Error*: No such catalog: ", filename[icat]);
126 2 bertin
 
127 170 bertin
/*---- Load the objects */
128 108 bertin
      e = 0;
129 183 bertin
      e = ext;
130 2 bertin
      tab = cat->tab;
131
      for (j=cat->ntab; j--; tab=tab->nexttab)
132
        if (!strcmp("LDAC_OBJECTS", tab->extname)
133
                ||  !strcmp("OBJECTS", tab->extname))
134 108 bertin
          {
135 114 bertin
          if (ext != ALL_EXTENSIONS)
136
            {
137 183 bertin
            if (e>0)
138 115 bertin
              {
139 183 bertin
              e--;
140 114 bertin
              continue;
141 115 bertin
              }
142 183 bertin
            else if (e<0)
143 114 bertin
              break;
144 183 bertin
            e--;
145 114 bertin
            }
146 170 bertin
 
147
/*-------- Read the data */
148
 
149
          read_keys(tab, pkeynames, NULL, nkeys, NULL);
150
 
151 184 bertin
          if ((key = name_to_key(tab, "ELONGATION")))
152
            elong = (float *)key->ptr;
153
          else
154
            {
155
            warning("ELONGATION parameter not found in catalog ",
156 119 bertin
                        filename[icat]);
157 184 bertin
            elong = NULL;
158
            }
159
          if ((key = name_to_key(tab, "FLAGS")))
160
            flags = (unsigned short *)key->ptr;
161
          else
162
            {
163
            warning("FLAGS parameter not found in catalog ", filename[icat]);
164
            flags = NULL;
165
            }
166
          if ((key = name_to_key(tab, "FLAGS_WEIGHT")))
167
            wflags = (unsigned short *)key->ptr;
168
          else
169
            wflags = NULL;
170
          if ((key = name_to_key(tab, "IMAFLAGS_ISO")))
171
            imaflags = (unsigned int *)key->ptr;
172
          else
173
            imaflags = NULL;
174 170 bertin
          if (!(key = name_to_key(tab, "FLUX_RADIUS")))
175
              {
176
              sprintf(str, "FLUS_RADIUS not found in catalog %s",
177
                        filename[icat]);
178
              error(EXIT_FAILURE, "*Error*: ", str);
179
              }
180
          hl = (float *)key->ptr;
181
          if (!(key = name_to_key(tab, "SNR_WIN")))
182
              {
183
              sprintf(str, "SNR_WIN not found in catalog %s",
184
                        filename[icat]);
185
              error(EXIT_FAILURE, "*Error*: ", str);
186
              }
187
          snr = (float *)key->ptr;
188 2 bertin
 
189 170 bertin
          for (n=tab->naxisn[1]; n--; hl++, snr++, flags++, elong++)
190 108 bertin
            {
191 170 bertin
            if (*snr>minsn
192 184 bertin
                && !(flags && (*flags&prefs.flag_mask))
193
                && !(wflags && (*wflags&prefs.wflag_mask))
194
                && !(imaflags && (*imaflags&prefs.imaflag_mask))
195
                && (!(elong && *elong>=maxelong))
196 79 bertin
                && (fval=2.0**hl)>=min
197
                && fval<max)
198 108 bertin
              {
199
              if (++nobj>nobjmax)
200
                {
201
                nobjmax += LSAMPLE_DEFSIZE;
202
                QREALLOC(fwhm, float, nobjmax);
203
                fwhmt=fwhm+nobj-1;
204
                }
205
              *(fwhmt++) = fval;
206
              }
207 2 bertin
            }
208
          }
209
      free_cat(&cat, 1);
210 79 bertin
      fwhmindex[i+1] = nobj;
211 2 bertin
      }
212
 
213 79 bertin
    if (prefs.var_type == VAR_NONE)
214 2 bertin
      {
215 79 bertin
      if (nobj)
216
        mode = compute_fwhmrange(fwhm, nobj, prefs.maxvar,
217
                prefs.fwhmrange[0],prefs.fwhmrange[1], &min, &max);
218
      else
219
        {
220
        warning("No source with appropriate FWHM found!!","");
221
        mode = min = max = 2.35/(1.0-1.0/INTERPFAC);
222
        }
223
      for (i=0; i<ncat; i++)
224
        {
225
        fwhmmin[i] = min;
226
        fwhmmax[i] = max;
227
        fwhmmode[i] = mode;
228
        }
229
      }
230
    else
231
      for (i=0; i<ncat; i++)
232
        {
233
        nobj = fwhmindex[i+1] - fwhmindex[i];
234
        if (nobj)
235 2 bertin
          {
236 79 bertin
          fwhmmode[i] = compute_fwhmrange(&fwhm[fwhmindex[i]],
237
                fwhmindex[i+1]-fwhmindex[i], prefs.maxvar,
238
                prefs.fwhmrange[0],prefs.fwhmrange[1], &fwhmmin[i],&fwhmmax[i]);
239 2 bertin
          }
240 79 bertin
        else
241
          {
242
          warning("No source with appropriate FWHM found!!","");
243
          fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0-1.0/INTERPFAC);
244
          }
245 2 bertin
        }
246 82 bertin
    free(fwhm);
247
    free(fwhmindex);
248 170 bertin
    free(pkeynames);
249 2 bertin
    }
250
  else
251 82 bertin
    for (i=0; i<ncat; i++)
252
      {
253
      fwhmmin[i] = (float)prefs.fwhmrange[0];
254
      fwhmmax[i] = (float)prefs.fwhmrange[1];
255
      fwhmmode[i] = (fwhmmin[i] + fwhmmax[i]) / 2.0;
256
      }
257 2 bertin
 
258 79 bertin
 
259 2 bertin
/* Load the samples */
260
  set = NULL;
261 79 bertin
  mode = BIG;
262 2 bertin
  for (i=0; i<ncat; i++)
263 79 bertin
    {
264 119 bertin
    icat = catindex + i;
265 110 bertin
    if (ext == ALL_EXTENSIONS)
266 170 bertin
      for (e=0; e<next; e++)
267 119 bertin
        set = read_samples(set, filename[icat], fwhmmin[i]/2.0, fwhmmax[i]/2.0,
268
                        e, next, icat, context, context->pc+i*context->npc);
269 110 bertin
    else
270 119 bertin
      set = read_samples(set, filename[icat], fwhmmin[i]/2.0, fwhmmax[i]/2.0,
271
                        ext, next, icat, context, context->pc+i*context->npc);
272 79 bertin
    if (fwhmmode[i]<mode)
273
      mode = fwhmmode[i];
274
    }
275 2 bertin
 
276 79 bertin
  set->fwhm = mode;
277
 
278 56 bertin
  sprintf(str, "%d samples loaded.", set->nsample);
279 112 bertin
//  NFPRINTF(OUTPUT, str);
280 79 bertin
 
281 56 bertin
  if (!set->nsample)
282
    warning("No appropriate source found!!","");
283 2 bertin
 
284 79 bertin
  free(fwhmmin);
285
  free(fwhmmax);
286
  free(fwhmmode);
287 150 bertin
/*
288
  if (set->badflags)
289
    printf("%d detections discarded with bad SExtractor flags\n",
290
        set->badflags);
291
  if (set->badsn)
292
    printf("%d detections discarded with S/N too low\n",
293
        set->badsn);
294
  if (set->badfrmin)
295
    printf("%d detections discarded with FWHM too small\n",
296
        set->badfrmin);
297
  if (set->badfrmax)
298
    printf("%d detections discarded with FWHM too large\n",
299
        set->badfrmax);
300
  if (set->badelong)
301
    printf("%d detections discarded with elongation too high\n",
302
        set->badelong);
303
  if (set->badpix)
304
    printf("%d detections discarded with too many bad pixels\n\n\n",
305
        set->badpix);
306
*/
307 2 bertin
  return set;
308
  }
309
 
310
 
311 79 bertin
/****** compute_fwhmrange *****************************************************
312
PROTO   float compute_fwhmrange(float *fwhm, int nfwhm,
313
                float minin, float maxin, float *minout, float *maxout)
314
PURPOSE Compute the FWHM range associated to a series of FWHM measurements.
315
INPUT   Pointer to an array of FWHMs,
316
        number of FWHMs,
317
        maximum allowed FWHM variation,
318
        minimum allowed FWHM,
319
        maximum allowed FWHM,
320
        pointer to the lower FWHM range (output),
321
        pointer to the upper FWHM range (output).
322
OUTPUT  FWHM mode.
323
NOTES   -.
324
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
325
VERSION 20/03/2008
326
*/
327
static float    compute_fwhmrange(float *fwhm, int nfwhm, float maxvar,
328
                float minin, float maxin, float *minout, float *maxout)
329
  {
330
   float        *fwhmt,*fwhmt2,
331
                df, dfmin,fmin;
332
   int          i, nw;
333
 
334
/* Sort FWHMs */
335 157 bertin
   fqmedian(fwhm, nfwhm);
336 79 bertin
 
337
/* Find the mode */
338
   nw = nfwhm/4;
339
   if (nw<4)
340
     nw = 1;
341
  dfmin = BIG;
342
  fmin = 0.0;
343
  fwhmt = fwhm;
344
  fwhmt2 = fwhm+nw;
345
  for (i=nfwhm-nw; i--; fwhmt++,fwhmt2++)
346
    {
347
    if ((df = *fwhmt2 - *fwhmt) < dfmin)
348
      {
349
      dfmin = df;
350
      fmin = (*fwhmt2 + *fwhmt)/2.0;
351
      }
352
    }
353
 
354
  if (nfwhm<2)
355
    fmin = *fwhm;
356
 
357
  dfmin = (float)pow((double)maxvar+1.0, 0.3333333);
358
  *minout = dfmin>0.0?fmin/dfmin:0.0;
359
  if (*minout<minin)
360
    *minout = minin;
361
  *maxout = fmin*dfmin*dfmin;
362
  if (*maxout>maxin)
363
    *maxout = maxin;
364
 
365
  return fmin;
366
  }
367
 
368
 
369 170 bertin
/****** read_samples *********************************************************
370
PROTO   setstruct *read_samples(setstruct *set, char *filename,
371
                        float frmin, float frmax,
372
                        int ext, int next, int catindex,
373
                        contextstruct *context, double *pcval)
374
 
375
PURPOSE Read point source data for a given set.
376
INPUT   Pointer to the data set,
377
        catalogue filename,
378
        minimum flux radius,
379
        maximum flux radius,
380
        current extension,
381
        number of extensions,
382
        catalogue index
383
        pointer to the context,
384
        pointer to the array of principal components, if available.
385
OUTPUT  Pointer to a set containing samples that match acceptance criteria.
386
NOTES   -.
387
AUTHOR  E. Bertin (IAP)
388 184 bertin
VERSION 19/07/2012
389 2 bertin
*/
390
setstruct *read_samples(setstruct *set, char *filename,
391 80 bertin
                        float frmin, float frmax,
392
                        int ext, int next, int catindex,
393 73 bertin
                        contextstruct *context, double *pcval)
394 2 bertin
 
395
  {
396
   catstruct            *cat;
397
   tabstruct            *tab, *keytab;
398
   keystruct            *key, *vigkey;
399
   samplestruct         *sample;
400
   t_type               contexttyp[MAXCONTEXT];
401 72 bertin
   void                 *contextvalp[MAXCONTEXT];
402 2 bertin
   static char          str[MAXCHAR], str2[MAXCHAR];
403
   char                 **kstr,
404
                        *head, *buf;
405 184 bertin
   unsigned int         *imaflags;
406
   unsigned short       *flags, *wflags;
407 147 bertin
   double               contextval[MAXCONTEXT],
408 170 bertin
                        *dxm,*dym, *cmin, *cmax, dval;
409
   float                *xm, *ym, *vignet,*vignett, *flux, *fluxrad, *elong,
410
                        *snr,
411 2 bertin
                        backnoise, backnoise2, gain, minsn,maxelong;
412
   static int           ncat;
413 147 bertin
   int                  *lxm,*lym,
414
                        i,j, n, nsample,nsamplemax,
415 109 bertin
                        vigw, vigh, vigsize, nobj, nt,
416 150 bertin
                        maxbad, maxbadflag, ldflag, ext2, pc, contflag;
417 147 bertin
   short                *sxm,*sym;
418 2 bertin
 
419
 
420
  maxbad = prefs.badpix_nmax;
421
  maxbadflag = prefs.badpix_flag;
422 86 bertin
  maxelong = (float)(prefs.maxellip < 1.0?
423
        (prefs.maxellip + 1.0)/(1.0 - prefs.maxellip)
424
        : 100.0);
425 2 bertin
  minsn = prefs.minsn;
426
 
427
/* If a NULL pointer is provided, we allocate a new set */
428
  if (!set)
429
    {
430 72 bertin
    set = init_set(context);
431 2 bertin
    nsample = nsamplemax = 0;
432
    ncat = 1;
433
    }
434
  else
435
    nsample = nsamplemax = set->nsample;
436
 
437
  cmin = cmax = (double *)NULL; /* To avoid gcc -Wall warnings */
438
  if (set->ncontext)
439
    {
440
    QMALLOC(cmin, double, set->ncontext);
441
    QMALLOC(cmax, double, set->ncontext);
442
    for (i=0; i<set->ncontext; i++)
443 59 bertin
      if (ncat>1 && set->nsample)
444
        {
445
        cmin[i] = set->contextoffset[i] - set->contextscale[i]/2.0;
446
        cmax[i] = cmin[i] + set->contextscale[i];
447
        }
448
      else
449
        {
450
        cmin[i] = BIG;
451
        cmax[i] = -BIG;
452
        }
453 2 bertin
    }
454
 
455
/*-- Read input catalog */
456
  if (!(cat = read_cat(filename)))
457
    error(EXIT_FAILURE, "*Error*: No such catalog: ", filename);
458
  head = (char *)NULL;  /* To avoid gcc -Wall warnings */
459
  ldflag = 1;
460
  ext2 = ext+1;
461
  tab = cat->tab;
462
  for (j=cat->ntab; j--; tab=tab->nexttab)
463
    if (!(ldflag = strcmp("LDAC_IMHEAD",tab->extname))
464
          || fitsread(head=tab->headbuf,"SEXBKDEV",&backnoise,H_FLOAT,T_FLOAT)
465
                ==RETURN_OK)
466
      if (!--ext2)
467
        break;
468
  if (j<0)
469
    error(EXIT_FAILURE, "*Error*: SExtractor table missing in ", filename);
470
  if (!ldflag)
471
    {
472
    key=read_key(tab, "Field Header Card");
473
    head = key->ptr;
474
    if (fitsread(head,"SEXBKDEV",&backnoise,H_FLOAT,T_FLOAT)==RETURN_ERROR)
475
      error(EXIT_FAILURE, "*Error*: Keyword not found:", "SEXBKDEV");
476
    }
477
  backnoise2 = backnoise*backnoise;
478
  if ((n=fitsfind(head, "END     ")) != RETURN_ERROR)
479
    {
480
    QCALLOC(set->head, char, ((n*80)/FBSIZE+1)*FBSIZE);
481
    memcpy(set->head, head, (n+1)*80);
482
    }
483
  if (fitsread(head, "SEXGAIN", &gain, H_FLOAT, T_FLOAT) == RETURN_ERROR)
484
    error(EXIT_FAILURE, "*Error*: Keyword not found:", "SEXGAIN");
485
 
486
  ext2 = ext+1;
487
  tab = cat->tab;
488
  for (j=cat->ntab; j--; tab=tab->nexttab)
489
    if (!strcmp("LDAC_OBJECTS", tab->extname)
490
                ||  !strcmp("OBJECTS", tab->extname))
491
      if (!--ext2)
492
        break;
493
  if (j<0)
494
    error(EXIT_FAILURE, "*Error*: OBJECTS table not found in catalog ",
495
                filename);
496
 
497
/* Init the single-row tab */
498 147 bertin
  dxm = dym = NULL;
499
  xm = ym = NULL;
500
  lxm = lym = NULL;
501
  sxm = sym = NULL;
502 2 bertin
  keytab = init_readobj(tab, &buf);
503
 
504 149 bertin
  if (!(key = name_to_key(keytab, prefs.center_key[0])))
505
    {
506
    sprintf(str, "*Error*: %s parameter not found in catalogue ",
507
        prefs.center_key[0]);
508
    error(EXIT_FAILURE, str, filename);
509
    }
510 147 bertin
  if (key->ttype == T_DOUBLE)
511
    dxm = (double *)key->ptr;
512
  else if (key->ttype == T_FLOAT)
513
    xm = (float *)key->ptr;
514
  else if (key->ttype == T_LONG)
515
    lxm = (int *)key->ptr;
516
  else
517
    sxm = (short *)key->ptr;
518
 
519 149 bertin
  if (!(key = name_to_key(keytab, prefs.center_key[1])))
520
    {
521
    sprintf(str, "*Error*: %s parameter not found in catalogue ",
522
        prefs.center_key[0]);
523
    error(EXIT_FAILURE, str, filename);
524
    }
525 147 bertin
   if (key->ttype == T_DOUBLE)
526
    dym = (double *)key->ptr;
527
  else if (key->ttype == T_FLOAT)
528
    ym = (float *)key->ptr;
529
  else if (key->ttype == T_LONG)
530
    lym = (int *)key->ptr;
531
  else
532
    sym = (short *)key->ptr;
533 2 bertin
 
534
  if (!(key = name_to_key(keytab, "FLUX_RADIUS")))
535
    error(EXIT_FAILURE, "*Error*: FLUX_RADIUS parameter not found in catalog ",
536
                filename);
537
  fluxrad = (float *)key->ptr;
538
 
539 149 bertin
  if (!(key = name_to_key(keytab, prefs.photflux_rkey)))
540
    {
541
    sprintf(str, "*Error*: %s parameter not found in catalogue ",
542
        prefs.photflux_rkey);
543
    error(EXIT_FAILURE, str, filename);
544
    }
545 2 bertin
  flux = (float *)key->ptr;
546 149 bertin
  n = prefs.photflux_num - 1;
547
  if (n)
548
    {
549
    if (key->naxis==1 && n<key->naxisn[0])
550
      flux += n;
551
    else
552
      {
553
      sprintf(str, "Not enough apertures for %s in catalogue %s: ",
554
        prefs.photflux_rkey,
555
        filename);
556
      warning(str, "using first aperture");
557
      }
558
    }
559 2 bertin
 
560 170 bertin
  if (!(key = name_to_key(keytab, "SNR_WIN")))
561 149 bertin
    {
562 170 bertin
    sprintf(str, "*Error*: SNR_WIN parameter not found in catalogue ");
563 149 bertin
    error(EXIT_FAILURE, str, filename);
564
    }
565 170 bertin
  snr = (float *)key->ptr;
566 2 bertin
 
567 184 bertin
  if ((key = name_to_key(keytab, "ELONGATION")))
568
    elong = (float *)key->ptr;
569
  else
570
    elong = NULL;
571 2 bertin
 
572 184 bertin
/* Load optional SExtractor FLAGS parameter */
573
  if ((key = name_to_key(keytab, "FLAGS")))
574
    flags = (unsigned short *)key->ptr;
575
  else
576
    flags = NULL;
577 2 bertin
 
578 184 bertin
/* Load optional SExtractor FLAGS_WEIGHT parameter */
579
  if ((key = name_to_key(keytab, "FLAGS_WEIGHT")))
580
    wflags = (unsigned short *)key->ptr;
581
  else
582
    wflags = NULL;
583
 
584
/* Load optional SExtractor IMAFLAGS_ISO parameter */
585
  if ((key = name_to_key(keytab, "IMAFLAGS_ISO")))
586
    imaflags = (unsigned int *)key->ptr;
587
  else
588
    imaflags = NULL;
589
 
590 2 bertin
  if (!(key = name_to_key(keytab, "VIGNET")))
591
    error(EXIT_FAILURE,
592
        "*Error*: VIGNET parameter not found in catalog ", filename);
593
  vignet = (float *)key->ptr;
594 184 bertin
  nobj = key->nobj;
595
 
596 2 bertin
  if (key->naxis != 2)
597
    error(EXIT_FAILURE, "*Error*: VIGNET should be a 2D vector", "");
598
  vigkey = key;
599
  vigw = *(vigkey->naxisn);
600
  vigh = *(vigkey->naxisn+1);
601
  vigsize = vigw*vigh;
602 189 rhl
  if (!set->nsample)
603 2 bertin
    {
604
    set->vigsize[0] = vigw;
605
    set->vigsize[1] = vigh;
606
    set->nvig = vigw*vigh;
607
    }
608
 
609
/* Try to load the set of context keys */
610 72 bertin
  kstr = context->name;
611 73 bertin
  pc = 0;
612 2 bertin
  for (i=0; i<set->ncontext; i++, kstr++)
613 77 bertin
    if (context->pcflag[i])
614 2 bertin
      {
615 77 bertin
      contextvalp[i] = &pcval[pc++];
616
      contexttyp[i] = T_DOUBLE;
617
      }
618
    else if (**kstr==(char)':')
619
      {
620 72 bertin
      contextvalp[i] = &contextval[i];
621 2 bertin
      contexttyp[i] = T_DOUBLE;
622 72 bertin
      if (fitsread(head, *kstr+1, contextvalp[i], H_FLOAT,T_DOUBLE)
623
                == RETURN_ERROR)
624 2 bertin
        {
625
        sprintf(str, "*Error*: %s parameter not found in the header of ",
626
                *kstr+1);
627
        error(EXIT_FAILURE, str, filename);
628
        }
629
      }
630
    else
631
      {
632
      if (!(key = name_to_key(keytab, *kstr)))
633
        {
634
        sprintf(str, "*Error*: %s parameter not found in catalog ", *kstr);
635
        error(EXIT_FAILURE, str, filename);
636
        }
637 72 bertin
      contextvalp[i] = key->ptr;
638 2 bertin
      contexttyp[i] = key->ttype;
639
      strcpy(set->contextname[i], key->name);
640
      }
641
  if (next>1)
642
    sprintf(str2, "[%d/%d]", ext+1, next);
643
  else
644
    strcpy(str2, "");
645
 
646
/* Now examine each vector of the shipment */
647
  nt = keytab->naxisn[1];
648
  for (n=0; nt; n++)
649
    {
650
    nt = read_obj(keytab,tab, buf);
651
    if (!(n%100))
652
      {
653
      sprintf(str,"Catalog #%d %s: Object #%d / %d samples stored",
654
        ncat, str2, n,nsample);
655 112 bertin
//      NFPRINTF(OUTPUT, str);
656 2 bertin
      }
657 170 bertin
 
658 2 bertin
/*---- Apply some selection over flags, fluxes... */
659 150 bertin
    contflag = 0;
660 184 bertin
    if (flags && (*flags&prefs.flag_mask))
661 2 bertin
      {
662 150 bertin
      contflag++;
663
      set->badflags++;
664
      }
665 184 bertin
    if (wflags && (*wflags&prefs.wflag_mask))
666
      {
667
      contflag++;
668
      set->badwflags++;
669
      }
670
    if (imaflags && (*imaflags&prefs.imaflag_mask))
671
      {
672
      contflag++;
673
      set->badwflags++;
674
      }
675 170 bertin
    if (*snr<minsn)
676 150 bertin
      {
677
      contflag++;
678
      set->badsn++;
679
      }
680
    if (*fluxrad<frmin)
681
      {
682
      contflag++;
683
      set->badfrmin++;
684
      }
685
    if (*fluxrad>frmax)
686
      {
687
      contflag++;
688
      set->badfrmax++;
689
      }
690 184 bertin
    if (elong && *elong>maxelong)
691 150 bertin
      {
692
      contflag++;
693
      set->badelong++;
694
      }
695
    if (contflag)
696
      continue;
697
/*-- ... and check the integrity of the sample */
698
    j = 0;
699
    vignett = vignet;
700
    for (i=vigsize; i--; vignett++)
701
      if (*vignett <= -BIG)
702
        j++;
703
    if (maxbadflag && j > maxbad)
704
      {
705
      set->badpix++;
706
      continue;
707
      }
708 2 bertin
 
709 150 bertin
/*-- Allocate memory for the first shipment */
710 189 rhl
    if (!set->nsample)
711 150 bertin
      {
712
      nsample = 0;
713
      nsamplemax = LSAMPLE_DEFSIZE;
714
      malloc_samples(set, nsamplemax);
715
      }
716
    else
717
      {
718
      if (set->vigsize[0] != vigw || set->vigsize[1] != vigh)
719
        error(EXIT_FAILURE, "*Error*: Incompatible VIGNET size found in ",
720 2 bertin
                filename);
721 150 bertin
      }
722 2 bertin
 
723 150 bertin
/*-- Increase storage space to receive new candidates if needed */
724
    if (nsample>=nsamplemax)
725
      {
726
       int      nadd=(int)(1.62*nsamplemax);
727
      nsamplemax = nadd>nsamplemax?nadd:nsamplemax+1;
728
      realloc_samples(set, nsamplemax);
729
      }
730 2 bertin
 
731 189 rhl
    sample = set->sample[nsample];
732 150 bertin
    sample->catindex = catindex;
733
    sample->extindex = ext;
734 2 bertin
 
735 150 bertin
/*-- Copy the vignet to the training set */
736
    memcpy(sample->vig, vignet, vigsize*sizeof(float));
737 2 bertin
 
738 150 bertin
    sample->norm = *flux;
739
    sample->backnoise2 = backnoise2;
740
    sample->gain = gain;
741
    if (dxm)
742
      sample->x = *dxm;
743
    else if (xm)
744 2 bertin
      sample->x = *xm;
745 150 bertin
    else if (lxm)
746
      sample->x = *lxm;
747
    else
748
      sample->x = *sxm;
749
    if (dym)
750
      sample->y = *dym;
751
    else if (ym)
752
      sample->y = *ym;
753
    else if (lym)
754
      sample->y = *lym;
755
    else
756
      sample->y = *sym;
757
    sample->dx = sample->x - (int)(sample->x+0.49999);
758
    sample->dy = sample->y - (int)(sample->y+0.49999);
759
    for (i=0; i<set->ncontext; i++)
760
      {
761
      dval = sample->context[i];
762
      ttypeconv(contextvalp[i], &dval, contexttyp[i], T_DOUBLE);
763
      sample->context[i] = dval;
764
/*---- Update min and max */
765
      if (dval<cmin[i])
766
        cmin[i] = dval;
767
      if (dval>cmax[i])
768
        cmax[i] = dval;
769 2 bertin
      }
770 150 bertin
    make_weights(set, sample);
771
    recenter_sample(sample, set, *fluxrad);
772 189 rhl
    nsample = ++set->nsample;
773 2 bertin
    }
774
 
775
/* Update the scaling */
776
  if (set->ncontext)
777
    {
778
    if (nsample)
779
      for (i=0; i<set->ncontext; i++)
780
        {
781
        set->contextscale[i] = cmax[i] - cmin[i];
782
        set->contextoffset[i] = (cmin[i] + cmax[i])/2.0;
783
        }
784
    free(cmin);
785
    free(cmax);
786
    }
787
  end_readobj(keytab,tab, buf);
788
  free_cat(&cat, 1);
789
 
790
/* Don't waste memory! */
791
  if (nsample)
792
    realloc_samples(set, nsample);
793
 
794
/* Increase the current catalog number */
795
  ncat++;
796
 
797
  return set;
798
  }