IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10897


Ignore:
Timestamp:
Jan 3, 2007, 12:29:46 PM (19 years ago)
Author:
eugene
Message:

working on MEF and glob files

Location:
trunk/Ohana/src/addstar/src
Files:
3 added
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/src/LoadStars.c

    r10889 r10897  
    11
    2 /* read an image header corresponding to a CMF data block */
     2/* read an image header corresponding to a CMF / CMP data block */
    33Image *ReadImageHeader (Header *header) {
    44

  • trunk/Ohana/src/addstar/src/gstars.c

    r10889 r10897  
    11# include "addstar.h"
    22
    3 /* break into sections, loop over images and load stars one at a time */
     3enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU};
     4/* note: MEF implies CMF */
    45
    5 Stars *gstars (char *file, int *NSTARS, int photcode, Image **imageList, int Nimages) {
     6Stars *gstars (char *filename, int *NSTARS, int photcode, Image **imageList, int Nimages) {
    67
     8  int i, Nfile, Nimage;
     9  char **file;
    710  FILE *f;
    8   int i, j, N, Nbytes, naxis;
    9   int itmp, hour, min, simple;
    10   char *name, *c, photname[64], line[80];
    11   double tmp, sec, dMs;
    12   Stars *stars, *rdstars;
    13   Header header;
    14   Image *images;
     11  glob_t globList;
     12  char **exthead, **extdata;
     13  int *extnum_head, *extnum_data, *extsize;
     14  Header *header, **headers;
     15  Image *image;
     16  Stars *inStars, *stars;
    1517
    16   /* XXX which of these values are required, which are just interesting? */
    17 
    18   # if (0)
    19   /* file may be a glob, in which case it should represent multiple chip files
    20      from a single mosaic image */
    21 
    22   // parse the word as a glob
    23   glob_t globList;
     18  // parse the filename as a glob
    2419  globList.gl_offs = 0;
    2520  glob (file, 0, NULL, &globList);
    … …  
    2823  // otherwise save all glob matches
    2924  if (globList.gl_pathc == 0) {
    30     if (!gfits_read_header (file, &header)) {
    31       fprintf (stderr, "ERROR: can't read header for %s\n", file);
    32       exit (1);
     25    Nfile = 1;
     26    ALLOCATE (file, char *, Nfile);
     27    file[0] = strcreate (filename);
     28  } else {
     29    Nfile = globList.gl_pathc;
     30    ALLOCATE (file, char *, Nfile);
     31    for (i = 0; i < Nfile; i++) {
     32      file[i] = strcreate (globList.gl_pathv[i]);
     33    }
     34  }
     35
     36  // open the first file, read the PHU header
     37  f = fopen (file[0], "r");
     38  if (f == NULL) {
     39    fprintf (stderr, "ERROR: can't read header for %s\n", file[0]);
     40    exit (1);
     41  }
     42  ALLOCATE (header, Header, 1);
     43  gfits_fread_header (f, header);
     44
     45  mode = GetFileMode (header);
     46
     47  stars = NULL;
     48  Nstars = 0;
     49
     50  // valid for all but MEF files
     51  if ((mode != SIMPLE_MEF) && (mode != MOSAIC_MEF)) {
     52    Nimage = Nfile;
     53    ALLOCATE (image, Image, Nimage);
     54    for (i = 0; i < Nfile; i++) {
     55      if (i > 0) {
     56        f = fopen (file[i], "r");
     57        if (f == NULL) {
     58          fprintf (stderr, "ERROR: can't read header for %s\n", file[i]);
     59          exit (1);
     60        }
     61        gfits_fread_header (f, header);
     62      }
     63
     64      ReadImageHeader (header, &image[i]);
     65   
     66      switch (mode) {
     67        case SIMPLE_CMP:
     68        case MOSAIC_CMP:
     69          inStars = ReadStarsTEXT (f, header, NULL, &image[i].nstars);
     70          inStars = FilterStars (inStars, &image[i]);
     71          stars = MergeStars (inStars, NinStars, stars, &Nstars);
     72          break;
     73
     74        case SIMPLE_CMF:
     75        case MOSAIC_CMF:
     76          inStars = ReadStarsFITS (f, header, NULL, &image[i].nstars);
     77          inStars = FilterStars (inStars, &image[i]);
     78          stars = MergeStars (inStars, NinStars, stars, &Nstars);
     79          break;
     80
     81        case MOSAIC_PHU:
     82          inStars = NULL;
     83          NinStars = 0;
     84          break;
     85      }
     86      fclose (f);
     87      gfits_free_header (header);
    3388    }
    3489  } else {
    35     for (i = 0; i < globList.gl_pathc; i++) {
    36       if (!gfits_read_header (globList.gl_pathv[i], &header)) {
    37         fprintf (stderr, "ERROR: can't read header for %s\n", globList.gl_pathv[i]);
    38         exit (1);
     90   
     91    /* we need to examine the extensions to determine the headers and the data */
     92    NHEADER = 10;
     93    ALLOCATE (headers, Header *, NHEADER);
     94
     95    // the first header is already loaded
     96    headers[0] = header;
     97    Nskip = gfits_matrix_size (header);
     98    fseek (f, Nskip, SEEK_CUR);
     99
     100    // load all headers into memory
     101    done = FALSE;
     102    for (i = 1; !done; i++) {
     103      ALLOCATE (header, Header, 1);
     104     
     105      status = gfits_fread_header (f, header);
     106      if (!status) {
     107        done = TRUE;
     108      } else {
     109        headers[i] = header;
     110        Nskip = gfits_matrix_size (header);
     111        fseek (f, Nskip, SEEK_CUR);
     112      }
     113      if (i == NHEADER - 1) {
     114        NHEADER += 10;
     115        REALLOCATE (headers, Header *, NHEADER);
    39116      }
    40117    }
    41   }
    42   # endif
     118    Nheader = i;
     119   
     120    // space to store the images, indexes to the matching headers
     121    Nimage = 0;
     122    NIMAGE = Nheader;
     123    ALLOCATE (image, Image, NIMAGE);
     124    ALLOCATE (exthead, char **, NIMAGE);
     125    ALLOCATE (extdata, char **, NIMAGE);
     126    ALLOCATE (extnum_head, int, NIMAGE);
     127    ALLOCATE (extnum_data, int, NIMAGE);
     128    ALLOCATE (extsize, int, Nheader);
    43129
    44   if (!gfits_read_header (file, &header)) {
    45     fprintf (stderr, "ERROR: can't read header for %s\n", file);
    46     exit (1);
    47   }
     130    if (mode == MOSAIC_MEF) {
     131      exthead[Nimage] = strcreate ("PHU");
     132      extdata[Nimage] = strcreate ("NONE");
     133      extnum_data[Nimage] = -1;
     134      extnum_head[Nimage] = 0;
     135      Nimage ++;
     136    }
    48137
    49   /* find image rootname */
    50   /* if we are reading a MEF mosaic image, attached [extname] to the name of the image */
    51   name = filebasename (file);
    52   strcpy (image[0].name, name);
    53   free (name);
     138    // now examine the headers, count the table entries, find corresponding headers
     139    for (i = 0; i < Nheader; i++) {
     140      extsize[i] = headers[i][0].size + gfits_matrix_size (headers[i]);
     141      gfits_scan (headers[i], "EXTTYPE", "%s", 1, tmpword);
     142      if (!strcmp (tmpword, "SMPDATA") || 
     143          !strcmp (tmpword, "PS1DATA")) {
     144        exttype[Nimage] = strcreate (tmpword);
     145        gfits_scan (headers[i], "EXTNAME", "%s", 1, tmpword);
     146        extdata[Nimage] = strcreate (tmpword);
     147        gfits_scan (header[i], "EXTHEAD", "%s", 1, tmpword);
     148        exthead[Nimage] = strcreate (tmpword);
     149        extnum_data[Nimage] = i;
     150        extnum_head[Nimage] = -1;
     151        // find the matching exthead entry
     152        for (j = 0; j < Nheader; j++) {
     153          gfits_scan (header[j], "EXTNAME", "%s", 1, tmpword);
     154          if (!strcmp (tmpword, exthead[Nimage])) {
     155            extnum_head[Nimage] = j;
     156          }
     157        }
     158        // skip or crash on table with missing matching header?
     159        if (extnum_head[Nimage] == -1) {
     160          fprintf (stderr, "ERROR: can't read header for %s\n", file[0]);
     161          exit (1);
     162        }
     163        Nimage ++;
     164      }
     165    }
    54166
    55   /* this is really lame.  get sky background from Flips data */
    56   /*** the problems are:
    57        1) need to define sky entry in image table
    58        2) need to have sky measurement (in reg.db)
    59        (could be derived from star data?)
    60        
    61     XXX EAM : make a different SKYPROBE format with sky in appropriate place?
    62   ***/
     167    // now run through the images, interpret the headers and read the stars
     168    for (i = 0; i < Nimage; i++) {
     169      Nhead = extnum_head[i];
     170      ReadImageHeader (header[Nhead], &image[i]);
    63171
    64   /** I'm essentially using subpix to be synonymous with SKYPROBE */
    65   if (SUBPIX) {
    66     char *p;
    67     int sky;
    68    
    69     sky = 0;
    70     gfits_scan (&header, "HISTORY", "%S", 10, line);
    71     p = strstr (line, "|I=");
    72     if (p != (char *) NULL) {
    73       p += 3;
    74       sky = atoi (p);
     172      if (!strcmp(exthead[i], "PHU")) continue;
     173
     174      // advance the pointer to the start of the corresponding table block
     175      Ndata = extnum_data[i];
     176      Nskip = 0;
     177      for (j = 0; j < Ndata; j++) {
     178        Nskip += extsize[j];
     179      }
     180      fseek (f, Nskip, SEEK_SET);
     181         
     182      inStars = ReadStarsFITS (f, header[Nhead], header[Ndata], &NinStars);
     183      inStars = FilterStars (inStars, &NinStars, &image[i]);
     184      stars = MergeStars (inStars, NinStars, stars, &Nstars);
    75185    }
    76     image[0].Myyyy = MIN (sky - 0x8000, 0xffff);
    77   }
    78 
    79   /* read in data section */
    80   f = fopen (file, "r");
    81   if (f == NULL) {
    82     fprintf (stderr, "ERROR: can't read data from %s\n", file);
    83     exit (1);
    84   }
    85   fseek (f, header.size, SEEK_SET);
    86 
    87   /* don't try to read stars for mosaic */
    88   if (!strcmp (&image[0].coords.ctype[4], "-DIS")) {
    89     rdstars = NULL;
    90     image[0].nstar = 0;
    91   } else {
    92     /* read from FITS table or from text table */
    93     naxis = 2;
    94     gfits_scan (&header, "NAXIS",  "%d", 1, &naxis);
    95     gfits_scan (&header, "SIMPLE", "%t", 1, &simple);
    96     if ((naxis == 0) && !TEXTMODE && simple) {
    97       Nbytes = gfits_matrix_size (&header);
    98       fseek (f, Nbytes, SEEK_CUR);
    99       rdstars = rfits (f, &image[0].nstar);
    100       if (rdstars == NULL) {
    101         fprintf (stderr, "ERROR: failed to read fits table\n");
    102         exit (1);
    103       }
    104     } else {
    105       /* allocate space for stars */
    106       if (!gfits_scan (&header, "NSTARS", "%d", 1, &image[0].nstar)) {
    107         fprintf (stderr, "ERROR: failed to find NSTARS\n");
    108         exit (1);
    109       }
    110       rdstars = rtext (f, &image[0].nstar);
    111     }
    112   }
    113   fclose (f);
    114 
    115   /* modify resulting star list */
    116   ALLOCATE (stars, Stars, image[0].nstar);
    117   for (N = j = 0; j < image[0].nstar; j++) {
    118     /* allow for some dynamic filtering of star list */
    119     if (SNLIMIT && rdstars[j].dM > SNLIMIT) continue;
    120     if (XMAX && (rdstars[j].X > XMAX)) continue;
    121     if (XMIN && (rdstars[j].X < XMIN)) continue;
    122     if (YMAX && (rdstars[j].Y > YMAX)) continue;
    123     if (YMIN && (rdstars[j].Y < YMIN)) continue;
    124     stars[N] = rdstars[j];
    125 
    126     XY_to_RD (&stars[N].R, &stars[N].D, stars[N].X, stars[N].Y, &image[0].coords);
    127     while (stars[N].R <    0.0) stars[N].R += 360.0;
    128     while (stars[N].R >= 360.0) stars[N].R -= 360.0;
    129     stars[N].found = -1;
    130     stars[N].code = photcode;
    131 
    132     if (SUBPIX) {
    133       dMs = get_subpix (stars[N].X, stars[N].Y);
    134       stars[N].M    -= dMs;
    135       stars[N].Mgal -= dMs;
    136       stars[N].Map  -= dMs;
    137       dMs = scat_subpix (stars[N].X, stars[N].Y);
    138       stars[N].dM = hypot (stars[N].dM, dMs);
    139     }
    140     N ++;
    141   }
    142   image[0].nstar = N;
    143   REALLOCATE (stars, Stars, image[0].nstar);
    144   free (rdstars);
    145 
    146   if (VERBOSE) fprintf (stderr, "read %d stars from target file\n", image[0].nstar);
    147   *NSTARS = image[0].nstar;
    148 
    149   return (stars);
    150 }
Note: See TracChangeset for help on using the changeset viewer.