Index: /trunk/psastro/src/Makefile.am =================================================================== --- /trunk/psastro/src/Makefile.am (revision 23809) +++ /trunk/psastro/src/Makefile.am (revision 23810) @@ -65,5 +65,4 @@ psastroExtractStars.c \ psastroExtractGhosts.c \ - psastroExtractFindChip.c \ psastroCleanup.c @@ -88,9 +87,12 @@ libpsastro_la_SOURCES = \ + psastroMaskUpdates.c \ + psastroMaskUtils.c \ + psastroLoadGhosts.c \ + psastroGhostUtils.c \ psastroErrorCodes.c \ psastroVersion.c \ psastroDefineFiles.c \ psastroAnalysis.c \ - psastroMaskUpdates.c \ psastroAstromGuess.c \ psastroLoadRefstars.c \ @@ -117,4 +119,5 @@ psastroMosaicSetAstrom.c \ psastroMosaicSetMatch.c \ + psastroFindChip.c \ psastroZeroPoint.c \ psastroDemoDump.c \ Index: /trunk/psastro/src/psastro.c =================================================================== --- /trunk/psastro/src/psastro.c (revision 23809) +++ /trunk/psastro/src/psastro.c (revision 23810) @@ -70,4 +70,5 @@ } + psFree (stats); psLogMsg("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete")); Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 23809) +++ /trunk/psastro/src/psastro.h (revision 23810) @@ -45,4 +45,19 @@ } pmLumFunc; #endif + +/** + * this structure defines the parameters to describe a ghost + */ +typedef struct { + psPlane *srcFP; ///< location in FPA coords of the source star + psPlane *FP; ///< location in FPA coords of the ghost center + psPlane *chip; ///< location in chip coords of the ghost center + double Mag; ///< instrumental magnitude of source star + psEllipseAxes inner; ///< inner elliptical annulus boundary + psEllipseAxes outer; ///< outer elliptical annulus boundary +} psastroGhost; + +psastroGhost *psastroGhostAlloc (void); +bool psastroLoadGhosts (pmConfig *config); bool psastroDataSave (pmConfig *config, psMetadata *stats); @@ -131,12 +146,27 @@ bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe); +// masking functions +pmCell *pmCellInChip (pmChip *chip, float x, float y); +bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip); +bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell); + +bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY); +bool psastroMaskEllipse (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes axes); +bool psastroMaskEllipticalAnnulus (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes eInner, psEllipseAxes eOuter); + +bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta); +void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW); +void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); +void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1); + +pmChip *psastroFindChip (double *xChip, double *yChip, pmFPA *fpa, double xFPA, double yFPA); +bool psastroChipBounds (pmFPA *fpa); + // psastroExtract functions 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); +psImage *psastroExtractStar (psImage *input, double x, double y, double dX, double dY); bool psastroExtractParseCamera (pmConfig *config); bool psastroExtractDataLoad (pmConfig *config); Index: unk/psastro/src/psastroExtractFindChip.c =================================================================== --- /trunk/psastro/src/psastroExtractFindChip.c (revision 23809) +++ (revision ) @@ -1,119 +1,0 @@ -/** @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) goto next_chip; - if (ptCH.x >= region->x1) goto next_chip; - if (ptCH.y < region->y0) goto next_chip; - if (ptCH.y >= region->y1) goto next_chip; - psFree (region); - - *xChip = ptCH.x; - *yChip = ptCH.y; - return chip; - - next_chip: - psFree (region); - } - - return NULL; -} - -bool psastroExtractFreeChipBounds () { - - psFree (chipXmin); - psFree (chipXmax); - psFree (chipYmin); - psFree (chipYmax); - return true; -} Index: /trunk/psastro/src/psastroExtractGhosts.c =================================================================== --- /trunk/psastro/src/psastroExtractGhosts.c (revision 23809) +++ /trunk/psastro/src/psastroExtractGhosts.c (revision 23810) @@ -89,9 +89,9 @@ // XXX can make a more clever model... - float xFPA = -1*ref->FP->x; - float yFPA = -1*ref->FP->y; + double xFPA = -1*ref->FP->x; + double yFPA = -1*ref->FP->y; - float xChip, yChip; - pmChip *ghostChip = psastroExtractFindChip (&xChip, &yChip, fpa, xFPA, yFPA); + double xChip, yChip; + pmChip *ghostChip = psastroFindChip (&xChip, &yChip, fpa, xFPA, yFPA); if (!ghostChip) continue; if (!ghostChip->cells) continue; Index: /trunk/psastro/src/psastroExtractStar.c =================================================================== --- /trunk/psastro/src/psastroExtractStar.c (revision 23809) +++ /trunk/psastro/src/psastroExtractStar.c (revision 23810) @@ -13,5 +13,5 @@ # include "psastroInternal.h" -psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY) { +psImage *psastroExtractStar (psImage *input, double x, double y, double dX, double dY) { // skip if no pixels are on the image Index: /trunk/psastro/src/psastroFindChip.c =================================================================== --- /trunk/psastro/src/psastroFindChip.c (revision 23810) +++ /trunk/psastro/src/psastroFindChip.c (revision 23810) @@ -0,0 +1,119 @@ +/** @file psastroFindChip.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 psastroChipBounds (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 *psastroFindChip (double *xChip, double *yChip, pmFPA *fpa, double xFPA, double yFPA) { + + *xChip = NAN; + *yChip = NAN; + + if (!chipXmin) { + psastroChipBounds (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) goto next_chip; + if (ptCH.x >= region->x1) goto next_chip; + if (ptCH.y < region->y0) goto next_chip; + if (ptCH.y >= region->y1) goto next_chip; + psFree (region); + + *xChip = ptCH.x; + *yChip = ptCH.y; + return chip; + + next_chip: + psFree (region); + } + + return NULL; +} + +bool psastroExtractFreeChipBounds () { + + psFree (chipXmin); + psFree (chipXmax); + psFree (chipYmin); + psFree (chipYmax); + return true; +} Index: /trunk/psastro/src/psastroGhostUtils.c =================================================================== --- /trunk/psastro/src/psastroGhostUtils.c (revision 23810) +++ /trunk/psastro/src/psastroGhostUtils.c (revision 23810) @@ -0,0 +1,47 @@ +/** @file psastroGhostUtils.c + * + * @brief Ghost support functions + * + * @ingroup libpsastro + * + * @author IfA + * @version $Revision: 1.19 $ + * @date $Date: 2009-02-09 21:25:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +static void psastroGhostFree (psastroGhost *ghost) { + + if (ghost == NULL) return; + + psFree (ghost->srcFP); + psFree (ghost->FP); + psFree (ghost->chip); + + return; +} + +psastroGhost *psastroGhostAlloc (void) { + + psastroGhost *ghost = (psastroGhost *) psAlloc(sizeof(psastroGhost)); + psMemSetDeallocator(ghost, (psFreeFunc) psastroGhostFree); + + ghost->srcFP = psPlaneAlloc(); + ghost->FP = psPlaneAlloc(); + ghost->chip = psPlaneAlloc(); + + ghost->Mag = 0.0; + + ghost->inner.major = 0.0; + ghost->inner.minor = 0.0; + ghost->inner.theta = 0.0; + + ghost->outer.major = 0.0; + ghost->outer.minor = 0.0; + ghost->outer.theta = 0.0; + + return ghost; +} + Index: /trunk/psastro/src/psastroLoadGhosts.c =================================================================== --- /trunk/psastro/src/psastroLoadGhosts.c (revision 23810) +++ /trunk/psastro/src/psastroLoadGhosts.c (revision 23810) @@ -0,0 +1,220 @@ +/** @file psastroLoadGhosts.c + * + * @brief calculate ghost FPA and Chip positions for the stars loaded on the FPA based on the model + * + * @ingroup psastro + * + * @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 GET_2D_POLY(NAME,OUT) \ + md = psMetadataLookupMetadata (&status, ghostModel, NAME); \ + if (!md) { \ + psError(PSASTRO_ERR_CONFIG, true, "Missing %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } \ + OUT = psPolynomial2DfromMetadata(md); \ + if (!OUT) { \ + psError(PSASTRO_ERR_CONFIG, true, "Trouble interpretting %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } + +# define GET_1D_POLY(NAME,OUT) \ + md = psMetadataLookupMetadata (&status, ghostModel, NAME); \ + if (!md) { \ + psError(PSASTRO_ERR_CONFIG, true, "Missing %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } \ + OUT = psPolynomial1DfromMetadata(md); \ + if (!OUT) { \ + psError(PSASTRO_ERR_CONFIG, true, "Trouble interpretting %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } + +/** + * calculate ghost FPA and Chip positions for the stars loaded on the FPA + */ +bool psastroLoadGhosts (pmConfig *config) { + + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + float zeropt, exptime; + psMetadata *md = NULL; + psPolynomial2D *centerX = NULL; + psPolynomial2D *centerY = NULL; + psPolynomial1D *outerMajor = NULL; + psPolynomial1D *outerMinor = NULL; + psPolynomial1D *innerMajor = NULL; + psPolynomial1D *innerMinor = NULL; + psMetadata *ghostModel = NULL; + + psLogMsg ("psastro", PS_LOG_INFO, "determine ghost positions"); + + // 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; + } + + bool REFSTAR_MASK_GHOST = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GHOST"); + if (!REFSTAR_MASK_GHOST) return true; + + char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL"); + if (!strcasecmp(ghostFile, "NONE")) return true; + + double MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG"); + + if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) { + psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + + GET_2D_POLY ("GHOST.CENTER.X", centerX); + GET_2D_POLY ("GHOST.CENTER.Y", centerY); + + GET_1D_POLY ("GHOST.OUTER.MAJOR", outerMajor); + GET_1D_POLY ("GHOST.OUTER.MINOR", outerMinor); + GET_1D_POLY ("GHOST.INNER.MAJOR", innerMajor); + GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor); + + // select the input astrometry data (also carries the refstars) + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!astrom) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + goto escape; + } + pmFPA *fpa = astrom->fpa; + + // 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"); + goto escape; + } + + // 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); + MAX_MAG += MagOffset; + + // 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; } + + 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 > MAX_MAG) continue; + + // Ghost model: + // X(ghost) = -X + dX(X,Y) --> GHOST.CENTER.X + // Y(ghost) = -Y + dY(X,Y) --> GHOST.CENTER.Y + // R1inner = R1inner_0 + r dR1inner + // R2inner = R2inner_0 + r dR1inner + // R1outer = R1outer_0 + r dR1inner + // R2outer = R2outer_0 + r dR1inner + + psastroGhost *ghost = psastroGhostAlloc (); + ghost->srcFP->x = ref->FP->x; + ghost->srcFP->y = ref->FP->y; + + // XXX it is stupid that this takes -X_fpa,-Y_fpa --> the analysis script is reporting the guess FPA position, and the fit is using that... + ghost->FP->x = -ref->FP->x + psPolynomial2DEval(centerX, -ref->FP->x, -ref->FP->y); + ghost->FP->y = -ref->FP->y + psPolynomial2DEval(centerY, -ref->FP->x, -ref->FP->y); + + float rSrc = hypot (ref->FP->x, ref->FP->y); + + ghost->inner.major = psPolynomial1DEval (innerMajor, rSrc); + ghost->inner.minor = psPolynomial1DEval (innerMinor, rSrc); + ghost->inner.theta = atan2(ref->FP->y, ref->FP->x); + + ghost->outer.major = psPolynomial1DEval (outerMajor, rSrc); + ghost->outer.minor = psPolynomial1DEval (outerMinor, rSrc); + ghost->outer.theta = atan2(ref->FP->y, ref->FP->x); + + // report instrumental ghost star mags + ghost->Mag = ref->Mag - MagOffset; + + // XXX this code yields a single chip: we need to provide results for any chips + // which encompass the full size of the ghost + pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, -ghost->srcFP->x, -ghost->srcFP->y); + fprintf (stderr, "raw chip position: %f, %f ", ghost->chip->x, ghost->chip->y); + + ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y); + fprintf (stderr, "-> model chip position: %f, %f\n", ghost->chip->x, ghost->chip->y); + + if (!ghostChip) goto skip; + if (!ghostChip->cells) goto skip; + if (!ghostChip->cells->n) goto skip; + pmCell *ghostCell = ghostChip->cells->data[0]; + if (!ghostCell) goto skip; + if (!ghostCell->readouts) goto skip; + if (!ghostCell->readouts->n) goto skip; + pmReadout *ghostReadout = ghostCell->readouts->data[0]; + if (!ghostReadout) goto skip; + + psArray *ghosts = psMetadataLookupPtr (&status, ghostReadout->analysis, "PSASTRO.GHOSTS"); + if (ghosts == NULL) { + ghosts = psArrayAllocEmpty (100); + if (!psMetadataAdd (ghostReadout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts)) { + psError(PSASTRO_ERR_CONFIG, false, "failure to add ghosts to readout"); + goto escape; + } + psFree (ghosts); + } + + psArrayAdd (ghosts, 100, ghost); + + skip: + psFree (ghost); + } + } + } + } + + psastroExtractFreeChipBounds(); + + psFree (centerX); + psFree (centerY); + psFree (innerMajor); + psFree (innerMinor); + psFree (outerMajor); + psFree (outerMinor); + psFree (ghostModel); + psFree (view); + return true; + +escape: + psFree (centerX); + psFree (centerY); + psFree (innerMajor); + psFree (innerMinor); + psFree (outerMajor); + psFree (outerMinor); + psFree (ghostModel); + psFree (view); + return false; +} Index: /trunk/psastro/src/psastroMaskMisc.c =================================================================== --- /trunk/psastro/src/psastroMaskMisc.c (revision 23810) +++ /trunk/psastro/src/psastroMaskMisc.c (revision 23810) @@ -0,0 +1,51 @@ +/** storage : I don't really need these print functions any more ***/ +# if (0) + bool REFSTAR_MASK_REGIONS = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_REGIONS"); + + // text region files for testing + FILE *f = NULL; + if (REFSTAR_MASK_REGIONS) { + char *filename = NULL; + char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME"); + psStringAppend (&filename, "refstars.mask.%s.dat", chipname); + FILE *f = fopen (filename, "w"); + if (!f) { + psWarning ("cannot create refstar mask file %s\n", filename); + continue; + } + psFree (filename); + } + + if (REFSTAR_MASK_REGIONS) { + fprintf (f, "CIRCLE %f %f %f %f\n", ref->chip->x, ref->chip->y, radius, radius); + } + + if (REFSTAR_MASK_REGIONS) { + // lower side + x0 = ref->chip->x + spikeWidth*sin(Theta); + y0 = ref->chip->y - spikeWidth*cos(Theta); + x1 = ref->chip->x + spikeLength*cos(Theta) + spikeWidth*sin(Theta); + y1 = ref->chip->y + spikeLength*sin(Theta) - spikeWidth*cos(Theta); + dx = x1 - x0; + dy = y1 - y0; + + fprintf (f, "LINE %f %f %f %f\n", x0, y0, dx, dy); + + // upper side + x0 = ref->chip->x - spikeWidth*sin(Theta); + y0 = ref->chip->y + spikeWidth*cos(Theta); + x1 = ref->chip->x + spikeLength*cos(Theta) - spikeWidth*sin(Theta); + y1 = ref->chip->y + spikeLength*sin(Theta) + spikeWidth*cos(Theta); + dx = x1 - x0; + dy = y1 - y0; + fprintf (f, "LINE %f %f %f %f\n", x0, y0, dx, dy); + } + + if (REFSTAR_MASK_REGIONS) { + fprintf (f, "LINE %f %f %f %f\n", ref->chip->x, ref->chip->y, 0.0, -100.0); + } + + if (REFSTAR_MASK_REGIONS) { + fclose (f); + } +# endif Index: /trunk/psastro/src/psastroMaskUpdates.c =================================================================== --- /trunk/psastro/src/psastroMaskUpdates.c (revision 23809) +++ /trunk/psastro/src/psastroMaskUpdates.c (revision 23810) @@ -19,14 +19,4 @@ } -pmCell *pmCellInChip (pmChip *chip, float x, float y); -bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip); -bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell); - -bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY); -bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta); -void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW); -void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); -void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1); - /** * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip @@ -40,5 +30,7 @@ float zeropt, exptime; - psImageMaskType maskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels + psImageMaskType ghostMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels + psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels + psImageMaskType starMaskValue = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels // psImageMaskType maskBlank = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels @@ -55,7 +47,12 @@ if (!REFSTAR_MASK) return true; + // convert star positions to ghost positions and add to the readout->analysis data + if (!psastroLoadGhosts (config)) { + psError(PSASTRO_ERR_CONFIG, false, "Error loading ghosts"); + return false; + } + psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask"); - bool REFSTAR_MASK_REGIONS = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_REGIONS"); bool REFSTAR_MASK_BLEED = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED"); @@ -134,18 +131,4 @@ if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE; - // text region files for testing - FILE *f = NULL; - if (REFSTAR_MASK_REGIONS) { - char *filename = NULL; - char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME"); - psStringAppend (&filename, "refstars.mask.%s.dat", chipname); - FILE *f = fopen (filename, "w"); - if (!f) { - psWarning ("cannot create refstar mask file %s\n", filename); - continue; - } - psFree (filename); - } - pmChip *refChip = pmFPAviewThisChip (view, refMask->fpa); @@ -173,9 +156,9 @@ // we mask pixels on the input mask image (chip-mosaic) - pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa); - pmReadout *readoutMask = NULL; - if (cellMask->readouts->n) { - readoutMask = cellMask->readouts->data[0]; - } + // pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa); + // pmReadout *readoutMask = NULL; + // if (cellMask->readouts->n) { + // readoutMask = cellMask->readouts->data[0]; + // } // process each of the readouts @@ -184,7 +167,11 @@ if (! readout->data_exists) { continue; } + // XXX why not do this? + pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa); + if (!readoutMask) continue; + // select the raw objects for this readout psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) { continue; } + if (!refstars) continue; // we need to generate the following masks regions: @@ -204,46 +191,21 @@ float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag); - if (REFSTAR_MASK_REGIONS) { - fprintf (f, "CIRCLE %f %f %f %f\n", ref->chip->x, ref->chip->y, radius, radius); - } - // XXX for now, assume cell binning is 1x1 relative to chip - if (readoutMask) { - psastroMaskCircle (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, radius, radius); - } - - // LINE for boundaries of the saturation spikes (scaled by magnitude) - float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (REFSTAR_MASK_SATSPIKE_MAG_MAX - ref->Mag); - float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH; + psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius); for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) { - float x0, y0, x1, y1, dx, dy; float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO; - if (REFSTAR_MASK_REGIONS) { - // lower side - x0 = ref->chip->x + spikeWidth*sin(Theta); - y0 = ref->chip->y - spikeWidth*cos(Theta); - x1 = ref->chip->x + spikeLength*cos(Theta) + spikeWidth*sin(Theta); - y1 = ref->chip->y + spikeLength*sin(Theta) - spikeWidth*cos(Theta); - dx = x1 - x0; - dy = y1 - y0; - - fprintf (f, "LINE %f %f %f %f\n", x0, y0, dx, dy); - - // upper side - x0 = ref->chip->x - spikeWidth*sin(Theta); - y0 = ref->chip->y + spikeWidth*cos(Theta); - x1 = ref->chip->x + spikeLength*cos(Theta) - spikeWidth*sin(Theta); - y1 = ref->chip->y + spikeLength*sin(Theta) + spikeWidth*cos(Theta); - dx = x1 - x0; - dy = y1 - y0; - fprintf (f, "LINE %f %f %f %f\n", x0, y0, dx, dy); - } - - if (readoutMask) { - psastroMaskBox (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta); - } + // LINE for boundaries of the saturation spikes (scaled by magnitude) + // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2; + + float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX; + float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag); + float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH; + // XXX we can make the width depend on the spike as well... + // The length should also be a function of the image background level + + psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta); } @@ -254,9 +216,5 @@ // LINE for boundaries of the bleed lines - if (REFSTAR_MASK_REGIONS) { - fprintf (f, "LINE %f %f %f %f\n", ref->chip->x, ref->chip->y, 0.0, -100.0); - } - - if (readoutMask && refCell) { + if (refCell) { float xCell = 0.0; float yCell = 0.0; @@ -271,12 +229,34 @@ float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); - psastroMaskRectangle (readoutMask->mask, maskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd); + psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd); } } } + + // select the raw objects for this readout (loaded in psastroExtract.c) + psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS"); + if (ghosts == NULL) { continue; } + + // identify the bright stars of interest + for (int i = 0; i < ghosts->n; i++) { + psastroGhost *ghost = ghosts->data[i]; + // XXX bright vs faint ghost bits? (OR with SUSPECT) + psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer); + } + + // select the raw objects for this readout, flag is they fall in a mask + psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) return false; + + // XXX finish this: + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[(int)(raw->chip->x)][(int)(raw->chip->y)]; + if (value) continue; + } } } - // output sequence for mask corresponding to this chip + // output sequence for mask corresponding to this chip (XXX this may not be needed...) *viewMask = *view; while ((cell = pmFPAviewNextCell (viewMask, outMask->fpa, 1)) != NULL) { @@ -290,9 +270,5 @@ if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; - if (REFSTAR_MASK_REGIONS) { - fclose (f); - } } if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; @@ -306,220 +282,2 @@ } -// XXX this is going to be very slow... -pmCell *pmCellInChip (pmChip *chip, float x, float y) { - -# if (0) - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - - // XXX fix the binning : currently not selected from concepts - // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y - // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y - int xBin = 1; - int yBin = 1; - - // Position on the cell - float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); - float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); -# endif - - for (int i = 0; i < chip->cells->n; i++) { - - pmCell *cell = chip->cells->data[i]; - psRegion *region = pmCellExtent (cell); - - if (x < region->x0) goto skip; - if (x > region->x1) goto skip; - if (y < region->y0) goto skip; - if (y > region->y1) goto skip; - - psFree (region); - return cell; - - skip: - psFree (region); - } - return NULL; -} - -/** - * convert chip coords to cell coords, given known cell - */ -bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip) { - - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - - // XXX fix the binning : currently not selected from concepts - // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y - // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y - int xBin = 1; - int yBin = 1; - - // Position on the cell - // ((pos)*(binning)*(cellParity) + (cell0)) - // XXX this is probably totally wrong now.... - // ((pos) - (cell0))*(cellParity)/(binning)) - *xCell = (xChip - x0Cell)*xParityCell/xBin; - *yCell = (yChip - y0Cell)*yParityCell/yBin; - - return true; -} - -/** - * convert chip coords to cell coords, given known cell - */ -bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell) { - - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - - // XXX fix the binning : currently not selected from concepts - // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y - // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y - int xBin = 1; - int yBin = 1; - - // Position on the cell - // ((pos)*(binning)*(cellParity) + (cell0)) - // XXX this is probably totally wrong now.... - // ((pos) - (cell0))*(cellParity)/(binning)) - *xChip = xCell*xBin*xParityCell + x0Cell; - *yChip = yCell*yBin*yParityCell + y0Cell; - - return true; -} - -// XXX should be doing an OR -bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY) { - - // XXX need to worry about row0, col0 - for (int ix = -dX; ix <= +dX; ix++) { - int jx = ix + x0; - if (jx < 0) continue; - if (jx >= mask->numCols) continue; - for (int iy = -dY; iy <= +dY; iy++) { - int jy = iy + y0; - if (jy < 0) continue; - if (jy >= mask->numRows) continue; - - double r2 = PS_SQR(ix/dX) + PS_SQR(iy/dY); - if (r2 > 1.0) continue; - - mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value; - } - } - return true; -} - -// XXX should be doing an OR -bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta) { - - // draw a series of lines (from -0.5*dW to +0.5*dW) of length dL, starting at x0, y0, angle theta - - float xs = x0; - float ys = y0; - - float xe = xs + dL*cos(theta); - float ye = ys + dL*sin(theta); - - psastroMaskLine (mask, value, xs, ys, xe, ye, (int) dW); - return true; -} - -/** - * identify the quadrant and draw the correct line - */ -void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW) { - - int FlipDirect, FlipCoords; - int X1, Y1, X2, Y2, dX, dY; - - /* rather than draw the line from float positions, we find the closest - integer end-points and draw the line between those pixels */ - - X1 = ROUND(x1); - Y1 = ROUND(y1); - X2 = ROUND(x2); - Y2 = ROUND(y2); - - dX = X2 - X1; - dY = Y2 - Y1; - - FlipCoords = (abs(dX) < abs(dY)); - FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2); - - if (!FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X1, Y1, X2, Y2, dW, FALSE); - if ( FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X2, Y2, X1, Y1, dW, FALSE); - if (!FlipDirect && FlipCoords) psastroMaskLineBresen (mask, value, Y1, X1, Y2, X2, dW, TRUE); - if ( FlipDirect && FlipCoords) psastroMaskLineBresen (mask, value, Y2, X2, Y1, X1, dW, TRUE); - - return; -} - -/** - * use the Bresenham line drawing technique - * integer-only Bresenham line-draw version which is fast - */ -void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) { - - int X, Y, dX, dY; - int e, e2; - - dX = X2 - X1; - dY = Y2 - Y1; - - Y = Y1; - e = 0; - for (X = X1; X <= X2; X++) { - if (X > 0) { - if (swapcoords) { - if (X >= mask->numRows) continue; - for (int y = Y - dW; y <= Y + dW; y++) { - if (y < 0) continue; - if (y >= mask->numCols) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[X][y] |= value; - } - } else { - if (X >= mask->numCols) continue; - for (int y = Y - dW; y <= Y + dW; y++) { - if (y < 0) continue; - if (y >= mask->numRows) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][X] |= value; - } - } - } - e += dY; - e2 = 2 * e; - if (e2 > dX) { - Y++; - e -= dX; - } - if (e2 < -dX) { - Y--; - e += dX; - } - } - return; -} - -void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1) { - - int xs = PS_MAX (0, PS_MIN (mask->numCols, PS_MIN (x0, x1))); - int xe = PS_MAX (0, PS_MIN (mask->numCols, PS_MAX (x0, x1))); - int ys = PS_MAX (0, PS_MIN (mask->numRows, PS_MIN (y0, y1))); - int ye = PS_MAX (0, PS_MIN (mask->numRows, PS_MAX (y0, y1))); - - for (int iy = ys; iy < ye; iy++) { - for (int ix = xs; ix < xe; ix++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= value; - } - } -} - Index: /trunk/psastro/src/psastroMaskUtils.c =================================================================== --- /trunk/psastro/src/psastroMaskUtils.c (revision 23810) +++ /trunk/psastro/src/psastroMaskUtils.c (revision 23810) @@ -0,0 +1,341 @@ +/** @file psastroMaskUpdates.c + * + * @brief + * + * @ingroup libpsastro + * + * @author IfA + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +# define ESCAPE { \ + psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate"); \ + psFree (view); \ + return false; \ +} + +// XXX this is going to be very slow... +pmCell *pmCellInChip (pmChip *chip, float x, float y) { + +# if (0) + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + // XXX fix the binning : currently not selected from concepts + // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y + // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y + int xBin = 1; + int yBin = 1; + + // Position on the cell + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); +# endif + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + psRegion *region = pmCellExtent (cell); + + if (x < region->x0) goto skip; + if (x > region->x1) goto skip; + if (y < region->y0) goto skip; + if (y > region->y1) goto skip; + + psFree (region); + return cell; + + skip: + psFree (region); + } + return NULL; +} + +/** + * convert chip coords to cell coords, given known cell + */ +bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip) { + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + // XXX fix the binning : currently not selected from concepts + // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y + // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y + int xBin = 1; + int yBin = 1; + + // Position on the cell + // ((pos)*(binning)*(cellParity) + (cell0)) + // XXX this is probably totally wrong now.... + // ((pos) - (cell0))*(cellParity)/(binning)) + *xCell = (xChip - x0Cell)*xParityCell/xBin; + *yCell = (yChip - y0Cell)*yParityCell/yBin; + + return true; +} + +/** + * convert chip coords to cell coords, given known cell + */ +bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell) { + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + // XXX fix the binning : currently not selected from concepts + // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y + // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y + int xBin = 1; + int yBin = 1; + + // Position on the cell + // ((pos)*(binning)*(cellParity) + (cell0)) + // XXX this is probably totally wrong now.... + // ((pos) - (cell0))*(cellParity)/(binning)) + *xChip = xCell*xBin*xParityCell + x0Cell; + *yChip = yCell*yBin*yParityCell + y0Cell; + + return true; +} + +bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY) { + + for (int ix = -dX; ix <= +dX; ix++) { + int jx = ix + x0; + if (jx < 0) continue; + if (jx >= mask->numCols) continue; + for (int iy = -dY; iy <= +dY; iy++) { + int jy = iy + y0; + if (jy < 0) continue; + if (jy >= mask->numRows) continue; + + double r2 = PS_SQR(ix/dX) + PS_SQR(iy/dY); + if (r2 > 1.0) continue; + + mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value; + } + } + return true; +} + + +void psastroMinMaxForShape (float *min, float *max, float y, psEllipseShape shape) { + + // XXX optimize this + float A = 1.0; + float B = shape.sxy*y*PS_SQR(shape.sx); + float C = PS_SQR(y*shape.sx/shape.sy) - PS_SQR(shape.sx); + + float T1 = PS_SQR(B); + float T2 = 4.0*A*C; + float R = T1 - T2; + if (R < 0) R = 0; + + *min = (-B - sqrt (R)) / (2.0*A); + *max = (-B + sqrt (R)) / (2.0*A); +} + +bool psastroMaskEllipse (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes axes) { + + psEllipseShape shape = psEllipseAxesToShape (axes); + + // phi is the coordinate along the elliptical path + // phiMin, phiMax are the points on the path for Ymin and Ymax + float phiMin = atan2 (-axes.minor*cos(axes.theta), axes.major*sin(axes.theta)); + float phiMax = phiMin + M_PI; + + float Ymin = -axes.major*cos(phiMin)*sin(axes.theta) + axes.minor*sin(phiMin)*cos(axes.theta); + float Ymax = -axes.major*cos(phiMax)*sin(axes.theta) + axes.minor*sin(phiMax)*cos(axes.theta); + if (Ymin > Ymax) PS_SWAP (Ymin, Ymax); + + for (int iy = Ymin; iy <= Ymax; iy++) { + int jy = iy + y0; + if (jy < 0) continue; + if (jy >= mask->numRows) continue; + + float Xmin, Xmax; + psastroMinMaxForShape (&Xmin, &Xmax, iy, shape); + fprintf (stderr, "mask %d : %f -> %f == %f -> %f (%x)\n", jy, Xmin, Xmax, Xmin + x0, Xmax + x0, value); + + for (int ix = Xmin; ix <= Xmax; ix++) { + int jx = ix + x0; + if (jx < 0) continue; + if (jx >= mask->numCols) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value; + } + } + return true; +} + +bool psastroMaskEllipticalAnnulus (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes eInner, psEllipseAxes eOuter) { + + psEllipseShape sInner = psEllipseAxesToShape (eInner); + psEllipseShape sOuter = psEllipseAxesToShape (eOuter); + + // phi is the coordinate along the elliptical path + // phiMin, phiMax are the points on the path for Ymin and Ymax + float phiMinInner = atan2 (-eInner.minor*cos(eInner.theta), eInner.major*sin(eInner.theta)); + float phiMaxInner = phiMinInner + M_PI; + float phiMinOuter = atan2 (-eOuter.minor*cos(eOuter.theta), eOuter.major*sin(eOuter.theta)); + float phiMaxOuter = phiMinOuter + M_PI; + + float YminInner = -eInner.major*cos(phiMinInner)*sin(eInner.theta) + eInner.minor*sin(phiMinInner)*cos(eInner.theta); + float YmaxInner = -eInner.major*cos(phiMaxInner)*sin(eInner.theta) + eInner.minor*sin(phiMaxInner)*cos(eInner.theta); + if (YminInner > YmaxInner) PS_SWAP (YminInner, YmaxInner); + + float YminOuter = -eOuter.major*cos(phiMinOuter)*sin(eOuter.theta) + eOuter.minor*sin(phiMinOuter)*cos(eOuter.theta); + float YmaxOuter = -eOuter.major*cos(phiMaxOuter)*sin(eOuter.theta) + eOuter.minor*sin(phiMaxOuter)*cos(eOuter.theta); + if (YminOuter > YmaxOuter) PS_SWAP (YminOuter, YmaxOuter); + + for (int iy = YminOuter; iy <= YmaxOuter; iy++) { + int jy = iy + y0; + if (jy < 0) continue; + if (jy >= mask->numRows) continue; + + float XminOuter, XmaxOuter; + psastroMinMaxForShape (&XminOuter, &XmaxOuter, iy, sOuter); + + if ((iy > YmaxInner) || (iy < YminInner)) { + for (int ix = XminOuter; ix <= XmaxOuter; ix++) { + int jx = ix + x0; + if (jx < 0) continue; + if (jx >= mask->numCols) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value; + } + } else { + float XminInner, XmaxInner; + psastroMinMaxForShape (&XminInner, &XmaxInner, iy, sInner); + + for (int ix = XminOuter; ix <= XminInner; ix++) { + int jx = ix + x0; + if (jx < 0) continue; + if (jx >= mask->numCols) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value; + } + for (int ix = XmaxInner; ix <= XmaxOuter; ix++) { + int jx = ix + x0; + if (jx < 0) continue; + if (jx >= mask->numCols) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value; + } + } + } + return true; +} + +bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta) { + + // draw a series of lines (from -0.5*dW to +0.5*dW) of length dL, starting at x0, y0, angle theta + + float xs = x0; + float ys = y0; + + float xe = xs + dL*cos(theta); + float ye = ys + dL*sin(theta); + + psastroMaskLine (mask, value, xs, ys, xe, ye, (int) dW); + return true; +} + +/** + * identify the quadrant and draw the correct line + */ +void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW) { + + int FlipDirect, FlipCoords; + int X1, Y1, X2, Y2, dX, dY; + + /* rather than draw the line from float positions, we find the closest + integer end-points and draw the line between those pixels */ + + X1 = ROUND(x1); + Y1 = ROUND(y1); + X2 = ROUND(x2); + Y2 = ROUND(y2); + + dX = X2 - X1; + dY = Y2 - Y1; + + FlipCoords = (abs(dX) < abs(dY)); + FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2); + + if (!FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X1, Y1, X2, Y2, dW, FALSE); + if ( FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X2, Y2, X1, Y1, dW, FALSE); + if (!FlipDirect && FlipCoords) psastroMaskLineBresen (mask, value, Y1, X1, Y2, X2, dW, TRUE); + if ( FlipDirect && FlipCoords) psastroMaskLineBresen (mask, value, Y2, X2, Y1, X1, dW, TRUE); + + return; +} + +/** + * use the Bresenham line drawing technique + * integer-only Bresenham line-draw version which is fast + */ +void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) { + + int X, Y, dX, dY; + int e, e2; + + dX = X2 - X1; + dY = Y2 - Y1; + + Y = Y1; + e = 0; + for (X = X1; X <= X2; X++) { + if (X > 0) { + if (swapcoords) { + if (X >= mask->numRows) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= mask->numCols) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[X][y] |= value; + } + } else { + if (X >= mask->numCols) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= mask->numRows) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][X] |= value; + } + } + } + e += dY; + e2 = 2 * e; + if (e2 > dX) { + Y++; + e -= dX; + } + if (e2 < -dX) { + Y--; + e += dX; + } + } + return; +} + +void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1) { + + int xs = PS_MAX (0, PS_MIN (mask->numCols, PS_MIN (x0, x1))); + int xe = PS_MAX (0, PS_MIN (mask->numCols, PS_MAX (x0, x1))); + int ys = PS_MAX (0, PS_MIN (mask->numRows, PS_MIN (y0, y1))); + int ye = PS_MAX (0, PS_MIN (mask->numRows, PS_MAX (y0, y1))); + + for (int iy = ys; iy < ye; iy++) { + for (int ix = xs; ix < xe; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= value; + } + } +} + +