public software.psfex

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

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

Line No. Rev Author Line
1 157 bertin
/*
2
*                               field.c
3 76 bertin
*
4 157 bertin
* Manage multiple PSFs.
5 76 bertin
*
6 157 bertin
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 76 bertin
*
8 157 bertin
*       This file part of:      PSFEx
9 76 bertin
*
10 177 bertin
*       Copyright:              (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
11 76 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 177 bertin
*       Last modified:          25/06/2012
26 157 bertin
*
27
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
28 76 bertin
 
29
#ifdef HAVE_CONFIG_H
30
#include        "config.h"
31
#endif
32
 
33
#include        <math.h>
34
#include        <stdio.h>
35
#include        <stdlib.h>
36
#include        <string.h>
37
 
38
#include        "define.h"
39
#include        "types.h"
40
#include        "globals.h"
41
#include        "fits/fitscat.h"
42 87 bertin
#include        "check.h"
43 76 bertin
#include        "fitswcs.h"
44 161 bertin
#include        "misc.h"
45 76 bertin
#include        "prefs.h"
46
#include        "psf.h"
47
#include        "field.h"
48
 
49 119 bertin
/****** field_init ************************************************************
50 150 bertin
PROTO   fieldstruct *field_init(char *catname)
51 76 bertin
PURPOSE Allocate and initialize a PSF MEF structure (groups of PSFs).
52
INPUT   Catalog filename.
53
OUTPUT  fieldstruct pointer.
54
NOTES   .
55
AUTHOR  E. Bertin (IAP)
56 150 bertin
VERSION 08/04/2010
57 76 bertin
 ***/
58
fieldstruct     *field_init(char *catname)
59
  {
60
   fieldstruct  *field;
61
   catstruct    *cat;
62
   tabstruct    *tab, *imatab;
63
   keystruct    *key;
64 112 bertin
   char         *pstr;
65 191 rhl
   int          next, next0, ntab;
66 76 bertin
 
67
  QCALLOC(field, fieldstruct, 1);
68
/* Compute the number of valid input extensions */
69
  if (!(cat = read_cat(catname)))
70
    error(EXIT_FAILURE, "*Error*: cannot open ", catname);
71
  tab = cat->tab;
72 84 bertin
  next0 = 0;
73 76 bertin
  for (ntab = 0 ; ntab<cat->ntab; ntab++, tab = tab->nexttab)
74
    {
75
/*--  Check for the next valid image extension */
76
    if ((tab->naxis != 2)
77
        || (strncmp(tab->xtension, "BINTABLE", 8)
78
                && strncmp(tab->xtension, "ASCTABLE", 8))
79
        || (strncmp(tab->extname, "LDAC_OBJECTS", 8)
80
                && strncmp(tab->extname, "OBJECTS", 8)))
81
      continue;
82 84 bertin
    next0++;
83 76 bertin
    }
84 84 bertin
  field->next = next0;
85
  QMALLOC(field->psf, psfstruct *, next0);
86 76 bertin
  strcpy(field->catname, catname);
87
/* A short, "relative" version of the filename */
88
  if (!(field->rcatname = strrchr(field->catname, '/')))
89
    field->rcatname = field->catname;
90
  else
91
    field->rcatname++;
92 112 bertin
  strcpy(field->rtcatname, field->rcatname);
93
  if ((pstr=strrchr(field->rtcatname, '.')))
94
    *pstr = '\0';
95 76 bertin
 
96 107 bertin
  if (!next0)
97
    {
98
    field_end(field);
99
    error(EXIT_FAILURE,"*Error*: No SExtractor FITS-LDAC catalog found in ",
100
        catname);
101
    }
102
 
103 84 bertin
  QMALLOC(field->wcs, wcsstruct *, next0);
104 76 bertin
/* Compute the number of valid input extensions */
105
  tab = cat->tab;
106
  next = 0;
107 114 bertin
  for (ntab = 0 ; ntab<cat->ntab; ntab++, tab = tab->nexttab)
108 76 bertin
/*--  Check for the next valid FITS extension */
109
    if ((!strcmp("LDAC_IMHEAD",tab->extname))
110
        && (key=read_key(tab, "Field Header Card")))
111
      {
112
/*---- Create a new table from scratch to hold the image header */
113
      imatab = new_tab("Image header");
114
      free(imatab->headbuf);
115
      imatab->headnblock = 1 + (key->nbytes-1)/FBSIZE;
116
      QCALLOC(imatab->headbuf, char, imatab->headnblock*FBSIZE);
117
      memcpy(imatab->headbuf, key->ptr, key->nbytes);
118
      imatab->cat = cat;
119
      readbasic_head(imatab);
120
      field->wcs[next++] = read_wcs(imatab);
121 84 bertin
      if (!imatab->headbuf
122
        || fitsread(imatab->headbuf, "OBJECT  ", field->ident,
123
        H_STRING,T_STRING)!= RETURN_OK)
124
        strcpy(field->ident, "no ident");
125 76 bertin
      free_tab(imatab);
126
      }
127 111 bertin
    else if ((!strcmp("LDAC_OBJECTS", tab->extname)
128 114 bertin
        ||  !strcmp("OBJECTS", tab->extname)) && tab->naxis == 2)
129 111 bertin
      field->ndet += tab->naxisn[1];
130 76 bertin
 
131
  free_cat(&cat, 1);
132
 
133 191 rhl
  field_init_finalize(field);
134 76 bertin
 
135
  return field;
136
  }
137
 
138 119 bertin
/****** field_end *************************************************************
139 76 bertin
PROTO   void field_end(fieldstruct *field)
140
PURPOSE Free a PSF MEF structure (groups of PSFs).
141
INPUT   Pointer to the fieldstruct.
142
OUTPUT  -.
143
NOTES   .
144
AUTHOR  E. Bertin (IAP)
145 150 bertin
VERSION 08/04/2010
146 76 bertin
 ***/
147
void    field_end(fieldstruct *field)
148
  {
149
   int  ext;
150
 
151
  for (ext=0; ext<field->next; ext++)
152
    {
153
    psf_end(field->psf[ext]);
154
    end_wcs(field->wcs[ext]);
155 119 bertin
    free(field->lcount[ext]);
156
    free(field->acount[ext]);
157 150 bertin
    free(field->count[ext]);
158
    free(field->modresi[ext]);
159
    free(field->modchi2[ext]);
160 76 bertin
    }
161
  free(field->psf);
162
  free(field->wcs);
163 87 bertin
  free(field->ccat);
164 119 bertin
  free(field->lcount);
165
  free(field->acount);
166 150 bertin
  free(field->count);
167
  free(field->modchi2);
168
  free(field->modresi);
169 76 bertin
  free(field);
170
 
171
  return;
172
  }
173
 
174
 
175
/****** field_locate *********************************************************
176
PROTO   void field_locate(fieldstruct *field)
177
PURPOSE Compute field position, scale and footprint.
178
INPUT   Pointer to field structure.
179
OUTPUT  A pointer to the created field structure.
180
NOTES   Global preferences are used.
181
AUTHOR  E. Bertin (IAP)
182 157 bertin
VERSION 05/10/2010
183 76 bertin
*/
184
void    field_locate(fieldstruct *field)
185
  {
186
   wcsstruct            *wcs;
187
   double               *scale[NAXIS],*scalet[NAXIS],
188
                        *wcsmean,
189
                        cosalpha,sinalpha, sindelta, dist, maxradius;
190
   int                  i, e, lat,lng, naxis;
191
 
192
/* Some initializations */
193
  cosalpha = sinalpha = sindelta = 0.0;
194
  wcs = field->wcs[0];
195
  naxis = wcs->naxis;
196
  wcsmean = field->meanwcspos;
197
  for (i=0; i<naxis; i++)
198
    {
199
    QMALLOC(scale[i], double, field->next);
200
    scalet[i] = scale[i];
201
    wcsmean[i] = 0.0;
202
    }
203
 
204
/* Go through each extension */
205
  for (e=0; e<field->next; e++)
206
    {
207
    wcs = field->wcs[e];
208
    lng = wcs->lng;
209
    lat = wcs->lat;
210
/*-- Locate set */
211
    if (lat != lng)
212
      {
213
      cosalpha += cos(wcs->wcsscalepos[lng]*DEG);
214
      sinalpha += sin(wcs->wcsscalepos[lng]*DEG);
215
      sindelta += sin(wcs->wcsscalepos[lat]*DEG);
216
      }
217
    for (i=0; i<naxis; i++)
218
      {
219
      if (lat==lng || (i!=lng && i!=lat))
220
        wcsmean[i] += wcs->wcsscalepos[i];
221
      *(scalet[i]++) = wcs->wcsscale[i];
222
      }
223
    }
224
 
225
/* Now make the stats on each axis */
226
  lng = field->wcs[0]->lng;
227
  lat = field->wcs[0]->lat;
228
  for (i=0; i<naxis; i++)
229
    {
230
    if (lat!=lng && (i==lng))
231
      {
232
      wcsmean[i] = atan2(sinalpha/field->next,cosalpha/field->next)/DEG;
233
      wcsmean[i] = fmod(wcsmean[i]+360.0, 360.0);
234
      }
235
    else if (lat!=lng && (i==lat))
236
      wcsmean[i] = asin(sindelta/field->next)/DEG;
237
    else
238
      wcsmean[i] /= field->next;
239 157 bertin
    field->meanwcsscale[i] = dqmedian(scale[i], field->next);
240 76 bertin
    }
241
 
242
/* Compute the radius of the field and mean airmass */
243
  maxradius = 0.0;
244
  for (e=0; e<field->next; e++)
245
    {
246
    wcs = field->wcs[e];
247
/*-- The distance is the distance to the center + the diagonal of the image */
248
    dist = wcs_dist(wcs, wcs->wcsscalepos, field->meanwcspos)
249
                + wcs->wcsmaxradius;
250
    if (dist>maxradius)
251
      maxradius = dist;
252
    }
253
 
254
  field->maxradius = maxradius;
255
 
256
/* Free memory */
257
  for (i=0; i<naxis; i++)
258
    free(scale[i]);
259
 
260
  return;
261
  }
262
 
263
 
264 119 bertin
/****** field_count ***********************************************************
265
PROTO   void field_count(fieldstruct **fields, setstruct *set, int counttype)
266
PURPOSE Count the number of sources (samples) per image area.
267
INPUT   Pointer to an array of fieldstruct pointers,
268 150 bertin
        Pointer to the set to be counted,
269 119 bertin
        Sample type.
270
OUTPUT  -.
271
NOTES   -.
272
AUTHOR  E. Bertin (IAP)
273
VERSION 30/03/2009
274
 ***/
275
void    field_count(fieldstruct **fields, setstruct *set, int counttype)
276
  {
277
   fieldstruct  *field;
278
   samplestruct *sample;
279
   int          c,e,n,s, w,h, x,y, size;
280
 
281
  size = (double)prefs.context_nsnap;
282 189 rhl
  for (s=0; s<set->nsample; ++s)
283 119 bertin
    {
284 189 rhl
    sample = set->sample[s];
285 119 bertin
    c = sample->catindex;
286
    e = sample->extindex;
287
    field = fields[c];
288
    w = field->wcs[e]->naxisn[0];
289
    h = field->wcs[e]->naxisn[1];
290
    x = (int)((sample->x-0.5)*size) / w;
291
    if (x<0)
292
      x = 0;
293
    else if (x>=size)
294
      x = size-1;
295
    y = (int)((sample->y-0.5)*size) / h;
296
    if (y<0)
297
      y = 0;
298
    else if (y>=size)
299
      y = size-1;
300
    n = y*size+x;
301
    if ((counttype & COUNT_LOADED))
302
      fields[c]->lcount[e][n]++;
303
    if ((counttype & COUNT_ACCEPTED))
304
      fields[c]->acount[e][n]++;
305
    }
306
 
307
  return;
308
  }
309
 
310
 
311 150 bertin
/****** field_stats **********************************************************
312
PROTO   void field_stats(fieldstruct **fields, setstruct *set)
313
PURPOSE Compute the average stats per image area.
314
INPUT   Pointer to an array of fieldstruct pointers,
315
        Pointer to the sample set.
316
OUTPUT  -.
317
NOTES   -.
318
AUTHOR  E. Bertin (IAP)
319 177 bertin
VERSION 25/06/2012
320 150 bertin
 ***/
321
void    field_stats(fieldstruct **fields, setstruct *set)
322
  {
323
   fieldstruct  *field;
324
   samplestruct *sample;
325
   int          c,e,n,s, w,h, x,y, size;
326
 
327
  size = (double)prefs.context_nsnap;
328
  for (s=set->nsample; s--; sample++)
329
    {
330 189 rhl
    sample = set->sample[s];
331 150 bertin
    c = sample->catindex;
332
    e = sample->extindex;
333
    field = fields[c];
334
    w = field->wcs[e]->naxisn[0];
335
    h = field->wcs[e]->naxisn[1];
336
    x = (int)((sample->x-0.5)*size) / w;
337
    if (x<0)
338
      x = 0;
339
    else if (x>=size)
340
      x = size-1;
341
    y = (int)((sample->y-0.5)*size) / h;
342
    if (y<0)
343
      y = 0;
344
    else if (y>=size)
345
      y = size-1;
346
    n = y*size+x;
347 177 bertin
    field->count[e][n]++;
348
    field->modchi2[e][n] += sample->chi2;
349
    field->modresi[e][n] += sample->modresi;
350 150 bertin
    }
351
 
352
  return;
353
  }
354
 
355
 
356 76 bertin
/****** field_psfsave *********************************************************
357
PROTO   void    field_psfsave(fieldstruct *field, char *filename)
358
PURPOSE Save PSF data as a Multi-extension FITS file.
359 176 bertin
INPUT   Pointer to the field structure,
360 76 bertin
        Filename,
361
        Extension number,
362
        Number of extensions.
363
OUTPUT  -.
364
NOTES   -.
365
AUTHOR  E. Bertin (IAP)
366 143 bertin
VERSION 30/10/2009
367 76 bertin
 ***/
368
void    field_psfsave(fieldstruct *field, char *filename)
369
  {
370
   catstruct    *cat;
371
   tabstruct    *tab;
372
   keystruct    *key;
373
   psfstruct    *psf;
374 142 bertin
   float        zero = 0.0;
375 76 bertin
   char         *head,
376
                str[80];
377
   int          i, ext, temp;
378
 
379
  cat = new_cat(1);
380
  init_cat(cat);
381
  sprintf(cat->filename, filename);
382
  if (open_cat(cat, WRITE_ONLY) != RETURN_OK)
383
    error(EXIT_FAILURE, "*Error*: cannot open for writing ", cat->filename);
384
/* Write primary HDU */
385
  save_tab(cat, cat->tab);
386
 
387
  for (ext=0; ext<field->next; ext++)
388
    {
389
    psf = field->psf[ext];
390
    tab = new_tab("PSF_DATA");
391
 
392
    head = tab->headbuf;
393
    addkeywordto_head(tab, "LOADED", "Number of loaded sources");
394
    fitswrite(head, "LOADED", &psf->samples_loaded, H_INT, T_LONG);
395
    addkeywordto_head(tab, "ACCEPTED", "Number of accepted sources");
396
    fitswrite(head, "ACCEPTED", &psf->samples_accepted, H_INT, T_LONG);
397
    addkeywordto_head(tab, "CHI2", "Final Chi2");
398
    fitswrite(head, "CHI2", &psf->chi2, H_FLOAT, T_DOUBLE);
399
    addkeywordto_head(tab, "POLNAXIS", "Number of context parameters");
400
    fitswrite(head, "POLNAXIS", &psf->poly->ndim, H_INT, T_LONG);
401
    for (i=0; i<psf->poly->ndim; i++)
402
      {
403
      sprintf(str, "POLGRP%1d", i+1);
404
      addkeywordto_head(tab, str, "Polynom group for this context parameter");
405
      temp = psf->poly->group[i]+1;
406
      fitswrite(head, str, &temp, H_INT, T_LONG);
407
      sprintf(str, "POLNAME%1d", i+1);
408
      addkeywordto_head(tab, str, "Name of this context parameter");
409
      fitswrite(head, str, psf->contextname[i], H_STRING, T_STRING);
410
      sprintf(str, "POLZERO%1d", i+1);
411
      addkeywordto_head(tab, str, "Offset value for this context parameter");
412
      fitswrite(head, str, &psf->contextoffset[i], H_EXPO, T_DOUBLE);
413
      sprintf(str, "POLSCAL%1d", i+1);
414
      addkeywordto_head(tab, str, "Scale value for this context parameter");
415
      fitswrite(head, str, &psf->contextscale[i], H_EXPO, T_DOUBLE);
416
      }
417
 
418
    addkeywordto_head(tab, "POLNGRP", "Number of context groups");
419
    fitswrite(head, "POLNGRP", &psf->poly->ngroup, H_INT, T_LONG);
420
    for (i=0; i<psf->poly->ngroup; i++)
421
      {
422
      sprintf(str, "POLDEG%1d", i+1);
423
      addkeywordto_head(tab, str, "Polynom degree for this context group");
424
      fitswrite(head, str, &psf->poly->degree[i], H_INT, T_LONG);
425
      }
426
 
427
/*-- Add and write important scalars as FITS keywords */
428
    addkeywordto_head(tab, "PSF_FWHM", "PSF FWHM");
429 142 bertin
    fitswrite(head, "PSF_FWHM", psf->samples_accepted? &psf->fwhm : &zero,
430
        H_FLOAT, T_FLOAT);
431 76 bertin
    addkeywordto_head(tab, "PSF_SAMP", "Sampling step of the PSF data");
432 142 bertin
    fitswrite(head, "PSF_SAMP", psf->samples_accepted? &psf->pixstep : &zero,
433
        H_FLOAT, T_FLOAT);
434 76 bertin
    addkeywordto_head(tab, "PSFNAXIS", "Dimensionality of the PSF data");
435
    fitswrite(head, "PSFNAXIS", &psf->dim, H_INT, T_LONG);
436
    for (i=0; i<psf->dim; i++)
437
      {
438
      sprintf(str, "PSFAXIS%1d", i+1);
439
      addkeywordto_head(tab, str, "Number of element along this axis");
440
      fitswrite(head, str, &psf->size[i], H_INT, T_LONG);
441
      }
442
 
443 143 bertin
/*-- PSF pixels */
444 76 bertin
    key = new_key("PSF_MASK");
445
    key->naxis = psf->dim;
446
    QMALLOC(key->naxisn, int, key->naxis);
447
    for (i=0; i<psf->dim; i++)
448
      key->naxisn[i] = psf->size[i];
449
    strcat(key->comment, "Tabulated PSF data");
450
    key->htype = H_FLOAT;
451
    key->ttype = T_FLOAT;
452
    key->nbytes = psf->npix*t_size[T_FLOAT];
453
    key->nobj = 1;
454
    key->ptr = psf->comp;
455
    add_key(key, tab, 0);
456
 
457 143 bertin
/*-- Basis coefficient (if applicable) */
458
    if (psf->basiscoeff)
459
      {
460
      key = new_key("PSF_COEFFS");
461
      key->naxis = psf->dim - 1;
462
      QMALLOC(key->naxisn, int, key->naxis);
463
      key->naxisn[0] = psf->nbasis;
464
      if (key->naxis>1)
465
        key->naxisn[1] = psf->size[2];
466
      strcat(key->comment, "PSF basis vector coefficients");
467
      key->htype = H_FLOAT;
468
      key->ttype = T_FLOAT;
469
      key->nbytes = psf->nbasis*psf->size[2]*t_size[T_FLOAT];
470
      key->nobj = 1;
471
      key->ptr = psf->basiscoeff;
472
      add_key(key, tab, 0);
473
      }
474 76 bertin
    save_tab(cat, tab);
475
/*-- But don't touch my arrays!! */
476
    blank_keys(tab);
477
    free_tab(tab);
478
    }
479
 
480
  free_cat(&cat , 1);
481
 
482
  return;
483
  }
484
 
485 111 bertin