Index: branches/cnb_branch_20090113/psastro/src/Makefile.am =================================================================== --- branches/cnb_branch_20090113/psastro/src/Makefile.am (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/Makefile.am (revision 21208) @@ -53,5 +53,4 @@ psastroErrorCodes.c \ psastroVersion.c \ - psastroVisual.c \ psastroDefineFiles.c \ psastroAnalysis.c \ Index: branches/cnb_branch_20090113/psastro/src/psastro.h =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastro.h (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastro.h (revision 21208) @@ -14,4 +14,5 @@ # define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) +#if 0 // this structure represents a fit to the logN / logS curve for a set of stars // logN = offset + slope * logS @@ -25,4 +26,5 @@ int sPeak; // sum of stars to peak bin } pmLumFunc; +#endif bool psastroDataSave (pmConfig *config); @@ -77,18 +79,4 @@ psString psastroVersionLong(void); -// psastroVisual functions -bool psastroSetVisual (bool mode); -bool psastroVisualClose(); -bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc); -bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); -bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe); -bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit); -bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld); -bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); -bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd); -bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); -bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale); -bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe); - // demo plots bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); Index: branches/cnb_branch_20090113/psastro/src/psastroArguments.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroArguments.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroArguments.c (revision 21208) @@ -88,5 +88,4 @@ if ((N = psArgumentGet (argc, argv, "-visual"))) { psArgumentRemove (N, &argc, argv); - psastroSetVisual (true); pmAstromSetVisual (true); } Index: branches/cnb_branch_20090113/psastro/src/psastroAstromGuess.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroAstromGuess.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroAstromGuess.c (revision 21208) @@ -29,6 +29,6 @@ psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE); if (!recipe) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); + return false; } @@ -36,5 +36,5 @@ bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); if (!status) { - useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); + useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); } @@ -42,6 +42,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -49,7 +49,7 @@ double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); - return false; - } + psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); + return false; + } psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -67,5 +67,5 @@ bool bilevelAstrometry = false; if (!useModel) { - psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); + psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); } @@ -75,35 +75,35 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - if (!useModel) { - if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; - } + if (!useModel) { + if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; + } if (newFPA) { newFPA = false; - while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; - while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; + while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; + while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; RAminSky = fpa->toSky->R - M_PI; RAmaxSky = fpa->toSky->R + M_PI; } - // report and save the current best guess for the chip 0,0 pixel coordinates - { - psPlane ptCH, ptFP, ptTP; - psSphere ptSky; - - ptCH.x = 0; - ptCH.y = 0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); - psDeproject (&ptSky, &ptTP, fpa->toSky); - psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); - - psVectorAppend (cornerL, ptFP.x); - psVectorAppend (cornerM, ptFP.y); - psVectorAppend (cornerP, ptTP.x); - psVectorAppend (cornerQ, ptTP.y); - psVectorAppend (cornerR, ptSky.r); - psVectorAppend (cornerD, ptSky.d); - } + // report and save the current best guess for the chip 0,0 pixel coordinates + { + psPlane ptCH, ptFP, ptTP; + psSphere ptSky; + + ptCH.x = 0; + ptCH.y = 0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + + psVectorAppend (cornerL, ptFP.x); + psVectorAppend (cornerM, ptFP.y); + psVectorAppend (cornerP, ptTP.x); + psVectorAppend (cornerQ, ptTP.y); + psVectorAppend (cornerR, ptSky.r); + psVectorAppend (cornerD, ptSky.d); + } // apply the new WCS guess data to all of the data in the readouts @@ -119,5 +119,5 @@ if (rawstars == NULL) { continue; } - *nStars += rawstars->n; + *nStars += rawstars->n; for (int i = 0; i < rawstars->n; i++) { pmAstromObj *raw = rawstars->data[i]; @@ -138,14 +138,14 @@ } - // dump or plot the resulting projected positions - if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - } - - psastroVisualPlotRawStars(rawstars, fpa, chip, recipe); - - if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotRawstars (rawstars, fpa, chip, recipe); - } + // dump or plot the resulting projected positions + if (psTraceGetLevel("psastro.dump") > 0) { + psastroDumpRawstars (rawstars, fpa, chip); + } + + pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe); + + if (psTraceGetLevel("psastro.plot") > 0) { + psastroPlotRawstars (rawstars, fpa, chip, recipe); + } } } @@ -157,12 +157,12 @@ psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR", PS_META_REPLACE, "", *nStars); if (*nStars == 0) { - psLogMsg ("psastro", 2, "no sources available for astrometry\n"); - psFree (view); - return true; - } - - psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", - DEG_RAD*RAmin, DEG_RAD*DECmin, - DEG_RAD*RAmax, DEG_RAD*DECmax); + psLogMsg ("psastro", 2, "no sources available for astrometry\n"); + psFree (view); + return true; + } + + psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", + DEG_RAD*RAmin, DEG_RAD*DECmin, + DEG_RAD*RAmax, DEG_RAD*DECmax); psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", RAmin); @@ -203,13 +203,13 @@ pmHDU *hdu = pmFPAviewThisHDU (view, fpa); if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (chip, hdu->header)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - return false; - } + if (!pmAstromReadBilevelChip (chip, hdu->header)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } } else { - if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - return false; - } + if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } } return true; @@ -225,12 +225,12 @@ // load mosaic-level astrometry? if (phu) { - char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); - if (ctype) { - *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } } if (*bilevelAstrometry) { - pmAstromReadBilevelMosaic (fpa, phu->header); - } + pmAstromReadBilevelMosaic (fpa, phu->header); + } psFree (view); return true; @@ -245,6 +245,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -277,55 +277,55 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - // XXX we are currently inconsistent with marking the good vs the bad data - // psastroChipAstrom sets data_exists to false if the fit is bad. this is - // probably wrong since it implies there is no data! - - // skip chips for which the astrometry failed (NASTRO == 0) - if (!chip->cells->n) goto skip_chip; - pmCell *cell = chip->cells->data[0]; - if (!cell) goto skip_chip; - - if (!cell->readouts->n) goto skip_chip; - pmReadout *readout = cell->readouts->data[0]; - if (!readout) goto skip_chip; - - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) goto skip_chip; - - int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); - if (!nAstro) goto skip_chip; - - float astError = psMetadataLookupF32 (&status, updates, "CERROR"); - if (fabs(astError) < 1e-6) goto skip_chip; - - psPlane ptCH, ptFP, ptTP; - psSphere ptSky; - - ptCH.x = 0; - ptCH.y = 0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); - psDeproject (&ptSky, &ptTP, fpa->toSky); - psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); - - // new corner locations based on the calibrated astrometry - psVectorAppend (cornerLn, ptFP.x); - psVectorAppend (cornerMn, ptFP.y); - psVectorAppend (cornerPn, ptTP.x); - psVectorAppend (cornerQn, ptTP.y); - psVectorAppend (cornerRn, ptSky.r); - psVectorAppend (cornerDn, ptSky.d); - psVectorAppend (cornerMK, 0); - continue; + // XXX we are currently inconsistent with marking the good vs the bad data + // psastroChipAstrom sets data_exists to false if the fit is bad. this is + // probably wrong since it implies there is no data! + + // skip chips for which the astrometry failed (NASTRO == 0) + if (!chip->cells->n) goto skip_chip; + pmCell *cell = chip->cells->data[0]; + if (!cell) goto skip_chip; + + if (!cell->readouts->n) goto skip_chip; + pmReadout *readout = cell->readouts->data[0]; + if (!readout) goto skip_chip; + + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) goto skip_chip; + + int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); + if (!nAstro) goto skip_chip; + + float astError = psMetadataLookupF32 (&status, updates, "CERROR"); + if (fabs(astError) < 1e-6) goto skip_chip; + + psPlane ptCH, ptFP, ptTP; + psSphere ptSky; + + ptCH.x = 0; + ptCH.y = 0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + + // new corner locations based on the calibrated astrometry + psVectorAppend (cornerLn, ptFP.x); + psVectorAppend (cornerMn, ptFP.y); + psVectorAppend (cornerPn, ptTP.x); + psVectorAppend (cornerQn, ptTP.y); + psVectorAppend (cornerRn, ptSky.r); + psVectorAppend (cornerDn, ptSky.d); + psVectorAppend (cornerMK, 0); + continue; skip_chip: - // new corner locations based on the calibrated astrometry - psVectorAppend (cornerLn, 0.0); - psVectorAppend (cornerMn, 0.0); - psVectorAppend (cornerPn, 0.0); - psVectorAppend (cornerQn, 0.0); - psVectorAppend (cornerRn, 0.0); - psVectorAppend (cornerDn, 0.0); - psVectorAppend (cornerMK, 1); + // new corner locations based on the calibrated astrometry + psVectorAppend (cornerLn, 0.0); + psVectorAppend (cornerMn, 0.0); + psVectorAppend (cornerPn, 0.0); + psVectorAppend (cornerQn, 0.0); + psVectorAppend (cornerRn, 0.0); + psVectorAppend (cornerDn, 0.0); + psVectorAppend (cornerMK, 1); } @@ -336,14 +336,14 @@ for (int i = 0; i < cornerRo->n; i++) { - - psPlane ptTP; - psSphere ptSky; - - ptSky.r = cornerRo->data.F32[i]; - ptSky.d = cornerDo->data.F32[i]; - - psProject (&ptTP, &ptSky, fpa->toSky); - psVectorAppend (cornerPs, ptTP.x); - psVectorAppend (cornerQs, ptTP.y); + + psPlane ptTP; + psSphere ptSky; + + ptSky.r = cornerRo->data.F32[i]; + ptSky.d = cornerDo->data.F32[i]; + + psProject (&ptTP, &ptSky, fpa->toSky); + psVectorAppend (cornerPs, ptTP.x); + psVectorAppend (cornerQs, ptTP.y); } @@ -351,9 +351,9 @@ map->x->coeffMask[1][1] = PS_POLY_MASK_SET; map->y->coeffMask[1][1] = PS_POLY_MASK_SET; - + // fit the valid chips, mask the invalid chips psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs); psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs); - + // apply the linear fit... psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); @@ -364,5 +364,5 @@ psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); - psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); + pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); @@ -374,5 +374,5 @@ float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]); float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]); - + psLogMsg ("psastro", 3, "boresite offset : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]); psLogMsg ("psastro", 3, "boresite angle : %f, scale: %f", angle*PS_DEG_RAD, scale); @@ -382,7 +382,7 @@ psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER"); if (!header) { - header = psMetadataAlloc(); - psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); - psFree (header); // drop this reference + header = psMetadataAlloc(); + psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); + psFree (header); // drop this reference } @@ -395,11 +395,11 @@ if (DEBUG) { - FILE *f = fopen ("corners.dat", "w"); - for (int i = 0; i < cornerRo->n; i++) { - fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", - cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], - cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); - } - fclose (f); + FILE *f = fopen ("corners.dat", "w"); + for (int i = 0; i < cornerRo->n; i++) { + fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", + cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], + cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); + } + fclose (f); } @@ -424,5 +424,5 @@ psFree (map); psFree (view); - + return true; Index: branches/cnb_branch_20090113/psastro/src/psastroCleanup.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroCleanup.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroCleanup.c (revision 21208) @@ -4,4 +4,5 @@ psFree (config); + pmAstromVisualClose (); psTimerStop (); @@ -11,6 +12,4 @@ pmConceptsDone (); pmConfigDone (); - psastroVisualClose (); - pmAstromVisualClose (); fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro"); // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro"); Index: branches/cnb_branch_20090113/psastro/src/psastroFixChips.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroFixChips.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroFixChips.c (revision 21208) @@ -13,5 +13,5 @@ bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); if (!status) { - fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); + fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); } if (!fixChips) return true; @@ -25,6 +25,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -48,6 +48,6 @@ // files associated with the science image if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError (PS_ERR_IO, false, "Can't load the astrometry model file"); - return false; + psError (PS_ERR_IO, false, "Can't load the astrometry model file"); + return false; } @@ -64,100 +64,100 @@ if (DEBUG) { - f = fopen ("corners.raw.dat", "w"); - chipName = NULL; + f = fopen ("corners.raw.dat", "w"); + chipName = NULL; } pmChip *obsChip = NULL; while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - // XXX we are currently inconsistent with marking the good vs the bad data - // psastroChipAstrom sets data_exists to false if the fit is bad. this is - // probably wrong since it implies there is no data! - - // skip chips for which the astrometry failed (NASTRO == 0) - if (!obsChip->cells->n) continue; - pmCell *cell = obsChip->cells->data[0]; - if (!cell) continue; - - if (!cell->readouts->n) continue; - pmReadout *readout = cell->readouts->data[0]; - if (!readout) continue; - - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) continue; - - int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); - if (!nAstro) continue; - - // set the chip astrometry using the astrom file - pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); - - psRegion *region = pmChipPixels (obsChip); - psPlane ptCP, ptFP; - - ptCP.x = region->x0; ptCP.y = region->y0; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - ptCP.x = region->x0; ptCP.y = region->y1; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - ptCP.x = region->x1; ptCP.y = region->y1; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - ptCP.x = region->x1; ptCP.y = region->y0; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - psFree (region); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // XXX we are currently inconsistent with marking the good vs the bad data + // psastroChipAstrom sets data_exists to false if the fit is bad. this is + // probably wrong since it implies there is no data! + + // skip chips for which the astrometry failed (NASTRO == 0) + if (!obsChip->cells->n) continue; + pmCell *cell = obsChip->cells->data[0]; + if (!cell) continue; + + if (!cell->readouts->n) continue; + pmReadout *readout = cell->readouts->data[0]; + if (!readout) continue; + + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) continue; + + int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); + if (!nAstro) continue; + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + psRegion *region = pmChipPixels (obsChip); + psPlane ptCP, ptFP; + + ptCP.x = region->x0; ptCP.y = region->y0; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + ptCP.x = region->x0; ptCP.y = region->y1; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + ptCP.x = region->x1; ptCP.y = region->y1; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + ptCP.x = region->x1; ptCP.y = region->y0; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + psFree (region); } xObs->n = yObs->n = xRef->n = yRef->n = nPts; if (DEBUG) fclose (f); - + psPlaneTransform *map = psPlaneTransformAlloc (1, 1); - + psVector *mask = psVectorAlloc (nPts, PS_TYPE_U8); psVectorInit (mask, 0); @@ -167,25 +167,25 @@ for (int i = 0; i < 3; i++) { - psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef); - psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n); - - psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef); - psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n); + psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef); + psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n); + + psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef); + psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n); } // loop over all chips, select the outliers, and replace the measured astrometry with the model - // the measured transformation above must be applied to make the comparison, and also then applied to the + // the measured transformation above must be applied to make the comparison, and also then applied to the // model transformation if (DEBUG) { - f = fopen ("corners.fit.dat", "w"); - for (int i = 0; i < xObs->n; i++) { - psPlane obsCoord, refCoord; - refCoord.x = xRef->data.F32[i]; - refCoord.y = yRef->data.F32[i]; - psPlaneTransformApply (&obsCoord, map, &refCoord); - fprintf (f, "%f %f %f %f %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y); - } - fclose (f); + f = fopen ("corners.fit.dat", "w"); + for (int i = 0; i < xObs->n; i++) { + psPlane obsCoord, refCoord; + refCoord.x = xRef->data.F32[i]; + refCoord.y = yRef->data.F32[i]; + psPlaneTransformApply (&obsCoord, map, &refCoord); + fprintf (f, "%f %f %f %f %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y); + } + fclose (f); } @@ -199,125 +199,125 @@ while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - // set the chip astrometry using the astrom file - pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); - - // bad Astrometry test: ref pixel or angle outside nominal - - psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; - psPlane obsCoord, refCoord, tmpCoord; - - // find location of 0,0 pixel in focal plane coords for this chip - psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); - - // find location of 0,0 pixel in focal plane coords for ref chip - // apply the global field rotation and offset before comparing - psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel); - psPlaneTransformApply (&refCoord, map, &tmpCoord); - - psPlane offPixel = {100.0, 0.0, 100.0, 0.0}; - psPlane obsOffPt, refOffPt; - - // find location of 0,0 pixel in focal plane coords for this chip - psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); - - // find location of 0,0 pixel in focal plane coords for ref chip - psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel); - psPlaneTransformApply (&refOffPt, map, &tmpCoord); - - double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); - double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); - - bool badAstrom = false; - badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; - badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; - badAstrom |= fabs(obsAngle - refAngle) > angleTol; - - fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); - - // XXX for now, just use first readout - pmCell *cell = obsChip->cells->data[0]; - pmReadout *readout = cell->readouts->data[0]; - - // update the header (pull or create local view to entry on readout->analysis) - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) { - updates = psMetadataAlloc (); - psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); - psFree (updates); - } - - psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x); - psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y); - psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle); - - // for successful chips, save the measured offsets in the header - if (!badAstrom) continue; - - // XXX for now, let's just fail on the bad chips. In the future, let's try to recover, but we still need to - // catch the failures relative to the model - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); - continue; - - psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n", - view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); - - psFree (obsChip->toFPA); - psFree (obsChip->fromFPA); - - // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter - // XXX this only works if toTPA is at most a linear transformation - psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY); - for (int i = 0; i <= refChip->toFPA->x->nX; i++) { - for (int j = 0; j <= refChip->toFPA->x->nY; j++) { - double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; - double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; - toFPA->x->coeff[i][j] = f1 + f2; - - double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; - double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; - toFPA->y->coeff[i][j] = g1 + g2; - } - } - toFPA->x->coeff[0][0] += map->x->coeff[0][0]; - toFPA->y->coeff[0][0] += map->y->coeff[0][0]; - - psRegion *region = pmChipPixels (obsChip); - obsChip->toFPA = toFPA; - obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50); - psFree (region); - - // use the new position to re-try the match fit - // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); - if (rawstars == NULL) { continue; } - - // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) { continue; } - - // the absolute minimum number of stars is 4 (for order = 1) - if ((rawstars->n < 4) || (refstars->n < 4)) { - readout->data_exists = false; - psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", - rawstars->n, refstars->n); - continue; - } - - psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars); - - // XXX update the header with info to reflect the failure - if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) { - readout->data_exists = false; - psLogMsg ("psastro", 3, "failed to find a solution\n"); - continue; - } - - pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL); - } - - psastroVisualPlotFixChips (input, xObs, yObs); + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + // bad Astrometry test: ref pixel or angle outside nominal + + psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; + psPlane obsCoord, refCoord, tmpCoord; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + // apply the global field rotation and offset before comparing + psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel); + psPlaneTransformApply (&refCoord, map, &tmpCoord); + + psPlane offPixel = {100.0, 0.0, 100.0, 0.0}; + psPlane obsOffPt, refOffPt; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel); + psPlaneTransformApply (&refOffPt, map, &tmpCoord); + + double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); + double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); + + bool badAstrom = false; + badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; + badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; + badAstrom |= fabs(obsAngle - refAngle) > angleTol; + + fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); + + // XXX for now, just use first readout + pmCell *cell = obsChip->cells->data[0]; + pmReadout *readout = cell->readouts->data[0]; + + // update the header (pull or create local view to entry on readout->analysis) + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) { + updates = psMetadataAlloc (); + psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); + psFree (updates); + } + + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle); + + // for successful chips, save the measured offsets in the header + if (!badAstrom) continue; + + // XXX for now, let's just fail on the bad chips. In the future, let's try to recover, but we still need to + // catch the failures relative to the model + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); + continue; + + psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n", + view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); + + psFree (obsChip->toFPA); + psFree (obsChip->fromFPA); + + // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter + // XXX this only works if toTPA is at most a linear transformation + psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY); + for (int i = 0; i <= refChip->toFPA->x->nX; i++) { + for (int j = 0; j <= refChip->toFPA->x->nY; j++) { + double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; + double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; + toFPA->x->coeff[i][j] = f1 + f2; + + double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; + double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; + toFPA->y->coeff[i][j] = g1 + g2; + } + } + toFPA->x->coeff[0][0] += map->x->coeff[0][0]; + toFPA->y->coeff[0][0] += map->y->coeff[0][0]; + + psRegion *region = pmChipPixels (obsChip); + obsChip->toFPA = toFPA; + obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50); + psFree (region); + + // use the new position to re-try the match fit + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) { continue; } + + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + // the absolute minimum number of stars is 4 (for order = 1) + if ((rawstars->n < 4) || (refstars->n < 4)) { + readout->data_exists = false; + psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", + rawstars->n, refstars->n); + continue; + } + + psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars); + + // XXX update the header with info to reflect the failure + if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) { + readout->data_exists = false; + psLogMsg ("psastro", 3, "failed to find a solution\n"); + continue; + } + + pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL); + } + + pmAstromVisualPlotFixChips (input, xObs, yObs); psFree (xObs); psFree (yObs); Index: branches/cnb_branch_20090113/psastro/src/psastroLoadRefstars.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroLoadRefstars.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroLoadRefstars.c (revision 21208) @@ -138,5 +138,5 @@ } - psastroVisualPlotRefStars (refstars, recipe); + pmAstromVisualPlotRefStars (refstars, recipe); if (psTraceGetLevel("psastro.plot") > 0) { @@ -239,7 +239,7 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data"); - photcode = psStringCopy ("NONE"); - return photcode; + psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data"); + photcode = psStringCopy ("NONE"); + return photcode; } assert (input->fpa); @@ -247,6 +247,6 @@ *maxRho = psMetadataLookupF32(&status, recipe, "DVO.GETSTAR.MAX.RHO"); if (!status) { - psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe"); - return NULL; + psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe"); + return NULL; } @@ -266,39 +266,39 @@ if (!status) ESCAPE ("missing DVO.GETSTAR.MIN.MAG.INST"); - // PHOTCODE.DATA is a multi of metadata items + // PHOTCODE.DATA is a multi of metadata items psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false); psMetadataItem *refItem = NULL; while ((refItem = psListGetAndIncrement (iter))) { - if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder"); - - char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); - if (!status) { - psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); - continue; - } - - // does this entry match the current filter? - if (strcmp (refFilter, filter)) continue; - - psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter); - - float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT"); - if (!status) { - psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); - continue; - } - photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE"); - if (!status) { - psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); - continue; - } - - // convert the minInst to a calibrated minimum magnitude - *minMag = minInst + 2.5*log10(exptime) + zeropt; - - psFree (iter); - return photcode; - } + if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder"); + + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + + // does this entry match the current filter? + if (strcmp (refFilter, filter)) continue; + + psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter); + + float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + + // convert the minInst to a calibrated minimum magnitude + *minMag = minInst + 2.5*log10(exptime) + zeropt; + + psFree (iter); + return photcode; + } psFree (iter); @@ -306,5 +306,5 @@ photcode = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.PHOTCODE"); PS_ASSERT (photcode, NULL); - + // give up and use fixed value *minMag = psMetadataLookupF32(NULL, recipe, "DVO.GETSTAR.MIN.MAG"); Index: branches/cnb_branch_20090113/psastro/src/psastroLuminosityFunction.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroLuminosityFunction.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroLuminosityFunction.c (revision 21208) @@ -129,5 +129,5 @@ lumFunc->sPeak = sPeak; - psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc); + pmAstromVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc); psFree (lnMag); Index: branches/cnb_branch_20090113/psastro/src/psastroMosaicOneChip.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroMosaicOneChip.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroMosaicOneChip.c (revision 21208) @@ -148,5 +148,5 @@ //plot results - psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe); + pmAstromVisualPlotMosaicOneChip(rawstars, refstars, match, recipe); psFree (fitStats); Index: branches/cnb_branch_20090113/psastro/src/psastroMosaicSetMatch.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroMosaicSetMatch.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroMosaicSetMatch.c (revision 21208) @@ -57,5 +57,5 @@ psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n); - psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe); + pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe); // XXX drop the old one Index: branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c (revision 21208) @@ -11,11 +11,11 @@ // default value for match/fit : radius is in pixels - REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); + REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); // run the match/fit sequence NITER times - REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); + REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); // correct radius to FP units (physical pixel scale in microns per pixel) - REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); + REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); RADIUS *= pixelScale; @@ -32,79 +32,79 @@ for (int iter = 0; iter < nIter; iter++) { - - char name[128]; - sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); - float radius = psMetadataLookupF32 (&status, recipe, name); - radius *= pixelScale; - if (!status || (radius == 0.0)) { - radius = RADIUS; - } + char name[128]; + + sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); + float radius = psMetadataLookupF32 (&status, recipe, name); + radius *= pixelScale; + if (!status || (radius == 0.0)) { + radius = RADIUS; + } - // use small radius to match stars - match = pmAstromRadiusMatchFP (rawstars, refstars, radius); - if (match == NULL) { - psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); - return false; - } + // use small radius to match stars + match = pmAstromRadiusMatchFP (rawstars, refstars, radius); + if (match == NULL) { + psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); + return false; + } - // modify the order to correspond to the actual number of matched stars: - int Ndof_min = 3; - int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1)); - order = PS_MIN (order, order_max); + // modify the order to correspond to the actual number of matched stars: + int Ndof_min = 3; + int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1)); + order = PS_MIN (order, order_max); - // if ((match->n < 11) && (order >= 3)) order = 2; - // if ((match->n < 7) && (order >= 2)) order = 1; - // if ((match->n < 4) && (order >= 1)) order = 0; + // if ((match->n < 11) && (order >= 3)) order = 2; + // if ((match->n < 7) && (order >= 2)) order = 1; + // if ((match->n < 4) && (order >= 1)) order = 0; - if (order < 1) { - psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); - psFree (match); - return false; - } + if (order < 1) { + psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); + psFree (match); + return false; + } - // create output toFPA; set masks appropriate to the Elixir DVO astrometry format - psFree (chip->toFPA); - chip->toFPA = psPlaneTransformAlloc (order, order); - for (int i = 0; i <= chip->toFPA->x->nX; i++) { - for (int j = 0; j <= chip->toFPA->x->nY; j++) { - if (i + j > order) { - chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; - chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; - } - } - } + // create output toFPA; set masks appropriate to the Elixir DVO astrometry format + psFree (chip->toFPA); + chip->toFPA = psPlaneTransformAlloc (order, order); + for (int i = 0; i <= chip->toFPA->x->nX; i++) { + for (int j = 0; j <= chip->toFPA->x->nY; j++) { + if (i + j > order) { + chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; + chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; + } + } + } - // XXX allow statitic to be set by the user - // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); - fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); + // XXX allow statitic to be set by the user + // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); - // improved fit for astrometric terms - results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); - if (!results) { - psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); - psFree (match); - psFree (fitStats); - return false; - } - - // determine fromFPA transformation and apply new transformation to raw & ref stars - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - - // toSky converts from FPA & TPA units (microns) to sky units (radians) - float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; - // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; - float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; - int astNstar = results->yStats->clippedNvalues; - psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); + // improved fit for astrometric terms + results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + if (!results) { + psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); + psFree (match); + psFree (fitStats); + return false; + } - if (iter < nIter - 1) { - psFree (fitStats); - psFree (results); - psFree (match); - } + // determine fromFPA transformation and apply new transformation to raw & ref stars + psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + + // toSky converts from FPA & TPA units (microns) to sky units (radians) + float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; + // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; + float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; + int astNstar = results->yStats->clippedNvalues; + psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); + + if (iter < nIter - 1) { + psFree (fitStats); + psFree (results); + psFree (match); + } } @@ -127,9 +127,9 @@ if (astError > maxError) { psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); - validSolution = false; + validSolution = false; } if (astNstar < minNstar) { psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); - validSolution = false; + validSolution = false; } @@ -138,9 +138,9 @@ psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "astrometry error (arcsec)", astError); if (validSolution) { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); } else { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); } psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars @@ -148,16 +148,16 @@ // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - + // XXX check if we correctly applied the new transformation: if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - psastroDumpMatchedStars ("match.dat", rawstars, refstars, match); - psastroDumpStars (refstars, "refstars.cal.dat"); + psastroDumpRawstars (rawstars, fpa, chip); + psastroDumpMatchedStars ("match.dat", rawstars, refstars, match); + psastroDumpStars (refstars, "refstars.cal.dat"); } - psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe); + pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe); if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotOneChipFit (rawstars, refstars, match, recipe); + psastroPlotOneChipFit (rawstars, refstars, match, recipe); } Index: branches/cnb_branch_20090113/psastro/src/psastroRemoveClumps.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroRemoveClumps.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroRemoveClumps.c (revision 21208) @@ -55,5 +55,5 @@ psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit); - psastroVisualPlotRemoveClumps (input, count, scale, limit); + pmAstromVisualPlotRemoveClumps (input, count, scale, limit); // find and exclude objects in bad pixels Index: branches/cnb_branch_20090113/psastro/src/psastroUtils.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroUtils.c (revision 21120) +++ branches/cnb_branch_20090113/psastro/src/psastroUtils.c (revision 21208) @@ -42,11 +42,11 @@ if (!chip->toFPA) { continue; } - if (chip->cells->n == 0) { continue; } - pmCell *cell = chip->cells->data[0]; + if (chip->cells->n == 0) { continue; } + pmCell *cell = chip->cells->data[0]; if (!cell->process || !cell->file_exists) { continue; } - if (cell->readouts->n == 0) { continue; } - pmReadout *readout = cell->readouts->data[0]; - if (! readout->data_exists) { continue; } + if (cell->readouts->n == 0) { continue; } + pmReadout *readout = cell->readouts->data[0]; + if (! readout->data_exists) { continue; } pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]); @@ -88,5 +88,5 @@ psastroMosaicSetAstrom (fpa); if (!useExternal) { - psastroVisualPlotCommonScale (fpa, oldScale); + pmAstromVisualPlotCommonScale (fpa, oldScale); } psFree (oldScale); Index: branches/cnb_branch_20090113/psastro/src/psastroVisual.c =================================================================== --- branches/cnb_branch_20090113/psastro/src/psastroVisual.c (revision 21120) +++ (revision ) @@ -1,1303 +1,0 @@ -/***********************************/ -/***Diagnostic plots for psastro****/ -/***********************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - -# include "psastroInternal.h" -# include -# include -# include - -# if (HAVE_KAPA) -# include - -# define KAPAX 700 -# define KAPAY 700 - -//variables to determine when things are plotted -static bool isVisual = false; -static bool plotRawStars = true; -static bool plotRefStars = false; -static bool plotLumFunc = true; -static bool plotRemoveClumps = true; -static bool plotOneChipFit = true; -static bool plotFixChips = true; -static bool plotAstromGuessCheck = true; -static bool plotMosaicMatches = true; -static bool plotCommonScale = true; -static bool plotMosaicOneChip = true; - -// variables to store plotting window indices -static int kapa = -1; -static int kapa2 = -1; - -// helper routine prototypes -bool psastroVisualInitGraph (int kapa, KapaSection *section, Graphdata *graphdata); -bool psastroVisualTriplePlot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing); -bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip); -bool psastroVisualTripleOverplot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing); -bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing); -bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title); - - -/*****************************/ -/** Initialization Routines **/ -/*****************************/ - -/** start or stop plotting */ -bool psastroSetVisual (bool mode) { - isVisual = mode; - return true; -} - - -/** open, name, and resize a window if necessary*/ -bool psastroVisualInitWindow( int *kapid, char *name ) { - if (*kapid == -1) { - *kapid = KapaOpenNamedSocket("kapa", name); - if (*kapid == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - KapaResize (*kapid, KAPAX, KAPAY); - } - return true; -} - - -/** ask the user how to proceed*/ -bool psastroVisualAskUser( bool *plotflag ) { - char key[10]; - fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?"); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - if (key[0] == 's') { - *plotflag = false; - } - if (key[0] == 'a') { - psastroSetVisual(false); - pmAstromSetVisual(false); - } - return true; -} - - -/** destroy windows at the end of a run*/ -bool psastroVisualClose() { - if(kapa != -1) - KiiClose(kapa); - if (kapa2 != -1) - KiiClose(kapa2); - return true; -} - - -/***********************/ -/** Plotting Routines **/ -/***********************/ - - -/** - * Plot raw stars as determined from first pass astrometry fit - * Called within psastroAstromGeuss - */ -bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) -{ - // make sure we want to plot this - if (!plotRawStars || !isVisual) return true; - - //set up plot region - if (!psastroVisualInitWindow (&kapa, "psastro:plots")) - return false; - Graphdata graphdata; - KapaSection section; - - KapaInitGraph (&graphdata); - KapaClearPlots (kapa); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 7; - graphdata.size = 0.5; - graphdata.style = 2; - - section.dx = 0.4; - section.dy = 0.4; - - //initialize and populate plotting vectors - bool status = false; - float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN"); - float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX"); - - psVector *xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); - psVector *yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); - psVector *zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); - - section.x = 0.0; - section.y = 0.5; - section.name = NULL; - psStringAppend (§ion.name, "a0"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - //Chip coordinates - int n = 0; - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - - xVec->data.F32[n] = raw->chip->x; - yVec->data.F32[n] = raw->chip->y; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel(kapa, "Chip", KAPA_LABEL_XP); - KapaSendLabel(kapa, "X", KAPA_LABEL_XM); - KapaSendLabel(kapa, "Y", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - //Focal Plane Coordinates - section.x = 0.5; - section.y = 0.5; - section.name = NULL; - psStringAppend (§ion.name, "a1"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - n = 0; - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - - xVec->data.F32[n] = raw->FP->x; - yVec->data.F32[n] = raw->FP->y; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "Focal Plane", KAPA_LABEL_XP); - KapaSendLabel (kapa, "L", KAPA_LABEL_XM); - KapaSendLabel (kapa, "M", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - // Tangent Plane Coordinates - section.x = 0.0; - section.y = 0.0; - section.name = NULL; - psStringAppend (§ion.name, "a2"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - n = 0; - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - - xVec->data.F32[n] = raw->TP->x; - yVec->data.F32[n] = raw->TP->y; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "Tangential Plane", KAPA_LABEL_XP); - KapaSendLabel (kapa, "P", KAPA_LABEL_XM); - KapaSendLabel (kapa, "Q", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - //sky coordinates - section.x = 0.5; - section.y = 0.0; - section.name = NULL; - psStringAppend (§ion.name, "a3"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - n = 0; - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - - xVec->data.F32[n] = DEG_RAD*raw->sky->r; - yVec->data.F32[n] = DEG_RAD*raw->sky->d; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "Sky", KAPA_LABEL_XP); - KapaSendLabel(kapa, "RA", KAPA_LABEL_XM); - KapaSendLabel(kapa, "Dec", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - // flip x (East increase to left) - SWAP (graphdata.xmin, graphdata.xmax); - KapaSetLimits (kapa, &graphdata); - - // plot label - section.x = 0.0; - section.y = 0.0; - section.dx = .95; - section.dy = .95; - section.name = NULL; - psStringAppend (§ion.name, "a5"); - KapaSetSection (kapa, §ion); - KapaSendLabel (kapa, "Raw Star Selection - Initial Astrometry", KAPA_LABEL_XP); - psFree (section.name); - - // pause and wait for user input: - psastroVisualAskUser( &plotRawStars ); - - psFree (xVec); - psFree (yVec); - psFree (zVec); - return true; -} - - -/** - * plot the location of references stars over the entire fpa - * invoked during psastroChooseRefStars - */ -bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe) -{ - //make sure we want to plot this - if (!isVisual || !plotRefStars) return true; - - //set up plotting variables - if (!psastroVisualInitWindow (&kapa, "psastro:plots")) - return false; - - Graphdata graphdata; - KapaInitGraph (&graphdata); - KapaClearSections (kapa); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 7; - graphdata.size = 0.5; - graphdata.style = 2; - - //initialize and populate plot vectors - bool status = false; - float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN"); - float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX"); - - psVector *xVec = psVectorAlloc (refstars->n, PS_TYPE_F32); - psVector *yVec = psVectorAlloc (refstars->n, PS_TYPE_F32); - psVector *zVec = psVectorAlloc (refstars->n, PS_TYPE_F32); - - int n = 0; - for (int i = 0; i < refstars->n; i++) { - pmAstromObj *ref = refstars->data[i]; - if (!isfinite(ref->Mag)) continue; - if (ref->Mag > rMagMax) continue; - if (ref->Mag < rMagMin) continue; - - xVec->data.F32[n] = DEG_RAD*ref->sky->r; - yVec->data.F32[n] = DEG_RAD*ref->sky->d; - zVec->data.F32[n] = ref->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "RA", KAPA_LABEL_XM); - KapaSendLabel (kapa, "Dec", KAPA_LABEL_YM); - KapaSendLabel (kapa, "Reference Stars", KAPA_LABEL_XP); - psastroVisualTriplePlot(kapa, &graphdata, xVec, yVec, zVec, false); - - // flip x (East increase to left) - SWAP (graphdata.xmin, graphdata.xmax); - KapaSetLimits (kapa, &graphdata); - - // pause and wait for user input: - psastroVisualAskUser(&plotRefStars); - - psFree (xVec); - psFree (yVec); - psFree (zVec); - return true; -} - - -/** - * Plot the two luminosity functions created within psastroRefStarSubset - * The luminosity functions are used to select a subset of reference stars, - * so we plot the cutoff that defines this subset - */ -bool psastroVisualPlotLuminosityFunction (psVector *lnMag, ///< Log(n) for each magnitude bin - psVector *Mag, ///< magnitude bins - pmLumFunc *lumFunc,///< Fit to the reference star luminosity function - pmLumFunc *rawFunc ///< Fit to the raw star luminoisty function - ) -{ - - // make sure we want to plot this - if ( !isVisual || !plotLumFunc ) return true; - - //set up plotting variables - if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) - return false; - - Graphdata graphdata; - KapaSection section1 = {"s1", .1, .1, .85, .35}; - KapaSection section2 = {"s2", .1, .5, .85, .35}; - KapaSection section; - if(rawFunc == NULL) - section = section1; - else - section = section2; - if (rawFunc == NULL) - KapaClearPlots(kapa); - KapaInitGraph(&graphdata); - - //Determine Plot Limits - psastroVisualScaleGraphdata(&graphdata, Mag, lnMag, false); - - //Make a line for the fit - float x[2] = {graphdata.xmin, graphdata.xmax}; - float y[2] = {lumFunc->offset + x[0] * lumFunc->slope, - lumFunc->offset + x[1] * lumFunc->slope}; - - //Plot Data - KapaSetSection(kapa, §ion); - KapaSetLimits(kapa, &graphdata); - - KapaSetFont (kapa, "helvetica", 14); - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "Magnitude", KAPA_LABEL_XM); - KapaSendLabel (kapa, "Log(N)", KAPA_LABEL_YM); - if (rawFunc == NULL) - KapaSendLabel (kapa, "Raw Star Luminosity Function", KAPA_LABEL_XP); - else - KapaSendLabel (kapa, - "Reference Star Luminosity Function, Shifted Raw Fit, and Cutoff", - KAPA_LABEL_XP); - graphdata.color=KapaColorByName("black"); - graphdata.style = 1; - KapaPrepPlot (kapa, lnMag->n, &graphdata); - KapaPlotVector(kapa, lnMag->n, Mag->data.F32, "x"); - KapaPlotVector(kapa, lnMag->n, lnMag->data.F32, "y"); - - //Overplot fit - graphdata.style=0; - KapaPrepPlot(kapa,2,&graphdata); - KapaPlotVector(kapa, 2, x, "x"); - KapaPlotVector(kapa, 2, y, "y"); - - //If rawFunc was supplied, overplot the raw star fit + cutoff - if( rawFunc != NULL) { - double mRef = 0.5*(lumFunc->mMin + lumFunc->mMax); - double logRho = mRef * lumFunc->slope + lumFunc->offset; - double mRaw = (logRho - rawFunc->offset) / rawFunc->slope; - double deltaM = mRef - mRaw; - double mRefMax = rawFunc->mMax + deltaM; - - float xraw[2] = {rawFunc->mMin + deltaM, rawFunc->mMax + deltaM}; - float yraw[2] = {rawFunc->offset + (rawFunc->slope) * rawFunc->mMin, - rawFunc->offset + (rawFunc->slope) * rawFunc->mMax}; - float x[2] = {mRefMax, mRefMax}; - float y[2] = {graphdata.ymin, graphdata.ymax}; - graphdata.color= KapaColorByName("red"); - KapaPrepPlot(kapa, 2, &graphdata); - KapaPlotVector(kapa, 2, x, "x"); - KapaPlotVector(kapa, 2, y, "y"); - KapaPrepPlot (kapa, 2, &graphdata); - KapaPlotVector (kapa, 2, xraw, "x"); - KapaPlotVector (kapa, 2, yraw, "y"); - - // pause and wait for user input: - psastroVisualAskUser (&plotLumFunc); - } - return true; -} // end of psastroVisualPlotLuminosityFunction - - -/** - * Plot the stars in a region, and indicate which stars are part of 'clumps' - * These stars are flagged during astrometric fitting, since dense regions are - * harder to cross-match than sparse ones. Called during psastroRemoveClumps. - */ -bool psastroVisualPlotRemoveClumps (psArray *input, ///< Array containing the field stars - psImage *count, ///< A 2D histogram of the field star distribution - int scale, ///< The pixel size of the histogram - float limit ///< The minimum numuber of stars in a bin flagged as a clump - ) -{ - - //make sure we want to plot this - if (!isVisual || !plotRemoveClumps) return true; - - //set up plot variables - if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) return false; - - KapaSection section; - Graphdata graphdata; - section.x = 0; - section.dx = .9; - section.y = 0; - section.dy = .9; - section.name = NULL; - psStringAppend( §ion.name, "a0"); - KapaInitGraph(&graphdata); - KapaSetSection(kapa, §ion); - psFree(section.name); - - graphdata.ptype = 7; - graphdata.size = 0.5; - graphdata.style = 2; - graphdata.color = KapaColorByName ("black"); - KapaClearPlots(kapa); - - //set up plot vectors - float Xmin = +FLT_MAX; - float Xmax = -FLT_MAX; - float Ymin = +FLT_MAX; - float Ymax = -FLT_MAX; - psVector *xVec = psVectorAlloc (input->n, PS_TYPE_F32); - psVector *yVec = psVectorAlloc (input->n, PS_TYPE_F32); - - //determine boundaries for histogram bin calculation - int n = 0; - for (int i=0; i< input->n; i++) { - pmAstromObj *obj = (pmAstromObj *)input->data[i]; - if (!isfinite(obj->FP->x)) continue; - if (!isfinite(obj->FP->y)) continue; - xVec->data.F32[n] = obj->FP->x; - yVec->data.F32[n] = obj->FP->y; - Xmin = PS_MIN (Xmin, xVec->data.F32[n]); - Xmax = PS_MAX (Xmax, xVec->data.F32[n]); - Ymin = PS_MIN (Ymin, yVec->data.F32[n]); - Ymax = PS_MAX (Ymax, yVec->data.F32[n]); - n++; - } - xVec->n = yVec->n = n; - - //plot stars - graphdata.xmax = Xmax; - graphdata.xmin = Xmin; - graphdata.ymax = Ymax; - graphdata.ymin = Ymin; - KapaSetLimits (kapa, &graphdata); - KapaSetFont (kapa, "helvetica", 14); - - KapaBox (kapa, &graphdata); - - KapaSendLabel (kapa, "L (pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa, "M (pixels)", KAPA_LABEL_YM); - KapaSendLabel (kapa, "Regions Flagged as Clumps (Red Boxes)", - KAPA_LABEL_XP); - - KapaPrepPlot (kapa, xVec->n, &graphdata); - KapaPlotVector (kapa, xVec->n, xVec->data.F32, "x"); - KapaPlotVector (kapa, yVec->n, yVec->data.F32, "y"); - - graphdata.color = KapaColorByName ("red"); - graphdata.style = 1; - - //overplot clumpy regions excluded from analysis - for(int i = 0; i < count->numCols; i++) { - for (int j = 0; j < count->numRows; j++) { - if(count->data.U32[j][i] <= limit) continue; //not a clump - float Xbot = (i - 5) * scale + Xmin; - float Ybot = (j - 5) * scale + Ymin; - if(Xbot < graphdata.xmin || Xbot > graphdata.xmax || - Ybot < graphdata.ymin || Ybot > graphdata.ymax) continue; - float x[5] = {Xbot, Xbot + scale, Xbot + scale, Xbot, Xbot}; - float y[5] = {Ybot, Ybot, Ybot + scale, Ybot + scale, Ybot}; - KapaPrepPlot (kapa, 5, &graphdata); - KapaPlotVector (kapa, 5, x, "x"); - KapaPlotVector (kapa, 5, y, "y"); - } - } - - //ask for user input and finish - psastroVisualAskUser (&plotRemoveClumps); - psFree (xVec); - psFree (yVec); - - return true; -} - - -/** - * Assess the goodness of fit for a signle chip by - * plotting the fit residuals - * invoked during psastroOneChipFit - */ -bool psastroVisualPlotOneChipFit (psArray *rawstars, ///< stars detected in the image - psArray *refstars, ///< reference stars over the same region - psArray *match, ///< contains which rawstars match to which refstars - psMetadata *recipe ///< data reduction recipe - ) -{ - - //make sure we want to plot this - if (!isVisual || !plotOneChipFit) - return true; - - //plot the residuals - if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)")) - return false; - - //ask for user input and finish - psastroVisualAskUser(&plotOneChipFit); - return true; -} - - -/** - * Plots the chip corners in the FP before and after chips with inconsistent solutions have been fixed. - * Invoked during psastroFixChips - */ -bool psastroVisualPlotFixChips (pmFPAfile *input, ///< focal plane array file - psVector *xOld, ///< old X location of chip cornerss - psVector *yOld ///< old Y location of chip corners - ) -{ - //make sure we want to plot this - if(!isVisual || !plotFixChips) return true; - - if(!psastroVisualInitWindow(&kapa, "psastro:plots")) return false; - - KapaSection section = {"s1", .05, .05, .9, .9}; - Graphdata graphdata; - KapaInitGraph( &graphdata); - KapaClearPlots (kapa); - graphdata.ptype = 2; - graphdata.style = 2; - - psVector *xNew = psVectorAllocEmpty (xOld->n, PS_TYPE_F32); - psVector *yNew = psVectorAllocEmpty (yOld->n, PS_TYPE_F32); - - // copy of the code in psastroFixChips that generated xOld, yOld, but for xNew, yNew - pmFPAview *view = pmFPAviewAlloc (0); - - pmChip *obsChip = NULL; - while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - psRegion *region = pmChipPixels (obsChip); - psPlane ptCP, ptFP; - - ptCP.x = region->x0; ptCP.y = region->y0; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - psVectorAppend (xNew, ptFP.x); - psVectorAppend (yNew, ptFP.y); - - ptCP.x = region->x0; ptCP.y = region->y1; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - psVectorAppend (xNew, ptFP.x); - psVectorAppend (yNew, ptFP.y); - - ptCP.x = region->x1; ptCP.y = region->y1; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - psVectorAppend (xNew, ptFP.x); - psVectorAppend (yNew, ptFP.y); - - ptCP.x = region->x1; ptCP.y = region->y0; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - psVectorAppend (xNew, ptFP.x); - psVectorAppend (yNew, ptFP.y); - - psFree (region); - } - - //set up graph - psastroVisualScaleGraphdata(&graphdata, xOld, yOld, true); - psastroVisualInitGraph(kapa, §ion, &graphdata); - KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM); - KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM); - KapaSendLabel (kapa, "Chip corners before (black) and after (red) FixChips", KAPA_LABEL_XP); - KapaPrepPlot (kapa, xOld->n, &graphdata); - KapaPlotVector (kapa, xOld->n, xOld->data.F32, "x"); - KapaPlotVector (kapa, xOld->n, yOld->data.F32, "y"); - - graphdata.ptype = 1; - graphdata.color = KapaColorByName("red"); - KapaPrepPlot (kapa, xNew->n, &graphdata); - KapaPlotVector (kapa, xNew->n, xNew->data.F32, "x"); - KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y"); - - psastroVisualAskUser(&plotFixChips); - psFree(xNew); - psFree(yNew); - psFree(view); - return true; -} - - -/** - * Plots the fpa chip corners projected on to the tangential plane before and after - * the astrometry solution has been applied. In psastroAstromGuessCheck, the old corners - * are then fit to the new corners to get a sense at how far off the initial WCS info was - * in offset, rotation, and scale. This procedure also plots the residuals of the fit from - * old to new coordinates - */ -bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, ///< P coordinates of chip corners before fitting - psVector *cornerQo, ///< Q coordinates of chip corners before fitting - psVector *cornerPn, ///< P coordinates of chip corners after fitting - psVector *cornerQn, ///< Q coordinates of chip corners after fitting - psVector *cornerPd, ///< P coordinate residuals of fit from old to new coordinates - psVector *cornerQd ///< Q coordinate residuals of fit from old to new coordinates - ) -{ - - //make sure we want to plot this - if (!isVisual || !plotAstromGuessCheck) return true; - - //set up graph window - if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) - return false; - Graphdata graphdata; - KapaSection section; - KapaInitGraph(&graphdata); - KapaClearPlots (kapa); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 2; - graphdata.size = 0.5; - graphdata.style = 2; - - section.dx = 0.4; - section.dy = 0.4; - - //Old Corners - section.x = 0.30; - section.y = 0.50; - section.name = NULL; - psStringAppend (§ion.name, "a0"); - KapaSetSection (kapa, §ion); - psFree(section.name); - - psastroVisualScaleGraphdata (&graphdata, cornerPo, cornerPo, true); - KapaSetLimits (kapa, &graphdata); - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM); - KapaSendLabel (kapa, - "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry", - KAPA_LABEL_XP); - KapaPrepPlot (kapa, cornerPo->n, &graphdata); - KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x"); - KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y"); - - // New Corners - graphdata.color = KapaColorByName("red"); - graphdata.ptype = 7; - graphdata.size = 1.5; - KapaPrepPlot (kapa, cornerPn->n, &graphdata); - KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x"); - KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y"); - - // Residuals - psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32); - psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32); - for(int i=0; i < cornerPn->n; i++) { - xResid->data.F32[i] = (cornerPd->data.F32[i]); - yResid->data.F32[i] = (cornerQd->data.F32[i]); - } - - graphdata.color = KapaColorByName("black"); - graphdata.size=0.5; - section.x = 0.3; - section.y = 0.0; - section.name = NULL; - psStringAppend (§ion.name, "a1"); - KapaSetSection (kapa, §ion); - psFree(section.name); - - psastroVisualScaleGraphdata (&graphdata, xResid, yResid, true); - KapaSetLimits (kapa, &graphdata); - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "dP", KAPA_LABEL_XM); - KapaSendLabel (kapa, "dQ", KAPA_LABEL_YM); - KapaSendLabel (kapa, - "Residual of the Fit from the Initial Astrometry to the Final Astrometry", - KAPA_LABEL_XP); - KapaPrepPlot (kapa, cornerPd->n, &graphdata); - KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x"); - KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y"); - - psFree(xResid); - psFree(yResid); - psastroVisualAskUser (&plotAstromGuessCheck); - return true; -} - - -/** - * Plots the pixel scales of the fpa before they are - * equalized in psastroMosaicCommonScale - */ -bool psastroVisualPlotCommonScale (pmFPA *fpa, ///< the fpa - psVector *oldScale ///< the old pixel scale of each chip in the fpa - ) -{ - //make sure we want to plot this - if (!isVisual || !plotCommonScale) return false; - - if (!psastroVisualInitWindow(&kapa, "psastro:plots")) return false; - - KapaSection section = {"s1", .05, .05, .9, .9}; - Graphdata graphdata; - psPlane ptCH, ptFP; - ptCH.x = 0; - ptCH.y = 0; - psVector *xVec = psVectorAlloc (oldScale->n, PS_TYPE_F32); - psVector *yVec = psVectorAlloc (oldScale->n, PS_TYPE_F32); - - int nobj = 0; - - //project each chip corner to the Focal Plane - for(int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip->process || !chip->file_exists) { continue; } - if (!chip->toFPA) { continue; } - - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - xVec->data.F32[nobj] = ptFP.x; - yVec->data.F32[nobj] = ptFP.y; - nobj++; - } - - //set up plot window - KapaInitGraph (&graphdata); - KapaClearPlots (kapa); - KapaSetSection (kapa, §ion); - KapaSetFont (kapa, "helvetica", 14); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, oldScale, false); - KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM); - KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM); - KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP); - - psastroVisualAskUser (&plotCommonScale); - return true; -} - - -/** - * plot the residuals between raw stars and ref stars after - * fitting in psastroMosaicOneChip - */ -bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, - psArray *match, psMetadata *recipe) -{ - - //make sure we want to plot this - if (!isVisual || !plotMosaicOneChip) return false; - - //plot the residuals - if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode")) - return false; - - //ask for user input and finish - psastroVisualAskUser(&plotMosaicOneChip); - - return true; -} - -/** psastroVisualPlotMosaicMatches - * Plot the matches between raw and reference stars during psastroVisualMosaicSetMatch - */ -bool psastroVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars, - psArray *match, int iteration, - psMetadata *recipe) -{ - //make sure we want to plot this - if (!isVisual || !plotMosaicMatches) return true; - - char title[60]; - sprintf(title, "Matches found during psastroMosaicSetMatch iteration %d", iteration); - - if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, title)) - return false; - - //ask for user input - psastroVisualAskUser (&plotMosaicMatches); - return true; -} - - -/*********************/ -/** Helper Routines **/ -/*********************/ - - -/** psastroVisualInitGraph (kapa, *section, *graphdata) - * Initializes graph, sets the section, sets the font, - * sets the limits, draws the box - */ -bool psastroVisualInitGraph (int kapa, KapaSection *section, Graphdata *graphdata) -{ - KapaSetSection (kapa, section); - KapaSetFont (kapa, "helvetica", 14); - KapaSetLimits (kapa, graphdata); - KapaBox (kapa, graphdata); - return true; -} - - -/** psastroVisualTriplePlot - * plot 2 vectors whose point sizes are scaled by a third vector - * autoscales the plot - */ -bool psastroVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) -{ - psastroVisualScaleGraphdata (graphdata, xVec, yVec, true); - //printf("%f %f %f %f \n",graphdata->xmin, graphdata->xmax, graphdata->ymin, graphdata->ymax); - // set the scale vector - psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); - psastroVisualCreateScaleVec (zVec, zScale, increasing); - - KapaSetFont (kapid, "helvetica", 14); - KapaSetLimits(kapid, graphdata); - KapaBox (kapid, graphdata); - - // the point size will be scaled from the z vector - graphdata->ptype = 7; - graphdata->style = 2; - graphdata->size = -1; - KapaPrepPlot (kapid, xVec->n, graphdata); - KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x"); - KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y"); - KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z"); - psFree (zScale); - return true; -} - - -/** psastroVisualTripleOverplot - * same as above, but does not rescale the plot or redraw the bounding box - */ -bool psastroVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) { - - // set the scale vector - psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); - psastroVisualCreateScaleVec (zVec, zScale, increasing); - - KapaSetFont (kapid, "helvetica", 14); - - // the point size will be scaled from the z vector - graphdata->size = -1; - KapaPrepPlot (kapid, xVec->n, graphdata); - KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x"); - KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y"); - KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z"); - psFree (zScale); - return true; -} - - -/** psastroVisualCreateScaleVec - * create a vector of plot point sizes from a vector - */ -bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing) { - // set limits based on data values - float zmin = +FLT_MAX; - float zmax = -FLT_MAX; - - for (int i = 0; i < zVec->n; i++) { - zmin = PS_MIN (zmin, zVec->data.F32[i]); - zmax = PS_MAX (zmax, zVec->data.F32[i]); - } - - float range = zmax - zmin; - if (range == 0.0) { - psVectorInit (zScale, 1.0); - } else { - for (int i = 0; i < zVec->n; i++) { - if (increasing) { - zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zVec->data.F32[i] - zmin)/range)); - } else { - zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zmax - zVec->data.F32[i])/range)); - } - } - } - return true; -} - - -/** psastroVisualScaleGraphdata - * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to - * prevent outliers from making te plot region too big. - */ -bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip) { - - graphdata->xmin = +FLT_MAX; - graphdata->xmax = -FLT_MAX; - graphdata->ymin = +FLT_MAX; - graphdata->ymax = -FLT_MAX; - - //determine standard deviation of xVec and yVec - psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psVectorStats (statsX, xVec, NULL, NULL, 0); - psVectorStats (statsY, yVec, NULL, NULL, 0); - - float xhi = +FLT_MAX, xlo = -FLT_MAX, yhi = +FLT_MAX, ylo = -FLT_MAX; - if (clip) { - xhi = statsX->sampleMedian + 3 *statsX->sampleStdev; - xlo = statsX->sampleMedian - 3 *statsX->sampleStdev; - yhi = statsY->sampleMedian + 3 *statsY->sampleStdev; - ylo = statsY->sampleMedian - 3 *statsY->sampleStdev; - } - - // abort if there is no good data - if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) { - graphdata->xmin = -1; - graphdata->ymin = -1; - graphdata->xmax = 1; - graphdata->ymax = 1; - psFree(statsX); - psFree(statsY); - return false; - } - - for(int i = 0; i < xVec->n; i++) { - if (!isfinite(xVec->data.F32[i])) continue; - if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue; - graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]); - graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]); - } - - for (int i = 0; i < yVec->n; i++) { - if (!isfinite(xVec->data.F32[i])) continue; - if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue; - graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]); - graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]); - } - - // add a whitespace border - float range = graphdata->xmax - graphdata->xmin; - if (range == 0) range = 1; - graphdata->xmin -= .05 * range; - graphdata->xmax += .05 * range; - - range = graphdata->ymax - graphdata->ymin; - if (range == 0) range = 1; - graphdata->ymin -= .05 * range; - graphdata->ymax += .05 * range; - - psFree (statsX); - psFree (statsY); - return true; -} - - -/** psastroVisualResdiPlot - * Plots the residuals between matched raw and reference stars. Used to assess the quality of an - * astrometry solution - */ -bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title) { - - //set up the first window - if (!psastroVisualInitWindow( &kapa, "psastro:plots")) return false; - - //initialize graph information - Graphdata graphdata; - KapaSection section; - - KapaInitGraph (&graphdata); - KapaClearPlots (kapa); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 7; - graphdata.size = 0.5; - graphdata.style = 2; - - section.dx = 0.4; - section.dy = 0.4; - - //initialize and populate the plotting vectors - bool status = false; - float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN"); - float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX"); - float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN"); - float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX"); - - psVector *xVec = psVectorAlloc (match->n, PS_TYPE_F32); - psVector *yVec = psVectorAlloc (match->n, PS_TYPE_F32); - psVector *zVec = psVectorAlloc (match->n, PS_TYPE_F32); - - // X vs dX - section.x = 0.0; - section.y = 0.5; - section.name = NULL; - psStringAppend (§ion.name, "a0"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - int n = 0; - for (int i = 0; i < match->n; i++) { - pmAstromMatch *pair = match->data[i]; - pmAstromObj *raw = rawstars->data[pair->raw]; - pmAstromObj *ref = refstars->data[pair->ref]; - - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - if (ref->Mag < rMagMin) continue; - if (ref->Mag > rMagMax) continue; - - xVec->data.F32[n] = raw->chip->x; - yVec->data.F32[n] = raw->chip->x - ref->chip->x; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "X", KAPA_LABEL_XM); - KapaSendLabel (kapa, "dX", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - // X vs dY - section.x = 0.5; - section.y = 0.5; - section.name = NULL; - psStringAppend (§ion.name, "a1"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - n = 0; - for (int i = 0; i < match->n; i++) { - pmAstromMatch *pair = match->data[i]; - pmAstromObj *raw = rawstars->data[pair->raw]; - pmAstromObj *ref = refstars->data[pair->ref]; - - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - if (ref->Mag < rMagMin) continue; - if (ref->Mag > rMagMax) continue; - - xVec->data.F32[n] = raw->chip->x; - yVec->data.F32[n] = raw->chip->y - ref->chip->y; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "X", KAPA_LABEL_XM); - KapaSendLabel (kapa, "dY", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - // Y vs dX - section.x = 0.0; - section.y = 0.0; - section.name = NULL; - psStringAppend (§ion.name, "a2"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - n = 0; - for (int i = 0; i < match->n; i++) { - pmAstromMatch *pair = match->data[i]; - pmAstromObj *raw = rawstars->data[pair->raw]; - pmAstromObj *ref = refstars->data[pair->ref]; - - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - if (ref->Mag < rMagMin) continue; - if (ref->Mag > rMagMax) continue; - - xVec->data.F32[n] = raw->chip->y; - yVec->data.F32[n] = raw->chip->x - ref->chip->x; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "Y", KAPA_LABEL_XM); - KapaSendLabel (kapa, "dX", KAPA_LABEL_YM); - psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); - - // Y vs dY - section.x = 0.5; - section.y = 0.0; - section.name = NULL; - psStringAppend (§ion.name, "a3"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - n = 0; - for (int i = 0; i < match->n; i++) { - pmAstromMatch *pair = match->data[i]; - pmAstromObj *raw = rawstars->data[pair->raw]; - pmAstromObj *ref = refstars->data[pair->ref]; - - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - if (ref->Mag < rMagMin) continue; - if (ref->Mag > rMagMax) continue; - - xVec->data.F32[n] = raw->chip->y; - yVec->data.F32[n] = raw->chip->y - ref->chip->y; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa, "Y", KAPA_LABEL_XM); - KapaSendLabel (kapa, "dY", KAPA_LABEL_YM); - pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false); - - section.x = 0.0; - section.y = 0.0; - section.dx = 0.95; - section.dy = 0.95; - section.name = NULL; - psStringAppend (§ion.name, "a5"); - KapaSetSection (kapa, §ion); - KapaSendLabel (kapa, title, KAPA_LABEL_XP); - psFree (section.name); - - - // X vs Y plot (different window) - if (!psastroVisualInitWindow( &kapa2, "psastro:plots")) - return false; - - KapaInitGraph (&graphdata); - KapaClearPlots (kapa2); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 2; - graphdata.style = 2; - - psFree (xVec); - psFree (yVec); - psFree (zVec); - - xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); - yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); - zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); - - // X vs Y by mag (raw) - n = 0; - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - if (!isfinite(raw->Mag)) continue; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - - xVec->data.F32[n] = raw->chip->x; - yVec->data.F32[n] = raw->chip->y; - zVec->data.F32[n] = raw->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - - KapaSendLabel (kapa2, "X", KAPA_LABEL_XM); - KapaSendLabel (kapa2, "Y", KAPA_LABEL_YM); - KapaSendLabel (kapa2, - "Chip Coordinates. Black = Raw Stars. Red = Ref Stars. Blue = Matched Stars" - , KAPA_LABEL_XP); - psastroVisualTriplePlot (kapa2, &graphdata, xVec, yVec, zVec, false); - - // X vs Y by mag (ref) - psFree (xVec); - psFree (yVec); - psFree (zVec); - - xVec = psVectorAlloc (refstars->n, PS_TYPE_F32); - yVec = psVectorAlloc (refstars->n, PS_TYPE_F32); - zVec = psVectorAlloc (refstars->n, PS_TYPE_F32); - - graphdata.color = KapaColorByName ("red"); - graphdata.ptype = 7; - graphdata.style = 2; - - n = 0; - for (int i = 0; i < refstars->n; i++) { - pmAstromObj *ref = refstars->data[i]; - if (!isfinite(ref->Mag)) continue; - if (ref->Mag < rMagMin) continue; - if (ref->Mag > rMagMax) continue; - - xVec->data.F32[n] = ref->chip->x; - yVec->data.F32[n] = ref->chip->y; - zVec->data.F32[n] = ref->Mag; - n++; - } - xVec->n = yVec->n = zVec->n = n; - psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false); - - //rescale the graph to include all points - float xmin = graphdata.xmin; - float ymin = graphdata.ymin; - float xmax = graphdata.xmax; - float ymax = graphdata.ymax; - psastroVisualScaleGraphdata(&graphdata, xVec, yVec, true); - graphdata.xmin = PS_MIN(xmin, graphdata.xmin); - graphdata.ymin = PS_MIN(ymin, graphdata.ymin); - graphdata.xmax = PS_MAX(xmax, graphdata.xmax); - graphdata.ymax = PS_MAX(ymax, graphdata.ymax); - KapaSetLimits (kapa2, &graphdata); - - //overplot matched stars in blue - psFree (xVec); - psFree (yVec); - psFree (zVec); - - xVec = psVectorAlloc (match->n, PS_TYPE_F32); - yVec = psVectorAlloc (match->n, PS_TYPE_F32); - zVec = psVectorAlloc (match->n, PS_TYPE_F32); - - graphdata.color = KapaColorByName ("blue"); - n = 0; - for (int i = 0; i < match->n; i++) { - pmAstromMatch *pair = match->data[i]; - pmAstromObj *raw = rawstars->data[pair->raw]; - pmAstromObj *ref = refstars->data[pair->ref]; - if (raw->Mag < iMagMin) continue; - if (raw->Mag > iMagMax) continue; - if (ref->Mag < rMagMin) continue; - if (ref->Mag > rMagMax) continue; - - xVec->data.F32[n] = raw->chip->x; - yVec->data.F32[n] = raw->chip->y; - zVec->data.F32[n] = iMagMin; - n++; - } - xVec->n = yVec->n = zVec->n = n; - psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false); - - psFree (xVec); - psFree (yVec); - psFree (zVec); - return true; -} - -// END OF PROGRAM -#else - -bool psastroSetVisual (bool mode) {return true}; -bool psastroVisualClose() {return true}; -bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) {return true}; -bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) {return true}; -bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe) {return true}; -bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {return true}; -bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld) {return true}; -bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true}; -bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd) {return true}; -bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true}; -bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale) {return true}; -bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe) {return true}; - -# endif -