IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24308 for trunk/Ohana


Ignore:
Timestamp:
Jun 3, 2009, 9:46:36 AM (17 years ago)
Author:
beaumont
Message:

merged relastro enhancements from cnb_branch_20090301

Location:
trunk/Ohana/src
Files:
14 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/libdvo/include/dvo.h

    r21508 r24308  
    108108# define ID_MEAS_POOR_ASTROM  0x0010  /* detection is astrometry outlier                                             */   
    109109# define ID_MEAS_SKIP_ASTROM  0x0020  /* detection was ignored for astrometry measurement                            */   
    110 # define ID_MEAS_UNDEF_1      0x0040  /* unused                                                                      */
    111 # define ID_MEAS_UNDEF_2      0x0080  /* unused                                                                      */
     110# define ID_MEAS_USED_OBJ     0x0040  /* detection was used during opdate objects                                    */
     111# define ID_MEAS_USED_CHIP    0x0080  /* detection was used during update chips                                      */
    112112# define ID_MEAS_BLEND_MEAS   0x0100  /* detection is within radius of multiple objects                              */
    113113# define ID_MEAS_BLEND_OBJ    0x0200  /* multiple detections within radius of object                                 */

  • trunk/Ohana/src/relastro

  • trunk/Ohana/src/relastro/Makefile

    r17243 r24308  
    5555$(SRC)/relastro_objects.$(ARCH).o    \
    5656$(SRC)/save_catalogs.$(ARCH).o       \
    57 $(SRC)/write_coords.$(ARCH).o
     57$(SRC)/write_coords.$(ARCH).o        \
     58$(SRC)/relastroVisual.$(ARCH).o
    5859
    5960$(RELASTRO): $(INC)/relastro.h $(KAPA_INCS)

  • trunk/Ohana/src/relastro/include/relastro.h

    r21508 r24308  
    4444  double Ro, dRo;
    4545  double Do, dDo;
    46  
     46
    4747  double uR, duR;
    4848  double uD, duD;
    … …  
    6363  char flags;
    6464  Coords coords;
    65 } Mosaic; 
     65} Mosaic;
    6666
    6767typedef struct {
    … …  
    8787
    8888double SIGMA_LIM;
     89int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
    8990double MIN_ERROR;
    9091
    … …  
    119120int TimeSelect;
    120121time_t TSTART, TSTOP;
     122
     123int FlagOutlier;
     124int    CLIP_THRESH;
    121125
    122126FitMode FIT_MODE;
    … …  
    145149int           corner_check        PROTO((double *x1, double *y1, double *x2, double *y2));
    146150void          dumpGrid            PROTO((void));
    147 void          dump_grid           PROTO((void)); 
     151void          dump_grid           PROTO((void));
    148152int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
    149153void          findImages          PROTO((Catalog *catalog, int Ncatalog));
    … …  
    220224int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
    221225void          set_ZP              PROTO((double ZERO));
    222 int           setrefcode          PROTO((Image *image, int Nimage)); 
     226int           setrefcode          PROTO((Image *image, int Nimage));
    223227void          skip_measurements   PROTO((Catalog *catalog, int pass));
    224228void          sortA               PROTO((double *X, int N));
    … …  
    265269int UpdateMeasures (Catalog *catalog, int Ncatalog);
    266270void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
     271void FlagOutliers(Catalog *catalog);
     272int MeasFilterTest(Measure *measure);
    267273
    268274int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);

  • trunk/Ohana/src/relastro/src/ConfigInit.c

    r17207 r24308  
    2020
    2121  GetConfig (config, "RELASTRO_SIGMA_LIM",     "%lf", 0, &SIGMA_LIM);
     22  GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW",   "%d",  0, &SRC_MEAS_TOOFEW);
    2223
    2324  // XXX these are used in relphot to identify poor stars / images -- define
    … …  
    5556    fprintf (stderr, "directory %s does not exist, giving up\n", CATDIR);
    5657    exit (1);
    57   } 
     58  }
    5859
    5960  /* update master photcode table if not defined */

  • trunk/Ohana/src/relastro/src/FitChip.c

    r17209 r24308  
    11# include "relastro.h"
     2# include "relastroVisual.h"
    23
    34// XXX make these user parameters
    … …  
    56# define FIT_CHIP_NITER     3
    67# define FIT_CHIP_NSIGMA    3.0
     8
    79
    810// XXX we should test if the fit is sufficiently constrained across the chip, or if the
    … …  
    2224
    2325  ALLOCATE (values, double, Nmatch);
    24 
    2526  for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) {
    2627
    … …  
    2829    for (i = Nscatter = 0; i < Nmatch; i++) {
    2930      if (raw[i].mask) continue;
    30       if (raw[i].dMag > FIT_CHIP_MAX_ERROR) continue;
     31      if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue;
    3132
    3233      dL = raw[i].L - ref[i].L;
    3334      dM = raw[i].M - ref[i].M;
    3435      dR = hypot (dL, dM);
    35    
     36
    3637      values[Nscatter] = dR;
    3738      Nscatter++;
    … …  
    4142    dsort (values, Nscatter);
    4243    dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];
     44    relastroVisualPlotScatter(values, dRmax, Nscatter);
     45    relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch);
    4346
    4447    // fit the requested order polynomial
    … …  
    5053    // generate the fit matches
    5154    for (i = 0; i < Nmatch; i++) {
    52       if (raw[i].mask) continue;
     55      if (raw[i].mask) {
     56        continue;
     57      }
     58      if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) {
     59        continue;
     60      }
    5361
    5462      // only keep objects within dRmax
    … …  
    5765      dR = hypot (dL, dM);
    5866      if (dR > dRmax) continue;
    59    
     67
    6068      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
    6169    }
    … …  
    6674      case 0:
    6775      case 1:
    68         skip = (fit[0].Npts < 8);
    69         break;
     76        skip = (fit[0].Npts < 8);
     77        break;
    7078      case 2:
    71         skip = (fit[0].Npts < 11);
    72         break;
     79        skip = (fit[0].Npts < 11);
     80        break;
    7381      case 3:
    74         skip = (fit[0].Npts < 15);
    75         break;
     82        skip = (fit[0].Npts < 15);
     83        break;
    7684      default:
    77         fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
    78         abort ();
     85        fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
     86        abort ();
    7987    }
    8088    if (skip) {
    … …  
    95103      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
    96104    }
     105
    97106  }
    98107

  • trunk/Ohana/src/relastro/src/GetAstromError.c

    r16810 r24308  
    3030 
    3131  dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
     32
     33  //XXX dXccd, dYccd are now working correctly
     34  //dPtotal = AS * dPobs;
     35
    3236  dPtotal = MAX (dPtotal, MIN_ERROR);
    3337

  • trunk/Ohana/src/relastro/src/ImageOps.c

    r17243 r24308  
    11# include "relastro.h"
     2#define TESTING
    23
    34static unsigned int *start;
    … …  
    9192  for (i = 0; VERBOSE && (i < Nimage); i++) {
    9293    name = GetPhotcodeNamebyCode (image[i].photcode);
    93     fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], 
    94              ohana_sec_to_date(image[i].tzero), name);
    95   } 
     94    fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i],
     95             ohana_sec_to_date(image[i].tzero), name);
     96  }
    9697}
    9798
    … …  
    101102  int i;
    102103  Measure *measure;
    103  
     104
    104105  measure = &catalog[cat].measure[meas];
    105106
    … …  
    110111    if (measure[0].t < start[i]) continue;
    111112    if (measure[0].t > stop[i]) continue;
    112    
     113
    113114    // index for (catalog, measure) -> image
    114115    bin[cat][meas] = i;
    … …  
    125126      REALLOCATE (clist[i], int, NLIST[i]);
    126127      REALLOCATE (mlist[i], int, NLIST[i]);
    127     }   
     128    }
    128129    return;
    129130  }
    … …  
    157158
    158159  plot_defaults (&graphdata);
    159   graphdata.ymin = PlotdMmin; 
     160  graphdata.ymin = PlotdMmin;
    160161  graphdata.ymax = PlotdMmax;
    161162  plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png");
    … …  
    163164  plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL);
    164165
    165 # define NBIN 200
     166# define NBIN 20
    166167  REALLOCATE (xlist, double, NBIN);
    167168  REALLOCATE (Mlist, double, NBIN);
    … …  
    191192  int i, m, c, n;
    192193  double X, Y, L, M, P, Q, R, D, dR, dD;
    193  
     194
    194195  Mosaic *mosaic;
    195196  Coords *moscoords, *imcoords;
    196  
     197
    197198  moscoords = NULL;
    198199  mosaic = getMosaicForImage (im);
    … …  
    224225    dR = 3600.0*(catalog[c].average[n].R - R);
    225226    dD = 3600.0*(catalog[c].average[n].D - D);
    226    
     227
    227228    if (fabs(catalog[c].measure[m].dR - dR) > 10.0) {
    228229      // XXXXX running into this still for last megacam exposure: wrong mosaic?
    … …  
    249250      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    250251    }
    251   } 
     252  }
    252253  return;
    253254}
    … …  
    259260
    260261  int i, m, c, n;
    261  
     262
    262263  Mosaic *mosaic;
    263264  Coords *moscoords;
    264265  StarData *raw;
    265  
     266
    266267  ALLOCATE (raw, StarData, Nlist[im]);
    267268
    … …  
    271272      mosaic = getMosaicForImage (im);
    272273      if (mosaic == NULL) {
    273         fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
    274         exit (1);
     274        fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
     275        exit (1);
    275276      }
    276277      moscoords = &mosaic[0].coords;
    … …  
    284285    raw[i].X = catalog[c].measure[m].Xccd;
    285286    raw[i].Y = catalog[c].measure[m].Yccd;
    286    
     287
    287288    raw[i].Mag  = catalog[c].measure[m].M;
    288289    raw[i].dMag = catalog[c].measure[m].dM;
    … …  
    292293
    293294    // an object with only one detection provides no information about the image calibration
     295    //XXX this is already taken care of in bcatalog
    294296    raw[i].mask = FALSE;
    295     if (catalog[c].average[n].Nmeasure < 2) {
    296       raw[i].mask = TRUE;
    297     }
     297    int mask = FALSE;
     298    if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
     299      mask = TRUE;
     300    }
     301    if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
     302      mask = TRUE;
     303    }
     304
     305    raw[i].mask = mask;
     306
    298307
    299308    switch (mode) {
    300309      case MODE_SIMPLE:
    301         /* note that for a Simple image, L,M = P,Q */
    302         XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
    303         raw[i].P = raw[i].L;
    304         raw[i].Q = raw[i].M;
    305         LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
    306         break;
     310        /* note that for a Simple image, L,M = P,Q */
     311        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     312        raw[i].P = raw[i].L;
     313        raw[i].Q = raw[i].M;
     314        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
     315        break;
    307316      case MODE_MOSAIC:
    308         XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
    309         XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
    310         LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
    311         break;
     317        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     318        XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
     319        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
     320        break;
    312321    default:
    313322      fprintf (stderr, "error: invalid mode in getImageRaw");
    314323      abort ();
    315324    }
    316   } 
     325  }
    317326
    318327  *Nstars = Nlist[im];
    … …  
    330339  Coords *moscoords;
    331340  StarData *ref;
    332  
     341
    333342  ALLOCATE (ref, StarData, Nlist[im]);
    334343
    … …  
    352361    ref[i].R = catalog[c].average[n].R;
    353362    ref[i].D = catalog[c].average[n].D;
    354    
     363
    355364    ref[i].Mag  = catalog[c].measure[m].M;
    356365    ref[i].dMag = catalog[c].measure[m].dM;
    … …  
    373382      break;
    374383      default:
    375         fprintf (stderr, "invalid case");
    376         abort();
    377     }
    378   }
    379  
     384        fprintf (stderr, "invalid case");
     385        abort();
     386    }
     387  }
     388
    380389  *Nstars = Nlist[im];
    381390  return (ref);
    382391}
     392
     393/** lifted from relphot/StarOps.clean_measures */
     394void FlagOutliers2D(Catalog *catalog);
     395
     396void FlagOutliers (Catalog *catalog) {
     397  // XXX FlagOutliers is just using FlagOutliers2D
     398  FlagOutliers2D(catalog);
     399  return;
     400
     401  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
     402  double Ns, theta, x, y;
     403  double *R, *D, *dR, *dD;
     404  StatType statsR, statsD;
     405
     406  Nsecfilt = GetPhotcodeNsecfilt();
     407  assert(catalog[0].Nsecfilt == Nsecfilt);
     408
     409  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
     410  Nmax = 0;
     411  for (i = 0; i < catalog[0].Naverage; i++) {
     412    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
     413  }
     414
     415  ALLOCATE (R, double, Nmax);
     416  ALLOCATE (D, double, Nmax);
     417  ALLOCATE (dR, double, Nmax);
     418  ALLOCATE (dD, double, Nmax);
     419
     420  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
     421  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
     422
     423  Ns = CLIP_THRESH;
     424  Ndel = Nave = 0;
     425     
     426  /* loop over each object in the catalog */
     427  for (j = 0; j < catalog[0].Naverage; j++) {
     428   
     429    /* accumulate list of valid measurements */
     430    m = catalog[0].average[j].measureOffset;
     431    N = 0;
     432    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     433     
     434      // skip measurements based on user selected criteria
     435      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
     436      R[N] = catalog[0].measure[m].dR;
     437      D[N] = catalog[0].measure[m].dD;
     438      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
     439      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
     440      if (isnan(R[N]) || isnan(D[N])) continue;
     441      N++;
     442    }
     443    if (N <= TOOFEW) continue;
     444   
     445    /* 3-sigma clip based on stats of inner 50% */
     446    initstats ("MEAN");
     447    liststats (R, dR, N, &statsR);
     448    liststats (D, dD, N, &statsD);
     449   
     450    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
     451    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
     452   
     453    /* compare per-object distance to this standard deviation, and flag outliers*/
     454    m = catalog[0].average[j].measureOffset;
     455    N = 0;
     456    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     457      //reset flag on each invocation
     458      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     459
     460      // skip measurements based on user selected criteria
     461      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
     462     
     463      x = catalog[0].measure[m].dR - statsR.median;
     464      y = catalog[0].measure[m].dD - statsD.median;
     465      theta = atan2(y,x);
     466      if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) +
     467                         SQR(statsD.sigma * Ns * sin(theta)))) {   
     468        catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
     469        Ndel++;
     470      }
     471      N++;
     472      Nave ++;
     473    }
     474
     475    //examine results
     476    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset,
     477                               catalog[0].average[j].Nmeasure,
     478                               statsR, statsD, Ns);
     479  }
     480 
     481  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
     482  free (R);
     483  free(dR);
     484  free(D);
     485  free(dD);
     486}
     487
     488
     489/** an alternative outlier rejection scheme */
     490void FlagOutliers2D (Catalog *catalog) {
     491
     492  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
     493  double *index;
     494  double Ns, theta, x, y;
     495  double *R, *D, *dR, *dD, *d2;
     496  StatType statsR, statsD;
     497
     498  Nsecfilt = GetPhotcodeNsecfilt();
     499  assert(catalog[0].Nsecfilt == Nsecfilt);
     500
     501  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
     502  Nmax = 0;
     503  for (i = 0; i < catalog[0].Naverage; i++) {
     504    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
     505  }
     506
     507  ALLOCATE (R, double, Nmax);
     508  ALLOCATE (D, double, Nmax);
     509  ALLOCATE (dR, double, Nmax);
     510  ALLOCATE (dD, double, Nmax);
     511  ALLOCATE (d2, double, Nmax);
     512  ALLOCATE (index, double, Nmax);
     513
     514  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
     515  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
     516
     517  Ns = CLIP_THRESH;
     518  Ndel = Nave = 0;
     519     
     520  /* loop over each object in the catalog */
     521  for (j = 0; j < catalog[0].Naverage; j++) {
     522   
     523    /* accumulate list of valid measurements */
     524    m = catalog[0].average[j].measureOffset;
     525    N = 0;
     526   
     527    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     528     
     529      //reset flag on each invocation
     530      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     531     
     532      // skip measurements based on user selected criteria
     533      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
     534      R[N] = catalog[0].measure[m].dR;
     535      D[N] = catalog[0].measure[m].dD;
     536      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
     537      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
     538      if (isnan(R[N]) || isnan(D[N])) continue;
     539      N++;
     540    }
     541    if (N <= TOOFEW) continue;
     542   
     543    /* calculate mean of all points*/
     544    initstats ("MEAN");
     545    liststats (R, dR, N, &statsR);
     546    liststats (D, dD, N, &statsD);
     547    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
     548    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
     549   
     550    /* calculate deviations of all points*/
     551    m = catalog[0].average[j].measureOffset;
     552    N = 0;
     553    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     554      //skip bad measurements
     555      if (!MeasFilterTest(&catalog[0].measure[m])) continue; 
     556      x = catalog[0].measure[m].dR - statsR.median;
     557      y = catalog[0].measure[m].dD - statsD.median;
     558      theta = atan2(y,x);
     559      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) +
     560                             SQR(statsD.sigma * Ns * sin(theta)));     
     561      index[N] = m;
     562      N++;
     563    }
     564   
     565    // sort d2
     566    dsortpair(d2, index, N);
     567    N = (N/2 > (N-1)) ? N/2 : N-1;
     568
     569    // recalculate image center, sigma based on closest 50% of points
     570    for(k = 0;  k < N; k++) {
     571      int ind = (int) index[k];
     572      R[k] = catalog[0].measure[ind].dR;
     573      D[k] = catalog[0].measure[ind].dD;
     574      dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA);
     575      dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC);
     576    }
     577    liststats (R, dR, N, &statsR);
     578    liststats (D, dD, N, &statsD);
     579    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
     580    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
     581   
     582    // use these new statistics to flag outliers
     583    m = catalog[0].average[j].measureOffset;
     584    N = 0;
     585    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     586      //skip bad measurements
     587      if (!MeasFilterTest(&catalog[0].measure[m])) continue; 
     588      x = catalog[0].measure[m].dR - statsR.median;
     589      y = catalog[0].measure[m].dD - statsD.median;
     590      theta = atan2(y,x);
     591      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) +
     592                             SQR(statsD.sigma * Ns * sin(theta)));     
     593      if ((d2[N]) > 1) {
     594        catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
     595        Ndel ++;
     596      }
     597      N++;
     598      Nave++;
     599    }  //done rejecting outliers
     600
     601    //examine results
     602    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset,
     603                               catalog[0].average[j].Nmeasure,
     604                               statsR, statsD, Ns);
     605   
     606  } //done looping over objects
     607 
     608  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
     609  free (R);
     610  free(dR);
     611  free(D);
     612  free(dD);
     613  free(d2);
     614  free(index);
     615}
     616
     617
     618
     619/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */
     620int MeasFilterTest(Measure *measure) {
     621  int found, k;
     622  long mask;
     623  PhotCode *code;
     624  float mag;
     625
     626  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
     627  if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
     628 
     629  /* select measurements by photcode, or equiv photcode, if specified */
     630  if (NphotcodesKeep > 0) {
     631    found = FALSE;
     632    for (k = 0; (k < NphotcodesKeep) && !found; k++) {
     633      if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
     634      if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     635    }
     636    if (!found) return FALSE;
     637  }
     638 
     639  if (NphotcodesSkip > 0) {
     640    found = FALSE;
     641    for (k = 0; (k < NphotcodesSkip) && !found; k++) {
     642      if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
     643      if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     644    }
     645    if (found) return FALSE;
     646  } 
     647 
     648  /* select measurements by time */
     649  if (TimeSelect) {
     650    if (measure[0].t < TSTART) return FALSE;
     651    if (measure[0].t > TSTOP) return FALSE;
     652  }
     653 
     654  /* select measurements by quality */
     655  if (PhotFlagSelect) {
     656    if (PhotFlagBad) {
     657      mask = PhotFlagBad;
     658    } else {
     659      code = GetPhotcodebyCode (measure[0].photcode);
     660      mask = code[0].astromBadMask;
     661    }
     662    if (mask & measure[0].photFlags) return FALSE;
     663  }
     664
     665  /* select measurements by measurement error */
     666  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
     667 
     668  /* select measurements by mag limit */
     669  if (ImagSelect) {
     670    mag = PhotInst (measure);
     671    if (mag < ImagMin || mag > ImagMax) return FALSE;
     672  }
     673 
     674  return TRUE;
     675}

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r23595 r24308  
    9797      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    9898
    99         //exclude measurements which have non-finite astrometry
    100         if (!finite(catalog[i].measure[m].dR) || !finite(catalog[i].measure[m].dD)) continue;
     99        //does the measurement pass the supplied filtering constraints?
     100        if (!MeasFilterTest(&catalog[i].measure[m])) {
     101          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
     102          continue;
     103        }
     104
     105        //outlier rejection
     106        if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
     107          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
     108          continue;
     109        }
    101110
    102111        // exclude measurements by previous outlier detection
    … …  
    109118        # endif
    110119
    111         /* exclude measurements by quality */
    112         if (PhotFlagSelect) {
    113           if (PhotFlagBad) {
    114             mask = PhotFlagBad;
    115           } else {
    116             code = GetPhotcodebyCode (catalog[i].measure[m].photcode);
    117             mask = code[0].astromBadMask;
    118           }
    119           if (mask & catalog[i].measure[m].photFlags) {
    120             catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
    121             continue;
    122           }
    123         }
    124        
    125         /* exclude measurements by mag limit */
    126         if (ImagSelect) {
    127           mag = PhotInst (&catalog[i].measure[m]);
    128           if (mag < ImagMin) {
    129             catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
    130             continue;
    131           }
    132           if (mag > ImagMax) {
    133             catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
    134             continue;
    135           }
    136         }
    137 
    138         /* select or exclude measurements by photcode, or equiv photcode, if specified */
    139         if (NphotcodesKeep > 0) {
    140           found = FALSE;
    141           for (kp = 0; (kp < NphotcodesKeep) && !found; kp++) {
    142             if (photcodesKeep[kp][0].code == catalog[i].measure[m].photcode) found = TRUE;
    143             if (photcodesKeep[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE;
    144           }
    145           if (!found) continue;
    146         }
    147         if (NphotcodesSkip > 0) {
    148           found = FALSE;
    149           for (kp = 0; (kp < NphotcodesSkip) && !found; kp++) {
    150             if (photcodesSkip[kp][0].code == catalog[i].measure[m].photcode) found = TRUE;
    151             if (photcodesSkip[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE;
    152           }
    153           if (found) continue;
    154         }
     120        catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
    155121
    156122        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    … …  
    179145
    180146      // too few measurements for average position (require 2 values)
    181       if (N < 2) {
     147      if (N < SRC_MEAS_TOOFEW) {
    182148        // XXX need to define PHOTOM and ASTROM object flags
    183         // catalog[i].average[j].code |= ID_STAR_FEW;
    184         continue;
     149        catalog[i].average[j].flags |= ID_STAR_FEW;
     150        if (N < 2) continue;
    185151      }
    186152
    … …  
    276242      }
    277243
     244      //make sure that the fit succeeded
     245      assert(finite(fit.Ro) && finite(fit.Do) &&
     246             finite(fit.dRo) && finite(fit.dDo) &&
     247             finite(fit.uR) && finite(fit.uD) &&
     248             finite(fit.duR) && finite(fit.duD) &&
     249             finite(fit.p) && finite(fit.dp));
     250
    278251      // the measure fields must be updated before the average fields
    279252      m = catalog[i].average[j].measureOffset;

  • trunk/Ohana/src/relastro/src/args.c

    r17205 r24308  
    4040    FIT_TARGET = TARGET_MOSAICS;
    4141  }
     42
     43  FlagOutlier = FALSE;
     44  if ((N = get_argument (argc, argv, "-clip"))) {
     45    remove_argument (N, &argc, argv);
     46    CLIP_THRESH = atof (argv[N]);
     47    remove_argument (N, &argc, argv);
     48    FlagOutlier = TRUE;
     49  }
     50
    4251  if (FIT_TARGET == TARGET_NONE) usage();
    4352

  • trunk/Ohana/src/relastro/src/bcatalog.c

    r21508 r24308  
    22
    33int bcatalog (Catalog *subcatalog, Catalog *catalog) {
    4  
     4
    55  int i, j, k, offset, found;
    66  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
    … …  
    2323  /* exclude stars not in range or with too few measurements */
    2424  for (i = 0; i < catalog[0].Naverage; i++) {
    25     if (catalog[0].average[i].Nmeasure < 2) continue;
    26 
     25    if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
     26      continue;
     27    }
     28   
    2729    /* start with all stars good */
    2830    subcatalog[0].average[Naverage] = catalog[0].average[i];
    … …  
    4446
    4547      offset = catalog[0].average[i].measureOffset + j;
    46 
    47       /* select measurements by photcode, or equiv photcode, if specified */
    48       if (NphotcodesKeep > 0) {
    49         found = FALSE;
    50         for (k = 0; (k < NphotcodesKeep) && !found; k++) {
    51           if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
    52           if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
    53         }
    54         if (!found) continue;
    55       }
    56       if (NphotcodesSkip > 0) {
    57         found = FALSE;
    58         for (k = 0; (k < NphotcodesSkip) && !found; k++) {
    59           if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
    60           if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
    61         }
    62         if (found) continue;
     48     
     49      //filter objects based on user supplied criteria
     50      if (!MeasFilterTest(&catalog[0].measure[offset])) {
     51        catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
     52        continue;
    6353      }
    6454
    65       /* select measurements by time */
    66       if (TimeSelect) {
    67         if (catalog[0].measure[offset].t < TSTART) continue;
    68         if (catalog[0].measure[offset].t > TSTOP) continue;
     55      //filter out outliers
     56      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
     57        catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
     58        continue;
    6959      }
    70 
    71       /* select measurements by quality */
    72       // XXX FIX THIS!!
    73       // if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
    74 
    75       /* select measurements by quality */
    76       if (PhotFlagSelect) {
    77         if (PhotFlagBad) {
    78           mask = PhotFlagBad;
    79         } else {
    80           code = GetPhotcodebyCode (catalog[0].measure[offset].photcode);
    81           mask = code[0].astromBadMask;
    82         }
    83         if (mask & catalog[0].measure[offset].photFlags) continue;
    84       }
    85 
    86       /* select measurements by measurement error */
    87       if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue;
    88 
    89       /* select measurements by mag limit */
    90       if (ImagSelect) {
    91         mag = PhotInst (&catalog[0].measure[offset]);
    92         if (mag < ImagMin) continue;
    93         if (mag > ImagMax) continue;
    94       }
     60      catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;
    9561
    9662      // re-assess on each run of relastro if a measurement should be used
    … …  
    11076      subcatalog[0].measure[Nmeasure]          = catalog[0].measure[offset];
    11177      subcatalog[0].measure[Nmeasure].averef   = Naverage;
    112       if (RESET) { 
    113         subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    114         subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
     78      if (RESET) {
     79        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     80        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
    11581      }
     82     
    11683      Nmeasure ++;
    11784      Nm ++;
    11885      if (Nmeasure == NMEASURE) {
    119         NMEASURE += 1000;
    120         REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
     86        NMEASURE += 1000;
     87        REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
    12188      }
    12289    }
    … …  
    139106
    140107  if (VERBOSE) {
    141     fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, 
    142              subcatalog[0].Naverage, subcatalog[0].Nmeasure);
     108    fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
     109             subcatalog[0].Naverage, subcatalog[0].Nmeasure);
    143110   }
    144111  return (TRUE);

  • trunk/Ohana/src/relastro/src/initialize.c

    r17213 r24308  
    2525      fprintf (stderr, "codename: %s\n", codename);
    2626      if ((photcodesKeep[NphotcodesKeep] = GetPhotcodebyName (codename)) == NULL) {
    27         fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
    28         exit (1);
     27        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
     28        exit (1);
    2929      }
    3030      NphotcodesKeep ++;
    … …  
    4646      fprintf (stderr, "codename: %s\n", codename);
    4747      if ((photcodesSkip[NphotcodesSkip] = GetPhotcodebyName (codename)) == NULL) {
    48         fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
    49         exit (1);
     48        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
     49        exit (1);
    5050      }
    5151      NphotcodesSkip ++;

  • trunk/Ohana/src/relastro/src/load_catalogs.c

    r21508 r24308  
    2929    }
    3030    if (VERBOSE && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
     31
     32    //outlier rejection
     33    if (FlagOutlier) {
     34      FlagOutliers(&tcatalog);
     35    }
    3136
    3237    // select only the brighter stars

  • trunk/Ohana/src/relastro/src/relastro_objects.c

    r21508 r24308  
    88  SkyList *skylist = NULL;
    99  Catalog catalog;
     10
    1011
    1112  // load the current sky table (layout of all SkyRegions)
    … …  
    3738    }
    3839
     40    if (FlagOutlier) {
     41      FlagOutliers(&catalog);
     42    }
     43
    3944    // XXX consider what gets reset (only ASTROM flags)
    4045    if (RESET) {
Note: See TracChangeset for help on using the changeset viewer.