Index: /branches/cnb_branch_20081011/psastro/src/Makefile.am =================================================================== --- /branches/cnb_branch_20081011/psastro/src/Makefile.am (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/Makefile.am (revision 20263) @@ -53,4 +53,5 @@ psastroErrorCodes.c \ psastroVersion.c \ + psastroVisual.c \ psastroDefineFiles.c \ psastroAnalysis.c \ Index: /branches/cnb_branch_20081011/psastro/src/psastro.h =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastro.h (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastro.h (revision 20263) @@ -75,4 +75,14 @@ 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 psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); +bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd); + // demo plots bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); Index: /branches/cnb_branch_20081011/psastro/src/psastroArguments.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroArguments.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroArguments.c (revision 20263) @@ -37,5 +37,5 @@ psArgumentRemove (N, &argc, argv); } - + // apply the chip correction based on the reference astrometry? if ((N = psArgumentGet (argc, argv, "-fixchips"))) { @@ -51,5 +51,5 @@ status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist"); if (status) { - // if supplied, assume -fixchips is desired + // if supplied, assume -fixchips is desired psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true); } @@ -81,4 +81,10 @@ } + // run in visual mode? + if ((N = psArgumentGet (argc, argv, "-visual"))) { + psArgumentRemove (N, &argc, argv); + psastroSetVisual (true); + } + // dump the configuration to a file? if ((N = psArgumentGet (argc, argv, "-dumpconfig"))) { Index: /branches/cnb_branch_20081011/psastro/src/psastroAstromGuess.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroAstromGuess.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroAstromGuess.c (revision 20263) @@ -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,12 +138,14 @@ } - // dump or plot the resulting projected positions - if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - } - - 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); + } + + psastroVisualPlotRawStars(rawstars, fpa, chip, recipe); + + if (psTraceGetLevel("psastro.plot") > 0) { + psastroPlotRawstars (rawstars, fpa, chip, recipe); + } } } @@ -155,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); @@ -201,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; @@ -223,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; @@ -243,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; } @@ -273,21 +275,21 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - 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); + 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); } @@ -298,14 +300,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); } @@ -313,8 +315,8 @@ map->x->coeffMask[1][1] = PS_POLY_MASK_SET; map->y->coeffMask[1][1] = PS_POLY_MASK_SET; - + psVectorFitPolynomial2D (map->x, NULL, 0, cornerPn, NULL, cornerPs, cornerQs); psVectorFitPolynomial2D (map->y, NULL, 0, cornerQn, NULL, cornerPs, cornerQs); - + // apply the linear fit... psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); @@ -325,4 +327,6 @@ psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); + psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); + psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); @@ -333,5 +337,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); @@ -341,7 +345,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 } @@ -354,11 +358,11 @@ if (0) { - 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); } @@ -381,5 +385,5 @@ psFree (map); psFree (view); - + return true; Index: /branches/cnb_branch_20081011/psastro/src/psastroCleanup.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroCleanup.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroCleanup.c (revision 20263) @@ -11,4 +11,5 @@ pmConceptsDone (); pmConfigDone (); + psastroVisualClose (); 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_20081011/psastro/src/psastroLoadRefstars.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroLoadRefstars.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroLoadRefstars.c (revision 20263) @@ -138,4 +138,6 @@ psastroDumpRefstars (refstars, "refstars.dat"); } + + psastroVisualPlotRefStars (refstars, recipe); if (psTraceGetLevel("psastro.plot") > 0) { Index: /branches/cnb_branch_20081011/psastro/src/psastroLuminosityFunction.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroLuminosityFunction.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroLuminosityFunction.c (revision 20263) @@ -129,7 +129,5 @@ lumFunc->sPeak = sPeak; -#if 0 - psastroLuminosityFunctionPlot(lnMag, Mag, lumFunc, rawFunc); -#endif + psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc); psFree (lnMag); Index: /branches/cnb_branch_20081011/psastro/src/psastroOneChipFit.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroOneChipFit.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroOneChipFit.c (revision 20263) @@ -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,74 +32,74 @@ 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: - 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; - } + // modify the order to correspond to the actual number of matched stars: + 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; + } - // 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); + } } @@ -122,9 +122,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; } @@ -133,9 +133,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 @@ -143,14 +143,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); + if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotOneChipFit (rawstars, refstars, match, recipe); + psastroPlotOneChipFit (rawstars, refstars, match, recipe); } Index: /branches/cnb_branch_20081011/psastro/src/psastroRemoveClumps.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroRemoveClumps.c (revision 20262) +++ /branches/cnb_branch_20081011/psastro/src/psastroRemoveClumps.c (revision 20263) @@ -11,11 +11,11 @@ float Ymax = obj->FP->y; for (int i = 0; i < input->n; i++) { - obj = (pmAstromObj *)input->data[i]; - if (!isfinite(obj->FP->x)) continue; - if (!isfinite(obj->FP->y)) continue; - Xmin = PS_MIN (Xmin, obj->FP->x); - Xmax = PS_MAX (Xmax, obj->FP->x); - Ymin = PS_MIN (Ymin, obj->FP->y); - Ymax = PS_MAX (Ymax, obj->FP->y); + obj = (pmAstromObj *)input->data[i]; + if (!isfinite(obj->FP->x)) continue; + if (!isfinite(obj->FP->y)) continue; + Xmin = PS_MIN (Xmin, obj->FP->x); + Xmax = PS_MAX (Xmax, obj->FP->x); + Ymin = PS_MIN (Ymin, obj->FP->y); + Ymax = PS_MAX (Ymax, obj->FP->y); } @@ -27,10 +27,10 @@ // accumulate 2D histogram in image for (int i = 0; i < input->n; i++) { - obj = (pmAstromObj *)input->data[i]; - if (!isfinite(obj->FP->x)) continue; - if (!isfinite(obj->FP->y)) continue; - int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols); - int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows); - count->data.U32[Yi][Xi] ++; + obj = (pmAstromObj *)input->data[i]; + if (!isfinite(obj->FP->x)) continue; + if (!isfinite(obj->FP->y)) continue; + int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols); + int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows); + count->data.U32[Yi][Xi] ++; } @@ -38,15 +38,15 @@ psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MAX | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); if (!psImageStats(stats, count, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); - psFree(stats); - psFree(count); - return NULL; + psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); + psFree(stats); + psFree(count); + return NULL; } if (stats->max < 1) { - psError(PS_ERR_UNKNOWN, false, "no valid sources in image\n"); - psFree(stats); - psFree(count); - return NULL; + psError(PS_ERR_UNKNOWN, false, "no valid sources in image\n"); + psFree(stats); + psFree(count); + return NULL; } @@ -55,14 +55,16 @@ psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit); + psastroVisualPlotRemoveClumps (input, count, scale, limit); + // find and exclude objects in bad pixels psArray *output = psArrayAllocEmpty (input->n); for (int i = 0; i < input->n; i++) { - obj = (pmAstromObj *)input->data[i]; - if (!isfinite(obj->FP->x)) continue; - if (!isfinite(obj->FP->y)) continue; - int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols); - int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows); - if (count->data.U32[Yi][Xi] > limit) continue; - psArrayAdd (output, 16, obj); + obj = (pmAstromObj *)input->data[i]; + if (!isfinite(obj->FP->x)) continue; + if (!isfinite(obj->FP->y)) continue; + int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols); + int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows); + if (count->data.U32[Yi][Xi] > limit) continue; + psArrayAdd (output, 16, obj); } @@ -77,11 +79,11 @@ for (int iy = 0; iy < count->numRows; iy++) { for (int ix = 0; ix < count->numCols; ix++) { - if (count->data.U32[iy][ix] > limit) { - psPlane *pixel = psPlaneAlloc(); - pixel->x = ix; - pixel->y = iy; - psArrayAdd (badpix, 16, pixel); - psFree (pixel); - } + if (count->data.U32[iy][ix] > limit) { + psPlane *pixel = psPlaneAlloc(); + pixel->x = ix; + pixel->y = iy; + psArrayAdd (badpix, 16, pixel); + psFree (pixel); + } } }