Index: /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtract.c =================================================================== --- /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtract.c (revision 23334) +++ /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtract.c (revision 23335) @@ -41,22 +41,4 @@ } - // load the reference stars overlapping the data stars - psArray *refs = psastroLoadRefstars(config); - if (!refs) { - psError (PSASTRO_ERR_UNKNOWN, false, "failed to load reference data\n"); - return false; - } - if (refs->n == 0) { - psError(PSASTRO_ERR_REFSTARS, true, "no reference stars found"); - psFree(refs); - return false; - } - - if (!psastroChooseRefstars (config, refs)) { - psError (PSASTRO_ERR_UNKNOWN, false, "failed to select reference data for chips\n"); - psFree(refs); - return false; - } - // run the full astrometry analysis (chip and/or mosaic) if (!psastroExtractAnalysis (config)) { @@ -65,19 +47,13 @@ } - // run the full astrometry analysis (chip and/or mosaic) - if (!psastroExtractAdjust (config)) { - psErrorStackPrint(stderr, "failure in psastro model adjust\n"); - exit (1); - } - - // save the model - if (!psastroExtractDataSave (config)) { - psErrorStackPrint(stderr, "error saving output data\n"); - exit (1); - } + // XXX save the results; is this needed? + // if (!psastroExtractDataSave (config)) { + // psErrorStackPrint(stderr, "error saving output data\n"); + // exit (1); + // } psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete")); - // psastroCleanup (config); + psastroCleanup (config); exit (EXIT_SUCCESS); } Index: /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c =================================================================== --- /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c (revision 23334) +++ /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c (revision 23335) @@ -20,5 +20,6 @@ # define TEST_OUTPUT 1 -static pmFits *outStars = NULL; + +psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY); /** @@ -32,4 +33,5 @@ pmReadout *readout = NULL; float zeropt, exptime; + char extname[81]; // select the current recipe @@ -40,5 +42,27 @@ } + // force this recipe value to be false (would cause an error in psastroChooseRefstars -- no rawstars!) + psMetadataAddBool (recipe, PS_LIST_TAIL, "PSASTRO.MATCH.LUMFUNC", PS_META_REPLACE, "do not match luminosity functions", false); + psLogMsg ("psastro", PS_LOG_INFO, "extracting bright-star subrasters"); + + // load the reference stars overlapping the data stars + psArray *refs = psastroLoadRefstars(config, "PSASTRO.EXTRACT.ASTROM"); + if (!refs) { + psError (PSASTRO_ERR_UNKNOWN, false, "failed to load reference data\n"); + return false; + } + if (refs->n == 0) { + psError(PSASTRO_ERR_REFSTARS, true, "no reference stars found"); + psFree(refs); + return false; + } + + if (!psastroChooseRefstars (config, refs, "PSASTRO.EXTRACT.ASTROM")) { + psError (PSASTRO_ERR_UNKNOWN, false, "failed to select reference data for chips\n"); + psFree(refs); + return false; + } + psFree(refs); double EXTRACT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "EXTRACT_MAX_MAG"); @@ -46,4 +70,12 @@ // we have two input pmFPAfiles: PSASTRO.EXTRACT.INPUT and PSASTRO.EXTRACT.ASTROM. both are in chip-mosaic or fpa format // we have already loaded the headers from PSASTRO.EXTRACT.ASTROM and determined the astrometry terms + + // select the input astrometry data (also carries the refstars) + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.ASTROM"); + if (!astrom) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + pmFPA *fpa = astrom->fpa; // select the input data sources @@ -53,12 +85,4 @@ return false; } - pmFPA *fpa = input->fpa; - - // select the input data sources - pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.ASTROM"); - if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; - } // only load data from INPUT @@ -77,18 +101,19 @@ EXTRACT_MAX_MAG += MagOffset; + int nExt = 0; + // open the output file handle: we are just saving a series of extensions to this file + char *filename = psMetadataLookupStr (&status, config->arguments, "OUTPUT"); + if (!filename) { + psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplieda"); + return false; + } + psFits *outStars = psFitsOpen (filename, "w"); if (!outStars) { - char *filename = psMetadataLookupStr (&stats, config->arguments, "OUTPUT"); - if (!filename) { - psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplieda"); - return false; - } - outStars = psFitsOpen (filename, "w"); - if (!outStars) { - psError(PS_ERR_IO, false, "error opening file %s\n", filename); - return false; - } - // call: psFitsWriteImage (fits, NULL, subraster, 0, "extname.NN"); - } + psError(PS_ERR_IO, false, "error opening file %s\n", filename); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); // open/load files as needed @@ -109,5 +134,5 @@ char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME"); psStringAppend (&filename, "bright.%s.dat", chipname); - FILE *f = fopen (filename, "w"); + f = fopen (filename, "w"); if (!f) { psWarning ("cannot create test output file %s\n", filename); @@ -132,5 +157,5 @@ for (int i = 0; i < refstars->n; i++) { pmAstromObj *ref = refstars->data[i]; - if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue; + if (ref->Mag > EXTRACT_MAX_MAG) continue; if (TEST_OUTPUT) { @@ -138,21 +163,49 @@ } - // XXX for now, assume cell binning is 1x1 relative to chip - psImage *subraster = psastroExtractStar (input->image, ref->chip->x, ref->chip->y, DX, DY); + pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa); + psImage *subraster = psastroExtractStar (inReadout->image, ref->chip->x, ref->chip->y, 100.0, 100.0); + if (!subraster) continue; + + snprintf (extname, 80, "extname.%05d", nExt); + psFitsWriteImage (outStars, NULL, subraster, 0, extname); + nExt ++; + + psFree (subraster); } - if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE; + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; } - if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE; + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; } if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + fclose (f); } if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; - // need to close the output fits file on some basis.. - fclose (f); - + psFitsClose (outStars); psFree (view); return true; } +psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY) { + + // skip if no pixels are on the image + // skip if center is not on the image + // if bounds fall off image, paste into a full-size image. + + if (x < 0) return NULL; + if (y < 0) return NULL; + if (x >= input->numCols) return NULL; + if (y >= input->numRows) return NULL; + + psRegion fullRegion = psRegionSet (x - dX, x + dX, y - dY, y + dY); + psRegion realRegion = psRegionForImage (input, fullRegion); + psImage *subset = psImageSubset (input, realRegion); + + // if ((subset->numCols < dX) || (subset->numRows < dY)) { + // psImage *fullsize = psImageAlloc (dX, dY, PS_TYPE_F32); + // // fill in the fullsize image with the values from the subset + // // I must already have done this elsewhere... + // } + return subset; +} Index: /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractArguments.c =================================================================== --- /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractArguments.c (revision 23334) +++ /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractArguments.c (revision 23335) @@ -37,4 +37,5 @@ // psMetadata *options = pmConfigRecipeOptions (config, PSASTRO_RECIPE); + // define the image pixel data if (!pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list")) { psError(PSASTRO_ERR_ARGUMENTS, true, "Missing -file (input) or -list (input)"); @@ -43,6 +44,10 @@ } - // XXX require this? - pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", "-astromlist"); + // define the astrometry data (XXX make this optional and use the image header?) + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", "-astromlist")) { + psError(PSASTRO_ERR_ARGUMENTS, true, "Missing -astrom (smf) or -astromlist (smf)"); + psErrorStackPrint(stderr, "%s", usage); + exit (1); + } // define the output filename Index: /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractDataLoad.c =================================================================== --- /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractDataLoad.c (revision 23334) +++ /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractDataLoad.c (revision 23335) @@ -29,4 +29,13 @@ pmChip *chip; + bool newFPA = true; + double RAmin = +FLT_MAX; + double RAmax = -FLT_MAX; + double DECmin = +FLT_MAX; + double DECmax = -FLT_MAX; + + double RAminSky = NAN; + double RAmaxSky = NAN; + psTimerStart ("psastro"); @@ -48,4 +57,5 @@ pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.ASTROM"); if (!astrom) psAbort ("PSASTRO.EXTRACT.ASTROM not listed in config->files"); + pmFPA *fpa = astrom->fpa; pmFPAfileActivate(config->files, false, NULL); @@ -76,7 +86,16 @@ } + // apply the header astrometry to the astrometry structures while ((chip = pmFPAviewNextChip (view, astrom->fpa, 1)) != NULL) { psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists || !chip->data_exists) { 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; + RAminSky = fpa->toSky->R - M_PI; + RAmaxSky = fpa->toSky->R + M_PI; + } // read WCS data from the corresponding header @@ -96,6 +115,49 @@ } } + + // determine RA,DEC of 4 corners, use to find RA_MIN,MAX, DEC_MIN,MAX + psRegion *region = pmChipPixels (chip); + psPlane ptCH[4], ptFP, ptTP; + psSphere ptSky; + ptCH[0].x = region->x0; + ptCH[0].y = region->y0; + ptCH[1].x = region->x1; + ptCH[1].y = region->y0; + ptCH[2].x = region->x1; + ptCH[2].y = region->y1; + ptCH[3].x = region->x0; + ptCH[3].y = region->y1; + psFree (region); + + for (int i = 0; i < 4; i++) { + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH[i]); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + + // rationalize ra to sky range centered on boresite + while (ptSky.r < RAminSky) ptSky.r += 2.0*M_PI; + while (ptSky.r > RAmaxSky) ptSky.r -= 2.0*M_PI; + + RAmin = PS_MIN (ptSky.r, RAmin); + RAmax = PS_MAX (ptSky.r, RAmax); + + DECmin = PS_MIN (ptSky.d, DECmin); + DECmax = PS_MAX (ptSky.d, DECmax); + + psLogMsg ("psastro", 2, "chip %d (corner %d) : %f %f -> %f %f -> %f %f -> %f %f\n", + view->chip, i, ptCH[i].x, ptCH[i].y, ptFP.x, ptFP.y, ptTP.x, ptTP.y, + DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + } } psLogMsg ("psastro", 3, "convert wcs terms to internal format : %f sec\n", psTimerMark ("psastro")); + + 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); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MAX", PS_META_REPLACE, "", RAmax); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", DECmin); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", DECmax); psFree (view); Index: /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractParseCamera.c =================================================================== --- /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractParseCamera.c (revision 23334) +++ /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractParseCamera.c (revision 23335) @@ -16,5 +16,4 @@ bool status = false; - pmFPAfile *input = NULL; // the input image(s) are required arguments; they define the camera @@ -29,12 +28,7 @@ status = false; pmFPAfile *astrom = pmFPAfileDefineFromArgs(&status, config, "PSASTRO.EXTRACT.ASTROM", "ASTROM"); - if (!status) { + if (!astrom || !status) { psError (PS_ERR_UNKNOWN, false, "failed to load astrometry definition"); return NULL; - } - if (astrom) { - psLogMsg ("psastro", 3, "using supplied astrometry\n"); - } else { - psLogMsg ("psastro", 3, "using header astrometry\n"); } @@ -54,8 +48,13 @@ psArray *chips = psStringSplitArray (chipLine, ",", false); if (chips->n > 0) { - pmFPASelectChip (output->fpa, -1, true); // deselect all chips + pmFPASelectChip (input->fpa, -1, true); // deselect all chips + pmFPASelectChip (astrom->fpa, -1, true); // deselect all chips for (int i = 0; i < chips->n; i++) { int chipNum = atoi(chips->data[i]); - if (! pmFPASelectChip(output->fpa, chipNum, false)) { + if (! pmFPASelectChip(input->fpa, chipNum, false)) { + psError(PSASTRO_ERR_CONFIG, true, "Chip number %d doesn't exist in camera.\n", chipNum); + return false; + } + if (! pmFPASelectChip(astrom->fpa, chipNum, false)) { psError(PSASTRO_ERR_CONFIG, true, "Chip number %d doesn't exist in camera.\n", chipNum); return false;