Index: branches/eam_branches/20090522/Ohana/src/addstar/include/addstar.h =================================================================== --- branches/eam_branches/20090522/Ohana/src/addstar/include/addstar.h (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/addstar/include/addstar.h (revision 24557) @@ -263,4 +263,6 @@ Stars *Convert_PS1_DEV_1 PROTO((FTable *table, int *nstars)); Stars *Convert_PS1_V1 PROTO((FTable *table, int *nstars)); +Stars *Convert_PS1_V1_Alt PROTO((FTable *table, int *nstars)); +Stars *Convert_PS1_V2 PROTO((FTable *table, int *nstars)); int InitStar PROTO((Stars *star)); Index: branches/eam_branches/20090522/Ohana/src/addstar/src/MatchHeaders.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/addstar/src/MatchHeaders.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/addstar/src/MatchHeaders.c (revision 24557) @@ -52,4 +52,5 @@ if (!strcmp (exttype, "PS1_DEV_1")) goto keep; if (!strcmp (exttype, "PS1_V1")) goto keep; + if (!strcmp (exttype, "PS1_V2")) goto keep; continue; Index: branches/eam_branches/20090522/Ohana/src/addstar/src/ReadStarsFITS.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/addstar/src/ReadStarsFITS.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/addstar/src/ReadStarsFITS.c (revision 24557) @@ -43,4 +43,7 @@ stars = Convert_PS1_V1 (&table, &Nstars); } + if (!strcmp (type, "PS1_V2")) { + stars = Convert_PS1_V2 (&table, &Nstars); + } if (stars == NULL) { fprintf (stderr, "ERROR: invalid table type %s\n", type); @@ -188,4 +191,17 @@ Stars *stars; CMF_PS1_V1 *ps1data; + + // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec. this was a + // mistake in two ways: a few format should have been defined (eg, CMF_PS1_V2), and the + // layout used did not have clean byte-boundaries for the corresponding structure. The + // former means we have two varieties of CMF_PS1_V1 out there; the latter means that the + // autocode tools do not work to read in the new version, even if we recognize it. Here we + // test for the existence of the broken version (table[0].headers[0].Naxis[0] == 136), and + // call a special conversion function if it is found. + + if (table[0].header[0].Naxis[0] == 136) { + stars = Convert_PS1_V1_Alt (table, nstars); + return (stars); + } ps1data = gfits_table_get_CMF_PS1_V1 (table, &Nstars, NULL); @@ -244,2 +260,136 @@ return (stars); } + +Stars *Convert_PS1_V1_Alt (FTable *table, int *nstars) { + + int i, Nstars; + double ZeroPt; + Stars *stars; + CMF_PS1_V2 *ps1data; + + // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec. this was a + // mistake in two ways: a few format should have been defined (eg, CMF_PS1_V2), and the + // layout used did not have clean byte-boundaries for the corresponding structure. The + // former means we have two varieties of CMF_PS1_V1 out there; the latter means that the + // autocode tools do not work to read in the new version, even if we recognize it. Here we + // test for the existence of the broken version (table[0].headers[0].Naxis[0] == 136), and + // call a special conversion function if it is found. + + ps1data = gfits_table_get_CMF_PS1_V1_Alt (table, &Nstars, NULL); + ZeroPt = GetZeroPoint(); + + ALLOCATE (stars, Stars, Nstars); + for (i = 0; i < Nstars; i++) { + InitStar (&stars[i]); + stars[i].measure.Xccd = ps1data[i].X; + stars[i].measure.Yccd = ps1data[i].Y; + stars[i].measure.dXccd = ShortPixels(ps1data[i].dX); + stars[i].measure.dYccd = ShortPixels(ps1data[i].dY); + + stars[i].measure.posangle = ShortDegree(ps1data[i].posangle); + stars[i].measure.pltscale = ps1data[i].pltscale; + + if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + } + stars[i].measure.dM = ps1data[i].dM; + stars[i].measure.dMcal = ps1data[i].dMcal; + stars[i].measure.Map = ps1data[i].Map + ZeroPt; + + stars[i].measure.Sky = ps1data[i].sky; + stars[i].measure.dSky = ps1data[i].dSky; + + stars[i].measure.psfChisq = ps1data[i].psfChisq; + stars[i].measure.psfQual = ps1data[i].psfQual; + stars[i].measure.psfNdof = ps1data[i].psfNdof; + stars[i].measure.psfNpix = ps1data[i].psfNpix; + stars[i].measure.crNsigma = ps1data[i].crNsigma; + stars[i].measure.extNsigma = ps1data[i].extNsigma; + + stars[i].measure.FWx = ShortPixels(ps1data[i].fx); + stars[i].measure.FWy = ShortPixels(ps1data[i].fy); + stars[i].measure.theta = ShortDegree(ps1data[i].df); + + stars[i].measure.Mxx = ShortPixels(ps1data[i].Mxx); + stars[i].measure.Mxy = ShortPixels(ps1data[i].Mxy); + stars[i].measure.Myy = ShortPixels(ps1data[i].Myy); + + stars[i].measure.photFlags = ps1data[i].flags; + + // this is may optionally be replaced by the internal sequence (see FilterStars.c) + stars[i].measure.detID = ps1data[i].detID; + + // the Average fields and the following Measure fields are set in FilterStars after + // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID, + // averef is set in find_matches, dbFlags is zero on ingest. + + // the following fields are currently not being set anywhere: t_msec + } + *nstars = Nstars; + return (stars); +} + +Stars *Convert_PS1_V2 (FTable *table, int *nstars) { + + int i, Nstars; + double ZeroPt; + Stars *stars; + CMF_PS1_V2 *ps1data; + + ps1data = gfits_table_get_CMF_PS1_V2 (table, &Nstars, NULL); + ZeroPt = GetZeroPoint(); + + ALLOCATE (stars, Stars, Nstars); + for (i = 0; i < Nstars; i++) { + InitStar (&stars[i]); + stars[i].measure.Xccd = ps1data[i].X; + stars[i].measure.Yccd = ps1data[i].Y; + stars[i].measure.dXccd = ShortPixels(ps1data[i].dX); + stars[i].measure.dYccd = ShortPixels(ps1data[i].dY); + + stars[i].measure.posangle = ShortDegree(ps1data[i].posangle); + stars[i].measure.pltscale = ps1data[i].pltscale; + + if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + } + stars[i].measure.dM = ps1data[i].dM; + stars[i].measure.dMcal = ps1data[i].dMcal; + stars[i].measure.Map = ps1data[i].Map + ZeroPt; + + stars[i].measure.Sky = ps1data[i].sky; + stars[i].measure.dSky = ps1data[i].dSky; + + stars[i].measure.psfChisq = ps1data[i].psfChisq; + stars[i].measure.psfQual = ps1data[i].psfQual; + stars[i].measure.psfNdof = ps1data[i].psfNdof; + stars[i].measure.psfNpix = ps1data[i].psfNpix; + stars[i].measure.crNsigma = ps1data[i].crNsigma; + stars[i].measure.extNsigma = ps1data[i].extNsigma; + + stars[i].measure.FWx = ShortPixels(ps1data[i].fx); + stars[i].measure.FWy = ShortPixels(ps1data[i].fy); + stars[i].measure.theta = ShortDegree(ps1data[i].df); + + stars[i].measure.Mxx = ShortPixels(ps1data[i].Mxx); + stars[i].measure.Mxy = ShortPixels(ps1data[i].Mxy); + stars[i].measure.Myy = ShortPixels(ps1data[i].Myy); + + stars[i].measure.photFlags = ps1data[i].flags; + + // this is may optionally be replaced by the internal sequence (see FilterStars.c) + stars[i].measure.detID = ps1data[i].detID; + + // the Average fields and the following Measure fields are set in FilterStars after + // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID, + // averef is set in find_matches, dbFlags is zero on ingest. + + // the following fields are currently not being set anywhere: t_msec + } + *nstars = Nstars; + return (stars); +} Index: branches/eam_branches/20090522/Ohana/src/getstar/src/ReadImageFiles.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/getstar/src/ReadImageFiles.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/getstar/src/ReadImageFiles.c (revision 24557) @@ -136,5 +136,6 @@ !strcmp (tmpword, "PS1_DEV_0") || !strcmp (tmpword, "PS1_DEV_1") || - !strcmp (tmpword, "PS1_V1")) { + !strcmp (tmpword, "PS1_V1") || + !strcmp (tmpword, "PS1_V2")) { exttype[Nimage] = strcreate (tmpword); Index: branches/eam_branches/20090522/Ohana/src/libautocode/Makefile.Targets =================================================================== --- branches/eam_branches/20090522/Ohana/src/libautocode/Makefile.Targets (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/libautocode/Makefile.Targets (revision 24557) @@ -48,4 +48,5 @@ $(ASRC)/cmf-ps1-dev-1.$(ARCH).o \ $(ASRC)/cmf-ps1-v1.$(ARCH).o \ +$(ASRC)/cmf-ps1-v2.$(ARCH).o \ $(ASRC)/cmf-smpdata.$(ARCH).o \ $(ASRC)/getstar-ps1-dev-0.$(ARCH).o \ @@ -107,4 +108,5 @@ $(AINC)/cmf-ps1-dev-1.h \ $(AINC)/cmf-ps1-v1.h \ +$(AINC)/cmf-ps1-v2.h \ $(AINC)/cmf-smpdata.h \ $(AINC)/getstar-ps1-dev-0.h \ Index: branches/eam_branches/20090522/Ohana/src/libautocode/def/autocode.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/libautocode/def/autocode.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/libautocode/def/autocode.c (revision 24557) @@ -7,7 +7,6 @@ if (size != $SIZE) { - fprintf (stderr, "ERROR: mismatch in data types $STRUCT: %d vs %d\n", - size, $SIZE); - exit (1); + fprintf (stderr, "WARNING: mismatch in data types $STRUCT: %d vs %d\n", size, $SIZE); + return (FALSE); } @@ -26,14 +25,50 @@ } +$STRUCT *gfits_downsize_and_convert_$STRUCT ($STRUCT *data, int size, int nitems) { + + int i; + unsigned char *byte, tmp; + $STRUCT *output; + + if ($SIZE > size) { + fprintf (stderr, "ERROR: uncorrectable mismatch in data types $STRUCT: %d vs %d\n", size, $SIZE); + exit (1); + } + + // allocate a new array + ALLOCATE (output, $STRUCT, nitems); + for (i = 0; i < nitems; i++) { + memcpy (&output[i], &data[i], $SIZE); + } + + /* provide initial values to avoid compiler warnings for non-BYTE_SWAP arch */ + i = tmp = 0; + byte = NULL; + +# ifdef BYTE_SWAP + byte = (unsigned char *) output; + for (i = 0; i < nitems; i++, byte += size) { + /** BYTE SWAP **/ + } +# endif + + return (output); +} + /*** add test of EXTNAME and header-defined columns? ***/ /* return internal structure representation */ $STRUCT *gfits_table_get_$STRUCT (FTable *ftable, int *Ndata, int *swapped) { - $STRUCT *data; + $STRUCT *data, *output; *Ndata = ftable[0].header[0].Naxis[1]; data = ($STRUCT *) ftable[0].buffer; if ((swapped == NULL) || (*swapped == FALSE)) { - gfits_convert_$STRUCT (data, sizeof ($STRUCT), *Ndata); + if (!gfits_convert_$STRUCT (data, sizeof ($STRUCT), *Ndata)) { + output = gfits_downsize_and_convert_$STRUCT (data, sizeof ($STRUCT), *Ndata); + free (ftable[0].buffer); + ftable[0].buffer = output; + // XXX do I need to change NX? + } gfits_table_scale_data (ftable); if (swapped != NULL) *swapped = TRUE; Index: branches/eam_branches/20090522/Ohana/src/libautocode/def/cmf-ps1-v2.d =================================================================== --- branches/eam_branches/20090522/Ohana/src/libautocode/def/cmf-ps1-v2.d (revision 24557) +++ branches/eam_branches/20090522/Ohana/src/libautocode/def/cmf-ps1-v2.d (revision 24557) @@ -0,0 +1,78 @@ +# name of structure type +STRUCT CMF_PS1_V2 +EXTNAME CMF_PS1_V2 +TYPE BINTABLE +SIZE 136 + +# elements of data structure / FITS table +FIELD detID, IPP_IDET, unsigned int, detection ID +FIELD X, X_PSF, float, x coord, pixels +FIELD Y, Y_PSF, float, y coord, pixels +FIELD dX, X_PSF_SIG, float, x coord error, pixels +FIELD dY, Y_PSF_SIG, float, y coord error, pixels +FIELD posangle, POSANGLE, float, Posangle at source, degrees +FIELD pltscale, PLTSCALE, float, Plate Scale at source, arcsec/pixel +FIELD M, PSF_INST_MAG, float, inst mags, mags +FIELD dM, PSF_INST_MAG_SIG, float, inst mag error, mags +FIELD Map, AP_MAG_STANDARD, float, standard aperture mag, mags +FIELD apRadius, AP_MAG_RADIUS, float, radius used for fit, pixels +FIELD Mpeak, PEAK_FLUX_AS_MAG, float, peak flux as a mag, mags +FIELD Mcalib, CAL_PSF_MAG, float, calibrated psf mag, mags +FIELD dMcal, CAL_PSF_MAG_SIG, float, zero point scatter, mags +FIELD RA, RA_PSF, double, PSF RA coord, degrees +FIELD DEC, DEC_PSF, double, PSF DEC coord, degrees +FIELD sky, SKY, float, sky flux, cnts/sec +FIELD dSky, SKY_SIG, float, sky flux error, cnts/sec +FIELD psfChisq, PSF_CHISQ, float, psf fit chisq +FIELD crNsigma, CR_NSIGMA, float, Nsigma deviations from PSF to CF +FIELD extNsigma, EXT_NSIGMA, float, Nsigma deviations from PSF to EXT +FIELD fx, PSF_MAJOR, float, psf fit major axis, pixels +FIELD fy, PSF_MINOR, float, psf fit minor axis, pixels +FIELD df, PSF_THETA, float, ellipse angle, degrees +FIELD psfQual, PSF_QF, float, quality factor +FIELD psfNdof, PSF_NDOF, int, psf degrees of freedom +FIELD psfNpix, PSF_NPIX, int, psf number of pixels +FIELD Mxx, MOMENTS_XX, float, second moment X, pixels^2 +FIELD Mxy, MOMENTS_XY, float, second moment Y, pixels^2 +FIELD Myy, MOMENTS_YY, float, second moment XY, pixels^2 +FIELD flags, FLAGS, int, analysis flags +FIELD nFrames, N_FRAMES, short, images overlapping peak +FIELD padding, PADDING, short, padding for 8byte records + +# for an object in an image, we have three triplets that tell us about the shape: +# second moments: Mxx, Mxy, Myy +# model shape parameters: F_major, F_minor, F_theta +# centroid errors: sigma_X, sigma_Y, sigma_XY + +# IPP_IDET 1J IPP_IDET, unsigned i +# X_PSF 1E X_PSF, float, +# Y_PSF 1E Y_PSF, float, +# X_PSF_SIG 1E X_PSF_SIG, float, +# Y_PSF_SIG 1E Y_PSF_SIG, float, +# RA_PSF 1E RA_PSF, float, +# DEC_PSF 1E DEC_PSF, float, +# POSANGLE 1E POSANGLE, float, +# PLTSCALE 1E PLTSCALE, float, +# PSF_INST_MAG 1E PSF_INST_MAG, float, +# PSF_INST_MAG_SIG 1E PSF_INST_MAG_SIG, float, +# AP_MAG_STANDARD 1E AP_MAG_STANDARD, float, +# AP_MAG_RADIUS 1E +# PEAK_FLUX_AS_MAG 1E PEAK_FLUX_AS_MAG, float, +# CAL_PSF_MAG 1E CAL_PSF_MAG, float, +# CAL_PSF_MAG_SIG 1E CAL_PSF_MAG_SIG, float, +# SKY 1E SKY, float, +# SKY_SIGMA 1E SKY_SIG, float, +# PSF_CHISQ 1E PSF_CHISQ, float, +# CR_NSIGMA 1E CR_NSIGMA, float, +# EXT_NSIGMA 1E EXT_NSIGMA, float, +# PSF_MAJOR 1E PSF_MAJOR, float, +# PSF_MINOR 1E PSF_MINOR, float, +# PSF_THETA 1E PSF_THETA, float, +# PSF_QF 1E PSF_QF, float, +# PSF_NDOF 1J PSF_NDOF, int, +# PSF_NPIX 1J PSF_NPIX, int, +# MOMENTS_XX 1E MOMENTS_XX, float, +# MOMENTS_XY 1E MOMENTS_XY, float, +# MOMENTS_YY 1E MOMENTS_YY, float, +# N_FRAMES 1I N_FRAMES, short, +# FLAGS 1J FLAGS, int, Index: branches/eam_branches/20090522/Ohana/src/libdvo/Makefile =================================================================== --- branches/eam_branches/20090522/Ohana/src/libdvo/Makefile (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/libdvo/Makefile (revision 24557) @@ -69,5 +69,6 @@ $(SRC)/skyregion_io.$(ARCH).o \ $(SRC)/skyregion_gsc.$(ARCH).o \ -$(SRC)/skyregion_ops.$(ARCH).o +$(SRC)/skyregion_ops.$(ARCH).o \ +$(SRC)/cmf-ps1-v1-alt.$(ARCH).o # $(SRC)/dvo_convert_panstarrs.$(ARCH).o \ @@ -118,2 +119,4 @@ cleandef: cd $(AUTO) && make clean + +.PRECIOUS: $(ASRC)/%.c Index: branches/eam_branches/20090522/Ohana/src/libdvo/include/dvo.h =================================================================== --- branches/eam_branches/20090522/Ohana/src/libdvo/include/dvo.h (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/libdvo/include/dvo.h (revision 24557) @@ -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 */ @@ -184,4 +184,7 @@ SkyRegion **regions; } SkyList; + +// special-case function: +CMF_PS1_V2 *gfits_table_get_CMF_PS1_V1_Alt (FTable *ftable, int *Ndata, int *swapped); # if (0) Index: branches/eam_branches/20090522/Ohana/src/libdvo/src/cmf-ps1-v1-alt.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/libdvo/src/cmf-ps1-v1-alt.c (revision 24557) +++ branches/eam_branches/20090522/Ohana/src/libdvo/src/cmf-ps1-v1-alt.c (revision 24557) @@ -0,0 +1,126 @@ +# include "dvo.h" + +CMF_PS1_V2 *gfits_table_get_CMF_PS1_V1_Alt (FTable *ftable, int *Ndata, int *swapped) { + + int i, nitems; + unsigned char *byte, *inbyte, *otbyte, tmp; + CMF_PS1_V2 *output; + + /* provide initial values to avoid compiler warnings for non-BYTE_SWAP arch */ + i = tmp = 0; + byte = NULL; + + // this function is a special case : it must have Nx = 136 + if (ftable[0].header[0].Naxis[0] != 136) { + fprintf (stderr, "ERROR: wrong format for CMF_PS1_V1_Alt: %d vs %d\n", ftable[0].header[0].Naxis[0], 136); + exit (2); + } + + *Ndata = ftable[0].header[0].Naxis[1]; + nitems = ftable[0].header[0].Naxis[1]; + + if ((swapped == NULL) || (*swapped == FALSE)) { + +# ifdef BYTE_SWAP + // we need to do the byte swap before applying the table scaling: + byte = (unsigned char *) ftable[0].buffer; + for (i = 0; i < nitems; i++, byte += 136) { + /** BYTE SWAP **/ + SWAP_WORD (0); // IPP_IDET + SWAP_WORD (4); // X_PSF + SWAP_WORD (8); // Y_PSF + SWAP_WORD (12); // X_PSF_SIG + SWAP_WORD (16); // Y_PSF_SIG + SWAP_DBLE (20); // RA_PSF + SWAP_DBLE (28); // DEC_PSF + SWAP_WORD (36); // POSANGLE + SWAP_WORD (40); // PLTSCALE + SWAP_WORD (44); // PSF_INST_MAG + SWAP_WORD (48); // PSF_INST_MAG_SIG + SWAP_WORD (52); // AP_MAG_STANDARD + SWAP_WORD (56); // AP_MAG_RADIUS + SWAP_WORD (60); // PEAK_FLUX_AS_MAG + SWAP_WORD (64); // CAL_PSF_MAG + SWAP_WORD (68); // CAL_PSF_MAG_SIG + SWAP_WORD (72); // SKY + SWAP_WORD (76); // SKY_SIG + SWAP_WORD (80); // PSF_CHISQ + SWAP_WORD (84); // CR_NSIGMA + SWAP_WORD (88); // EXT_NSIGMA + SWAP_WORD (92); // PSF_MAJOR + SWAP_WORD (96); // PSF_MINOR + SWAP_WORD (100); // PSF_THETA + SWAP_WORD (104); // PSF_QF + SWAP_WORD (108); // PSF_NDOF + SWAP_WORD (112); // PSF_NPIX + SWAP_WORD (116); // MOMENTS_XX + SWAP_WORD (120); // MOMENTS_XY + SWAP_WORD (124); // MOMENTS_YY + SWAP_WORD (128); // FLAGS + SWAP_BYTE (132); // N_FRAMES + SWAP_BYTE (134); // PADDING + } +# endif + + gfits_table_scale_data (ftable); + if (swapped != NULL) *swapped = TRUE; + } + + byte = (unsigned char *) ftable[0].buffer; + + // allocate a new output data buffer + ALLOCATE (output, CMF_PS1_V2, nitems); + inbyte = (unsigned char *) byte; + otbyte = (unsigned char *) output; + + // the data in the input table does not line up with the output structure: copy carefully. + for (i = 0; i < nitems; i++, inbyte += 136, otbyte += 136) { + memcpy (&otbyte[0], &inbyte[0], 20); + memcpy (&otbyte[20], &inbyte[36], 36); + memcpy (&otbyte[56], &inbyte[20], 16); + memcpy (&otbyte[72], &inbyte[72], 62); + } + + free (ftable[0].buffer); + ftable[0].buffer = output; + + // XXX other mods to make ftable consistent with CMF_PS1_V2? (Nx, EXTNAME?) + + return (output); +} + +// data organization (input vs output) +// FITS struct +// WORD 0 IPP_IDET 0 IPP_IDET +// WORD 4 X_PSF 4 X_PSF +// WORD 8 Y_PSF 8 Y_PSF +// WORD 12 X_PSF_SIG 12 X_PSF_SIG +// WORD 16 Y_PSF_SIG 16 Y_PSF_SIG +// DBLE 20 RA_PSF 20 POSANGLE +// DBLE 28 DEC_PSF 24 PLTSCALE +// WORD 36 POSANGLE 28 PSF_INST_MAG +// WORD 40 PLTSCALE 32 PSF_INST_MAG_SIG +// WORD 44 PSF_INST_MAG 36 AP_MAG_STANDARD +// WORD 48 PSF_INST_MAG_SIG 40 AP_MAG_RADIUS +// WORD 52 AP_MAG_STANDARD 44 PEAK_FLUX_AS_MAG +// WORD 56 AP_MAG_RADIUS 48 CAL_PSF_MAG +// WORD 60 PEAK_FLUX_AS_MAG 52 CAL_PSF_MAG_SIG +// WORD 64 CAL_PSF_MAG 56 RA_PSF +// WORD 68 CAL_PSF_MAG_SIG 64 DEC_PSF +// WORD 72 SKY 72 SKY +// WORD 76 SKY_SIG 76 SKY_SIG +// WORD 80 PSF_CHISQ 80 PSF_CHISQ +// WORD 84 CR_NSIGMA 84 CR_NSIGMA +// WORD 88 EXT_NSIGMA 88 EXT_NSIGMA +// WORD 92 PSF_MAJOR 92 PSF_MAJOR +// WORD 96 PSF_MINOR 96 PSF_MINOR +// WORD 100 PSF_THETA 100 PSF_THETA +// WORD 104 PSF_QF 104 PSF_QF +// WORD 108 PSF_NDOF 108 PSF_NDOF +// WORD 112 PSF_NPIX 112 PSF_NPIX +// WORD 116 MOMENTS_XX 116 MOMENTS_XX +// WORD 120 MOMENTS_XY 120 MOMENTS_XY +// WORD 124 MOMENTS_YY 124 MOMENTS_YY +// WORD 128 FLAGS 128 FLAGS +// BYTE 132 N_FRAMES 132 N_FRAMES +// BYTE 134 PADDING 134 PADDING Index: branches/eam_branches/20090522/Ohana/src/relastro/Makefile =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/Makefile (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/Makefile (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/include/relastro.h =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/include/relastro.h (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/include/relastro.h (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/ConfigInit.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/ConfigInit.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/FitChip.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/FitChip.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/FitChip.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/GetAstromError.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/GetAstromError.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/ImageOps.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/ImageOps.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/ImageOps.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/UpdateObjects.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/UpdateObjects.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/args.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/args.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/args.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/bcatalog.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/bcatalog.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/bcatalog.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/initialize.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/initialize.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/initialize.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/load_catalogs.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/load_catalogs.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/relastroVisual.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/relastroVisual.c (revision 24557) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/relastroVisual.c (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/relastroVisual.h =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/relastroVisual.h (revision 24557) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/relastroVisual.h (revision 24557) @@ -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: branches/eam_branches/20090522/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- branches/eam_branches/20090522/Ohana/src/relastro/src/relastro_objects.c (revision 24238) +++ branches/eam_branches/20090522/Ohana/src/relastro/src/relastro_objects.c (revision 24557) @@ -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) {