IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 23649


Ignore:
Timestamp:
Mar 31, 2009, 8:37:27 PM (17 years ago)
Author:
beaumont
Message:

added outlier detection to relastro

Location:
branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h

    r23344 r23649  
    120120int TimeSelect;
    121121time_t TSTART, TSTOP;
     122
     123int FlagOutlier;
    122124
    123125FitMode FIT_MODE;
    … …  
    266268int UpdateMeasures (Catalog *catalog, int Ncatalog);
    267269void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
     270void FlagOutliers(Catalog *catalog, int Ncatalog);
     271int MeasFilterTest(Measure *measure);
    268272
    269273int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c

    r23594 r23649  
    386386  return (ref);
    387387}
     388
     389/** lifted from relphot/StarOps.clean_measures */
     390void FlagOutliers (Catalog *catalog, int Ncatalog) {
     391
     392  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
     393  double Ns;
     394  double *R, *D, *dR, *dD;
     395  StatType statsR, statsD;
     396
     397  Nsecfilt = GetPhotcodeNsecfilt();
     398  assert(catalog[0].Nsecfilt == Nsecfilt);
     399
     400  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
     401  /* Nmeasure is now different, need to reallocate */
     402  Nmax = 0;
     403  for (i = 0; i < Ncatalog; i++) {
     404    for (j = 0; j < catalog[i].Naverage; j++) {
     405      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
     406    }
     407  }
     408  ALLOCATE (R, double, Nmax);
     409  ALLOCATE (D, double, Nmax);
     410  ALLOCATE (dR, double, Nmax);
     411  ALLOCATE (dD, double, Nmax);
     412
     413  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
     414  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
     415
     416  Ns = 3;
     417  Ndel = Nave = 0;
     418  for (i = 0; i < Ncatalog; i++) {
     419    for (j = 0; j < catalog[i].Naverage; j++) {
     420     
     421      /* accumulate list of valid measurements */
     422      m = catalog[i].average[j].measureOffset;
     423      N = 0;
     424      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
     425
     426        // skip measurements based on user selected criteria
     427        if (!MeasFilterTest(&catalog[i].measure[m])) continue;
     428
     429        R[N] = getMeanR(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     430        D[N] = getMeanD(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     431        dR[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_RA);
     432        dD[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_DEC);
     433        if (isnan(R[N]) || isnan(D[N])) continue;
     434        N++;
     435      }
     436      if (N <= TOOFEW) continue;
     437
     438      /* 3-sigma clip based on stats of inner 50% */
     439      initstats ("INNER_MEAN");
     440      liststats (R, dR, N, &statsR);
     441      liststats (R, dR, N, &statsD);
     442
     443      statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
     444      statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
     445
     446      for (k = m = 0; k < N; k++) {
     447        if ((fabs (R[k] - statsR.median) < Ns*statsR.sigma) &&
     448            (fabs (D[k] - statsD.median) < Ns*statsD.sigma)) {
     449          R[m] = R[k];
     450          dR[m] = dR[k];
     451          D[m] = D[k];
     452          dD[m] = dD[k];
     453          m++;
     454        }
     455      }
     456      initstats ("MEAN");
     457      liststats (R, dR, m, &statsR);
     458      liststats (D, dD, m, &statsD);
     459      statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
     460      statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
     461
     462      /* apply to list of all relevant measurements*/
     463      m = catalog[i].average[j].measureOffset;
     464      N = 0;
     465      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
     466        // skip measurements based on user selected criteria
     467        if (!MeasFilterTest(&catalog[i].measure[m])) continue;
     468
     469        R[N] = getMeanR(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     470        D[N] = getMeanD(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     471        dR[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_RA);
     472        dD[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_DEC);
     473       
     474        if((fabs(R[N] - statsR.median) > Ns * statsR.sigma) ||
     475           (fabs(D[N] - statsR.median) > Ns * statsD.sigma)) {
     476          catalog[i].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
     477          Ndel++;
     478        }
     479        N++;
     480        Nave ++;
     481      }
     482      initstats (STATMODE);
     483   
     484      if (TRUE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
     485      free (R);
     486      free(dR);
     487      free(D);
     488      free(dD);
     489    }
     490  }
     491}
     492
     493 
     494/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */
     495int MeasFilterTest(Measure *measure) {
     496  int found, k;
     497  long mask;
     498  PhotCode *code;
     499  float mag;
     500
     501  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
     502  if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
     503
     504  /* select measurements by photcode, or equiv photcode, if specified */
     505  if (NphotcodesKeep > 0) {
     506    found = FALSE;
     507    for (k = 0; (k < NphotcodesKeep) && !found; k++) {
     508      if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
     509      if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     510    }
     511    if (!found) return FALSE;
     512  }
     513 
     514  if (NphotcodesSkip > 0) {
     515    found = FALSE;
     516    for (k = 0; (k < NphotcodesSkip) && !found; k++) {
     517      if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
     518      if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     519    }
     520    if (found) return FALSE;
     521  } 
     522 
     523  /* select measurements by time */
     524  if (TimeSelect) {
     525    if (measure[0].t < TSTART) return FALSE;
     526    if (measure[0].t > TSTOP) return FALSE;
     527  }
     528 
     529  /* select measurements by quality */
     530  if (PhotFlagSelect) {
     531    if (PhotFlagBad) {
     532      mask = PhotFlagBad;
     533    } else {
     534      code = GetPhotcodebyCode (measure[0].photcode);
     535      mask = code[0].astromBadMask;
     536    }
     537    if (mask & measure[0].photFlags) return FALSE;
     538  }
     539
     540  /* select measurements by measurement error */
     541  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
     542 
     543  /* select measurements by mag limit */
     544  if (ImagSelect) {
     545    mag = PhotInst (measure);
     546    if (mag < ImagMin || mag > ImagMax) return FALSE;
     547  }
     548 
     549  return TRUE;
     550}

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/UpdateObjects.c

    r23594 r23649  
    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])) continue;
     101
     102        //outlier rejection
     103        if (FlagOutliers && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
     104          continue;
     105        }
    101106
    102107        // exclude measurements by previous outlier detection
    … …  
    109114        # endif
    110115
    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         }
    155 
    156116        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    157117        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    … …  
    179139
    180140      // too few measurements for average position (require 2 values)
    181       if (N < 2) {
     141      if (N < SRC_MEAS_TOOFEW) {
    182142        // XXX need to define PHOTOM and ASTROM object flags
    183143        // catalog[i].average[j].code |= ID_STAR_FEW;

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/args.c

    r17205 r23649  
    3939    remove_argument (N, &argc, argv);
    4040    FIT_TARGET = TARGET_MOSAICS;
     41  }
     42  if ((N = get_argument (argc, argv, "-clip"))) {
     43    remove_argument (N, &argc, argv);
     44    FlagOutlier = TRUE;
    4145  }
    4246  if (FIT_TARGET == TARGET_NONE) usage();

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c

    r23594 r23649  
    88  int mask;
    99  PhotCode *code;
    10   int skipFew = skipPhotCodeKeep = skipPhotCodeSkip = skipTime = skipPhotFlag = skipSigmaLim = skipImagSelect = 0;
    1110
    1211  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
    … …  
    2524  for (i = 0; i < catalog[0].Naverage; i++) {
    2625    if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
    27       skipFew++;
    2826      continue;
    2927    }
    … …  
    4846
    4947      offset = catalog[0].average[i].measureOffset + j;
     48     
     49      //filter objects based on user supplied criteria
     50      if (!MeasFilterTest(&catalog[0].measure[offset])) continue;
    5051
    51       /* select measurements by photcode, or equiv photcode, if specified */
    52       if (NphotcodesKeep > 0) {
    53         found = FALSE;
    54         for (k = 0; (k < NphotcodesKeep) && !found; k++) {
    55           if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
    56           if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
    57         }
    58         if (!found) {
    59           skipPhotCodeKeep++;
    60           continue;
    61         }
    62       }
    63       if (NphotcodesSkip > 0) {
    64         found = FALSE;
    65         for (k = 0; (k < NphotcodesSkip) && !found; k++) {
    66           if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
    67           if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
    68         }
    69         if (found) {
    70           skipPhotCodeSkip++;
    71           continue;
    72         }
    73       }
    74 
    75       /* select measurements by time */
    76       if (TimeSelect) {
    77         if (catalog[0].measure[offset].t < TSTART) {
    78           skipTime++;
    79           continue;
    80         }
    81         if (catalog[0].measure[offset].t > TSTOP) {
    82           skipTime++;
    83           continue;
    84         }
    85       }
     52      //filter out outliers
     53      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) continue;
    8654     
    87       /* select measurements by quality */
    88       // XXX FIX THIS!!
    89       // if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
    90      
    91       /* select measurements by quality */
    92       if (PhotFlagSelect) {
    93         if (PhotFlagBad) {
    94           mask = PhotFlagBad;
    95         } else {
    96           code = GetPhotcodebyCode (catalog[0].measure[offset].photcode);
    97           mask = code[0].astromBadMask;
    98         }
    99         if (mask & catalog[0].measure[offset].photFlags) {
    100           skipPhotFlag++;
    101           continue;
    102         }
    103       }
    104      
    105       /* select measurements by measurement error */
    106       if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) {
    107         skipSigmaLim++;
    108         continue;
    109       }
    110      
    111       /* select measurements by mag limit */
    112       if (ImagSelect) {
    113         mag = PhotInst (&catalog[0].measure[offset]);
    114         if (mag < ImagMin || mag > ImagMax) {
    115           skipImagSelect++;
    116           continue;
    117         }
    118       }
    11955      // re-assess on each run of relastro if a measurement should be used
    12056
    … …  
    16298
    16399  if (VERBOSE) {
    164     fprintf(stderr, "Reasons for exclusion in bcatalog:\n");
    165     fprintf(stderr, "\ntoo few: %d \nphotCodeSkip: %d \nphotCodeKeep: %d \ntime: %d\n",
    166             skipfew, skipPhotCodeKeep, skipPhotCodeSkip, skipTime);
    167     fprintf(stderr, "photFlag: %d \nmagnitude error: %d \nImag: %d\n",
    168             skipPhotFlag, skipSigmaLim, skipImagSelect);
    169   }
    170 
    171   if (VERBOSE) {
    172100    fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
    173101             subcatalog[0].Naverage, subcatalog[0].Nmeasure);

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/load_catalogs.c

    r21508 r23649  
    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, 1);
     35    }
    3136
    3237    // select only the brighter stars

  • branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastro_objects.c

    r21508 r23649  
    3737    }
    3838
     39    if (FlagOutliers) {
     40      FlagOutliers(&catalog, 1);
     41    }
     42
    3943    // XXX consider what gets reset (only ASTROM flags)
    4044    if (RESET) {
Note: See TracChangeset for help on using the changeset viewer.