Index: /trunk/Ohana/src/photdbc/include/photdbc.h =================================================================== --- /trunk/Ohana/src/photdbc/include/photdbc.h (revision 30615) +++ /trunk/Ohana/src/photdbc/include/photdbc.h (revision 30616) @@ -51,8 +51,13 @@ double DMGAIN; double CHISQ_MAX; +double SIGMA_MIN_KEEP; double SIGMA_MAX; double AVE_SIGMA_LIM; int NMEAS_MIN; +int NMEAS_MIN_FILTERED; +int NCODE_MIN; double ZERO_POINT; + +int ExcludeByMinSigma; int ExcludeByInstMag; @@ -65,4 +70,8 @@ SkyRegion REGION; PhotCodeData photcodes; + +char *PHOTCODE_DROP_LIST, *PHOTCODE_SKIP_LIST; +int NphotcodesDrop, NphotcodesSkip; +PhotCode **photcodesDrop, **photcodesSkip; # define FLAG_AREA 0X0001 Index: /trunk/Ohana/src/photdbc/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/photdbc/src/ConfigInit.c (revision 30615) +++ /trunk/Ohana/src/photdbc/src/ConfigInit.c (revision 30616) @@ -26,7 +26,11 @@ if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); - WarnConfig (config, "PHOTDBC_JOIN_RADIUS", "%lf", 0, &JOIN_RADIUS); + // XXX join_stars (in photdbc.c) is currently disabled + // WarnConfig (config, "PHOTDBC_JOIN_RADIUS", "%lf", 0, &JOIN_RADIUS); + + // XXX unique_measures (in photdbc.c) is currently disabled // WarnConfig (config, "UNIQ_RADIUS", "%lf", 0, &UNIQ_RADIUS); + // XXX flag_measures (in photdbc.c) is currently disabled // WarnConfig (config, "XMIN", "%lf", 0, &XMIN); // WarnConfig (config, "XMAX", "%lf", 0, &XMAX); @@ -37,11 +41,14 @@ // WarnConfig (config, "MMIN", "%lf", 0, &tmp); MMIN = 1000*tmp; // WarnConfig (config, "MMAX", "%lf", 0, &tmp); MMAX = 1000*tmp; + // WarnConfig (config, "CHISQ_MAX", "%lf", 0, &CHISQ_MAX); + + // XXX get_mags (in photdbc.c) is currently disabled // WarnConfig (config, "DMSYS", "%lf", 0, &tmp); DMSYS = SQ(1000*tmp); // WarnConfig (config, "DMGAIN", "%lf", 0, &DMGAIN); - // WarnConfig (config, "CHISQ_MAX", "%lf", 0, &CHISQ_MAX); ScanConfig (config, "SIGMA_MAX", "%lf", 0, &SIGMA_MAX); ScanConfig (config, "AVE_SIGMA_LIM", "%lf", 0, &AVE_SIGMA_LIM); ScanConfig (config, "NMEAS_MIN", "%d", 0, &NMEAS_MIN); + ScanConfig (config, "NMEAS_MIN_FILTERED", "%d", 0, &NMEAS_MIN_FILTERED); WarnConfig (config, "GSCFILE", "%s", 0, GSCFILE); Index: /trunk/Ohana/src/photdbc/src/args.c =================================================================== --- /trunk/Ohana/src/photdbc/src/args.c (revision 30615) +++ /trunk/Ohana/src/photdbc/src/args.c (revision 30616) @@ -19,5 +19,4 @@ ExcludeByInstMag = FALSE; if ((N = get_argument (argc, argv, "-instmag"))) { - remove_argument (N, &argc, argv); ExcludeByInstMag = TRUE; remove_argument (N, &argc, argv); @@ -25,12 +24,21 @@ remove_argument (N, &argc, argv); INST_MAG_MAX = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + ExcludeByMinSigma = FALSE; + if ((N = get_argument (argc, argv, "-min-sigma"))) { + ExcludeByMinSigma = TRUE; + remove_argument (N, &argc, argv); + SIGMA_MIN_KEEP = atof(argv[N]); + remove_argument (N, &argc, argv); } ExcludeByMaxMinMag = FALSE; if ((N = get_argument (argc, argv, "-maxminmag"))) { - remove_argument (N, &argc, argv); ExcludeByMaxMinMag = TRUE; remove_argument (N, &argc, argv); MAX_MIN_MAG = atof(argv[N]); + remove_argument (N, &argc, argv); } @@ -61,4 +69,18 @@ } + PHOTCODE_DROP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode-drop"))) { + remove_argument (N, &argc, argv); + PHOTCODE_DROP_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_SKIP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode-skip"))) { + remove_argument (N, &argc, argv); + PHOTCODE_SKIP_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + if (argc != 2) usage(); @@ -81,20 +103,35 @@ fprintf (stderr, "USAGE: photdbc (output)\n\n"); fprintf (stderr, " this program takes an existing DVO database and makes a copy, applying a number of optional\n"); - fprintf (stderr, " filters in the process. \n"); - fprintf (stderr, "\n"); - fprintf (stderr, " photdbc (output)\n"); - fprintf (stderr, "\n"); - fprintf (stderr, " allowed filters / restrictions include:\n"); - fprintf (stderr, "\n"); + fprintf (stderr, " filters in the process. \n\n"); + + fprintf (stderr, " photdbc (output)\n\n"); + + fprintf (stderr, " allowed filters / restrictions include:\n\n"); + fprintf (stderr, " -region Rmin Rmax Dmin Dmax : limit operation to the specified region \n"); - fprintf (stderr, " -join : join measurements between stars using JOIN_RADIUS\n"); - fprintf (stderr, " -ccdregion X Y X Y : only keep detections within the specified detector region\n"); - fprintf (stderr, " (can this be limited to specific photcodes? cameras? detectors?)\n"); - fprintf (stderr, " -photcode_limits code Mmin Mmax : allow multiples of these\n"); - fprintf (stderr, "\n"); - fprintf (stderr, " SIGMA_MAX\n"); - fprintf (stderr, " NMEAS_MIN\n"); - fprintf (stderr, " INST_MAG_MAX\n"); - fprintf (stderr, " INST_MAG_MIN\n"); + + fprintf (stderr, " -photcode-drop : remove these photcodes from the output (REF or DEP only)\n"); + fprintf (stderr, " -photcode-skip : ignore these photcodes when assessing the validity (keep unless object is dropped)\n"); + + fprintf (stderr, " -instmag (min) (max) : range of valid instrumental magnitudes (or measurements are dropped)\n"); + fprintf (stderr, " -min-sigma (sigma) : object must have one measurement error less than sigma or object is dropped\n"); + + fprintf (stderr, " -maxminmag (mag) : object must have one magnitude less than this or object is dropped\n"); + + fprintf (stderr, " option options:\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -params : list the current parameters\n\n"); + + fprintf (stderr, "ptolemy.rc config values used by this program:\n"); + fprintf (stderr, " SIGMA_MAX : drop measurements with errors greater than this\n"); + fprintf (stderr, " NMEAS_MIN : drop objects with fewer measurements than this\n"); + fprintf (stderr, " NMEAS_MIN_FILTERED : drop objects with fewer measurements than this after filtering above\n"); + fprintf (stderr, " AVE_SIGMA_LIM : drop objects if all average mags are greater than this\n"); + exit (2); } + +// fprintf (stderr, " -join : join measurements between stars using JOIN_RADIUS\n"); +// fprintf (stderr, " -ccdregion X Y X Y : only keep detections within the specified detector region\n"); +// fprintf (stderr, " (can this be limited to specific photcodes? cameras? detectors?)\n"); +// fprintf (stderr, " -photcode_limits code Mmin Mmax : allow multiples of these\n"); Index: /trunk/Ohana/src/photdbc/src/find_images.c =================================================================== --- /trunk/Ohana/src/photdbc/src/find_images.c (revision 30615) +++ /trunk/Ohana/src/photdbc/src/find_images.c (revision 30616) @@ -122,5 +122,6 @@ these are measurements from external sources, like USNO */ - assignMcal (&image[nimage], (double *) NULL, -1); + // assignMcal (&image[nimage], (double *) NULL, -1); + image[nimage].Mcal = 0; image[nimage].code = ID_IMAGE_NEW; Index: /trunk/Ohana/src/photdbc/src/initialize.c =================================================================== --- /trunk/Ohana/src/photdbc/src/initialize.c (revision 30615) +++ /trunk/Ohana/src/photdbc/src/initialize.c (revision 30616) @@ -2,4 +2,7 @@ void initialize (int argc, char **argv) { + + int NPHOTCODES; + char *codename, *ptr, *list; /* are these set correctly? */ @@ -11,4 +14,46 @@ ConfigInit (&argc, argv); args (argc, argv); + + NphotcodesDrop = 0; + photcodesDrop = NULL; + if (PHOTCODE_DROP_LIST != NULL) { + NPHOTCODES = 10; + ALLOCATE (photcodesDrop, PhotCode *, NPHOTCODES); + + /* parse the comma-separated list of photcodesDrop */ + list = PHOTCODE_DROP_LIST; + while ((codename = strtok_r (list, ",", &ptr)) != NULL) { + list = NULL; // pass NULL on successive strtok_r calls + fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_DROP_LIST); + fprintf (stderr, "codename: %s\n", codename); + if ((photcodesDrop[NphotcodesDrop] = GetPhotcodebyName (codename)) == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); + exit (1); + } + NphotcodesDrop ++; + CHECK_REALLOCATE (photcodesDrop, PhotCode *, NPHOTCODES, NphotcodesDrop, 10); + } + } + + NphotcodesSkip = 0; + photcodesSkip = NULL; + if (PHOTCODE_SKIP_LIST != NULL) { + NPHOTCODES = 10; + ALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES); + + /* parse the comma-separated list of photcodesSkip */ + list = PHOTCODE_SKIP_LIST; + while ((codename = strtok_r (list, ",", &ptr)) != NULL) { + list = NULL; // pass NULL on successive strtok_r calls + fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_SKIP_LIST); + 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); + } + NphotcodesSkip ++; + CHECK_REALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES, NphotcodesSkip, 10); + } + } if (SHOW_PARAMS) { Index: /trunk/Ohana/src/photdbc/src/make_subcatalog.c =================================================================== --- /trunk/Ohana/src/photdbc/src/make_subcatalog.c (revision 30615) +++ /trunk/Ohana/src/photdbc/src/make_subcatalog.c (revision 30616) @@ -5,11 +5,24 @@ int make_subcatalog (Catalog *subcatalog, Catalog *catalog) { - off_t i, j, offset; + int found; + off_t i, j, k, offset; off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; - double mag, minMag; - int keep; + double mag, minMag, minSigma; + int keep, *secSkip; + PhotCode *photcode; Nsecfilt = GetPhotcodeNsecfilt (); assert (catalog[0].Nsecfilt == Nsecfilt); + + // set up a list of SEC entries to ignore when evaluating a source + ALLOCATE (secSkip, int, Nsecfilt); + for (i = 0; i < Nsecfilt; i++) { + secSkip[i] = FALSE; + photcode = GetPhotcodebyNsec(i); + for (k = 0; k < NphotcodesSkip; k++) { + if (photcodesSkip[k][0].code != photcode[0].code) continue; + secSkip[i] = TRUE; + } + } /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ @@ -38,4 +51,17 @@ } + // exclude stars with too few measurements + if (NCODE_MIN) { + // drop if all of the allowed average photcodes have ncode < NCODE_MIN values + keep = FALSE; + for (j = 0; !keep && (j < Nsecfilt); j++) { + if (secSkip[j]) continue; + if (catalog[0].secfilt[Nsecfilt*i+j].Ncode >= NCODE_MIN) { + keep = TRUE; + } + } + if (!keep) continue; + } + /* assign average and secfilt values */ subcatalog[0].average[Naverage] = catalog[0].average[i]; @@ -45,5 +71,6 @@ } - minMag = 32; + minMag = 32; + minSigma = 32; Nm = 0; for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { @@ -52,15 +79,31 @@ # if 0 - if (DropPhotcode) { - ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode); - if (ecode == photcode[0].code) continue; - } if (DropNonStellar && catalog[0].measure[offset].dophot != 1) continue; # endif - /* exclude measurements by measurement error */ + // remove certain photcodes from the output measurements + if (NphotcodesDrop > 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; + } + + // ignore certain photcodes to assess the measurements + 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) goto keep; + } + + // exclude measurements by measurement error -- drop exactly this measurement if (SIGMA_MAX && (catalog[0].measure[offset].dM > SIGMA_MAX)) continue; - /* select measurements by mag limit */ + // select measurements by mag limit -- drop exactly this measurement if (ExcludeByInstMag) { mag = PhotInst (&catalog[0].measure[offset]); @@ -69,8 +112,16 @@ } + // check measurements for this object -- drop object if no measurements pass + if (ExcludeByMinSigma) { + minSigma = MIN (minSigma, catalog[0].measure[offset].dM); + } + + // check measurements for this object -- drop object if no measurements pass if (ExcludeByMaxMinMag) { mag = PhotSys (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); minMag = MIN (minMag, mag); } + + keep: subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; @@ -90,9 +141,16 @@ } - // after measurement exclusion, exclude stars with too few measurements - if (Nm < NMEAS_MIN) { + // exclude faint objects + if (ExcludeByMinSigma && (minSigma > SIGMA_MIN_KEEP)) { Nmeasure -= Nm; continue; } + + // after measurement exclusion, exclude stars with too few measurements + if (NMEAS_MIN_FILTERED && (Nm < NMEAS_MIN_FILTERED)) { + Nmeasure -= Nm; + continue; + } + subcatalog[0].average[Naverage].Nmissing = 0; subcatalog[0].average[Naverage].Nmeasure = Nm; @@ -122,2 +180,24 @@ return (TRUE); } + +/** + the purpose of this function is to create a subset database. there are several ways this could be done: + + * reject all measurements from all objects that fail to meet some criteria (objects may lose measurements and/or be dropped) + + * only reject an object if all measurements fail to meet some criteria (ie, keep all measurements if any measurement passes) + + * only apply the keep / reject criteria to certain photcodes + + + *** + + * options: + + minimum value for nmeas : NMEAS_MIN + minimum value for nmeas : NMEAS_MIN_FILTERED + minimum value for ncode : NCODE_MIN (drop if all ncode < NCODE_MIN) (respect photcodes) + reject faint meas : SIGMA_MAX + reject faint source : SIGMA_MIN_KEEP (keep source source if its minimum dMag < this limit) + + **/ Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 30615) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 30616) @@ -26,4 +26,5 @@ double dPos; int mask; + int Nmeas; } StarData; @@ -94,6 +95,10 @@ double SIGMA_LIM; -int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored +int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored double MIN_ERROR; + +int IMFIT_CLIP_NITER; // number of clipping iterations to perform in FitChip +double IMFIT_CLIP_NSIGMA; // number of sigma to clip in FitChip +double IMFIT_SYS_SIGMA_LIM; // max dMag for objects used to measure systematic scatter double RADIUS; // match radius for high-speed objects @@ -141,4 +146,5 @@ int FlagOutlier; int CLIP_THRESH; +int USE_BASIC_CHECK; FitMode FIT_MODE; @@ -296,5 +302,5 @@ void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im); void FlagOutliers(Catalog *catalog); -int MeasFilterTest(Measure *measure); +int MeasFilterTest(Measure *measure, int applySigmaLim); int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); @@ -325,10 +331,10 @@ int UpdateObjectOffsets (SkyList *skylist); -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); - - +int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj); +// 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); +void relastroSetVisual(int state); +int relastroGetVisual(void); int FixProblemImages (SkyList *skylist); @@ -352,2 +358,4 @@ int createStarMapPoints(); int checkStarMap(int N); + +int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit); Index: /trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 30616) @@ -19,6 +19,10 @@ if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); - GetConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + GetConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM); // exclude measurements on this basis GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW); + + if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NITER", "%d", 0, &IMFIT_CLIP_NITER)) IMFIT_CLIP_NITER = 3; + if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NSIGMA", "%lf", 0, &IMFIT_CLIP_NSIGMA)) IMFIT_CLIP_NSIGMA = 3.0; + if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01; // XXX these are used in relphot to identify poor stars / images -- define Index: /trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitChip.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/FitChip.c (revision 30616) @@ -2,44 +2,51 @@ # include "relastroVisual.h" -// XXX make these user parameters -# define FIT_CHIP_NITER 3 -# define FIT_CHIP_NSIGMA 3.0 - -// XXX save the fit[0].Npts value in the image table? - -// XXX save measurements of the fit quality (scatter, chisq) in the image table - int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) { - int i, Nscatter, Niter, skip; - CoordFit *fit; - double dL, dM, dR, dRmax, *values; - - ALLOCATE (values, double, Nmatch); - for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) { + int i, NstatFull, Nstat, Niter, skip; + float dLsig, dMsig, dRsig; + float dLsigFull, dMsigFull, dRsigFull; + float dL, dM, dR, dRmax; + + CoordFit *fit = NULL; + + dRmax = 0.0; + + for (Niter = 0; Niter < IMFIT_CLIP_NITER; Niter ++) { + + // measure the scatter for the unmarked (good) measurements with dM < limit + // SIGMA_LIM here is redundant with ImageOps.c:915 + GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, SIGMA_LIM); + dRmax = IMFIT_CLIP_NSIGMA*dRsig; + + // (a) skip unfittable points + // (b) mark good and bad points for fit (in and outliers) // measure the scatter distribution (use only the bright end detections) - for (i = Nscatter = 0; i < Nmatch; i++) { + for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) continue; + if (isnan(raw[i].dMag)) { + raw[i].mask |= 0x0004; + continue; + } + if (raw[i].dMag > SIGMA_LIM) { + raw[i].mask |= 0x0008; + continue; + } // is this redundant with ImageOps.c:915? dL = raw[i].L - ref[i].L; dM = raw[i].M - ref[i].M; dR = hypot (dL, dM); - - values[Nscatter] = dR; - Nscatter++; - } - - if (Nscatter > 5) { - // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma - // XXX this test is not sensible for Nscatter < XXX (5?) - dsort (values, Nscatter); - dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)]; - relastroVisualPlotScatter(values, dRmax, Nscatter); - relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch); - } else { - dRmax = 0; - } + if (dR > dRmax) { + raw[i].mask |= 0x0010; + continue; + } + } + + // figures to assess the fitting process: + // x vs dx, x vs dy, y vs dx, y vs dy : + // residual vector field + // dx vs mag, dy vs mag + relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch); // fit the requested order polynomial @@ -47,24 +54,10 @@ image[0].coords.Npolyterms = CHIPORDER; } + if (fit) fit_free (fit); fit = fit_init (image[0].coords.Npolyterms); // generate the fit matches for (i = 0; i < Nmatch; i++) { - if (raw[i].mask) { - continue; - } - if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) { - continue; - } - - // only keep objects within dRmax - dL = raw[i].L - ref[i].L; - dM = raw[i].M - ref[i].M; - dR = hypot (dL, dM); - - // fprintf (stderr, "fit %f %f -> %f %f : %f %f (%f vs %f) wt: %f\n", raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].L, raw[i].M, dR, dRmax, raw[i].dPos); - - if ((dRmax > 0.0) && (dR > dRmax)) continue; - + if (raw[i].mask) continue; fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos); } @@ -88,12 +81,9 @@ } if (skip) { - if (VERBOSE) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms); + if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms); fit_free (fit); - free (values); image[0].flags |= ID_IMAGE_ASTROM_FEW; return FALSE; } - - // fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch); // measure the fit, update the coords & object coordinates @@ -110,6 +100,4 @@ } - fit_free (fit); - for (i = 0; i < Nmatch; i++) { XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords); @@ -118,6 +106,85 @@ } - free (values); + // count mask classes + int nMask1 = 0; + int nMask2 = 0; + int nMask3 = 0; + int nMask4 = 0; + int nMask5 = 0; + for (i = 0; i < Nmatch; i++) { + if (!raw[i].mask) continue; + if (raw[i].mask & 0x01) nMask1 ++; + if (raw[i].mask & 0x02) nMask2 ++; + if (raw[i].mask & 0x04) nMask3 ++; + if (raw[i].mask & 0x08) nMask4 ++; + if (raw[i].mask & 0x10) nMask5 ++; + } + + GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM); + GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM); + + int Nm = 0; + int Ns = 0; + for (i = 0; i < Nmatch; i++) { + if (raw[i].mask) continue; + Nm += raw[i].Nmeas; + Ns++; + } + image[0].nLinkAstrom = (Nm / Ns); + + if (VERBOSE) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5); + + image[0].dXpixSys = dLsig; + image[0].dYpixSys = dMsig; + image[0].nFitAstrom = fit[0].Npts; + + if (fit) fit_free (fit); + return TRUE; +} + +// measure the scatter distribution (limit selected objects by dMag) +int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit) { + + int i, Ns; + float dL, dM, dLsum, dLsum2, dMsum, dMsum2, *dR; + + ALLOCATE (dR, float, Nstars); + + Ns = 0; + dLsum = dLsum2 = dMsum = dMsum2 = 0.0; + for (i = 0; i < Nstars; i++) { + if (raw[i].mask) continue; + if (isnan(raw[i].dMag)) continue; + if ((SigmaLimit > 0.0) && (raw[i].dMag > SigmaLimit)) continue; + + dL = raw[i].L - ref[i].L; + dM = raw[i].M - ref[i].M; + + dLsum += dL; + dLsum2 += dL*dL; + dMsum += dM; + dMsum2 += dM*dM; + + dR[Ns] = hypot (dL, dM); + Ns++; + } + + *dLsig = sqrt((dLsum2 / Ns) - SQ(dLsum/Ns)); + *dMsig = sqrt((dMsum2 / Ns) - SQ(dMsum/Ns)); + *nKeep = Ns; + + if (Ns < 5) { + *dRsig = NAN; + free (dR); + return (FALSE); + } + + // for a 2D Gaussian, 40% of the points are within R = 1 sigma + fsort (dR, Ns); + *dRsig = dR[(int)(0.40*Ns)]; + + free (dR); + return (TRUE); } @@ -162,2 +229,3 @@ fclose (f); */ + Index: /trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- /trunk/Ohana/src/relastro/src/GetAstromError.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/GetAstromError.c (revision 30616) @@ -1,24 +1,31 @@ # include "relastro.h" +# define WEIGHTED_ERRORS 1 float GetAstromError (Measure *measure, int mode) { - //BIG HACKXXXXXXXX - return 0.1; + PhotCode *code; - float dPobs, dPsys, dPtotal, dM, AS, MS; + float dPobs, dPsys, dPtotal, dM, AS, MS, dX, dY; + + if (!WEIGHTED_ERRORS) { + // if we don't understand the errors at all, this at least lets us get things roughly + // right: + return 0.1; + } + switch (mode) { case ERROR_MODE_RA: - dPobs = measure[0].dXccd / 100.0; // need to redefine this as RAerr + dPobs = FromShortPixels(measure[0].dXccd); // dXccd is a value in pixels break; case ERROR_MODE_DEC: - dPobs = measure[0].dYccd / 100.0; // need to redefine this as RAerr + dPobs = FromShortPixels(measure[0].dYccd); // dYccd is a value in pixels break; case ERROR_MODE_POS: - dPobs = hypot (measure[0].dXccd, measure[0].dYccd) / 100.0; // need to redefine this as RAerr + dX = FromShortPixels(measure[0].dXccd); // dXccd is a value in pixels + dY = FromShortPixels(measure[0].dYccd); // dYccd is a value in pixels + dPobs = hypot (dX, dY); break; default: abort(); } - /* the astrometric errors are not being carried yet (but should be!) */ - /* we use the photometric mag error as a weighting term */ code = GetPhotcodebyCode (measure[0].photcode); @@ -27,8 +34,5 @@ dPsys = code[0].astromErrSys; dM = measure[0].dM; - dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM)); - - //XXX dXccd, dYccd are now working correctly - //dPtotal = AS * dPobs; + dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM)); dPtotal = MAX (dPtotal, MIN_ERROR); @@ -36,2 +40,17 @@ } +/* for a long time, psphot was either not reported position errors, or was reporting + * completely wrong astrometry errors. This function lets us handle, in the + * configuration, different strategies to generating a position error + * + * astrometry systematic error : this is the minimum expected per-position error. You + * should probably measure this from you data. watch out for the chicken and egg problem! + * + * astrometry error scale : this field lets you accept position errors in (say) pixels and + * convert them with this term to arcsec. AS : pixel scale in arcsec + * + * astrometry mag scale : this field lets you define position errors based on the + * photometry error. the scale factor should be something like a typical seeing number + * (in arcsec) for the given instrument + * + */ Index: /trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relastro/src/ImageOps.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/ImageOps.c (revision 30616) @@ -146,4 +146,6 @@ free (clist); free (mlist); + free (Nlist); + free (NLIST); } @@ -325,5 +327,5 @@ for (j = 0; j < average[0].Nmeasure; j++) { off = average[0].measureOffset + j; - fprintf (stderr, "%f, %f\n", measure[off].dR, measure[off].dD); + fprintf (stderr, "dR, dD, mag, dMag: %f, %f, %f, %f\n", measure[off].dR, measure[off].dD, measure[off].M, measure[off].dM); } return; @@ -336,5 +338,5 @@ off_t i, m, c, n, nPos; double X, Y, L, M, P, Q, R, D, dR, dD; - double dPos, DPOS_MAX_ASEC; + double dPos, dPosSys, DPOS_MAX_ASEC; Mosaic *mosaic; @@ -359,9 +361,32 @@ } - // accumulate the rms position offsets. if this value, or any specific entry, is too - // large, we will reset the image to the original coords at the end of the analysis - + // these are used to accumulate the rms position offsets. if this value, or any + // specific entry, is too large, we will reset the image to the original coords at the + // end of the analysis dPos = 0.0; nPos = 0; + + // convert the image systematic error in pixels to a value in arcsec + { + double dLsig, dMsig; + double Ro, Do, Rx, Dx, dP0, dP1; + Coords *coords; + + // these values are in pixels, but we to convert to arcsec + dLsig = image[0].dXpixSys; + dMsig = image[0].dYpixSys; + + if (moscoords == NULL) { + coords = imcoords; + } else { + coords = moscoords; + } + XY_to_LM (&Ro, &Do, 0.0, 0.0, coords); + XY_to_LM (&Rx, &Dx, dLsig, 0.0, coords); + dP0 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec + XY_to_LM (&Rx, &Dx, 0.0, dLsig, coords); + dP1 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec + dPosSys = 0.5 * (dP0 + dP1); + } for (i = 0; i < Nlist[im]; i++) { @@ -390,10 +415,10 @@ // complain if the new location is far from the average location // NOTE: This should never happen, or our StarMap tests are not working - if (fabs(dR) > 1.5*ADDSTAR_RADIUS) { + if (fabs(dR) > 3.0*ADDSTAR_RADIUS) { fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); dump_measures (&catalog[c].average[n], catalog[c].measure); // abort (); } - if (fabs(dD) > 1.5*ADDSTAR_RADIUS) { + if (fabs(dD) > 3.0*ADDSTAR_RADIUS) { fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); dump_measures (&catalog[c].average[n], catalog[c].measure); @@ -418,5 +443,5 @@ catalog[c].measure[m].dR = dR; catalog[c].measure[m].dD = dD; - + if (catalog[c].measure[m].dR > +180.0*3600.0) { // average on high end of boundary, move star up @@ -429,4 +454,7 @@ catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); } + + // set the systematic error for this image: + catalog[c].measure[m].dRsys = ToShortPixels(dPosSys); } @@ -544,16 +572,13 @@ // 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; - int mask = FALSE; + // XXX this is already taken care of in bcatalog + raw[i].mask = 0x0000; if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) { - mask = TRUE; + raw[i].mask |= 0x0001; } if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) { - mask = TRUE; - } - - raw[i].mask = mask; - + raw[i].mask |= 0x0002; + } + raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is switch (mode) { @@ -685,5 +710,5 @@ // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m])) continue; + if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; R[N] = catalog[0].measure[m].dR; D[N] = catalog[0].measure[m].dD; @@ -711,5 +736,5 @@ // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m])) continue; + if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; x = catalog[0].measure[m].dR - statsR.median; @@ -725,8 +750,6 @@ } - //examine results - relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, - catalog[0].average[j].Nmeasure, - statsR, statsD, Ns); + // examine results + // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns); } @@ -784,5 +807,5 @@ // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m])) continue; + if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; R[N] = catalog[0].measure[m].dR; D[N] = catalog[0].measure[m].dD; @@ -806,5 +829,5 @@ for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { //skip bad measurements - if (!MeasFilterTest(&catalog[0].measure[m])) continue; + if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; x = catalog[0].measure[m].dR - statsR.median; y = catalog[0].measure[m].dD - statsD.median; @@ -838,5 +861,5 @@ for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { //skip bad measurements - if (!MeasFilterTest(&catalog[0].measure[m])) continue; + if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; x = catalog[0].measure[m].dR - statsR.median; y = catalog[0].measure[m].dD - statsD.median; @@ -852,8 +875,6 @@ } //done rejecting outliers - //examine results - relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, - catalog[0].average[j].Nmeasure, - statsR, statsD, Ns); + // examine results + // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns); } //done looping over objects @@ -871,5 +892,6 @@ /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ -int MeasFilterTest(Measure *measure) { +// we only optionally apply the sigma limit: for object averages, this should not be used (should it?) +int MeasFilterTest(Measure *measure, int applySigmaLim) { int found, k; long mask; @@ -917,5 +939,7 @@ /* select measurements by measurement error */ - if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE; + if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) { + return FALSE; + } /* select measurements by mag limit */ Index: /trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- /trunk/Ohana/src/relastro/src/StarMaps.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/StarMaps.c (revision 30616) @@ -79,5 +79,5 @@ starmap[N].stars[ybin*NX_MAP + xbin] ++; } - MARKTIME("assign stars to starmap bins: %f sec\n", dtime); + if (VERBOSE2) { MARKTIME("assign stars to starmap bins: %f sec\n", dtime); } return (TRUE); @@ -117,5 +117,5 @@ } - if (VERBOSE) fprintf (stderr, "starmap: %d points for image %s\n", starmap[i].Npoints, images[i].name); + if (VERBOSE2) fprintf (stderr, "starmap: %d points for image %s\n", starmap[i].Npoints, images[i].name); } @@ -145,8 +145,14 @@ } - if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + if (VERBOSE2) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); - if (dLmax > DPOS_MAX) return (FALSE); - if (dMmax > DPOS_MAX) return (FALSE); + if (dLmax > DPOS_MAX) { + if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + return (FALSE); + } + if (dMmax > DPOS_MAX) { + if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + return (FALSE); + } return (TRUE); } Index: /trunk/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateChips.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/UpdateChips.c (revision 30616) @@ -1,31 +1,61 @@ # include "relastro.h" +int plotChipFits (double *Ro, double *Do, char *mode, int Nimage); +int saveCenter (Image *image, double *Ro, double *Do, int im); int UpdateChips (Catalog *catalog, int Ncatalog) { + int Nskip, Nmosaic, NnewFit, NoldFit; + /* we can measure new image parameters for each non-mosaic chip independently */ - off_t i, Nimage, Nraw, Nref; + off_t i, Nimage, Nraw, Nref, nFitAstr; Image *image; StarData *raw, *ref; Coords *oldCoords; + float dXpixSys, dYpixSys; + double *Ro, *Do; + char *mode; + + Nskip = Nmosaic = NnewFit = NoldFit = 0; image = getimages (&Nimage); + // save fit results for summary plot + ALLOCATE (Ro, double, Nimage); + ALLOCATE (Do, double, Nimage); + ALLOCATE (mode, char, Nimage); + for (i = 0; i < Nimage; i++) { /* skip all except WRP images */ - if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue; + if (strcmp(&image[i].coords.ctype[4], "-WRP")) { + Nmosaic ++; + mode[i] = 0; + continue; + } /* convert measure coordinates to raw entries */ raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC); - if (!raw) continue; + if (!raw) { + Nskip ++; + mode[i] = 0; + continue; + } /* convert average coordinates to ref entries */ ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC); - if (!ref) continue; + if (!ref) { + Nskip ++; + mode[i] = 0; + continue; + } // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions. assert (Nraw == Nref); + // save these in case of failure saveCoords (&image[i].coords, i); + dXpixSys = image[i].dXpixSys; + dYpixSys = image[i].dYpixSys; + nFitAstr = image[i].nFitAstrom; // FitChip does iterative, clipped fitting @@ -33,6 +63,15 @@ if (!FitChip (raw, ref, Nraw, &image[i])) { if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + + // restore status quo ante oldCoords = getCoords (i); memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + image[i].dXpixSys = dXpixSys; + image[i].dYpixSys = dYpixSys; + image[i].nFitAstrom = nFitAstr; + + saveCenter (image, &Ro[i], &Do[i], i); + mode[i] = 1; + NoldFit ++; free (raw); free (ref); @@ -42,14 +81,178 @@ if (!checkStarMap (i)) { if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + // restore status quo ante oldCoords = getCoords (i); memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + image[i].dXpixSys = dXpixSys; + image[i].dYpixSys = dYpixSys; + image[i].nFitAstrom = nFitAstr; + + saveCenter (image, &Ro[i], &Do[i], i); + mode[i] = 2; image[i].flags |= ID_IMAGE_ASTROM_POOR; - } - + NoldFit ++; + free (raw); + free (ref); + continue; + } + + saveCenter (image, &Ro[i], &Do[i], i); + mode[i] = 3; + NnewFit ++; free (raw); free (ref); } + plotChipFits (Ro, Do, mode, Nimage); + + fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic); return (TRUE); } +int saveCenter (Image *image, double *Ro, double *Do, int im) { + + Mosaic *mosaic; + Coords *moscoords, *imcoords; + double X, Y, L, M, P, Q, R, D; + + moscoords = NULL; + if (!strcmp(&image[im].coords.ctype[4], "-WRP")) { + mosaic = getMosaicForImage (im); + if (mosaic == NULL) return FALSE; // if we cannot find the associated image, skip it + moscoords = &mosaic[0].coords; + } + imcoords = &image[im].coords; + + if (!strcmp(&image[im].coords.ctype[4], "-WRP")) { + X = 0.5*image[im].NX; + Y = 0.5*image[im].NY; + } else { + X = 0.0; + Y = 0.0; + } + + if (moscoords == NULL) { + // this is a Simple image (not a mosaic) + // note that for a Simple image, L,M = P,Q + XY_to_LM (&L, &M, X, Y, imcoords); + LM_to_RD (&R, &D, L, M, imcoords); + } else { + XY_to_LM (&L, &M, X, Y, imcoords); + XY_to_LM (&P, &Q, L, M, moscoords); + LM_to_RD (&R, &D, P, Q, moscoords); + } + + double Rmid; + if (UserCatalog) { + Rmid = UserCatalogRA; + } else { + Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax); + } + + R = ohana_normalize_angle_to_midpoint (R, Rmid); + + *Ro = R; + *Do = D; + + return (TRUE); +} + +int plotChipFits (double *Ro, double *Do, char *mode, int Nimage) { + + static int kapa = -1; + + int i, N; + double Rmin, Rmax, Dmin, Dmax; + float *xvec, *yvec; + Graphdata graphdata; + + if (!relastroGetVisual()) return (TRUE); + + if (kapa == -1) { + kapa = KapaOpenNamedSocket("kapa", "relastro"); + if (kapa == -1) { + fprintf (stderr, "can't open kapa window\n"); + return FALSE; + } + } + + Rmin = +720; + Rmax = -720; + Dmin = +90; + Dmax = -90; + + // find the R, D range + for (i = 0; i < Nimage; i++) { + if (!mode[i]) continue; + + Rmin = MIN(Rmin, Ro[i]); + Rmax = MAX(Rmax, Ro[i]); + Dmin = MIN(Dmin, Do[i]); + Dmax = MAX(Dmax, Do[i]); + } + + ALLOCATE (xvec, float, Nimage); + ALLOCATE (yvec, float, Nimage); + + bzero (&graphdata, sizeof(Graphdata)); + plot_defaults (&graphdata); + graphdata.xmin = Rmin; + graphdata.xmax = Rmax; + graphdata.ymin = Dmin; + graphdata.ymax = Dmax; + graphdata.style = 2; + graphdata.size = 1; + + KapaSetFont (kapa, "helvetica", 14); + KapaSetLimits (kapa, &graphdata); + KapaBox (kapa, &graphdata); + + // *** good images *** + N = 0; + for (i = 0; i < Nimage; i++) { + if (mode[i] != 3) continue; + xvec[N] = Ro[i]; + yvec[N] = Do[i]; + N++; + } + graphdata.ptype = 7; + graphdata.color = KapaColorByName("black"); + KapaPrepPlot (kapa, N, &graphdata); + KapaPlotVector (kapa, N, xvec, "x"); + KapaPlotVector (kapa, N, yvec, "y"); + + // *** reject fit *** + N = 0; + for (i = 0; i < Nimage; i++) { + if (mode[i] != 1) continue; + xvec[N] = Ro[i]; + yvec[N] = Do[i]; + N++; + } + graphdata.ptype = 3; + graphdata.color = KapaColorByName("red"); + KapaPrepPlot (kapa, N, &graphdata); + KapaPlotVector (kapa, N, xvec, "x"); + KapaPlotVector (kapa, N, yvec, "y"); + + // *** divergent fit *** + N = 0; + for (i = 0; i < Nimage; i++) { + if (mode[i] != 2) continue; + xvec[N] = Ro[i]; + yvec[N] = Do[i]; + N++; + } + graphdata.ptype = 2; + graphdata.color = KapaColorByName("blue"); + KapaPrepPlot (kapa, N, &graphdata); + KapaPlotVector (kapa, N, xvec, "x"); + KapaPlotVector (kapa, N, yvec, "y"); + + free (xvec); + free (yvec); + + return (TRUE); +} + +// XXX if (!FindMosaicForImage (image, Nimage, i)) { } Index: /trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 30616) @@ -40,4 +40,6 @@ UpdateMeasures (&catalog, 1); + UpdateObjects (&catalog, 1); + freeImageBins (1); Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 30616) @@ -103,5 +103,5 @@ //does the measurement pass the supplied filtering constraints? - if (!MeasFilterTest(&catalog[i].measure[m])) { + if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) { catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; continue; @@ -134,9 +134,18 @@ // dX, dY : error in arcsec -- - // dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); - // dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC); - - dX[N] = 0.1; - dY[N] = 0.1; + dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); + dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC); + + // add systematic error in quadrature, if desired + // only do this after the fit has converged (or you will never improve the poor images) + // if (INCLUDE_SYS_ERR) { + // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys); + // dX[N] = hypot(dX[N], dRsys); + // dY[N] = hypot(dY[N], dRsys); + // } + + // dX[N] = 0.1; + // dY[N] = 0.1; + dT[N] = catalog[i].measure[m].dt; Index: /trunk/Ohana/src/relastro/src/args.c =================================================================== --- /trunk/Ohana/src/relastro/src/args.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/args.c (revision 30616) @@ -13,19 +13,20 @@ remove_argument (N, &argc, argv); FIT_TARGET = TARGET_OBJECTS; - - // check for object fitting modes (not valid for images) - if ((N = get_argument (argc, argv, "-pm"))) { - remove_argument (N, &argc, argv); - FIT_MODE = FIT_PM_ONLY; - } - if ((N = get_argument (argc, argv, "-par"))) { - remove_argument (N, &argc, argv); - FIT_MODE = FIT_PAR_ONLY; - } - if ((N = get_argument (argc, argv, "-pmpar"))) { - remove_argument (N, &argc, argv); - FIT_MODE = FIT_PM_AND_PAR; - } - } + } + + // check for object fitting modes + if ((N = get_argument (argc, argv, "-pm"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PM_ONLY; + } + if ((N = get_argument (argc, argv, "-par"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PAR_ONLY; + } + if ((N = get_argument (argc, argv, "-pmpar"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PM_AND_PAR; + } + if ((N = get_argument (argc, argv, "-high-speed"))) { // XXX include a parallax / no-parallax option @@ -93,4 +94,10 @@ usage (); } + } + + USE_BASIC_CHECK = FALSE; + if ((N = get_argument (argc, argv, "-basic-image-search"))) { + remove_argument (N, &argc, argv); + USE_BASIC_CHECK = TRUE; } @@ -144,4 +151,9 @@ VERBOSE = VERBOSE2 = TRUE; remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-visual"))) { + remove_argument (N, &argc, argv); + relastroSetVisual(TRUE); } Index: /trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relastro/src/bcatalog.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/bcatalog.c (revision 30616) @@ -45,6 +45,6 @@ offset = catalog[0].average[i].measureOffset + j; - // filter objects based on user supplied criteria - if (!MeasFilterTest(&catalog[0].measure[offset])) { + // filter objects based on user supplied criteria, including SIGMA_LIM + if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; continue; @@ -78,5 +78,5 @@ subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; } - + Nmeasure ++; Nm ++; Index: /trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/load_catalogs.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/load_catalogs.c (revision 30616) @@ -26,9 +26,9 @@ pcatalog[0].Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) { + if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE2, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename); exit (1); } - if (VERBOSE && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename); + if (VERBOSE2 && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename); //outlier rejection Index: /trunk/Ohana/src/relastro/src/plotstuff.c =================================================================== --- /trunk/Ohana/src/relastro/src/plotstuff.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/plotstuff.c (revision 30616) @@ -146,3 +146,12 @@ graphdata[0].ymax = dUNDEF; + graphdata[0].ticktextPad = NAN; + graphdata[0].labelPadXm = NAN; + graphdata[0].labelPadXp = NAN; + graphdata[0].labelPadYm = NAN; + graphdata[0].labelPadYp = NAN; + graphdata[0].padXm = NAN; + graphdata[0].padXp = NAN; + graphdata[0].padYm = NAN; + graphdata[0].padYp = NAN; } Index: /trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/relastro.c (revision 30616) @@ -78,4 +78,6 @@ MARKTIME("update chips: %f sec\n", dtime); } + // create summary plots of the process + // relastroVisualSummaryChips(); break; Index: /trunk/Ohana/src/relastro/src/relastroVisual.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastroVisual.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/relastroVisual.c (revision 30616) @@ -11,13 +11,14 @@ #define KAPAY 700 -static int kapa = -1; +static int kapa1 = -1; static int kapa2 = -1; static int kapa3 = -1; +static int kapa4 = -1; static int isVisual = FALSE; -static int plotRawRef = FALSE; -static int plotScatter = FALSE; -static int plotResid = FALSE; -static int plotVector = FALSE; +static int plotRawRef = TRUE; +// static int plotScatter = TRUE; +// static int plotResid = TRUE; +// static int plotVector = TRUE; static int plotOutliers = TRUE; @@ -29,9 +30,8 @@ fprintf(stderr, "Failure to open kapa.\n"); isVisual = 0; - return 0; + return FALSE; } - // KapaResize (*kapid, KAPAX, KAPAY); - } - return 1; + } + return TRUE; } @@ -49,5 +49,13 @@ isVisual = 0; } - return 1; + return TRUE; +} + +void relastroSetVisual(int state) { + isVisual = state; +} + +int relastroGetVisual(void) { + return isVisual; } @@ -70,11 +78,10 @@ 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; + return TRUE; +} + +/** 4-panel plot of x vs dx, y vs dx, x vs dy, y vs dy **/ +int relastroVisualRawRef(int *kapaID, float *x, float *y, float *dx, float *dy, float *dPos, int npts) { + Graphdata graphdata; KapaSection section; @@ -83,5 +90,5 @@ KapaInitGraph(&graphdata); - KapaClearPlots(*kapaID); + KapaClearSections(*kapaID); KapaSetFont(*kapaID, "helvetica", 14); @@ -89,9 +96,11 @@ section.x = 0.0; section.y = 0.0; section.dx = .45, section.dy = .45; + section.bg = KapaColorByName("white"); graphdata.ptype = 7; graphdata.style = 2; + graphdata.etype |= 0x01; KapaSetSection(*kapaID, §ion); - if(!scaleGraphdata(x, xVec, &graphdata, npts)) return 0; + if(!scaleGraphdata(x, dx, &graphdata, npts)) return 0; KapaSetLimits(*kapaID, &graphdata); KapaBox(*kapaID, &graphdata); @@ -100,9 +109,11 @@ KapaPrepPlot(*kapaID, npts, &graphdata); KapaPlotVector(*kapaID, npts, x, "x"); - KapaPlotVector(*kapaID, npts, xVec, "y"); + KapaPlotVector(*kapaID, npts, dx, "y"); + KapaPlotVector(*kapaID, npts, dPos, "dym"); + KapaPlotVector(*kapaID, npts, dPos, "dyp"); section.x = .5; section.y = 0; section.name="1"; KapaSetSection(*kapaID, §ion); - if(!scaleGraphdata(x, yVec, &graphdata, npts)) return 0; + if(!scaleGraphdata(x, dy, &graphdata, npts)) return 0; KapaSetLimits(*kapaID, &graphdata); KapaBox(*kapaID, &graphdata); @@ -111,9 +122,11 @@ KapaPrepPlot(*kapaID, npts, &graphdata); KapaPlotVector(*kapaID, npts, x, "x"); - KapaPlotVector(*kapaID, npts, yVec, "y"); + KapaPlotVector(*kapaID, npts, dy, "y"); + KapaPlotVector(*kapaID, npts, dPos, "dym"); + KapaPlotVector(*kapaID, npts, dPos, "dyp"); section.x = .0; section.y = .5; section.name="2"; KapaSetSection(*kapaID, §ion); - if(!scaleGraphdata(y, xVec, &graphdata, npts)) return 0;; + if(!scaleGraphdata(y, dx, &graphdata, npts)) return 0;; KapaSetLimits(*kapaID, &graphdata); KapaBox(*kapaID, &graphdata); @@ -122,9 +135,11 @@ KapaPrepPlot(*kapaID, npts, &graphdata); KapaPlotVector(*kapaID, npts, y, "x"); - KapaPlotVector(*kapaID, npts, xVec, "y"); + KapaPlotVector(*kapaID, npts, dx, "y"); + KapaPlotVector(*kapaID, npts, dPos, "dym"); + KapaPlotVector(*kapaID, npts, dPos, "dyp"); section.x = .5; section.y = .5; section.name="3"; KapaSetSection(*kapaID, §ion); - if(!scaleGraphdata(y, yVec, &graphdata, npts)) return 0; + if(!scaleGraphdata(y, dy, &graphdata, npts)) return 0; KapaSetLimits(*kapaID, &graphdata); KapaBox(*kapaID, &graphdata); @@ -133,30 +148,32 @@ 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; + KapaPlotVector(*kapaID, npts, dy, "y"); + KapaPlotVector(*kapaID, npts, dPos, "dym"); + KapaPlotVector(*kapaID, npts, dPos, "dyp"); + + return TRUE; +} + + +/** Plot a vector field (circles at vector origin, scaled lines giving vector directions */ +int relastroVisualVectorField(int *kapaID, float *x, float *y, float *dx, float *dy, int npts, double maxVecLength) { Graphdata graphdata; - float singleX[2], singleY[2]; + float *xVec, *yVec; float vecScaleFactor; float graphSize; int i; char plotTitle[50]; - sprintf(plotTitle, "Maximum Vector Size = %5.1e", maxVecLength); - + + if (!npts) return FALSE; if (!initWindow(kapaID)) return 0; + snprintf(plotTitle, 50, "Maximum Vector Size = %5.1e", maxVecLength); + KapaInitGraph(&graphdata); - KapaClearPlots(*kapaID); - if(!scaleGraphdata(x, y, &graphdata, npts)) return 0; + KapaClearSections(*kapaID); + KapaSetFont(*kapaID, "helvetica", 14); + + if (!scaleGraphdata(x, y, &graphdata, npts)) return 0; graphSize = graphdata.xmax - graphdata.xmin; @@ -164,11 +181,7 @@ 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 - + + // fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n", graphSize, maxVecLength, vecScaleFactor); KapaSetFont (*kapaID, "helvetica", 14); @@ -179,93 +192,142 @@ graphdata.ptype = 7; graphdata.style = 2; + graphdata.color = KapaColorByName("black"); KapaPrepPlot(*kapaID, npts, &graphdata); KapaPlotVector(*kapaID, npts, x, "x"); KapaPlotVector(*kapaID, npts, y, "y"); - //plot each vector individually - graphdata.ptype = 0; - graphdata.style = 0; + ALLOCATE (xVec, float, 2*npts); + ALLOCATE (yVec, float, 2*npts); + for (i = 0; i < npts; i++) { + xVec[2*i + 0] = x[i]; + yVec[2*i + 0] = y[i]; + xVec[2*i + 1] = x[i] + dx[i] * vecScaleFactor; + yVec[2*i + 1] = y[i] + dy[i] * vecScaleFactor; + } + + graphdata.ptype = 100; // line segements by point pair + graphdata.style = 2; 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]; - } + KapaPrepPlot (*kapaID, 2*npts, &graphdata); + KapaPlotVector(*kapaID, 2*npts, xVec, "x"); + KapaPlotVector(*kapaID, 2*npts, yVec, "y"); + + free (xVec); + free (yVec); + return TRUE; +} + +int relastroVisualPlotScatter(int *kapaID, float *dXfit, float *dYfit, float *mag, int npts) { Graphdata graphdata; KapaSection section; - section.x = 0; section.y = 0; section.dx = 1; section.dy = 1; - section.name = "junk"; + + if (!initWindow(kapaID)) return 0; 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; -} - - + KapaClearSections(*kapaID); + KapaSetFont(*kapaID, "helvetica", 14); + + section.name = "a"; + section.x = 0.0; section.y = 0.0; section.dx = 1.0; section.dy = 0.5; + section.bg = KapaColorByName("white"); + KapaSetSection(*kapaID, §ion); + + graphdata.ptype = 2; + graphdata.style = 2; + + if(!scaleGraphdata(mag, dXfit, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "mag", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dXfit", KAPA_LABEL_YM); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, mag, "x"); + KapaPlotVector(*kapaID, npts, dXfit, "y"); + + section.name = "b"; + section.x = 0.0; section.y = 0.5; section.dx = 1.0; section.dy = 0.5; + section.bg = KapaColorByName("white"); + KapaSetSection(*kapaID, §ion); + + graphdata.ptype = 2; + graphdata.style = 2; + + if(!scaleGraphdata(mag, dYfit, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "mag", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dYfit", KAPA_LABEL_YM); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, mag, "x"); + KapaPlotVector(*kapaID, npts, dYfit, "y"); + + return TRUE; +} /** 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; +int relastroVisualPlotFittedStars(int *kapaID, float *rawX, float *rawY, float *refX, float *refY, int numNoFit, float *rawXfit, float *rawYfit, float *refXfit, float *refYfit, int numFit) { + + Graphdata graphdata; + + if (!initWindow(kapaID)) return 0; + + KapaInitGraph(&graphdata); + KapaClearPlots(*kapaID); + + graphdata.ptype = 7; + graphdata.style = 2; + + if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) { + fprintf(stderr, "Not enough finite points for plotting"); + return 0; + } + + KapaSetFont(*kapaID, "helvetica", 14); + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel( *kapaID, "X", KAPA_LABEL_XM); + KapaSendLabel( *kapaID, "Y", KAPA_LABEL_YM); + KapaSendLabel( *kapaID, "orange, red, green, blue: (raw, ref), (nofit, fit)", + KAPA_LABEL_XP); + + graphdata.color = KapaColorByName("orange"); + graphdata.size = 1; + KapaPrepPlot(*kapaID, numNoFit, &graphdata); + KapaPlotVector(*kapaID, numNoFit, rawX, "x"); + KapaPlotVector(*kapaID, numNoFit, rawY, "y"); + + graphdata.color = KapaColorByName("red"); + graphdata.size = 2; + KapaPrepPlot(*kapaID, numNoFit, &graphdata); + KapaPlotVector(*kapaID, numNoFit, refX, "x"); + KapaPlotVector(*kapaID, numNoFit, refY, "y"); + + graphdata.color = KapaColorByName("green"); + graphdata.size = 1; + KapaPrepPlot(*kapaID, numFit, &graphdata); + KapaPlotVector(*kapaID, numFit, rawXfit, "x"); + KapaPlotVector(*kapaID, numFit, rawYfit, "y"); + + graphdata.color = KapaColorByName("blue"); + graphdata.size = 2; + KapaPrepPlot(*kapaID, numFit, &graphdata); + KapaPlotVector(*kapaID, numFit, refXfit, "x"); + KapaPlotVector(*kapaID, numFit, refYfit, "y"); + + return TRUE; +} + +/** create various plots using the data in raw and ref **/ +int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj) { float *rawX, *rawY, *refX, *refY; float *rawXfit, *rawYfit, *refXfit, *refYfit; float *magRaw, *magRef, *magRawfit, *magReffit; - float *xVec, *yVec; + float *dXfit, *dYfit, *dPos; int numFit = 0, numNoFit = 0; double dL, dM, dR; + + if (!isVisual) return TRUE; ALLOCATE(rawX, float, numObj); @@ -273,21 +335,28 @@ ALLOCATE(refX, float, numObj); ALLOCATE(refY, float, numObj); + ALLOCATE(magRaw, float, numObj); + ALLOCATE(magRef, 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); + ALLOCATE(dXfit, float, numObj); + ALLOCATE(dYfit, float, numObj); + ALLOCATE(dPos, float, numObj); int i; for(i = 0; i < numObj; i++) { - if (raw[i].mask) continue; + if (raw[i].mask) continue; // XXX dL = raw[i].L - ref[i].L; dM = raw[i].M - ref[i].M; dR = hypot (dL, dM); + + // XXX change the selection to a mask-based thing if (dR > dRmax) { + // UNFITTED values rawX[numNoFit] = raw[i].X; rawY[numNoFit] = raw[i].Y; @@ -298,10 +367,14 @@ 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; + // FITTED values + rawXfit[numFit] = raw[i].X; + rawYfit[numFit] = raw[i].Y; + refXfit[numFit] = ref[i].X; + refYfit[numFit] = ref[i].Y; + magRawfit[numFit] = raw[i].Mag; + magReffit[numFit] = ref[i].Mag; + dPos[numFit] = ref[i].dPos; + dXfit[numFit] = raw[i].X - ref[i].X; + dYfit[numFit] = raw[i].Y - ref[i].Y; numFit++; } @@ -310,69 +383,21 @@ 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); + // 4-panel plot of x vs dx, y vs dx, x vs dy, y vs dy + relastroVisualRawRef(&kapa1, rawXfit, rawYfit, dXfit, dYfit, dPos, numFit); + + // vector line plot for the fit + relastroVisualVectorField(&kapa2, rawXfit, rawYfit, dXfit, dYfit, numFit, dRmax); + + // dXfit, dYfit vs mag + relastroVisualPlotScatter(&kapa3, dXfit, dYfit, magRawfit, numFit); + + // plot the positions of the fitted stars on the chip + relastroVisualPlotFittedStars(&kapa4, rawX, rawY, refX, refY, numNoFit, rawXfit, rawYfit, refXfit, refYfit, numFit); askUser(&plotRawRef); + FREE(dXfit); + FREE(dYfit); + FREE(dPos); FREE(rawX); FREE(rawY); @@ -388,5 +413,5 @@ FREE(magReffit); - return 1; + return TRUE; } @@ -403,6 +428,6 @@ KapaSection section; - if (!isVisual || !plotOutliers) return 1; - if (!initWindow(&kapa)) return 0; + if (!isVisual || !plotOutliers) return TRUE; + if (!initWindow(&kapa1)) return 0; // populate vectors @@ -426,5 +451,5 @@ for(i = 0; i < Nmeasure; i++, m++) { meas = catalog[0].measure[m]; - if (!MeasFilterTest(&meas)) continue; + if (!MeasFilterTest(&meas, FALSE)) continue; xmin = MIN(xmin, meas.dR); xmax = MAX(xmax, meas.dR); @@ -459,8 +484,9 @@ section.x = 0; section.y = 0; section.dx = 1; section.dy = 1; section.name = "junk"; + section.bg = KapaColorByName("white"); KapaInitGraph(&graphdata); - KapaClearPlots(kapa); - KapaSetFont(kapa, "helvetica", 14); + KapaClearPlots(kapa1); + KapaSetFont(kapa1, "helvetica", 14); graphdata.ptype = 7; @@ -472,32 +498,30 @@ 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)", + KapaSetSection(kapa1, §ion); + KapaSetLimits(kapa1, &graphdata); + KapaBox(kapa1, &graphdata); + + KapaSendLabel( kapa1, "RA (arcsec)", KAPA_LABEL_XM); + KapaSendLabel( kapa1, "Dec (arcsec)", KAPA_LABEL_YM); + KapaSendLabel( kapa1, "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"); + KapaPrepPlot(kapa1, Nin, &graphdata); + KapaPlotVector(kapa1, Nin, Rin, "x"); + KapaPlotVector(kapa1, Nin, Din, "y"); graphdata.color = KapaColorByName("red"); - KapaPrepPlot(kapa, Nout, &graphdata); - KapaPlotVector(kapa, Nout, Rout, "x"); - KapaPlotVector(kapa, Nout, Dout, "y"); + KapaPrepPlot(kapa1, Nout, &graphdata); + KapaPlotVector(kapa1, Nout, Rout, "x"); + KapaPlotVector(kapa1, 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"); - - - + KapaPrepPlot(kapa1, 100, &graphdata); + KapaPlotVector(kapa1, 100, xCirc, "x"); + KapaPlotVector(kapa1, 100, yCirc, "y"); + askUser(&plotOutliers); @@ -510,3 +534,14 @@ } - +# if (0) +int relastroVisualSummaryChips() { + + // plot the dXsys, dYsys histograms + + // plot x vs dx, y vs dy, etc for all mosaics + + // plot a map of median star scatter + +} +# endif + Index: /trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- /trunk/Ohana/src/relastro/src/select_images.c (revision 30615) +++ /trunk/Ohana/src/relastro/src/select_images.c (revision 30616) @@ -15,4 +15,10 @@ void dsortindex (double *X, off_t *Y, int N); off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) { @@ -25,5 +31,8 @@ Coords tcoords; SkyRegionCoords *skycoords; - + struct timeval start, stop; + + double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; + double *RmaxSky; off_t *index; @@ -40,4 +49,6 @@ return NULL; } + + gettimeofday (&start, (void *) NULL); // the comparison is made in the catalog local projection. below we set crval1,2 @@ -52,4 +63,9 @@ ALLOCATE (RmaxSky, double, skylist[0].Nregions); ALLOCATE (index, off_t, skylist[0].Nregions); + + RminSkyRegion = +360.0; + RmaxSkyRegion = -360.0; + DminSkyRegion = +90.0; + DmaxSkyRegion = -90.0; /* compare with each region file */ @@ -80,10 +96,19 @@ skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; - } + + RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin); + RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax); + DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin); + DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax); + } + RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); + MARKTIME("create sky region coords: %f sec\n", dtime); dsortindex (RmaxSky, index, skylist[0].Nregions); + MARKTIME("sort sky coords: %f sec\n", dtime); if (VERBOSE) fprintf (stderr, "finding images\n"); BuildChipMatch (timage, Ntimage); + MARKTIME("build chip match: %f sec\n", dtime); nimage = 0; @@ -124,18 +149,51 @@ } - /* define image corners */ - Xi[0] = 0; Yi[0] = 0; - Xi[1] = timage[i].NX; Yi[1] = 0; - Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; - Xi[3] = 0; Yi[3] = timage[i].NY; - Xi[4] = 0; Yi[4] = 0; + /* define image corners - note the DIS images (mosaic phu) are special */ + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) { + Xi[0] = -0.5*timage[i].NX; Yi[0] = -0.5*timage[i].NY; + Xi[1] = +0.5*timage[i].NX; Yi[1] = -0.5*timage[i].NY; + Xi[2] = +0.5*timage[i].NX; Yi[2] = +0.5*timage[i].NY; + Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY; + Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY; + } else { + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; + } found = FALSE; /* transform corners to ra,dec */ double RminImage = 360.0; + double RmaxImage = 0.0; + double DminImage = +90.0; + double DmaxImage = -90.0; + // int leftside = FALSE; for (j = 0; j < 5; j++) { XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion); + RminImage = MIN(RminImage, Ri[j]); - } + RmaxImage = MAX(RmaxImage, Ri[j]); + DminImage = MIN(DminImage, Di[j]); + DmaxImage = MAX(DmaxImage, Di[j]); + } + if (RmaxImage - RminImage > 180.0) { + double tmp = RminImage; + RmaxImage = RminImage; + RminImage = tmp - 360.0; + } + + // check that this image is even in range of the searched region + if (DminImage > DmaxSkyRegion) continue; + if (DmaxImage < DminSkyRegion) continue; + + // the sky region RA is defined to be 0 - 360.0 + if (RminImage > RmaxSkyRegion) continue; + if (RmaxImage < RminSkyRegion) continue; + + // image overlaps region, keep it + if (USE_BASIC_CHECK) goto found_it; // RA(nStart) is guaranteed to be < RminImage: @@ -161,8 +219,10 @@ for (j = 0; (j < 4) && !found; j++) { found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]); + if (found) goto found_it; } /* check if catalog corner inside image */ for (j = 0; (j < 4) && !found; j++) { found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]); + if (found) goto found_it; } /* check if edges cross */ @@ -170,29 +230,33 @@ for (k = 0; (k < 4) && !found; k++) { found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]); + if (found) goto found_it; } } - if (!found) continue; - - image[nimage] = timage[i]; - /* always allow 'few' images to succeed, if possible */ - if (image[nimage].flags & ID_IMAGE_ASTROM_FEW) { - image[nimage].flags &= ~ID_IMAGE_ASTROM_FEW; - } - if (RESET) { - // XXX do we need / want to do this in relastro? - assignMcal (&image[nimage], (double *) NULL, -1); - image[nimage].dMcal = NAN; - image[nimage].flags &= ~badImage; - } - line_number[nimage] = i; - nimage ++; - if (nimage == NIMAGE) { - NIMAGE += 100; - REALLOCATE (image, Image, NIMAGE); - REALLOCATE (line_number, off_t, NIMAGE); - } - } - } - + } + if (!found) continue; + + found_it: + image[nimage] = timage[i]; + /* always allow 'few' images to succeed, if possible */ + if (image[nimage].flags & ID_IMAGE_ASTROM_FEW) { + image[nimage].flags &= ~ID_IMAGE_ASTROM_FEW; + } + if (RESET) { + // XXX do we need / want to do this in relastro? + // assignMcal (&image[nimage], (double *) NULL, -1); + // image[nimage].Mcal = NAN; + // image[nimage].dMcal = NAN; + image[nimage].flags &= ~badImage; + } + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += 100; + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, off_t, NIMAGE); + } + } + MARKTIME("finish image selection: %f sec\n", dtime); + if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); Index: /trunk/Ohana/src/relphot/Makefile =================================================================== --- /trunk/Ohana/src/relphot/Makefile (revision 30615) +++ /trunk/Ohana/src/relphot/Makefile (revision 30616) @@ -29,4 +29,5 @@ $(SRC)/StarOps.$(ARCH).o \ $(SRC)/args.$(ARCH).o \ +$(SRC)/help.$(ARCH).o \ $(SRC)/bcatalog.$(ARCH).o \ $(SRC)/global_stats.$(ARCH).o \ Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 30615) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 30616) @@ -13,4 +13,6 @@ float Mcal; float dMcal; + float dMsys; + short nFitPhotom; short Xm; float secz; @@ -50,4 +52,5 @@ double STAR_CHISQ; double MIN_ERROR; +double IMFIT_SYS_SIGMA_LIM; int VERBOSE; @@ -103,4 +106,6 @@ SkyRegion UserPatch; int UserPatchSelect; + +int USE_BASIC_CHECK; # ifdef GRID_V1 @@ -215,2 +220,5 @@ void write_coords PROTO((Header *header, Coords *coords)); int relphot_objects (void); + +void relphot_usage (void); +void relphot_help (int argc, char **argv); Index: /trunk/Ohana/src/relphot/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relphot/src/ConfigInit.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/ConfigInit.c (revision 30616) @@ -20,4 +20,7 @@ GetConfig (config, "MAG_LIM", "%lf", 0, &MAG_LIM); GetConfig (config, "SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + + if (!ScanConfig (config, "RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01; + GetConfig (config, "STAR_SCATTER", "%lf", 0, &STAR_SCATTER); GetConfig (config, "IMAGE_SCATTER", "%lf", 0, &IMAGE_SCATTER); Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 30616) @@ -124,4 +124,6 @@ free (clist); free (mlist); + free (Nlist); + free (NLIST); } @@ -299,6 +301,6 @@ off_t i, j, m, c, n, N, Nmax, mark, bad; - int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys; - float Msys, Mrel, Mmos, Mgrid; + int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright; + float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2; double *list, *dlist; StatType stats; @@ -318,4 +320,8 @@ Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0; + + // counters to measure the bright-end scatter + McalBright = McalBright2 = 0.0; + Nbright = 0; for (i = 0; i < Nimage; i++) { @@ -361,4 +367,9 @@ list[N] = Msys - Mrel - Mmos - Mgrid; dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR); + if (catalog[c].measure[m].dM < IMFIT_SYS_SIGMA_LIM) { + McalBright += list[N]; + McalBright2 += SQ(list[N]); + Nbright ++; + } N++; } @@ -378,5 +389,7 @@ liststats (list, dlist, N, &stats); image[i].Mcal = stats.mean; - image[i].dMcal = stats.sigma; + image[i].dMcal = stats.error; + image[i].dMagSys = stats.sigma; + image[i].nFitPhotom = N; image[i].Xm = 100.0*log10(stats.chisq); } Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 30616) @@ -73,4 +73,5 @@ mosaic[Nmosaic].Mcal = 0.0; mosaic[Nmosaic].dMcal = 0.0; + mosaic[Nmosaic].dMsys = 0.0; mosaic[Nmosaic].Xm = 0.0; mosaic[Nmosaic].flags = image[i].flags; @@ -189,4 +190,6 @@ image[im].dMcal = mosaic[i].dMcal; image[im].Xm = mosaic[i].Xm; + image[im].dMagSys = mosaic[i].dMsys; + image[im].nFitPhotom = mosaic[i].nFitPhotom; image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW); image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR); @@ -400,5 +403,5 @@ mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]); if (mark) { - fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]); + if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]); } mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW; Nfew ++; @@ -410,5 +413,7 @@ if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas, N); mosaic[i].Mcal = stats.mean; - mosaic[i].dMcal = stats.sigma; + mosaic[i].dMcal = stats.error; + mosaic[i].dMsys = stats.sigma; + mosaic[i].nFitPhotom = N; mosaic[i].Xm = 100.0*log10(stats.chisq); } Index: /trunk/Ohana/src/relphot/src/args.c =================================================================== --- /trunk/Ohana/src/relphot/src/args.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/args.c (revision 30616) @@ -1,4 +1,3 @@ # include "relphot.h" -void usage (void); int args (int argc, char **argv) { @@ -53,4 +52,10 @@ } + USE_BASIC_CHECK = FALSE; + if ((N = get_argument (argc, argv, "-basic-image-search"))) { + remove_argument (N, &argc, argv); + USE_BASIC_CHECK = TRUE; + } + VERBOSE = VERBOSE2 = FALSE; if ((N = get_argument (argc, argv, "-v"))) { @@ -228,32 +233,6 @@ if (UpdateAverages && (argc == 1)) return TRUE; if (UserPatchSelect && (argc == 2)) return TRUE; - if (argc != 3) usage (); + if (argc != 3) relphot_usage (); return TRUE; } - -void usage () { - fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); - fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\n"); - fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); - fprintf (stderr, " options: \n"); - fprintf (stderr, " -time (start) (stop)\n"); - fprintf (stderr, " -v : verbose output\n"); - fprintf (stderr, " -vv : more verbose output\n"); - fprintf (stderr, " -outroot (outroot)\n"); - fprintf (stderr, " -plot\n"); - fprintf (stderr, " -plotdelay (seconds)\n"); - fprintf (stderr, " -statmode (mode)\n"); - fprintf (stderr, " -refcode (name) : give extra weight to this photcode\n"); - fprintf (stderr, " -n (nloop)\n"); - fprintf (stderr, " -reset\n"); - fprintf (stderr, " -update\n"); - fprintf (stderr, " -params\n"); - fprintf (stderr, " -mosaic (mosaic)\n"); - fprintf (stderr, " -imfreeze\n"); - fprintf (stderr, " -grid\n"); - fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); - fprintf (stderr, " -instmag min max\n"); - fprintf (stderr, " \n"); - exit (2); -} Index: /trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/bcatalog.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/bcatalog.c (revision 30616) @@ -119,6 +119,6 @@ if (VERBOSE) { - fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog\n", - subcatalog[0].Naverage, subcatalog[0].Nmeasure, i); + fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", + subcatalog[0].Naverage, subcatalog[0].Nmeasure, i, catalog[0].filename); fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few\n", Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew); Index: /trunk/Ohana/src/relphot/src/help.c =================================================================== --- /trunk/Ohana/src/relphot/src/help.c (revision 30616) +++ /trunk/Ohana/src/relphot/src/help.c (revision 30616) @@ -0,0 +1,44 @@ +# include "relphot.h" + +void relphot_usage (void) { + fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); + fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\n"); + fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); + fprintf (stderr, " use -h for more usage information\n"); + exit (2); +} + +void relphot_help (int argc, char **argv) { + + /* check for help request */ + if (get_argument (argc, argv, "-help")) goto show_help; + if (get_argument (argc, argv, "-h")) goto show_help; + if (argc == 1) relphot_usage(); + return; + +show_help: + fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); + fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\n"); + fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); + fprintf (stderr, " options: \n"); + fprintf (stderr, " -time (start) (stop)\n"); + fprintf (stderr, " -v : verbose output\n"); + fprintf (stderr, " -vv : more verbose output\n"); + fprintf (stderr, " -outroot (outroot)\n"); + fprintf (stderr, " -plot\n"); + fprintf (stderr, " -plotdelay (seconds)\n"); + fprintf (stderr, " -statmode (mode)\n"); + fprintf (stderr, " -refcode (name) : give extra weight to this photcode\n"); + fprintf (stderr, " -n (nloop)\n"); + fprintf (stderr, " -reset\n"); + fprintf (stderr, " -update\n"); + fprintf (stderr, " -params\n"); + fprintf (stderr, " -mosaic (mosaic)\n"); + fprintf (stderr, " -imfreeze\n"); + fprintf (stderr, " -grid\n"); + fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); + fprintf (stderr, " -instmag min max\n"); + fprintf (stderr, " \n"); + exit (2); +} + Index: /trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- /trunk/Ohana/src/relphot/src/initialize.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/initialize.c (revision 30616) @@ -5,4 +5,5 @@ int N; + relphot_help (argc, argv); ConfigInit (&argc, argv); args (argc, argv); Index: /trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- /trunk/Ohana/src/relphot/src/load_catalogs.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/load_catalogs.c (revision 30616) @@ -6,5 +6,5 @@ Catalog *catalog, tcatalog; - if (VERBOSE) fprintf (stderr, "loading catalog data\n"); + if (VERBOSE2) fprintf (stderr, "loading catalog data\n"); ALLOCATE (catalog, Catalog, skylist[0].Nregions); @@ -22,9 +22,15 @@ tcatalog.Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog - if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) { + if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE2, "r")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename); exit (1); } - if (VERBOSE && !tcatalog.Naves_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); + if (!tcatalog.Naves_disk) { + if (VERBOSE2) { fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); } + dvo_catalog_unlock (&tcatalog); + dvo_catalog_free (&tcatalog); + continue; + } + Nstar_total += tcatalog.Naverage; Nmeas_total += tcatalog.Nmeasure; Index: /trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/load_images.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/load_images.c (revision 30616) @@ -1,3 +1,9 @@ # include "relphot.h" + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) { @@ -6,7 +12,10 @@ off_t Nimage, Nsubset, Nchar; off_t *LineNumber; + struct timeval start, stop; SkyTable *sky = NULL; SkyList *skylist = NULL; + + gettimeofday (&start, (void *) NULL); // load the current sky table (layout of all SkyRegions) @@ -29,12 +38,18 @@ exit (2); } + MARKTIME("read image table: %f sec\n", dtime); // select the images which overlap the selected sky regions subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + MARKTIME("selected images: %f sec\n", dtime); gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); + MARKTIME("converted ftable to vtable: %f sec\n", dtime); initImages (subset, Nsubset); + MARKTIME("init images: %f sec\n", dtime); + initMosaics (subset, Nsubset); + MARKTIME("init mosaics: %f sec\n", dtime); return (skylist); Index: /trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/relphot.c (revision 30616) @@ -85,5 +85,5 @@ if (PLOTSTUFF) { plot_star_coords (catalog, Ncatalog); - plot_mosaic_fields (catalog); + // plot_mosaic_fields (catalog); } Index: /trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/select_images.c (revision 30615) +++ /trunk/Ohana/src/relphot/src/select_images.c (revision 30616) @@ -15,4 +15,10 @@ void dsortindex (double *X, off_t *Y, int N); off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) { @@ -25,5 +31,8 @@ Coords tcoords; SkyRegionCoords *skycoords; - + struct timeval start, stop; + + double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion, RmidSkyRegion; + double *RmaxSky; off_t *index; @@ -35,4 +44,6 @@ return NULL; } + + gettimeofday (&start, (void *) NULL); // the comparison is made in the catalog local projection. below we set crval1,2 @@ -47,4 +58,9 @@ ALLOCATE (RmaxSky, double, skylist[0].Nregions); ALLOCATE (index, off_t, skylist[0].Nregions); + + RminSkyRegion = +360.0; + RmaxSkyRegion = -360.0; + DminSkyRegion = +90.0; + DmaxSkyRegion = -90.0; /* compare with each region file */ @@ -75,10 +91,19 @@ skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; - } + + RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin); + RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax); + DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin); + DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax); + } + RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); + MARKTIME("create sky region coords: %f sec\n", dtime); dsortindex (RmaxSky, index, skylist[0].Nregions); + MARKTIME("sort sky coords: %f sec\n", dtime); if (VERBOSE) fprintf (stderr, "finding images\n"); BuildChipMatch (timage, Ntimage); + MARKTIME("build chip match: %f sec\n", dtime); nimage = 0; @@ -100,4 +125,5 @@ } + // this adds 1.3 sec for 3M images if (!FindMosaicForImage (timage, Ntimage, i)) { fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); @@ -113,16 +139,40 @@ found = FALSE; - /* transform corners to ra,dec */ + /* transform corners to ra,dec -- costs ~3sec for 3M images */ double RminImage = 360.0; + double RmaxImage = 0.0; + double DminImage = +90.0; + double DmaxImage = -90.0; for (j = 0; j < 5; j++) { XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion); + RminImage = MIN(RminImage, Ri[j]); - } - - // RA(nStart) is guaranteed to be < RminImage: + RmaxImage = MAX(RmaxImage, Ri[j]); + DminImage = MIN(DminImage, Di[j]); + DmaxImage = MAX(DmaxImage, Di[j]); + } + if (RmaxImage - RminImage > 180.0) { + double tmp = RminImage; + RmaxImage = RminImage; + RminImage = tmp - 360.0; + } + + // check that this image is even in range of the searched region + if (DminImage > DmaxSkyRegion) continue; + if (DmaxImage < DminSkyRegion) continue; + + // the sky region RA is defined to be 0 - 360.0 + if (RminImage > RmaxSkyRegion) continue; + if (RmaxImage < RminSkyRegion) continue; + + // image overlaps region, keep it + if (USE_BASIC_CHECK) goto found_it; + + // RA(nStart) is guaranteed to be < RminImage: -- costs 0.5sec for 3M images nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions); /* compare with each region file */ - for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) { + for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { m = index[iSky]; @@ -142,8 +192,10 @@ for (j = 0; (j < 4) && !found; j++) { found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]); + if (found) goto found_it; } /* check if catalog corner inside image */ for (j = 0; (j < 4) && !found; j++) { found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]); + if (found) goto found_it; } /* check if edges cross */ @@ -151,28 +203,33 @@ for (k = 0; (k < 4) && !found; k++) { found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]); + if (found) goto found_it; } } - if (!found) continue; - - image[nimage] = timage[i]; - /* always allow 'few' images to succeed, if possible */ - if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { - image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR); - } - if (RESET) { - assignMcal (&image[nimage], (double *) NULL, -1); - image[nimage].dMcal = NAN; - image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR; - } - line_number[nimage] = i; - nimage ++; - if (nimage == NIMAGE) { - NIMAGE += 100; - REALLOCATE (image, Image, NIMAGE); - REALLOCATE (line_number, off_t, NIMAGE); - } - } - } - + } + if (!found) continue; + + found_it: + image[nimage] = timage[i]; + /* always allow 'few' images to succeed, if possible */ + if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { + image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR); + } + if (RESET) { + // XXX assignMcal (&image[nimage], (double *) NULL, -1); + // XXX this needs to be thought through a bit more + image[nimage].Mcal = 0.0; + image[nimage].dMcal = NAN; + image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR; + } + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += 100; + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, off_t, NIMAGE); + } + } + MARKTIME("finish image selection: %f sec\n", dtime); + if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); @@ -295,2 +352,21 @@ return (Nlo); } + +off_t getRegionStopByRA (double R, double *Rref, off_t Nregions) { + + // use bisection to find the overlapping mosaic + + off_t Nlo, Nhi, N; + + // find the last mosaic before start + Nlo = 0; Nhi = Nregions; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (Rref[N] < R) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Nregions); + } + } + return (Nlo); +}