Index: trunk/psastro/src/Makefile.am =================================================================== --- trunk/psastro/src/Makefile.am (revision 23590) +++ trunk/psastro/src/Makefile.am (revision 23591) @@ -47,4 +47,8 @@ psastroExtractDataLoad.c \ psastroExtractAnalysis.c \ + psastroExtractStar.c \ + psastroExtractStars.c \ + psastroExtractGhosts.c \ + psastroExtractFindChip.c \ psastroCleanup.c Index: trunk/psastro/src/psastro.h =================================================================== --- trunk/psastro/src/psastro.h (revision 23590) +++ trunk/psastro/src/psastro.h (revision 23591) @@ -127,14 +127,19 @@ bool psastroMetadataStats (pmConfig *config); -bool psastroZeroPoint (pmConfig *config); -bool psastroZeroPointReadout(pmReadout *readout, float zeropt, float exptime); -bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe); +bool psastroZeroPoint (pmConfig *config); +bool psastroZeroPointReadout(pmReadout *readout, float zeropt, float exptime); +bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe); // psastroExtract functions -bool psastroExtractAnalysis (pmConfig *config); -psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY); -bool psastroExtractParseCamera (pmConfig *config); -bool psastroExtractDataLoad (pmConfig *config); -pmConfig *psastroExtractArguments (int argc, char **argv); +bool psastroExtractAnalysis (pmConfig *config); +bool psastroExtractStars (pmConfig *config); +bool psastroExtractGhosts (pmConfig *config); +pmChip *psastroExtractFindChip (float *xChip, float *yChip, pmFPA *fpa, float xFPA, float yFPA); +bool psastroExtractChipBounds (pmFPA *fpa); + +psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY); +bool psastroExtractParseCamera (pmConfig *config); +bool psastroExtractDataLoad (pmConfig *config); +pmConfig *psastroExtractArguments (int argc, char **argv); ///@} Index: trunk/psastro/src/psastroExtractAnalysis.c =================================================================== --- trunk/psastro/src/psastroExtractAnalysis.c (revision 23590) +++ trunk/psastro/src/psastroExtractAnalysis.c (revision 23591) @@ -6,5 +6,5 @@ * * @author IfA - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @version $Revision: 1.7 $ * @date $Date: 2009-02-07 02:03:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii @@ -12,14 +12,4 @@ # include "psastroInternal.h" - -# define ESCAPE(MSG) { \ - psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate: %s", MSG); \ - psFree (view); \ - return false; \ - } - -# define TEST_OUTPUT 0 - -psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY); /** @@ -29,9 +19,4 @@ bool status; - pmChip *chip = NULL; - pmCell *cell = NULL; - pmReadout *readout = NULL; - float zeropt, exptime; - char extname[81]; // select the current recipe @@ -45,5 +30,5 @@ 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"); + psLogMsg ("psastro", PS_LOG_INFO, "loading reference stars"); // load the reference stars overlapping the data stars @@ -66,185 +51,6 @@ psFree(refs); - double EXTRACT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "EXTRACT_MAX_MAG"); - - // 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 - pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.INPUT"); - if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; - } - - // only load data from INPUT - pmFPAfileActivate(config->files, false, NULL); - pmFPAfileActivate(config->files, true, "PSASTRO.EXTRACT.INPUT"); - - // really error-out here? or just skip? - if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) { - psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); - return false; - } - - // recipe values are given in instrumental magnitudes - // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) - float MagOffset = zeropt + 2.5*log10(exptime); - 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) { - psError(PS_ERR_IO, false, "error opening file %s\n", filename); - return false; - } - - pmFPAview *view = pmFPAviewAlloc (0); - - // generate a (very) basic header: - psMetadata *phu = psMetadataAlloc (); - // select the filter; default to fixed photcode and mag limit otherwise - char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); - if (!status) ESCAPE("missing FPA.FILTER in concepts"); - - psMetadataAddStr (phu, PS_LIST_TAIL, "FILTER", 0, "filter", filter); - psMetadataAddF32 (phu, PS_LIST_TAIL, "EXPTIME", 0, "exptime", exptime); - psMetadataAddF32 (phu, PS_LIST_TAIL, "ZEROPT", 0, "zeropt", zeropt); - psFitsWriteBlank (outStars, phu, NULL); - psFree (phu); - - // open/load files as needed - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE("failed on FPA BEFORE"); - - // this loop selects the matched stars for all chips - while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); - if (!chip->process || !chip->file_exists) { continue; } - if (!chip->fromFPA) { continue; } - - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE("failed on Chip BEFORE"); - - // text output file for testing - FILE *f = NULL; - if (TEST_OUTPUT) { - char *filename = NULL; - char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME"); - psStringAppend (&filename, "bright.%s.dat", chipname); - f = fopen (filename, "w"); - if (!f) { - psWarning ("cannot create test output file %s\n", filename); - continue; - } - psFree (filename); - } - - while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); - if (!cell->process || !cell->file_exists) { continue; } - - // process each of the readouts - while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { - if (! readout->data_exists) { continue; } - - // select the raw objects for this readout (loaded in psastroExtract.c) - psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) { continue; } - - // identify the bright stars of interest - for (int i = 0; i < refstars->n; i++) { - pmAstromObj *ref = refstars->data[i]; - if (ref->Mag > EXTRACT_MAX_MAG) continue; - - if (TEST_OUTPUT) { - fprintf (f, "stars %f %f %f %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->Mag); - } - - pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa); - psImage *subraster = psastroExtractStar (inReadout->image, ref->chip->x, ref->chip->y, 200.0, 200.0); - if (!subraster) continue; - - // generate a (very) basic header: - psMetadata *header = psMetadataAlloc (); - psMetadataAddF32 (header, PS_LIST_TAIL, "CHIP_X", 0, "chip coordinate", ref->chip->x); - psMetadataAddF32 (header, PS_LIST_TAIL, "CHIP_Y", 0, "chip coordinate", ref->chip->y); - psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_X", 0, "focal-plane coordinate", ref->FP->x); - psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_Y", 0, "focal-plane coordinate", ref->FP->y); - psMetadataAddF32 (header, PS_LIST_TAIL, "MAG", 0, "magnitude", ref->Mag); - - snprintf (extname, 80, "extname.%05d", nExt); - psFitsWriteImage (outStars, header, subraster, 0, extname); - nExt ++; - - psFree (header); - psFree (subraster); - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Readout AFTER");; - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Cell AFTER");; - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Chip AFTER");; - - if (f) fclose (f); - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on FPA AFTER");; - - psFitsClose (outStars); - psFree (view); + psastroExtractStars (config); + psastroExtractGhosts (config); 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; - - int xo = x; // index of center pixel (eg, 6.00 to 6.99 -> 6) - int yo = y; - - // build an output image of size 2dX + 1, 2dY + 1 - // psRegion fullRegion = psRegionSet (xo - dX, x + dX + 1, y - dY, y + dY + 1); - // psRegion realRegion = psRegionForImage (input, fullRegion); - // psImage *subset = psImageSubset (input, realRegion); - - psImage *subset = psImageAlloc (2*dX+1, 2*dY+1, PS_TYPE_F32); - psImageInit (subset, 0.0); - - // fill in the subset image with the values from the subset - // I must already have done this elsewhere... - - for (int iy = yo - dY; iy < yo + dY + 1; iy++) { - for (int ix = xo - dX; ix < xo + dX + 1; ix++) { - - if (ix < input->col0) continue; - if (iy < input->row0) continue; - if (ix >= input->numCols) continue; - if (iy >= input->numRows) continue; - - int jx = ix + dX - xo; // runs from 0 to 2dX - int jy = iy + dY - yo; - subset->data.F32[jy][jx] = input->data.F32[iy][ix]; - } - } - return subset; -} Index: trunk/psastro/src/psastroExtractArguments.c =================================================================== --- trunk/psastro/src/psastroExtractArguments.c (revision 23590) +++ trunk/psastro/src/psastroExtractArguments.c (revision 23591) @@ -13,5 +13,5 @@ # include "psastroStandAlone.h" -static char *usage = "USAGE: psastroExtract [-outroot root] -astrom (cmffiles)"; +static char *usage = "USAGE: psastroExtract -file (input) -astrom (cmffiles) -outroot (outroot) [-chip]"; pmConfig *psastroExtractArguments (int argc, char **argv) { Index: trunk/psastro/src/psastroExtractFindChip.c =================================================================== --- trunk/psastro/src/psastroExtractFindChip.c (revision 23591) +++ trunk/psastro/src/psastroExtractFindChip.c (revision 23591) @@ -0,0 +1,106 @@ +/** @file psastroExtractFindChip.c + * + * @brief calculate chip position for the given FPA coords + * + * @ingroup psastroExtract + * + * @author IfA + * @version $Revision: 1.7 $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +static psVector *chipXmin = NULL; +static psVector *chipXmax = NULL; +static psVector *chipYmin = NULL; +static psVector *chipYmax = NULL; + +bool psastroExtractChipBounds (pmFPA *fpa) { + + chipXmin = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + chipXmax = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + chipYmin = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + chipYmax = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + + // this loop selects the matched stars for all chips + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + // determine RA,DEC of 4 corners, use to find RA_MIN,MAX, DEC_MIN,MAX + psRegion *region = pmChipPixels (chip); + psPlane ptCH[4], ptFP; + + 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); + + double Xmin = +FLT_MAX; + double Xmax = -FLT_MAX; + double Ymin = +FLT_MAX; + double Ymax = -FLT_MAX; + + for (int j = 0; j < 4; j++) { + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH[j]); + Xmin = PS_MIN (ptFP.x, Xmin); + Xmax = PS_MAX (ptFP.x, Xmax); + Ymin = PS_MIN (ptFP.y, Ymin); + Ymax = PS_MAX (ptFP.y, Ymax); + } + + // fpa-range for the given chip + chipXmin->data.F32[i] = Xmin; + chipXmax->data.F32[i] = Xmax; + chipYmin->data.F32[i] = Ymin; + chipYmax->data.F32[i] = Ymax; + } + + return true; +} + +pmChip *psastroExtractFindChip (float *xChip, float *yChip, pmFPA *fpa, float xFPA, float yFPA) { + + *xChip = NAN; + *yChip = NAN; + + if (!chipXmin) { + psastroExtractChipBounds (fpa); + } + + for (int i = 0; i < fpa->chips->n; i++) { + + if (xFPA < chipXmin->data.F32[i]) continue; + if (xFPA >= chipXmax->data.F32[i]) continue; + if (yFPA < chipYmin->data.F32[i]) continue; + if (yFPA >= chipYmax->data.F32[i]) continue; + + pmChip *chip = fpa->chips->data[i]; + psRegion *region = pmChipPixels (chip); + + psPlane ptCH, ptFP; + ptFP.x = xFPA; + ptFP.y = yFPA; + psPlaneTransformApply (&ptCH, chip->fromFPA, &ptFP); + + if (ptCH.x < region->x0) continue; + if (ptCH.x >= region->x1) continue; + if (ptCH.y < region->y0) continue; + if (ptCH.y >= region->y1) continue; + + *xChip = ptCH.x; + *yChip = ptCH.y; + return chip; + } + + return NULL; +} Index: trunk/psastro/src/psastroExtractGhosts.c =================================================================== --- trunk/psastro/src/psastroExtractGhosts.c (revision 23591) +++ trunk/psastro/src/psastroExtractGhosts.c (revision 23591) @@ -0,0 +1,183 @@ +/** @file psastroExtractGhosts.c + * + * @brief calculate ghost FPA and Chip positions for the stars loaded on the FPA + * + * @ingroup psastroExtract + * + * @author IfA + * @version $Revision: 1.7 $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate: %s", MSG); \ + psFree (view); \ + return false; \ + } + +/** + * calculate ghost FPA and Chip positions for the stars loaded on the FPA + */ +bool psastroExtractGhosts (pmConfig *config) { + + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + float zeropt, exptime; + char extname[81]; + + psLogMsg ("psastro", PS_LOG_INFO, "extracting bright-star subrasters"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + double EXTRACT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "EXTRACT_MAX_MAG"); + + // 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 + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.INPUT"); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + + // only load data from INPUT + pmFPAfileActivate(config->files, false, NULL); + pmFPAfileActivate(config->files, true, "PSASTRO.EXTRACT.INPUT"); + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + return false; + } + + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + float MagOffset = zeropt + 2.5*log10(exptime); + EXTRACT_MAX_MAG += MagOffset; + + int nExt = 0; + + // open the output file handle: we are just saving a series of extensions to this file + char *fileroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT"); + if (!fileroot) { + psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplied"); + return false; + } + char *filename = NULL; + psStringAppend (&filename, "%s.gh.fits", fileroot); + + psFits *outStars = psFitsOpen (filename, "w"); + if (!outStars) { + psError(PS_ERR_IO, false, "error opening file %s\n", filename); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + + // generate a (very) basic header: + psMetadata *phu = psMetadataAlloc (); + // select the filter; default to fixed photcode and mag limit otherwise + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + if (!status) ESCAPE("missing FPA.FILTER in concepts"); + + psMetadataAddStr (phu, PS_LIST_TAIL, "FILTER", 0, "filter", filter); + psMetadataAddF32 (phu, PS_LIST_TAIL, "EXPTIME", 0, "exptime", exptime); + psMetadataAddF32 (phu, PS_LIST_TAIL, "ZEROPT", 0, "zeropt", zeropt); + psFitsWriteBlank (outStars, phu, NULL); + psFree (phu); + + // open/load files as needed + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE("failed on FPA BEFORE"); + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE("failed on Chip BEFORE"); + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout (loaded in psastroExtract.c) + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + // identify the bright stars of interest + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + if (ref->Mag > EXTRACT_MAX_MAG) continue; + + // XXX can make a more clever model... + float xFPA = -1*ref->FP->x; + float yFPA = -1*ref->FP->y; + + float xChip, yChip; + pmChip *ghostChip = psastroExtractFindChip (&xChip, &yChip, fpa, xFPA, yFPA); + if (!ghostChip) continue; + if (!ghostChip->cells) continue; + if (!ghostChip->cells->n) continue; + pmCell *ghostCell = ghostChip->cells->data[0]; + if (!ghostCell) continue; + if (!ghostCell->readouts) continue; + if (!ghostCell->readouts->n) continue; + pmReadout *ghostReadout = ghostCell->readouts->data[0]; + if (!ghostReadout) continue; + + // XXX box size needs to be larger... + psImage *subraster = psastroExtractStar (ghostReadout->image, xChip, yChip, 200.0, 200.0); + if (!subraster) continue; + + // generate a (very) basic header: + psMetadata *header = psMetadataAlloc (); + psMetadataAddF32 (header, PS_LIST_TAIL, "CHIP_X", 0, "chip coordinate", ref->chip->x); + psMetadataAddF32 (header, PS_LIST_TAIL, "CHIP_Y", 0, "chip coordinate", ref->chip->y); + psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_X", 0, "focal-plane coordinate", ref->FP->x); + psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_Y", 0, "focal-plane coordinate", ref->FP->y); + psMetadataAddF32 (header, PS_LIST_TAIL, "MAG", 0, "magnitude", ref->Mag); + + snprintf (extname, 80, "extname.%05d", nExt); + psFitsWriteImage (outStars, header, subraster, 0, extname); + nExt ++; + + psFree (header); + psFree (subraster); + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Readout AFTER");; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Cell AFTER");; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Chip AFTER");; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on FPA AFTER");; + + psFitsClose (outStars); + psFree (filename); + psFree (view); + return true; +} Index: trunk/psastro/src/psastroExtractStar.c =================================================================== --- trunk/psastro/src/psastroExtractStar.c (revision 23591) +++ trunk/psastro/src/psastroExtractStar.c (revision 23591) @@ -0,0 +1,54 @@ +/** @file psastroExtractStar.c + * + * @brief + * + * @ingroup psastroExtract + * + * @author IfA + * @version $Revision: 1.7 $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +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; + + int xo = x; // index of center pixel (eg, 6.00 to 6.99 -> 6) + int yo = y; + + // build an output image of size 2dX + 1, 2dY + 1 + // psRegion fullRegion = psRegionSet (xo - dX, x + dX + 1, y - dY, y + dY + 1); + // psRegion realRegion = psRegionForImage (input, fullRegion); + // psImage *subset = psImageSubset (input, realRegion); + + psImage *subset = psImageAlloc (2*dX+1, 2*dY+1, PS_TYPE_F32); + psImageInit (subset, 0.0); + + // fill in the subset image with the values from the subset + // I must already have done this elsewhere... + + for (int iy = yo - dY; iy < yo + dY + 1; iy++) { + for (int ix = xo - dX; ix < xo + dX + 1; ix++) { + + if (ix < input->col0) continue; + if (iy < input->row0) continue; + if (ix >= input->numCols) continue; + if (iy >= input->numRows) continue; + + int jx = ix + dX - xo; // runs from 0 to 2dX + int jy = iy + dY - yo; + subset->data.F32[jy][jx] = input->data.F32[iy][ix]; + } + } + return subset; +} Index: trunk/psastro/src/psastroExtractStars.c =================================================================== --- trunk/psastro/src/psastroExtractStars.c (revision 23591) +++ trunk/psastro/src/psastroExtractStars.c (revision 23591) @@ -0,0 +1,166 @@ +/** @file psastroExtractAnalysis.c + * + * @brief + * + * @ingroup psastroExtract + * + * @author IfA + * @version $Revision: 1.7 $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate: %s", MSG); \ + psFree (view); \ + return false; \ + } + +/** + * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip + */ +bool psastroExtractStars (pmConfig *config) { + + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + float zeropt, exptime; + char extname[81]; + + psLogMsg ("psastro", PS_LOG_INFO, "extracting bright-star subrasters"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + double EXTRACT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "EXTRACT_MAX_MAG"); + + // 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 + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.INPUT"); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + + // only load data from INPUT + pmFPAfileActivate(config->files, false, NULL); + pmFPAfileActivate(config->files, true, "PSASTRO.EXTRACT.INPUT"); + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + return false; + } + + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + float MagOffset = zeropt + 2.5*log10(exptime); + EXTRACT_MAX_MAG += MagOffset; + + int nExt = 0; + + // open the output file handle: we are just saving a series of extensions to this file + char *fileroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT"); + if (!fileroot) { + psLogMsg ("psastro", PS_LOG_INFO, "output fileroot not supplied"); + return false; + } + char *filename = NULL; + psStringAppend (&filename, "%s.st.fits", fileroot); + + psFits *outStars = psFitsOpen (filename, "w"); + if (!outStars) { + psError(PS_ERR_IO, false, "error opening file %s\n", filename); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + + // generate a (very) basic header: + psMetadata *phu = psMetadataAlloc (); + // select the filter; default to fixed photcode and mag limit otherwise + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + if (!status) ESCAPE("missing FPA.FILTER in concepts"); + + psMetadataAddStr (phu, PS_LIST_TAIL, "FILTER", 0, "filter", filter); + psMetadataAddF32 (phu, PS_LIST_TAIL, "EXPTIME", 0, "exptime", exptime); + psMetadataAddF32 (phu, PS_LIST_TAIL, "ZEROPT", 0, "zeropt", zeropt); + psFitsWriteBlank (outStars, phu, NULL); + psFree (phu); + + // open/load files as needed + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE("failed on FPA BEFORE"); + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE("failed on Chip BEFORE"); + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout (loaded in psastroExtract.c) + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + // identify the bright stars of interest + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + if (ref->Mag > EXTRACT_MAX_MAG) continue; + + pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa); + psImage *subraster = psastroExtractStar (inReadout->image, ref->chip->x, ref->chip->y, 200.0, 200.0); + if (!subraster) continue; + + // generate a (very) basic header: + psMetadata *header = psMetadataAlloc (); + psMetadataAddF32 (header, PS_LIST_TAIL, "CHIP_X", 0, "chip coordinate", ref->chip->x); + psMetadataAddF32 (header, PS_LIST_TAIL, "CHIP_Y", 0, "chip coordinate", ref->chip->y); + psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_X", 0, "focal-plane coordinate", ref->FP->x); + psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_Y", 0, "focal-plane coordinate", ref->FP->y); + psMetadataAddF32 (header, PS_LIST_TAIL, "MAG", 0, "magnitude", ref->Mag); + + snprintf (extname, 80, "extname.%05d", nExt); + psFitsWriteImage (outStars, header, subraster, 0, extname); + nExt ++; + + psFree (header); + psFree (subraster); + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Readout AFTER");; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Cell AFTER");; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on Chip AFTER");; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE("failed on FPA AFTER");; + + psFitsClose (outStars); + psFree (view); + return true; +}