Index: /trunk/psphot/src/LocalSky.c =================================================================== --- /trunk/psphot/src/LocalSky.c (revision 4215) +++ /trunk/psphot/src/LocalSky.c (revision 4215) @@ -0,0 +1,112 @@ +# include "psphot.h" + +psSource *pmSourceDefinePixels(const psImageData *imdata, + const psPeak *peak, + psF32 Radius) +{ + // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) + psS32 SubImageStartRow = PS_MAX (0, peak->y - Radius); + psS32 SubImageEndRow = PS_MIN (image->numRows, peak->y + Radius + 1); + psS32 SubImageStartCol = PS_MAX (0, peak->x - Radius); + psS32 SubImageEndCol = PS_MIN (image->numCols, peak->x + Radius + 1); + + // Grab a subimage of the original image of size (2 * outerRadius). + psImage *subImage = psImageSubset(imdata->image, + SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow); + + psImage *subImageMask = psImageSubset(imdata->mask, + SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow); + + psImage *subImageNoise = psImageSubset(imdata->noise, + SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow); + + psSource *mySource = pmSourceAlloc(); + mySource->peak = peak; + mySource->pixels = subImage; + mySource->mask = subImageMask; + mySource->noise = subImageNoise; + + return(mySource); +} + +// this sets and clears bit 0x80 +bool pmSourceLocalSky_EAM (psSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + + psImage *image = source->pixels; + psImage *mask = source->mask; + psImage *peak = source->peak; + + // Loop through the subimage, mask off pixels in the inner square. + // this uses a static mask value of 0x08 + psS32 StartCol = PS_MAX (0, peak->x - mask->col0 - Radius); + psS32 StartRow = PS_MAX (0, peak->y - mask->row0 - Radius); + psS32 EndCol = PS_MIN (mask->numCols, peak->x - mask->col0 + Radius + 1); + psS32 EndRow = PS_MIN (mask->numRows, peak->y - mask->row0 + Radius + 1); + for (psS32 row = StartRow; row < EndRow; row++) { + for (psS32 col = StartCol; col < EndCol; col++) { + mask->data.U8[row][col] |= 0x80; + } + } + + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, image, mask, 0xff); + + // clear the mask on the inner pixels + for (psS32 row = StartRow; row < EndRow; row++) { + for (psS32 col = StartCol; col < EndCol; col++) { + mask->data.U8[row][col] &= 0x7f; + } + } + + mySource->moments = pmMomentsAlloc(); + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + mySource->moments->Sky = (psF32) tmpF64; + + return (true); +} + +// set mask pixels within the circle +bool psImageMaskCircle(psImage *mask, + psF32 x, + psF32 y, + psF32 radius, + int maskValue) +{ + // XXX EAM : for the circle to stay on the image + // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) + psS32 SubImageStartRow = PS_MAX (0, x + mask->row0 - radiusS32); + psS32 SubImageEndRow = PS_MIN (mask->numRows, SubImageCenterRow + radiusS32 + 1); + psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); + psS32 SubImageEndCol = PS_MIN (mask->numCols, SubImageCenterCol + radiusS32 + 1); + + // Loop through the subimage mask, set maskValue bits + for (psS32 row = 0 ; row < source->mask->numRows; row++) { + for (psS32 col = 0 ; col < source->mask->numCols; col++) { + + if (CheckRadius2((psF32) radiusS32, + (psF32) radiusS32, + radius, + (psF32) col, + (psF32) row)) { + source->mask->data.U8[row][col] = 0; + } else { + source->mask->data.U8[row][col] = 1; + } + } + } + return(true); +} + Index: /trunk/psphot/src/choose_psf_model.c =================================================================== --- /trunk/psphot/src/choose_psf_model.c (revision 4214) +++ /trunk/psphot/src/choose_psf_model.c (revision 4215) @@ -26,15 +26,7 @@ float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS"); - float XBORDER = psMetadataLookupF32 (&status, config, "XBORDER"); - float YBORDER = psMetadataLookupF32 (&status, config, "YBORDER"); - psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER); - keep = psRegionForImage (keep, image, keep); - for (int i = 0; i < stars->n; i++) { psSource *source = stars->data[i]; pmSourceSetPixelCircle (source, image, RADIUS); - pmSourceMaskRegion (source, keep); - // XXX EAM - this is silly: reconsider how we use the masks - // XXX - we are needlessly re-defining the pixels here } Index: /trunk/psphot/src/find_peaks.c =================================================================== --- /trunk/psphot/src/find_peaks.c (revision 4214) +++ /trunk/psphot/src/find_peaks.c (revision 4215) @@ -1,5 +1,5 @@ # include "psphot.h" -psArray *find_peaks (psImage *image, psMetadata *config, psStats *sky) +psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky) { bool status = false; @@ -16,4 +16,5 @@ NSIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_NSIGMA"); + psImage *image = imdata->image; psImage *smooth = psImageCopy (NULL, image, PS_TYPE_F32); psImageSmooth (smooth, SIGMA, NSIGMA); Index: /trunk/psphot/src/image_stats.c =================================================================== --- /trunk/psphot/src/image_stats.c (revision 4214) +++ /trunk/psphot/src/image_stats.c (revision 4215) @@ -1,8 +1,10 @@ # include "psphot.h" -psStats *image_stats (psImage *image, psMetadata *config) +psStats *image_stats (psImageData *imdata, psMetadata *config) { psStats *stats = NULL; psStats *sky = NULL; + + psImage *image = imdata->image; // get image stats on a subset of the image (random 1e5 pts) Index: /trunk/psphot/src/psphot-utils.c =================================================================== --- /trunk/psphot/src/psphot-utils.c (revision 4214) +++ /trunk/psphot/src/psphot-utils.c (revision 4215) @@ -160,4 +160,34 @@ } +// mask the area contained by the region +void psImageMaskRegion (psImage *image, psRegion *region, int maskValue) { + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + if (ix + image->col0 < region->x0) continue; + if (ix + image->col0 >= region->x1) continue; + if (iy + image->row0 < region->y0) continue; + if (iy + image->row0 >= region->y1) continue; + image->data.U8[iy][ix] |= maskValue; + } + } +} + +// mask the area not contained by the region +void psImageKeepRegion (psImage *image, psRegion *region, int maskValue) { + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + if (ix + image->col0 < region->x0) goto maskit; + if (ix + image->col0 >= region->x1) goto maskit; + if (iy + image->row0 < region->y0) goto maskit; + if (iy + image->row0 >= region->y1) goto maskit; + continue; + maskit: + image->data.U8[iy][ix] |= maskValue; + } + } +} + // create a subset of sources for the PS_SOURCE_PSFSTAR entries // XXX deprecated Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 4214) +++ /trunk/psphot/src/psphot.c (revision 4215) @@ -4,13 +4,10 @@ int main (int argc, char **argv) { - psMetadata *config = NULL; // user-provided configuration information - psMetadata *header = NULL; // input image header - psImage *image = NULL; // input image data - psImage *noise = NULL; // input image data - psImage *mask = NULL; // input image data - psArray *sources = NULL; - psArray *peaks = NULL; // a list of pmPeaks - pmPSF *psf = NULL; - psStats *sky = NULL; + psMetadata *config = NULL; // user-provided configuration information + psImageData *imdata = NULL; + psArray *sources = NULL; + psArray *peaks = NULL; // a list of pmPeaks + pmPSF *psf = NULL; + psStats *sky = NULL; psLogArguments (&argc, argv); @@ -20,14 +17,14 @@ // load input data (image and config) // create or load mask and noise images - setup (&image, &mask, &noise, &header, config); + imdata = setup (config); // measure image stats for initial guess - sky = image_stats (image, config); + sky = image_stats (imdata, config); // find the peaks in the image - peaks = find_peaks (image, config, sky); + peaks = find_peaks (imdata, config, sky); // construct sources and measure basic stats - sources = source_moments (image, config, peaks); + sources = source_moments (imdata, config, peaks); // source analysis is done in S/N order (brightest first) Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 4214) +++ /trunk/psphot/src/psphot.h (revision 4215) @@ -5,4 +5,11 @@ # include # include + +typedef struct { + psImage *image; + psImage *mask; + psImage *noise; + psMetadata *header; +} psImageData; // data to test a given PSF model type Index: /trunk/psphot/src/setup.c =================================================================== --- /trunk/psphot/src/setup.c (revision 4214) +++ /trunk/psphot/src/setup.c (revision 4215) @@ -1,53 +1,62 @@ # include "psphot.h" -bool setup (psImage **image, psImage **mask, psImage **noise, psMetadata **header, psMetadata **config) +psImageData *setup (psMetadata *config) { - int Nfail; - psFits *file = NULL; // FITS pointer to input image psRegion region = {0,0,0,0}; // a region representing the entire array psTimerStart ("psphot"); - char *input = pmConfigLookupSTR (&status, *config, *header, "INPUT"); - file = psFitsAlloc (input); + char *input = psMetadataLookupSTR (&status, config, "INPUT"); + psFits *file = psFitsAlloc (input); - *header = psMetadataAlloc (); - psMetadataAdd (*header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL); - psMetadataAdd (*header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL); - *header = psFitsReadHeader (*header, file); - *image = psFitsReadImage (*image, file, region, 0); + psImageData *imdata = psAlloc(sizeof(psImageData)); + psMetadata *header = psMetadataAlloc (); + psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL); + psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL); + header = psFitsReadHeader (header, file); + psImage *image = psFitsReadImage (NULL, file, region, 0); psFree (file); // grab these values from the approrpiate location (image header if necessary) bool status = false; - float RDNOISE = pmConfigLookupF32 (&status, *config, *header, "RDNOISE"); - float GAIN = pmConfigLookupF32 (&status, *config, *header, "GAIN"); - psMetadataAdd (*config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE); - psMetadataAdd (*config, PS_LIST_TAIL, "GAIN", PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN); + float RDNOISE = pmConfigLookupF32 (&status, config, header, "RDNOISE"); + float GAIN = pmConfigLookupF32 (&status, config, header, "GAIN"); + psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE); + psMetadataAdd (config, PS_LIST_TAIL, "GAIN", PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN); - // load the noise image if it exists, otherwise build from input - char *noiseName = pmConfigLookupSTR (&status, *config, *header, "NOISE"); + imdata->image = image; + imdata->header = header; + + // load the noise image if it is supplied, otherwise build from input + char *noiseName = psMetadataLookupSTR (&status, config, "NOISE"); if (status == true) { file = psFitsAlloc (noiseName); - *noise = psFitsReadImage (*noise, file, region, 0); + imdata->noise = psFitsReadImage (NULL, file, region, 0); psFree (file); } else { - *value = psScalar (1.0 / GAIN); - *noise = BinaryOp (NULL, *image, '/', value); + psScalar *value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64); + imdata->noise = BinaryOp (NULL, imdata->image, '/', value); psFree (value); - *value = psScalar (PS_SQR(RDNOISE/GAIN)); - *noise = BinaryOp (*noise, *noise, '+', value); + psScalar *value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64); + imdata->noise = BinaryOp (imdata->noise, imdata->noise, '+', value); psFree (value); } - char *maskName = pmConfigLookupSTR (&status, *config, *header, "MASK"); + char *maskName = psMetadataLookupSTR (&status, config, "MASK"); if (status == true) { file = psFitsAlloc (maskName); - *mask = psFitsReadImage (*mask, file, region, 0); + imdata->mask = psFitsReadImage (NULL, file, region, 0); + // check if U8 & force? psFree (file); } else { - *mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); + imdata->mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); } + + float XBORDER = psMetadataLookupF32 (&status, config, "XBORDER"); + float YBORDER = psMetadataLookupF32 (&status, config, "YBORDER"); + psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER); + keep = psRegionForImage (keep, image, keep); + psImageKeepRegion (imdata->mask, keep, 0x01); psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); Index: /trunk/psphot/src/source_moments.c =================================================================== --- /trunk/psphot/src/source_moments.c (revision 4214) +++ /trunk/psphot/src/source_moments.c (revision 4215) @@ -1,5 +1,5 @@ # include "psphot.h" -psArray *source_moments (psImage *image, psMetadata *config, psArray *allpeaks) +psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks) { bool status = false; @@ -8,16 +8,4 @@ psTimerStart ("psphot"); - - // reduce the initial peaks to a subset of the better peaks - // is this the appropriate place for this culling? - - float XBORDER = psMetadataLookupF32 (&status, config, "XBORDER"); - float YBORDER = psMetadataLookupF32 (&status, config, "YBORDER"); - float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE"); - - psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER); - keep = psRegionForImage (keep, image, keep); - // I need to exclude the saturated objects here, but include them again later...? - peaks = pmPeaksSubset (allpeaks, 100000, keep); // determine properties (sky, moments) of initial sources @@ -32,7 +20,8 @@ for (int i = 0; i < peaks->n; i++) { // provide the user arguments as metadata? - psSource *source = pmSourceLocalSky (image, (psPeak *)peaks->data[i], PS_STAT_SAMPLE_MEDIAN, INNER, OUTER); + psSource *source = pmSourceDefinePixels (imdata, (psPeak *)peaks->data[i], OUTER); if (source == NULL) continue; - pmSourceMaskRegion (source, keep); + + pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER); pmSourceMoments (source, RADIUS); psArrayAdd (sources, 100, source);