IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27489


Ignore:
Timestamp:
Mar 27, 2010, 10:39:57 AM (16 years ago)
Author:
eugene
Message:

adding checks for bad fit and repairs

Location:
branches/eam_branches/relastro.20100326
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relastro.20100326/Makefile

    r24308 r27489  
    5656$(SRC)/save_catalogs.$(ARCH).o       \
    5757$(SRC)/write_coords.$(ARCH).o        \
     58$(SRC)/CoordOps.$(ARCH).o        \
     59$(SRC)/FixProblemImages.$(ARCH).o        \
    5860$(SRC)/relastroVisual.$(ARCH).o
    5961

  • branches/eam_branches/relastro.20100326/include/relastro.h

    r27435 r27489  
    151151void          dump_grid           PROTO((void));
    152152int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
    153 void          findImages          PROTO((Catalog *catalog, int Ncatalog));
     153void          findImages          PROTO((Catalog *catalog, int Ncatalog, int MATCHCAT));
    154154int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
    155155Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
    … …  
    166166void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
    167167int           gcatalog            PROTO((Catalog *catalog, int FINAL));
    168 Coords       *getCoords           PROTO((off_t meas, int cat));
     168// Coords       *getCoords           PROTO((off_t meas, int cat));
    169169float         getMcal             PROTO((off_t meas, int cat));
    170170float         getMgrid            PROTO((off_t meas, int cat));
    … …  
    178178void          initGrid            PROTO((int dX, int dY));
    179179void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
    180 void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
     180void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int FULLINIT));
    181181void          initImages          PROTO((Image *input, off_t N));
    182182void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
    … …  
    199199int           main                PROTO((int argc, char **argv));
    200200void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
    201 void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat));
     201void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat, int MATCHCAT));
    202202void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
    203203GSCRegion    *name_region         PROTO((char *name, off_t *Nregions));
    … …  
    250250void fit_free (CoordFit *fit);
    251251void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
    252 void fit_eval (CoordFit *fit);
     252int fit_eval (CoordFit *fit);
    253253void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1);
    254254double **poly2d_dx (double **poly, int Nx, int Ny);
    … …  
    256256double **poly2d_copy (double **poly, int Nx, int Ny);
    257257double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
    258 CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords);
     258int fit_apply_coords (CoordFit *fit, Coords *coords);
    259259int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
    260260CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
    261 void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
     261int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
    262262void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
    263263void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
    … …  
    303303                               StatType statsR, StatType statsD, double thresh);
    304304
     305
     306
     307int FixProblemImages (SkyList *skylist);
     308int *getCatlist (int *N, off_t im);
     309
     310void initCoords (void);
     311void getOffsets (double *dPos, off_t *nPos, off_t N);
     312void saveOffsets (double dPos, off_t nPos, off_t N);
     313void setBadCoords (off_t N);
     314int badCoords (off_t N);
     315Coords *getCoords (off_t N);
     316int saveCoords (Coords *coords, off_t N);
     317void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im);

  • branches/eam_branches/relastro.20100326/src/CoordOps.c

    r27488 r27489  
    1010void initCoords (void) {
    1111
    12   images = getImages (&N);
     12  off_t N;
     13  Image *images;
     14
     15  images = getimages (&N);
    1316
    1417  NoldCoords = N;
    … …  
    3235}
    3336
    34 Coords getCoords (off_t N) {
     37Coords *getCoords (off_t N) {
    3538
    3639  if (N < 0) return NULL;
    … …  
    6871}
    6972 
     73void getOffsets (double *dPos, off_t *nPos, off_t N) {
     74
     75  if (N < 0) return;
     76  if (N >= NoldCoords) return;
     77
     78  *dPos = dPosSum[N];
     79  *nPos = nPosSum[N];
     80 
     81  return;
     82}
     83 

  • branches/eam_branches/relastro.20100326/src/FitChip.c

    r27488 r27489  
    1717// XXX save measurements of the fit quality (scatter, chisq) in the image table
    1818
    19 void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
     19int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
    2020
    2121  int i, Nscatter, Niter, skip;
    2222  CoordFit *fit;
    2323  double dL, dM, dR, dRmax, *values;
    24   Coords oldCoords;
    2524
    2625  ALLOCATE (values, double, Nmatch);
    … …  
    9190      fit_free (fit);
    9291      free (values);
    93       return;
     92      return FALSE;
    9493    }
    9594

  • branches/eam_branches/relastro.20100326/src/FixProblemImages.c

    r27488 r27489  
    44// original coordinates and recalculate the positions
    55
    6 int FixProblemImages () {
     6int FixProblemImages (SkyList *skylist) {
    77
    8   off_t im, Nimage;
     8  off_t i, Nimage;
    99  Image *image;
     10  SkyList sublist;
     11
     12  // allocate so we can reallocate below
     13  ALLOCATE (sublist.regions, SkyRegion *, 1);
     14  ALLOCATE (sublist.filename, char *, 1);
    1015
    1116  image = getimages (&Nimage);
    1217
     18  // first check on the dPos reported for each image
    1319  for (i = 0; i < Nimage; i++) {
    14    
     20    double dPosSum, dPos;
     21    off_t nPos;
     22    getOffsets (&dPosSum, &nPos, i);
     23    dPos = sqrt(dPosSum / nPos);
     24    if (dPos > 4.0) {
     25      setBadCoords (i);
     26    }
     27  }
     28
     29  for (i = 0; i < Nimage; i++) {
     30    int j, cat, Ncat, *catlist, Ncatlist;
     31    Catalog *catalog;
     32
    1533    // check if this image should be fixed
    1634    if (!badCoords(i)) continue;
    1735
    18     // get list of catalogs containing measurements for the bad images:
    19     catlist = getCatList(im, &Ncat);
    20     measlist = getMeasList(im, &Nmeas);
     36    // I need a list of the catalogs for this image
     37    catlist = getCatlist(&Ncatlist, i);
    2138
    22     for (i = 0; i < Nlist[im]; i++) {
    23       m = mlist[im][i];
    24       c = clist[im][i];
    25      
     39    // allocate Ncatlist skylist regions
     40    REALLOCATE (sublist.regions, SkyRegion *, Ncatlist);
     41    REALLOCATE (sublist.filename, char *, Ncatlist);
     42    sublist.Nregions = Ncatlist;
     43    sublist.ownElements = FALSE; // this list is only holding a view to the elements
    2644
    27 
    28     // set up the basic catalog info
    29     catalog.filename  = skylist[0].filename[i];
    30     catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    31     catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    32     catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    33     catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    34 
    35     if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
    36       fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
    37       exit (1);
    38     }
    39     if (!catalog.Naves_disk) {
    40       if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
    41       dvo_catalog_unlock (&catalog);
    42       dvo_catalog_free (&catalog);
    43       continue;
     45    // copy the desired catalogs from skylist to skylistSubset
     46    for (j = 0; j < Ncatlist; j++) {
     47      cat = catlist[j];
     48      sublist.filename[j] = skylist[0].filename[cat];
     49      sublist.regions[j] = skylist[0].regions[cat];
    4450    }
    4551
     52    catalog = load_catalogs (&sublist, &Ncat, FALSE);
     53    assert (Ncat == Ncatlist);
     54
    4655    // match measurements with images
    47     initImageBins (&catalog, 1);
    48     findImages (&catalog, 1);
     56    initImageBins (catalog, Ncat, FALSE);
     57    findImages (catalog, Ncat, FALSE);
    4958
    5059    // update the detection coordinates using the new image parameters
    51     UpdateMeasures (&catalog, 1);
     60    resetImageRaw (catalog, Ncat, i);
    5261
    53     freeImageBins (1);
     62    freeImageBins (Ncat);
    5463
    5564    // write the updated detections to disk
    56     save_catalogs (&catalog, 1);
     65    save_catalogs (catalog, Ncat);
    5766  }
    5867 

  • branches/eam_branches/relastro.20100326/src/ImageOps.c

    r27488 r27489  
    1111static Image        *image;   // list of available images
    1212static off_t        Nimage;   // number of available images
     13
     14static int         *Ncatlist;  // catalogs associated with each image
     15static int         *NCATLIST;  // catalogs associated with each image
     16static int         **catlist;  // catalogs associated with each image
    1317
    1418// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
    … …  
    3236}
    3337
     38int *getCatlist (int *N, off_t im) {
     39
     40  *N = Ncatlist[im];
     41  return (catlist[im]);
     42}
     43
    3444void initImages (Image *input, off_t N) {
    3545
    … …  
    8696}
    8797
    88 void initImageBins (Catalog *catalog, int Ncatalog) {
     98// these are really image & catalog indexes
     99void initImageBins (Catalog *catalog, int Ncatalog, int FULLINIT) {
    89100
    90101  off_t i, j;
    … …  
    102113
    103114  for (i = 0; i < Nimage; i++) {
    104     Nlist[i] = 0;
     115    Nlist[i] =   0;
    105116    NLIST[i] = 100;
    106117    ALLOCATE (clist[i], int, NLIST[i]);
    107118    ALLOCATE (mlist[i], off_t, NLIST[i]);
     119  }
     120
     121  if (FULLINIT) {
     122    ALLOCATE (Ncatlist, int,  Nimage);
     123    ALLOCATE (NCATLIST, int,  Nimage);
     124    ALLOCATE (catlist, int *, Nimage);
     125
     126    for (i = 0; i < Nimage; i++) {
     127      Ncatlist[i] =  0;
     128      NCATLIST[i] = 32;
     129      ALLOCATE (catlist[i], int, NCATLIST[i]);
     130    }
    108131  }
    109132}
    … …  
    126149
    127150/* match measurements to images */
    128 void findImages (Catalog *catalog, int Ncatalog) {
     151void findImages (Catalog *catalog, int Ncatalog, int MATCHCAT) {
    129152
    130153  off_t i, j;
    … …  
    136159      // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
    137160      // if (photcode[0].code != ecode) continue;
    138       matchImage (catalog, j, i);
     161      matchImage (catalog, j, i, MATCHCAT);
    139162    }
    140163  }
    … …  
    149172# if USE_IMAGE_ID
    150173// this is the imageID-based match
    151 void matchImage (Catalog *catalog, off_t meas, int cat) {
     174void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {
    152175
    153176  off_t idx, ID;
    154177  Measure *measure;
     178  int i, found;
    155179
    156180  measure = &catalog[cat].measure[meas];
    … …  
    178202    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
    179203  }
     204
     205  if (MATCHCAT) {
     206    // index for image -> catalog list
     207    found = FALSE;
     208    for (i = 0; !found && (i < Ncatlist[idx]); i++) {
     209      if (catlist[idx][i] == cat) found = TRUE;
     210    }
     211    if (!found) {
     212      catlist[idx][Ncatlist[idx]] = cat;
     213      Ncatlist[idx] ++;
     214      if (Ncatlist[idx] == NCATLIST[idx]) {
     215        NCATLIST[idx] += 32;
     216        REALLOCATE (catlist[idx], int, NCATLIST[idx]);
     217      }
     218    }
     219  }
     220
    180221  return;
    181222}
    … …  
    183224# else
    184225// this is the time-based match
    185 void matchImage (Catalog *catalog, off_t meas, int cat) {
     226void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {
    186227
    187228  off_t i;
    … …  
    221262# endif
    222263
     264/*
    223265Coords *getCoords (off_t meas, int cat) {
    224266
    … …  
    229271  return (&image[i].coords);
    230272}
     273*/
    231274
    232275void plot_images () {
    … …  
    280323void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
    281324
    282   off_t i, m, c, n;
     325  off_t i, m, c, n, nPos;
    283326  double X, Y, L, M, P, Q, R, D, dR, dD;
     327  double dPos;
    284328
    285329  Mosaic *mosaic;
    … …  
    325369    dD = 3600.0*(catalog[c].average[n].D - D);
    326370
    327     if (fabs(catalog[c].measure[m].dR - dR) > 10.0) {
     371    if (fabs(catalog[c].measure[m].dR - dR) > 2.0) {
    328372      fprintf (stderr, "!");
    329373      setBadCoords (im); // report a failure for this image
    330374      return;
    331375    }
    332     if (fabs(catalog[c].measure[m].dD - dD) > 10.0) {
     376    if (fabs(catalog[c].measure[m].dD - dD) > 2.0) {
    333377      fprintf (stderr, "*");
    334378      setBadCoords (im); // report a failure for this image
    … …  
    355399
    356400  saveOffsets (dPos, nPos, im);
     401
     402  return;
     403}
     404
     405// return StarData values for detections in the specified image, converting coordinates from the
     406// chip positions: X,Y -> L,M -> P,Q -> R,D
     407void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
     408
     409  off_t i, m, c, n;
     410  double X, Y, L, M, P, Q, R, D, dR, dD;
     411
     412  Mosaic *mosaic;
     413  Coords *moscoords, *imcoords, *oldcoords;
     414
     415  // check if this image is bad and should be skipped
     416  if (!badCoords(im)) {
     417    fprintf (stderr, "ERROR: inconsistent result?");
     418    exit (1);
     419  }
     420
     421  // replace the current coords with the old coords:
     422  oldcoords = getCoords (im);
     423  memcpy (&image[im].coords, oldcoords, sizeof(Coords));
     424
     425  moscoords = NULL;
     426  mosaic = getMosaicForImage (im);
     427  if (mosaic != NULL) {
     428    moscoords = &mosaic[0].coords;
     429  }
     430  imcoords = &image[im].coords;
     431
     432  for (i = 0; i < Nlist[im]; i++) {
     433    m = mlist[im][i];
     434    c = clist[im][i];
     435
     436    X = catalog[c].measure[m].Xccd;
     437    Y = catalog[c].measure[m].Yccd;
     438    n = catalog[c].measure[m].averef;
     439
     440    dR = dD = 0.0;
     441    if (moscoords == NULL) {
     442      // this is a Simple image (not a mosaic)
     443      // note that for a Simple image, L,M = P,Q
     444      XY_to_LM (&L, &M, X, Y, imcoords);
     445      LM_to_RD (&R, &D, L, M, imcoords);
     446    } else {
     447      XY_to_LM (&L, &M, X, Y, imcoords);
     448      XY_to_LM (&P, &Q, L, M, moscoords);
     449      LM_to_RD (&R, &D, P, Q, moscoords);
     450    }
     451
     452    catalog[c].measure[m].dR = dR;
     453    catalog[c].measure[m].dD = dD;
     454
     455    if (catalog[c].measure[m].dR > +180.0*3600.0) {
     456      // average on high end of boundary, move star up
     457      R += 360.0;
     458      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     459    }
     460    if (catalog[c].measure[m].dR < -180.0*3600.0) {
     461      // average on low end of boundary, move star down
     462      R -= 360.0;
     463      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     464    }
     465  }
    357466
    358467  return;

  • branches/eam_branches/relastro.20100326/src/UpdateObjectOffsets.c

    r27488 r27489  
    3434
    3535    // match measurements with images
    36     initImageBins (&catalog, 1);
    37     findImages (&catalog, 1);
     36    initImageBins (&catalog, 1, FALSE);
     37    findImages (&catalog, 1, FALSE);
    3838
    3939    // update the detection coordinates using the new image parameters

  • branches/eam_branches/relastro.20100326/src/fitpoly.c

    r27488 r27489  
    146146    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
    147147  }     
    148 
    149   status = dgaussjordan (matrix, vector, fit[0].Nelems, 2);
    150   if (!status) {
     148 
     149  if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) {
    151150    return (FALSE);
    152151  }

  • branches/eam_branches/relastro.20100326/src/relastro.c

    r27488 r27489  
    4747
    4848  /* match measurements with images */
    49   initImageBins (catalog, Ncatalog);
     49  initImageBins (catalog, Ncatalog, TRUE);
    5050  MARKTIME("make image bins: %f sec\n", dtime);
    5151
    52   findImages (catalog, Ncatalog);
     52  findImages (catalog, Ncatalog, TRUE);
    5353  MARKTIME("set up image indexes: %f sec\n", dtime);
    5454
    … …  
    8686
    8787  // iterate over catalogs to make detection coordinates consistant
    88   FixProblemImages ();
     88  FixProblemImages (skylist);
    8989
    9090  // save the updated image parameters
Note: See TracChangeset for help on using the changeset viewer.