Index: unk/psphot/src/modelTestArguments.c =================================================================== --- /trunk/psphot/src/modelTestArguments.c (revision 5985) +++ (revision ) @@ -1,73 +1,0 @@ -# include "psphot.h" -static int usage (); - -psMetadata *modelTestArguments (int *argc, char **argv) { - - int N; - unsigned int Nfail; - int mode = PS_DATA_STRING | PS_META_REPLACE; - psF32 UseCoords_X, UseCoords_Y; - - // basic pslib options - psLogSetFormat ("M"); - psArgumentVerbosity (argc, argv); - - // identify options in args list - these go on config - bool UseCoords = false; - if ((N = psArgumentGet (*argc, argv, "-coords"))) { - UseCoords = true; - psArgumentRemove (N, argc, argv); - UseCoords_X = atof (argv[N]); - psArgumentRemove (N, argc, argv); - UseCoords_Y = atof (argv[N]); - psArgumentRemove (N, argc, argv); - } - - char *model = NULL; - if ((N = psArgumentGet (*argc, argv, "-model"))) { - psArgumentRemove (N, argc, argv); - model = psStringCopy (argv[N]); - psArgumentRemove (N, argc, argv); - } - - char *fitmode = NULL; - if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { - psArgumentRemove (N, argc, argv); - fitmode = psStringCopy (argv[N]); - psArgumentRemove (N, argc, argv); - } - - if (*argc != 3) usage (); - - // load config information - psMetadata *config = psMetadataAlloc (); - psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); - config = psMetadataConfigParse (config, &Nfail, argv[2], FALSE); - - // identify input image - psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); - - if (model != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODEL", mode, "", model); - } - if (fitmode != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODE", mode, "", fitmode); - } - if (UseCoords) { - mode = PS_DATA_F32 | PS_META_REPLACE; - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", mode, "", UseCoords_X); - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", mode, "", UseCoords_Y); - } - - return (config); -} - -static int usage () { - - fprintf (stderr, "USAGE: psphotModelTest (image.fits) (config)\n"); - fprintf (stderr, "options: \n"); - fprintf (stderr, " -coords x y : specify object center\n"); - fprintf (stderr, " -model (name) : specify model of interest\n"); - fprintf (stderr, " -fitmode (flt|psf) : specify how to fit model\n"); - exit (2); -} Index: unk/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5985) +++ (revision ) @@ -1,48 +1,0 @@ -# include "psphot.h" - -psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky) -{ - bool status = false; - float NSIGMA; - float SIGMA; - float threshold; - - // smooth the image - - psTimerStart ("psphot"); - - SIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_SIGMA"); - NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA"); - - psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); - psImageSmooth (smooth, SIGMA, NSIGMA); - psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); - - psTimerStart ("psphot"); - - // set peak threshold - NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT"); - - threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; - // threshold = NSIGMA*sky->sampleStdev; - psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold); - - // find the peaks in the smoothed image - psArray *peaks = pmFindImagePeaks (smooth, threshold); - if (peaks == NULL) psAbort ("find peaks", "no peaks found"); - psFree (smooth); - - // optional dump of all peak data - char *output = psMetadataLookupPtr (&status, config, "PEAKS_OUTPUT_FILE"); - if (status && (output != NULL) && (output[0])) { - pmPeaksWriteText (peaks, output); - psFree (output); - } - - psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); - return (peaks); -} - -// In this function, we smooth the image, then search for the peaks -// Should we also subtract a super-binned image? (as an option?) -// We need to gracefully handle no source detections Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5985) +++ /trunk/psphot/src/psphot.c (revision 5986) @@ -1,4 +1,5 @@ # include "psphot.h" +// XXX need a better structure for handling optional sequence int main (int argc, char **argv) { @@ -9,28 +10,16 @@ pmPSF *psf = NULL; psStats *sky = NULL; - pmPSFClump psfClump; bool status; psTimerStart ("complete"); - psphotModelGroupInit (); - + // load command-line arguments and options config = psphotArguments (&argc, argv); // load input data (config and images (signal, noise, mask) - // XXX we have memory leaks here -- may be from psMetadata imdata = psphotSetup (config); - char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); - if (!status) breakPt = psStringCopy ("NONE"); - - int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE"); - if (!status) FITMODE = 2; - - // run model fitting tests on a single source - if (psMetadataLookupBool (&status, config, "TEST_FIT")) { - psphotModelTest (imdata, config); - exit (0); - } + // run a single-model test if desired + psphotModelTest (imdata, config); // measure image stats for initial guess @@ -42,31 +31,20 @@ // find the peaks in the image - peaks = pmPeaksSigmaLimit (imdata, config, sky); - if (!strcasecmp (breakPt, "PEAKS")) exit (0); + peaks = psphotFindPeaks (imdata, config, sky); // construct sources and measure basic stats sources = psphotSourceStats (imdata, config, peaks); - if (!strcasecmp (breakPt, "STATS")) exit (0); // classify sources based on moments, brightness - // XXX - this has 1 leak - psfClump = pmSourcePSFClump (sources, config); - - // group into STAR, COSMIC, GALAXY, SATURATED, etc. - pmSourceRoughClass (sources, config, psfClump); + psphotRoughClass (sources, config); // mark blended peaks PS_SOURCE_BLEND psphotBasicDeblend (sources, config, sky); - if (!strcasecmp (breakPt, "DEBLEND")) exit (0); - - // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, psphotSortBySN); - psphotDumpMoments (config, sources); - if (!strcasecmp (breakPt, "CLASS")) exit (0); // use bright stellar objects to measure PSF psf = psphotChoosePSF (config, sources, sky); - if (!strcasecmp (breakPt, "PSFFIT")) exit (0); + // XXX change FITMODE to a string + int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE"); switch (FITMODE) { case 0: @@ -101,5 +79,4 @@ // write out data in appropriate format - psphotSamplePSFs (config, psf, imdata->image); psphotOutput (imdata, config, sources, psf, sky); psLogMsg ("psphot", 3, "wrote output: %f sec\n", psTimerMark ("psphot")); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5985) +++ /trunk/psphot/src/psphot.h (revision 5986) @@ -25,5 +25,7 @@ eamReadout *psphotSetup (psMetadata *config); psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); +psArray *psphotFindPeaks (eamReadout *imdata, psMetadata *config, psStats *sky); psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); +bool psphotRoughClass (psArray *sources, psMetadata *config); pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); @@ -42,5 +44,4 @@ bool pmSourceFitFixed (pmSource *source, pmModel *model); -psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky); pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); pmModel *pmSourceSelectModel (pmSource *source); Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 5985) +++ /trunk/psphot/src/psphotArguments.c (revision 5986) @@ -117,5 +117,10 @@ psMetadata *config = psMetadataAlloc (); + // add default values + psMetadataAddStr (config, PS_LIST_TAIL, "BREAK_POINT", 0, "", "NONE"); + psMetadataAddS32 (config, PS_LIST_TAIL, "FIT_MODE", 0, "", 2); + psMetadataAdd (config, PS_LIST_TAIL, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); + config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE); Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 5985) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 5986) @@ -132,4 +132,8 @@ } psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f sec)\n", Nblend, psTimerMark ("psphot")); + + char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); + if (!strcasecmp (breakPt, "DEBLEND")) exit (0); + return true; } Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 5985) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 5986) @@ -12,4 +12,7 @@ psTimerStart ("psphot"); + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); // array to store candidate PSF stars @@ -88,4 +91,7 @@ psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); + if (!strcasecmp (breakPt, "PSFFIT")) exit (0); + return (psf); } Index: /trunk/psphot/src/psphotFindPeaks.c =================================================================== --- /trunk/psphot/src/psphotFindPeaks.c (revision 5986) +++ /trunk/psphot/src/psphotFindPeaks.c (revision 5986) @@ -0,0 +1,52 @@ +# include "psphot.h" + +psArray *psphotFindPeaks (eamReadout *imdata, psMetadata *config, psStats *sky) { + + bool status = false; + float NSIGMA; + float SIGMA; + float threshold; + + // smooth the image + + psTimerStart ("psphot"); + + SIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_SIGMA"); + NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA"); + + psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); + psImageSmooth (smooth, SIGMA, NSIGMA); + psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); + + psTimerStart ("psphot"); + + // set peak threshold + NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT"); + + threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + // threshold = NSIGMA*sky->sampleStdev; + psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold); + + // find the peaks in the smoothed image + psArray *peaks = pmFindImagePeaks (smooth, threshold); + if (peaks == NULL) psAbort ("find peaks", "no peaks found"); + psFree (smooth); + + // optional dump of all peak data + char *output = psMetadataLookupPtr (&status, config, "PEAKS_OUTPUT_FILE"); + if (status && (output != NULL) && (output[0])) { + pmPeaksWriteText (peaks, output); + psFree (output); + } + + psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); + + char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); + if (!strcasecmp (breakPt, "PEAKS")) exit (0); + + return (peaks); +} + +// In this function, we smooth the image, then search for the peaks +// Should we also subtract a super-binned image? (as an option?) +// We need to gracefully handle no source detections Index: /trunk/psphot/src/psphotGrowthCurve.c =================================================================== --- /trunk/psphot/src/psphotGrowthCurve.c (revision 5986) +++ /trunk/psphot/src/psphotGrowthCurve.c (revision 5986) @@ -0,0 +1,65 @@ +# include "psphot.h" + +// XXX we can probably move this into pmGrowthCurve.c +// XXX need to change the way we grab an image, or else use +// the 'center' option + +bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf) { + + // bool status; + float xc, yc, dx, dy; + float fitMag, apMag; + float radius; + + // create template model + pmModel *modelRef = pmModelAlloc(psf->type); + + // use this image region + xc = 0.5*imdata->image->numCols; + yc = 0.5*imdata->image->numRows; + dx = psf->growth->maxRadius + 1; + dy = psf->growth->maxRadius + 1; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[0] = 0; + modelRef->params->data.F32[1] = 1000; + modelRef->params->data.F32[2] = xc; + modelRef->params->data.F32[3] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + + // generate working image for this source + psRegion region = {xc - dx, xc + dx, yc - dy, yc + dy}; + psImage *view = psImageSubset (imdata->image, region); + psImage *image = psImageCopy (NULL, view, PS_TYPE_F32); + psImage *mask = psImageCopy (NULL, view, PS_TYPE_U8); + psImageInit (mask, 0); + + // loop over a range of source fluxes + for (int i = 0; i < psf->growth->radius->n; i++) { + + psImageInit (image, 0.0); + + radius = psf->growth->radius->data.F32[i]; + + psImageKeepCircle (mask, xc, yc, radius, "OR", PSPHOT_MASK_MARKED); + + pmSourceAddModel (image, mask, model, false, false); + + pmSourcePhotometry (&fitMag, &apMag, model, image, mask); + + psImageKeepCircle (mask, xc, yc, radius, "AND", ~PSPHOT_MASK_MARKED); + + psf->growth->apMag->data.F32[i] = apMag; + } + psf->growth->fitMag = fitMag; + + psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius); + psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef; + + psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f\n", psf->growth->apLoss); + + return true; +} Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 5985) +++ /trunk/psphot/src/psphotModelTest.c (revision 5986) @@ -11,4 +11,7 @@ psMetadataItem *item = NULL; + + // run model fitting tests on a single source + if (!psMetadataLookupBool (&status, config, "TEST_FIT")) return false; // what fitting mode to use? @@ -166,4 +169,5 @@ psphotSaveImage (NULL, source->pixels, "resid.fits"); psphotSaveImage (NULL, source->mask, "mask.fits"); - return true; + + exit (0); } Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5985) +++ /trunk/psphot/src/psphotOutput.c (revision 5986) @@ -15,4 +15,5 @@ if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage); if (psfFile != NULL) psphotWritePSF (psf, psfFile); + psphotSamplePSFs (config, psf, imdata->image); if (outputFile == NULL) { Index: /trunk/psphot/src/psphotRoughClass.c =================================================================== --- /trunk/psphot/src/psphotRoughClass.c (revision 5986) +++ /trunk/psphot/src/psphotRoughClass.c (revision 5986) @@ -0,0 +1,21 @@ +# include "psphot.h" + +bool psphotRoughClass (psArray *sources, psMetadata *config) { + + bool status; + pmPSFClump psfClump; + + // XXX - this has 1 leak + psfClump = pmSourcePSFClump (sources, config); + + // group into STAR, COSMIC, GALAXY, SATURATED, etc. + pmSourceRoughClass (sources, config, psfClump); + + // optional printout of source moments only + psphotDumpMoments (config, sources); + + char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); + if (!strcasecmp (breakPt, "CLASS")) exit (0); + + return true; +} Index: /trunk/psphot/src/psphotSetup.c =================================================================== --- /trunk/psphot/src/psphotSetup.c (revision 5985) +++ /trunk/psphot/src/psphotSetup.c (revision 5986) @@ -17,4 +17,7 @@ psTimerStart ("psphot"); + + // load implementation-specific models + psphotModelGroupInit (); // setup header template, specify COMMENT and HISTORY as MULTI Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 5985) +++ /trunk/psphot/src/psphotSourceStats.c (revision 5986) @@ -61,4 +61,8 @@ psLogMsg ("psphot", 3, "%d moments: %f sec\n", sources->n, psTimerMark ("psphot")); + + char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); + if (!strcasecmp (breakPt, "STATS")) exit (0); + return (sources); }