Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 24233) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 24308) @@ -108,6 +108,6 @@ # define ID_MEAS_POOR_ASTROM 0x0010 /* detection is astrometry outlier */ # define ID_MEAS_SKIP_ASTROM 0x0020 /* detection was ignored for astrometry measurement */ -# define ID_MEAS_UNDEF_1 0x0040 /* unused */ -# define ID_MEAS_UNDEF_2 0x0080 /* unused */ +# define ID_MEAS_USED_OBJ 0x0040 /* detection was used during opdate objects */ +# define ID_MEAS_USED_CHIP 0x0080 /* detection was used during update chips */ # define ID_MEAS_BLEND_MEAS 0x0100 /* detection is within radius of multiple objects */ # define ID_MEAS_BLEND_OBJ 0x0200 /* multiple detections within radius of object */ Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 24233) +++ trunk/Ohana/src/relastro/Makefile (revision 24308) @@ -55,5 +55,6 @@ $(SRC)/relastro_objects.$(ARCH).o \ $(SRC)/save_catalogs.$(ARCH).o \ -$(SRC)/write_coords.$(ARCH).o +$(SRC)/write_coords.$(ARCH).o \ +$(SRC)/relastroVisual.$(ARCH).o $(RELASTRO): $(INC)/relastro.h $(KAPA_INCS) Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 24233) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 24308) @@ -44,5 +44,5 @@ double Ro, dRo; double Do, dDo; - + double uR, duR; double uD, duD; @@ -63,5 +63,5 @@ char flags; Coords coords; -} Mosaic; +} Mosaic; typedef struct { @@ -87,4 +87,5 @@ double SIGMA_LIM; +int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored double MIN_ERROR; @@ -119,4 +120,7 @@ int TimeSelect; time_t TSTART, TSTOP; + +int FlagOutlier; +int CLIP_THRESH; FitMode FIT_MODE; @@ -145,5 +149,5 @@ int corner_check PROTO((double *x1, double *y1, double *x2, double *y2)); void dumpGrid PROTO((void)); -void dump_grid PROTO((void)); +void dump_grid PROTO((void)); int edge_check PROTO((double *x1, double *y1, double *x2, double *y2)); void findImages PROTO((Catalog *catalog, int Ncatalog)); @@ -220,5 +224,5 @@ int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int mark)); void set_ZP PROTO((double ZERO)); -int setrefcode PROTO((Image *image, int Nimage)); +int setrefcode PROTO((Image *image, int Nimage)); void skip_measurements PROTO((Catalog *catalog, int pass)); void sortA PROTO((double *X, int N)); @@ -265,4 +269,6 @@ int UpdateMeasures (Catalog *catalog, int Ncatalog); void fixImageRaw (Catalog *catalog, int Ncatalog, int im); +void FlagOutliers(Catalog *catalog); +int MeasFilterTest(Measure *measure); int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 24233) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 24308) @@ -20,4 +20,5 @@ GetConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW); // XXX these are used in relphot to identify poor stars / images -- define @@ -55,5 +56,5 @@ fprintf (stderr, "directory %s does not exist, giving up\n", CATDIR); exit (1); - } + } /* update master photcode table if not defined */ Index: trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- trunk/Ohana/src/relastro/src/FitChip.c (revision 24233) +++ trunk/Ohana/src/relastro/src/FitChip.c (revision 24308) @@ -1,3 +1,4 @@ # include "relastro.h" +# include "relastroVisual.h" // XXX make these user parameters @@ -5,4 +6,5 @@ # define FIT_CHIP_NITER 3 # define FIT_CHIP_NSIGMA 3.0 + // XXX we should test if the fit is sufficiently constrained across the chip, or if the @@ -22,5 +24,4 @@ ALLOCATE (values, double, Nmatch); - for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) { @@ -28,10 +29,10 @@ for (i = Nscatter = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - if (raw[i].dMag > FIT_CHIP_MAX_ERROR) continue; + if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue; dL = raw[i].L - ref[i].L; dM = raw[i].M - ref[i].M; dR = hypot (dL, dM); - + values[Nscatter] = dR; Nscatter++; @@ -41,4 +42,6 @@ dsort (values, Nscatter); dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)]; + relastroVisualPlotScatter(values, dRmax, Nscatter); + relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch); // fit the requested order polynomial @@ -50,5 +53,10 @@ // generate the fit matches for (i = 0; i < Nmatch; i++) { - if (raw[i].mask) continue; + if (raw[i].mask) { + continue; + } + if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) { + continue; + } // only keep objects within dRmax @@ -57,5 +65,5 @@ dR = hypot (dL, dM); if (dR > dRmax) continue; - + fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos); } @@ -66,15 +74,15 @@ case 0: case 1: - skip = (fit[0].Npts < 8); - break; + skip = (fit[0].Npts < 8); + break; case 2: - skip = (fit[0].Npts < 11); - break; + skip = (fit[0].Npts < 11); + break; case 3: - skip = (fit[0].Npts < 15); - break; + skip = (fit[0].Npts < 15); + break; default: - fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms); - abort (); + fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms); + abort (); } if (skip) { @@ -95,4 +103,5 @@ XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords); } + } Index: trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- trunk/Ohana/src/relastro/src/GetAstromError.c (revision 24233) +++ trunk/Ohana/src/relastro/src/GetAstromError.c (revision 24308) @@ -30,4 +30,8 @@ dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM)); + + //XXX dXccd, dYccd are now working correctly + //dPtotal = AS * dPobs; + dPtotal = MAX (dPtotal, MIN_ERROR); Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 24233) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 24308) @@ -1,3 +1,4 @@ # include "relastro.h" +#define TESTING static unsigned int *start; @@ -91,7 +92,7 @@ for (i = 0; VERBOSE && (i < Nimage); i++) { name = GetPhotcodeNamebyCode (image[i].photcode); - fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], - ohana_sec_to_date(image[i].tzero), name); - } + fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], + ohana_sec_to_date(image[i].tzero), name); + } } @@ -101,5 +102,5 @@ int i; Measure *measure; - + measure = &catalog[cat].measure[meas]; @@ -110,5 +111,5 @@ if (measure[0].t < start[i]) continue; if (measure[0].t > stop[i]) continue; - + // index for (catalog, measure) -> image bin[cat][meas] = i; @@ -125,5 +126,5 @@ REALLOCATE (clist[i], int, NLIST[i]); REALLOCATE (mlist[i], int, NLIST[i]); - } + } return; } @@ -157,5 +158,5 @@ plot_defaults (&graphdata); - graphdata.ymin = PlotdMmin; + graphdata.ymin = PlotdMmin; graphdata.ymax = PlotdMmax; plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png"); @@ -163,5 +164,5 @@ plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL); -# define NBIN 200 +# define NBIN 20 REALLOCATE (xlist, double, NBIN); REALLOCATE (Mlist, double, NBIN); @@ -191,8 +192,8 @@ int i, m, c, n; double X, Y, L, M, P, Q, R, D, dR, dD; - + Mosaic *mosaic; Coords *moscoords, *imcoords; - + moscoords = NULL; mosaic = getMosaicForImage (im); @@ -224,5 +225,5 @@ dR = 3600.0*(catalog[c].average[n].R - R); dD = 3600.0*(catalog[c].average[n].D - D); - + if (fabs(catalog[c].measure[m].dR - dR) > 10.0) { // XXXXX running into this still for last megacam exposure: wrong mosaic? @@ -249,5 +250,5 @@ catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); } - } + } return; } @@ -259,9 +260,9 @@ int i, m, c, n; - + Mosaic *mosaic; Coords *moscoords; StarData *raw; - + ALLOCATE (raw, StarData, Nlist[im]); @@ -271,6 +272,6 @@ mosaic = getMosaicForImage (im); if (mosaic == NULL) { - fprintf (stderr, "mosaic not found for image %s\n", image[im].name); - exit (1); + fprintf (stderr, "mosaic not found for image %s\n", image[im].name); + exit (1); } moscoords = &mosaic[0].coords; @@ -284,5 +285,5 @@ raw[i].X = catalog[c].measure[m].Xccd; raw[i].Y = catalog[c].measure[m].Yccd; - + raw[i].Mag = catalog[c].measure[m].M; raw[i].dMag = catalog[c].measure[m].dM; @@ -292,27 +293,35 @@ // an object with only one detection provides no information about the image calibration + //XXX this is already taken care of in bcatalog raw[i].mask = FALSE; - if (catalog[c].average[n].Nmeasure < 2) { - raw[i].mask = TRUE; - } + int mask = FALSE; + if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) { + mask = TRUE; + } + if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) { + mask = TRUE; + } + + raw[i].mask = mask; + switch (mode) { case MODE_SIMPLE: - /* note that for a Simple image, L,M = P,Q */ - XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); - raw[i].P = raw[i].L; - raw[i].Q = raw[i].M; - LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords); - break; + /* note that for a Simple image, L,M = P,Q */ + XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); + raw[i].P = raw[i].L; + raw[i].Q = raw[i].M; + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords); + break; case MODE_MOSAIC: - XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); - XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords); - LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords); - break; + XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); + XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords); + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords); + break; default: fprintf (stderr, "error: invalid mode in getImageRaw"); abort (); } - } + } *Nstars = Nlist[im]; @@ -330,5 +339,5 @@ Coords *moscoords; StarData *ref; - + ALLOCATE (ref, StarData, Nlist[im]); @@ -352,5 +361,5 @@ ref[i].R = catalog[c].average[n].R; ref[i].D = catalog[c].average[n].D; - + ref[i].Mag = catalog[c].measure[m].M; ref[i].dMag = catalog[c].measure[m].dM; @@ -373,10 +382,294 @@ break; default: - fprintf (stderr, "invalid case"); - abort(); - } - } - + fprintf (stderr, "invalid case"); + abort(); + } + } + *Nstars = Nlist[im]; return (ref); } + +/** lifted from relphot/StarOps.clean_measures */ +void FlagOutliers2D(Catalog *catalog); + +void FlagOutliers (Catalog *catalog) { + // XXX FlagOutliers is just using FlagOutliers2D + FlagOutliers2D(catalog); + return; + + int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt; + double Ns, theta, x, y; + double *R, *D, *dR, *dD; + StatType statsR, statsD; + + Nsecfilt = GetPhotcodeNsecfilt(); + assert(catalog[0].Nsecfilt == Nsecfilt); + + if (VERBOSE) fprintf (stderr, "marking poor measures\n"); + Nmax = 0; + for (i = 0; i < catalog[0].Naverage; i++) { + Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure); + } + + ALLOCATE (R, double, Nmax); + ALLOCATE (D, double, Nmax); + ALLOCATE (dR, double, Nmax); + ALLOCATE (dD, double, Nmax); + + /* it makes no sense to mark 3-sigma outliers with <5 measurements */ + TOOFEW = MAX (5, SRC_MEAS_TOOFEW); + + Ns = CLIP_THRESH; + Ndel = Nave = 0; + + /* loop over each object in the catalog */ + for (j = 0; j < catalog[0].Naverage; j++) { + + /* accumulate list of valid measurements */ + m = catalog[0].average[j].measureOffset; + N = 0; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + + // skip measurements based on user selected criteria + if (!MeasFilterTest(&catalog[0].measure[m])) continue; + R[N] = catalog[0].measure[m].dR; + D[N] = catalog[0].measure[m].dD; + dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA); + dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC); + if (isnan(R[N]) || isnan(D[N])) continue; + N++; + } + if (N <= TOOFEW) continue; + + /* 3-sigma clip based on stats of inner 50% */ + initstats ("MEAN"); + liststats (R, dR, N, &statsR); + liststats (D, dD, N, &statsD); + + statsR.sigma = MAX (MIN_ERROR, statsR.sigma); + statsD.sigma = MAX (MIN_ERROR, statsD.sigma); + + /* compare per-object distance to this standard deviation, and flag outliers*/ + m = catalog[0].average[j].measureOffset; + N = 0; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + //reset flag on each invocation + catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM; + + // skip measurements based on user selected criteria + if (!MeasFilterTest(&catalog[0].measure[m])) continue; + + x = catalog[0].measure[m].dR - statsR.median; + y = catalog[0].measure[m].dD - statsD.median; + theta = atan2(y,x); + if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + + SQR(statsD.sigma * Ns * sin(theta)))) { + catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM; + Ndel++; + } + N++; + Nave ++; + } + + //examine results + relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, + catalog[0].average[j].Nmeasure, + statsR, statsD, Ns); + } + + if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); + free (R); + free(dR); + free(D); + free(dD); +} + + +/** an alternative outlier rejection scheme */ +void FlagOutliers2D (Catalog *catalog) { + + int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt; + double *index; + double Ns, theta, x, y; + double *R, *D, *dR, *dD, *d2; + StatType statsR, statsD; + + Nsecfilt = GetPhotcodeNsecfilt(); + assert(catalog[0].Nsecfilt == Nsecfilt); + + if (VERBOSE) fprintf (stderr, "marking poor measures\n"); + Nmax = 0; + for (i = 0; i < catalog[0].Naverage; i++) { + Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure); + } + + ALLOCATE (R, double, Nmax); + ALLOCATE (D, double, Nmax); + ALLOCATE (dR, double, Nmax); + ALLOCATE (dD, double, Nmax); + ALLOCATE (d2, double, Nmax); + ALLOCATE (index, double, Nmax); + + /* it makes no sense to mark 3-sigma outliers with <5 measurements */ + TOOFEW = MAX (5, SRC_MEAS_TOOFEW); + + Ns = CLIP_THRESH; + Ndel = Nave = 0; + + /* loop over each object in the catalog */ + for (j = 0; j < catalog[0].Naverage; j++) { + + /* accumulate list of valid measurements */ + m = catalog[0].average[j].measureOffset; + N = 0; + + for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + + //reset flag on each invocation + catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM; + + // skip measurements based on user selected criteria + if (!MeasFilterTest(&catalog[0].measure[m])) continue; + R[N] = catalog[0].measure[m].dR; + D[N] = catalog[0].measure[m].dD; + dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA); + dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC); + if (isnan(R[N]) || isnan(D[N])) continue; + N++; + } + if (N <= TOOFEW) continue; + + /* calculate mean of all points*/ + initstats ("MEAN"); + liststats (R, dR, N, &statsR); + liststats (D, dD, N, &statsD); + statsR.sigma = MAX (MIN_ERROR, statsR.sigma); + statsD.sigma = MAX (MIN_ERROR, statsD.sigma); + + /* calculate deviations of all points*/ + m = catalog[0].average[j].measureOffset; + N = 0; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + //skip bad measurements + if (!MeasFilterTest(&catalog[0].measure[m])) continue; + x = catalog[0].measure[m].dR - statsR.median; + y = catalog[0].measure[m].dD - statsD.median; + theta = atan2(y,x); + d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + + SQR(statsD.sigma * Ns * sin(theta))); + index[N] = m; + N++; + } + + // sort d2 + dsortpair(d2, index, N); + N = (N/2 > (N-1)) ? N/2 : N-1; + + // recalculate image center, sigma based on closest 50% of points + for(k = 0; k < N; k++) { + int ind = (int) index[k]; + R[k] = catalog[0].measure[ind].dR; + D[k] = catalog[0].measure[ind].dD; + dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA); + dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC); + } + liststats (R, dR, N, &statsR); + liststats (D, dD, N, &statsD); + statsR.sigma = MAX (MIN_ERROR, statsR.sigma); + statsD.sigma = MAX (MIN_ERROR, statsD.sigma); + + // use these new statistics to flag outliers + m = catalog[0].average[j].measureOffset; + N = 0; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + //skip bad measurements + if (!MeasFilterTest(&catalog[0].measure[m])) continue; + x = catalog[0].measure[m].dR - statsR.median; + y = catalog[0].measure[m].dD - statsD.median; + theta = atan2(y,x); + d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + + SQR(statsD.sigma * Ns * sin(theta))); + if ((d2[N]) > 1) { + catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM; + Ndel ++; + } + N++; + Nave++; + } //done rejecting outliers + + //examine results + relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, + catalog[0].average[j].Nmeasure, + statsR, statsD, Ns); + + } //done looping over objects + + if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); + free (R); + free(dR); + free(D); + free(dD); + free(d2); + free(index); +} + + + +/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ +int MeasFilterTest(Measure *measure) { + int found, k; + long mask; + PhotCode *code; + float mag; + + if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE; + if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks? + + /* select measurements by photcode, or equiv photcode, if specified */ + if (NphotcodesKeep > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + if (!found) return FALSE; + } + + if (NphotcodesSkip > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesSkip) && !found; k++) { + if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE; + if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + if (found) return FALSE; + } + + /* select measurements by time */ + if (TimeSelect) { + if (measure[0].t < TSTART) return FALSE; + if (measure[0].t > TSTOP) return FALSE; + } + + /* select measurements by quality */ + if (PhotFlagSelect) { + if (PhotFlagBad) { + mask = PhotFlagBad; + } else { + code = GetPhotcodebyCode (measure[0].photcode); + mask = code[0].astromBadMask; + } + if (mask & measure[0].photFlags) return FALSE; + } + + /* select measurements by measurement error */ + if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE; + + /* select measurements by mag limit */ + if (ImagSelect) { + mag = PhotInst (measure); + if (mag < ImagMin || mag > ImagMax) return FALSE; + } + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 24233) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 24308) @@ -97,6 +97,15 @@ for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - //exclude measurements which have non-finite astrometry - if (!finite(catalog[i].measure[m].dR) || !finite(catalog[i].measure[m].dD)) continue; + //does the measurement pass the supplied filtering constraints? + if (!MeasFilterTest(&catalog[i].measure[m])) { + catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; + continue; + } + + //outlier rejection + if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) { + catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; + continue; + } // exclude measurements by previous outlier detection @@ -109,48 +118,5 @@ # endif - /* exclude measurements by quality */ - if (PhotFlagSelect) { - if (PhotFlagBad) { - mask = PhotFlagBad; - } else { - code = GetPhotcodebyCode (catalog[i].measure[m].photcode); - mask = code[0].astromBadMask; - } - if (mask & catalog[i].measure[m].photFlags) { - catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM; - continue; - } - } - - /* exclude measurements by mag limit */ - if (ImagSelect) { - mag = PhotInst (&catalog[i].measure[m]); - if (mag < ImagMin) { - catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM; - continue; - } - if (mag > ImagMax) { - catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM; - continue; - } - } - - /* select or exclude measurements by photcode, or equiv photcode, if specified */ - if (NphotcodesKeep > 0) { - found = FALSE; - for (kp = 0; (kp < NphotcodesKeep) && !found; kp++) { - if (photcodesKeep[kp][0].code == catalog[i].measure[m].photcode) found = TRUE; - if (photcodesKeep[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE; - } - if (!found) continue; - } - if (NphotcodesSkip > 0) { - found = FALSE; - for (kp = 0; (kp < NphotcodesSkip) && !found; kp++) { - if (photcodesSkip[kp][0].code == catalog[i].measure[m].photcode) found = TRUE; - if (photcodesSkip[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE; - } - if (found) continue; - } + catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ; R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); @@ -179,8 +145,8 @@ // too few measurements for average position (require 2 values) - if (N < 2) { + if (N < SRC_MEAS_TOOFEW) { // XXX need to define PHOTOM and ASTROM object flags - // catalog[i].average[j].code |= ID_STAR_FEW; - continue; + catalog[i].average[j].flags |= ID_STAR_FEW; + if (N < 2) continue; } @@ -276,4 +242,11 @@ } + //make sure that the fit succeeded + assert(finite(fit.Ro) && finite(fit.Do) && + finite(fit.dRo) && finite(fit.dDo) && + finite(fit.uR) && finite(fit.uD) && + finite(fit.duR) && finite(fit.duD) && + finite(fit.p) && finite(fit.dp)); + // the measure fields must be updated before the average fields m = catalog[i].average[j].measureOffset; Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 24233) +++ trunk/Ohana/src/relastro/src/args.c (revision 24308) @@ -40,4 +40,13 @@ FIT_TARGET = TARGET_MOSAICS; } + + FlagOutlier = FALSE; + if ((N = get_argument (argc, argv, "-clip"))) { + remove_argument (N, &argc, argv); + CLIP_THRESH = atof (argv[N]); + remove_argument (N, &argc, argv); + FlagOutlier = TRUE; + } + if (FIT_TARGET == TARGET_NONE) usage(); Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 24233) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 24308) @@ -2,5 +2,5 @@ int bcatalog (Catalog *subcatalog, Catalog *catalog) { - + int i, j, k, offset, found; int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; @@ -23,6 +23,8 @@ /* exclude stars not in range or with too few measurements */ for (i = 0; i < catalog[0].Naverage; i++) { - if (catalog[0].average[i].Nmeasure < 2) continue; - + if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) { + continue; + } + /* start with all stars good */ subcatalog[0].average[Naverage] = catalog[0].average[i]; @@ -44,53 +46,17 @@ offset = catalog[0].average[i].measureOffset + j; - - /* select measurements by photcode, or equiv photcode, if specified */ - if (NphotcodesKeep > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesKeep) && !found; k++) { - if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; - } - if (!found) continue; - } - if (NphotcodesSkip > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; - } - if (found) continue; + + //filter objects based on user supplied criteria + if (!MeasFilterTest(&catalog[0].measure[offset])) { + catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; + continue; } - /* select measurements by time */ - if (TimeSelect) { - if (catalog[0].measure[offset].t < TSTART) continue; - if (catalog[0].measure[offset].t > TSTOP) continue; + //filter out outliers + if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) { + catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; + continue; } - - /* select measurements by quality */ - // XXX FIX THIS!! - // if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue; - - /* select measurements by quality */ - if (PhotFlagSelect) { - if (PhotFlagBad) { - mask = PhotFlagBad; - } else { - code = GetPhotcodebyCode (catalog[0].measure[offset].photcode); - mask = code[0].astromBadMask; - } - if (mask & catalog[0].measure[offset].photFlags) continue; - } - - /* select measurements by measurement error */ - if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue; - - /* select measurements by mag limit */ - if (ImagSelect) { - mag = PhotInst (&catalog[0].measure[offset]); - if (mag < ImagMin) continue; - if (mag > ImagMax) continue; - } + catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP; // re-assess on each run of relastro if a measurement should be used @@ -110,13 +76,14 @@ subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; subcatalog[0].measure[Nmeasure].averef = Naverage; - if (RESET) { - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + if (RESET) { + subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; + subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; } + Nmeasure ++; Nm ++; if (Nmeasure == NMEASURE) { - NMEASURE += 1000; - REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); } } @@ -139,6 +106,6 @@ if (VERBOSE) { - fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, - subcatalog[0].Naverage, subcatalog[0].Nmeasure); + fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, + subcatalog[0].Naverage, subcatalog[0].Nmeasure); } return (TRUE); Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 24233) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 24308) @@ -25,6 +25,6 @@ fprintf (stderr, "codename: %s\n", codename); if ((photcodesKeep[NphotcodesKeep] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); + exit (1); } NphotcodesKeep ++; @@ -46,6 +46,6 @@ fprintf (stderr, "codename: %s\n", codename); if ((photcodesSkip[NphotcodesSkip] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); + exit (1); } NphotcodesSkip ++; Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 24233) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 24308) @@ -29,4 +29,9 @@ } if (VERBOSE && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename); + + //outlier rejection + if (FlagOutlier) { + FlagOutliers(&tcatalog); + } // select only the brighter stars Index: trunk/Ohana/src/relastro/src/relastroVisual.c =================================================================== --- trunk/Ohana/src/relastro/src/relastroVisual.c (revision 24308) +++ trunk/Ohana/src/relastro/src/relastroVisual.c (revision 24308) @@ -0,0 +1,512 @@ +/** Diagnostic plots for relastro + * @author Chris Beaumont, IfA + */ + +#include "kapa.h" +#include "relastro.h" + +#define TESTING + +#define KAPAX 700 +#define KAPAY 700 + +static int kapa = -1; +static int kapa2 = -1; +static int kapa3 = -1; + +static int isVisual = FALSE; +static int plotRawRef = FALSE; +static int plotScatter = FALSE; +static int plotResid = FALSE; +static int plotVector = FALSE; +static int plotOutliers = TRUE; + +/** Spawn a kapa window */ +static int initWindow(int *kapid) { + if (*kapid == -1) { + *kapid = KapaOpenNamedSocket("kapa", "relastro_plots"); + if (*kapid == -1) { + fprintf(stderr, "Failure to open kapa.\n"); + isVisual = 0; + return 0; + } + // KapaResize (*kapid, KAPAX, KAPAY); + } + return 1; +} + +/** Ask the user how to proceed */ +static int askUser(int *flag) { + char key[10]; + fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?"); + if (!fgets(key, 8, stdin)) { + fprintf(stderr, "warning: Unable to read option"); + } + if (key[0] == 's') { + *flag = 0; + } + if (key[0] == 'a') { + isVisual = 0; + } + return 1; +} + +/** Size graphdata to encompass all points */ +static int scaleGraphdata(float x[], float y[], Graphdata *graphdata, int n) { + float xlo = FLT_MAX, xhi = -FLT_MAX, ylo = FLT_MAX, yhi = -FLT_MAX; + int i; + int goodData = 0; + for(i = 0; i < n; i++) { + goodData++; + if(!finite(x[i]) || !finite(y[i])) continue; + if(x[i] < xlo) xlo = x[i]; + if(x[i] > xhi) xhi = x[i]; + if(y[i] < ylo) ylo = y[i]; + if(y[i] > yhi) yhi = y[i]; + } + if (goodData < 2) return 0; + graphdata->xmin = xlo; + graphdata->ymin = ylo; + graphdata->xmax = xhi; + graphdata->ymax = yhi; + return 1; +} + +static int residPlot(float x[], float y[], + float xVec[], float yVec[], + int npts, int *kapaID) { + if (!isVisual || !plotResid) return TRUE; + Graphdata graphdata; + KapaSection section; + + if (!initWindow(kapaID)) return FALSE; + + KapaInitGraph(&graphdata); + KapaClearPlots(*kapaID); + KapaSetFont(*kapaID, "helvetica", 14); + + section.name = "0"; + section.x = 0.0; section.y = 0.0; + section.dx = .45, section.dy = .45; + graphdata.ptype = 7; + graphdata.style = 2; + + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(x, xVec, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "x", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dx", KAPA_LABEL_YM); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, x, "x"); + KapaPlotVector(*kapaID, npts, xVec, "y"); + + section.x = .5; section.y = 0; section.name="1"; + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(x, yVec, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "x", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dy", KAPA_LABEL_YP); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, x, "x"); + KapaPlotVector(*kapaID, npts, yVec, "y"); + + section.x = .0; section.y = .5; section.name="2"; + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(y, xVec, &graphdata, npts)) return 0;; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "y", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dx", KAPA_LABEL_YM); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, y, "x"); + KapaPlotVector(*kapaID, npts, xVec, "y"); + + section.x = .5; section.y = .5; section.name="3"; + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(y, yVec, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "y", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dy", KAPA_LABEL_YP); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, y, "x"); + KapaPlotVector(*kapaID, npts, yVec, "y"); + + return 1; +} + + +/**Plot a vector field*/ +static int plotVectorField(float x[], float y[], + float xVec[], float yVec[], + int npts, int *kapaID, double maxVecLength) { + + if(!plotVector) return 1; + + Graphdata graphdata; + float singleX[2], singleY[2]; + float vecScaleFactor; + float graphSize; + int i; + char plotTitle[50]; + sprintf(plotTitle, "Maximum Vector Size = %5.1e", maxVecLength); + + if (!initWindow(kapaID)) return 0; + + KapaInitGraph(&graphdata); + KapaClearPlots(*kapaID); + if(!scaleGraphdata(x, y, &graphdata, npts)) return 0; + + graphSize = graphdata.xmax - graphdata.xmin; + if ((graphdata.ymax - graphdata.ymin) > graphSize) { + graphSize = graphdata.ymax - graphdata.ymin; + } + + vecScaleFactor = graphSize * 0.02 / (float)maxVecLength; +#ifdef TESTING + fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n", + graphSize, maxVecLength, vecScaleFactor); +#endif + + + KapaSetFont (*kapaID, "helvetica", 14); + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, plotTitle, KAPA_LABEL_XP); + + graphdata.ptype = 7; + graphdata.style = 2; + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, x, "x"); + KapaPlotVector(*kapaID, npts, y, "y"); + + //plot each vector individually + graphdata.ptype = 0; + graphdata.style = 0; + graphdata.color = KapaColorByName("blue"); + for(i = 0; i < npts; i++) { + singleX[0] = x[i]; + singleY[0] = y[i]; + singleX[1] = x[i] + xVec[i] * vecScaleFactor; + singleY[1] = y[i] + yVec[i] * vecScaleFactor; + KapaPrepPlot(*kapaID, 2, &graphdata); + KapaPlotVector(*kapaID, 2, singleX, "x"); + KapaPlotVector(*kapaID, 2, singleY, "y"); + } + return 1; +} + +int relastroVisualPlotScatter(double values[], double thresh, int npts) { + float *x, *data; + int i; + float xline[2], yline[2]; + if (!isVisual || !plotScatter) return 1; + if (!initWindow(&kapa2)) return 0; + + ALLOCATE(x, float, npts); + ALLOCATE(data, float, npts); + + for(i = 0; i < npts; i++) { + x[i] = i; + data[i] = (float) values[i]; + } + + Graphdata graphdata; + KapaSection section; + section.x = 0; section.y = 0; section.dx = 1; section.dy = 1; + section.name = "junk"; + + KapaInitGraph(&graphdata); + KapaClearPlots(kapa2); + KapaSetSection(kapa2, §ion); + + graphdata.ptype = 0; + graphdata.style = 0; + graphdata.xmin = 0; + graphdata.xmax = npts; + graphdata.ymin = 0; + graphdata.ymax = data[npts-1]; + + KapaSetFont(kapa2, "helvetica", 14); + KapaSetLimits(kapa2, &graphdata); + KapaBox(kapa2, &graphdata); + KapaSendLabel( kapa2, "Object", KAPA_LABEL_XM); + KapaSendLabel( kapa2, "Offset(pixels)", KAPA_LABEL_YM); + KapaSendLabel( kapa2, "Astrometric Offset with fit cutoff", + KAPA_LABEL_XP); + KapaPrepPlot(kapa2, npts, &graphdata); + KapaPlotVector(kapa2, npts, x, "x"); + KapaPlotVector(kapa2, npts, data, "y"); + + graphdata.color = KapaColorByName("red"); + KapaPrepPlot(kapa2, 2, &graphdata); + yline[0] = (float) thresh; + yline[1] = (float) thresh; + xline[0] = graphdata.xmin; + xline[1] = graphdata.xmax; + KapaPlotVector(kapa2, 2, xline, "x"); + KapaPlotVector(kapa2, 2, yline, "y"); + + askUser(&plotScatter); + return 1; +} + + + +/** plot raw vs ref (L, M). Only those whose distance is < drMax are used in fit*/ +int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj) { + + if( !isVisual || !plotRawRef) return 1; + if( !initWindow(&kapa)) return 0; + + float *rawX, *rawY, *refX, *refY; + float *rawXfit, *rawYfit, *refXfit, *refYfit; + float *magRaw, *magRef, *magRawfit, *magReffit; + float *xVec, *yVec; + int numFit = 0, numNoFit = 0; + double dL, dM, dR; + + ALLOCATE(rawX, float, numObj); + ALLOCATE(rawY, float, numObj); + ALLOCATE(refX, float, numObj); + ALLOCATE(refY, float, numObj); + ALLOCATE(rawXfit, float, numObj); + ALLOCATE(rawYfit, float, numObj); + ALLOCATE(refXfit, float, numObj); + ALLOCATE(refYfit, float, numObj); + ALLOCATE(magRaw, float, numObj); + ALLOCATE(magRef, float, numObj); + ALLOCATE(magRawfit, float, numObj); + ALLOCATE(magReffit, float, numObj); + + int i; + for(i = 0; i < numObj; i++) { + if (raw[i].mask) continue; + + dL = raw[i].L - ref[i].L; + dM = raw[i].M - ref[i].M; + dR = hypot (dL, dM); + if (dR > dRmax) { + rawX[numNoFit] = raw[i].X; + rawY[numNoFit] = raw[i].Y; + refX[numNoFit] = ref[i].X; + refY[numNoFit] = ref[i].Y; + magRaw[numNoFit] = raw[i].Mag; + magRef[numNoFit] = ref[i].Mag; + numNoFit++; + } else { + rawXfit[numFit] = raw[i].X; + rawYfit[numFit] = raw[i].Y; + refXfit[numFit] = ref[i].X; + refYfit[numFit] = ref[i].Y; + magRaw[numFit] = raw[i].Mag; + magRef[numFit] = ref[i].Mag; + numFit++; + } + } + + if (numFit == 0) return 0; + + Graphdata graphdata; + + KapaInitGraph(&graphdata); + KapaClearPlots(kapa); + + graphdata.ptype = 7; + graphdata.style = 2; + + if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) { + fprintf(stderr, "Not enough finite points for plotting"); + return 0; + } + + KapaSetFont(kapa, "helvetica", 14); + KapaSetLimits(kapa, &graphdata); + KapaBox(kapa, &graphdata); + KapaSendLabel( kapa, "X", KAPA_LABEL_XM); + KapaSendLabel( kapa, "Y", KAPA_LABEL_YM); + KapaSendLabel( kapa, "orange, red, green, blue: (raw, ref), (nofit, fit)", + KAPA_LABEL_XP); + + graphdata.color = KapaColorByName("orange"); + graphdata.size = 1; + KapaPrepPlot(kapa, numNoFit, &graphdata); + KapaPlotVector(kapa, numNoFit, rawX, "x"); + KapaPlotVector(kapa, numNoFit, rawY, "y"); + + graphdata.color = KapaColorByName("red"); + graphdata.size = 2; + KapaPrepPlot(kapa, numNoFit, &graphdata); + KapaPlotVector(kapa, numNoFit, refX, "x"); + KapaPlotVector(kapa, numNoFit, refY, "y"); + + graphdata.color = KapaColorByName("green"); + graphdata.size = 1; + KapaPrepPlot(kapa, numFit, &graphdata); + KapaPlotVector(kapa, numFit, rawXfit, "x"); + KapaPlotVector(kapa, numFit, rawYfit, "y"); + + graphdata.color = KapaColorByName("blue"); + graphdata.size = 2; + KapaPrepPlot(kapa, numFit, &graphdata); + KapaPlotVector(kapa, numFit, refXfit, "x"); + KapaPlotVector(kapa, numFit, refYfit, "y"); + + ALLOCATE(xVec, float, numFit); + ALLOCATE(yVec, float, numFit); + + //plot the fitted objects as vectors + for(i = 0; i < numFit; i++) { + xVec[i] = rawXfit[i] - refXfit[i]; + yVec[i] = rawYfit[i] - refYfit[i]; + } + + plotVectorField(rawXfit, rawYfit, xVec, yVec, numFit, &kapa3, dRmax); + if(!residPlot(rawXfit, rawYfit, xVec, yVec, numFit, &kapa2)) { + fprintf(stderr, "Unable to plot residuals"); + return 0; + } + + FREE(xVec); + FREE(yVec); + + askUser(&plotRawRef); + + FREE(rawX); + FREE(rawY); + FREE(refX); + FREE(refY); + FREE(rawXfit); + FREE(rawYfit); + FREE(refXfit); + FREE(refYfit); + FREE(magRaw); + FREE(magRef); + FREE(magRawfit); + FREE(magReffit); + + return 1; +} + +int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, + StatType statsR, StatType statsD, double thresh) { + + float *Din, *Rin, *Dout, *Rout; + double xmin, xmax, ymin, ymax, range; + float xCirc[100], yCirc[100]; + int m, i; + int Nin, Nout; + Measure meas; + Graphdata graphdata; + KapaSection section; + + if (!isVisual || !plotOutliers) return 1; + if (!initWindow(&kapa)) return 0; + + // populate vectors + ALLOCATE(Din, float, Nmeasure); + ALLOCATE(Rin, float, Nmeasure); + ALLOCATE(Dout, float, Nmeasure); + ALLOCATE(Rout, float, Nmeasure); + + //create the threshhold ellipse + for(i = 0; i < 100; i++) { + xCirc[i] = statsR.median + thresh * statsR.sigma * cos(2 * 3.14 / 99. * i); + yCirc[i] = statsD.median + thresh * statsD.sigma * sin(2 * 3.14 / 99. * i); + } + + m = offset; + Nin = Nout = 0; + xmin = +FLT_MAX; + xmax = -FLT_MAX; + ymax = -FLT_MAX; + ymin = +FLT_MAX; + for(i = 0; i < Nmeasure; i++, m++) { + meas = catalog[0].measure[m]; + if (!MeasFilterTest(&meas)) continue; + xmin = MIN(xmin, meas.dR); + xmax = MAX(xmax, meas.dR); + ymin = MIN(ymin, meas.dD); + ymax = MAX(ymax, meas.dD); + + if (meas.dbFlags & ID_MEAS_POOR_ASTROM) { + Rout[Nout] = (float)(meas.dR); + Dout[Nout] = (float)(meas.dD); + fprintf(stderr, "r: %f\td: %f\t outlier: 1\n", Rout[Nout], Dout[Nout]); + Nout++; + } else { + Rin[Nin] = (float)(meas.dR); + Din[Nin] = (float)(meas.dD); + fprintf(stderr, "r: %f\td: %f\t outlier: 0\n", Rin[Nin], Din[Nin]); + Nin++; + } + } + + range = (xmax - xmin); + xmin -= .1 * range; + xmax += .1 * range; + range = (ymax - ymin); + ymax += .1 * range; + ymin -= .1 * range; + + //temporary fix + xmin = -1; xmax = 1; ymin = -1; ymax = 1; + + + //initialize graph info + section.x = 0; section.y = 0; section.dx = 1; section.dy = 1; + section.name = "junk"; + + KapaInitGraph(&graphdata); + KapaClearPlots(kapa); + KapaSetFont(kapa, "helvetica", 14); + + graphdata.ptype = 7; + graphdata.style = 2; + graphdata.size = 3; + graphdata.xmin = xmin; + graphdata.xmax = xmax; + graphdata.ymin = ymin; + graphdata.ymax = ymax; + + KapaSetSection(kapa, §ion); + KapaSetLimits(kapa, &graphdata); + KapaBox(kapa, &graphdata); + + KapaSendLabel( kapa, "RA (arcsec)", KAPA_LABEL_XM); + KapaSendLabel( kapa, "Dec (arcsec)", KAPA_LABEL_YM); + KapaSendLabel( kapa, "Points flagged as outliers (red)", + KAPA_LABEL_XP); + + graphdata.color = KapaColorByName("green"); + KapaPrepPlot(kapa, Nin, &graphdata); + KapaPlotVector(kapa, Nin, Rin, "x"); + KapaPlotVector(kapa, Nin, Din, "y"); + + graphdata.color = KapaColorByName("red"); + KapaPrepPlot(kapa, Nout, &graphdata); + KapaPlotVector(kapa, Nout, Rout, "x"); + KapaPlotVector(kapa, Nout, Dout, "y"); + + graphdata.color = KapaColorByName("black"); + graphdata.ptype = 0; + graphdata.style = 0; + KapaPrepPlot(kapa, 100, &graphdata); + KapaPlotVector(kapa, 100, xCirc, "x"); + KapaPlotVector(kapa, 100, yCirc, "y"); + + + + askUser(&plotOutliers); + + FREE(Rout); + FREE(Dout); + FREE(Rin); + FREE(Din); + + return TRUE; +} + + Index: trunk/Ohana/src/relastro/src/relastroVisual.h =================================================================== --- trunk/Ohana/src/relastro/src/relastroVisual.h (revision 24308) +++ trunk/Ohana/src/relastro/src/relastroVisual.h (revision 24308) @@ -0,0 +1,9 @@ +#ifndef RELASTRO_VISUAL_H +#define RELASTROPVISUAL_H + +int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj); +int relastroVisualPlotScatter(double values[], double thresh, int npts); +int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, + StatType statsR, StatType statsD, double thresh); + +#endif Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 24233) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 24308) @@ -8,4 +8,5 @@ SkyList *skylist = NULL; Catalog catalog; + // load the current sky table (layout of all SkyRegions) @@ -37,4 +38,8 @@ } + if (FlagOutlier) { + FlagOutliers(&catalog); + } + // XXX consider what gets reset (only ASTROM flags) if (RESET) {