Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5057) +++ /trunk/psphot/Makefile (revision 5058) @@ -3,6 +3,5 @@ @echo "USAGE: make psphot" -# CC = gcc -g -Wall -O2 -CC = gcc -g -Wall +CC = gcc SRC = src BIN = bin @@ -16,5 +15,5 @@ # LIBS = -lpsmodule $(LPSLIB) LIBS = $(LPSLIB) -CFLAGS = $(INCS) -std=c99 -Wall -Werror +CFLAGS = $(INCS) -std=c99 -Wall -Werror -g LFLAGS = $(LIBS) Index: /trunk/psphot/src/pmPSF.c =================================================================== --- /trunk/psphot/src/pmPSF.c (revision 5057) +++ /trunk/psphot/src/pmPSF.c (revision 5058) @@ -24,6 +24,9 @@ pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF)); - psf->type = type; - psf->chisq = 0.0; + psf->type = type; + psf->chisq = 0.0; + psf->ApResid = 0.0; + psf->dApResid = 0.0; + psf->skyBias = 0.0; Nparams = pmModelParameterCount (type); Index: /trunk/psphot/src/pmPSF.h =================================================================== --- /trunk/psphot/src/pmPSF.h (revision 5057) +++ /trunk/psphot/src/pmPSF.h (revision 5058) @@ -13,4 +13,7 @@ psArray *params; ///< Model parameters (psPolynomial2D) float chisq; ///< PSF goodness statistic + float ApResid; + float dApResid; + float skyBias; int nPSFstars; ///< number of stars used to measure PSF } Index: /trunk/psphot/src/pmPSFtry.c =================================================================== --- /trunk/psphot/src/pmPSFtry.c (revision 5057) +++ /trunk/psphot/src/pmPSFtry.c (revision 5058) @@ -31,18 +31,17 @@ pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) { + pmModelType type; + pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); // XXX probably need to increment ref counter - test->modelType = pmModelSetType (modelName); - test->psf = pmPSFAlloc (test->modelType); - test->sources = psMemCopy(sources); - test->modelFLT = psArrayAlloc (sources->n); - test->modelPSF = psArrayAlloc (sources->n); - test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); - test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F64); - test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); - test->ApResid = 0; - test->dApResid = 0; - test->skyBias = 0; + type = pmModelSetType (modelName); + test->psf = pmPSFAlloc (type); + test->sources = psMemCopy(sources); + test->modelFLT = psArrayAlloc (sources->n); + test->modelPSF = psArrayAlloc (sources->n); + test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); + test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F64); + test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); for (int i = 0; i < test->modelFLT->n; i++) { @@ -84,5 +83,5 @@ pmSource *source = try->sources->data[i]; - pmModel *model = pmSourceModelGuess (source, try->modelType); + pmModel *model = pmSourceModelGuess (source, try->psf->type); x = source->peak->x; y = source->peak->y; @@ -164,5 +163,5 @@ pmPSFtryMetric (try, RADIUS); psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", - modelName, try->ApResid, try->dApResid, try->skyBias); + modelName, try->psf->ApResid, try->psf->dApResid, try->psf->skyBias); return (try); @@ -279,7 +278,7 @@ stats = psVectorStats (stats, daResid, NULL, maskB, 1); - try->ApResid = poly->coeff[0]; - try->dApResid = stats->clippedStdev; - try->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + try->psf->ApResid = poly->coeff[0]; + try->psf->dApResid = stats->clippedStdev; + try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); psFree (rflux); Index: /trunk/psphot/src/pmPSFtry.h =================================================================== --- /trunk/psphot/src/pmPSFtry.h (revision 5057) +++ /trunk/psphot/src/pmPSFtry.h (revision 5058) @@ -5,5 +5,4 @@ // data to try a given PSF model type typedef struct { - pmModelType modelType; pmPSF *psf; psArray *sources; // pointers to the original sources @@ -13,7 +12,4 @@ psVector *metric; psVector *fitMag; - float ApResid; - float dApResid; - float skyBias; } pmPSFtry; Index: /trunk/psphot/src/psLibUtils.h =================================================================== --- /trunk/psphot/src/psLibUtils.h (revision 5057) +++ /trunk/psphot/src/psLibUtils.h (revision 5058) @@ -15,13 +15,5 @@ # define psMemCopy(A)(psMemIncrRefCounter((A))) -// XXX EAM : bug : missing from psLib/*.h -bool psMinimizeGaussNewtonDelta (psVector *delta, - const psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yErr, - psMinimizeLMChi2Func func); - +// XXX EAM : my version using varience instead of stdev bool psMinimizeGaussNewtonDelta_EAM (psVector *delta, const psVector *params, @@ -32,5 +24,5 @@ psMinimizeLMChi2Func func); -// minimize +// minimize : using varience vs sigma and parameter limits psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, psVector *Beta, psVector *Params, const psImage *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda); psBool psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func); Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5057) +++ /trunk/psphot/src/psphot.c (revision 5058) @@ -47,5 +47,5 @@ // write out data in appropriate format - psphotOutput (imdata, config, sources); + psphotOutput (imdata, config, sources, psf, sky); exit (0); } Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5057) +++ /trunk/psphot/src/psphot.h (revision 5058) @@ -26,5 +26,5 @@ bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); -void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources); +void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); @@ -38,9 +38,9 @@ // output functions -bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources); -bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources); -bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources); -bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources); -bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources); +bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); int pmSourcesDophotType (pmSource *source); bool pmPeaksWriteText (psArray *sources, char *filename); Index: /trunk/psphot/src/psphotApplyPSF.c =================================================================== --- /trunk/psphot/src/psphotApplyPSF.c (revision 5057) +++ /trunk/psphot/src/psphotApplyPSF.c (revision 5058) @@ -28,5 +28,5 @@ // set the object surface-brightness limit for fitted pixels float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; - psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); // this function specifies the radius at this the model hits the given flux @@ -72,5 +72,5 @@ psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); if (!status || (model->params->data.F32[1] < 0)) { - psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); + psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? psFree (model); @@ -91,5 +91,5 @@ psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); - psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub); + psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub); return (true); } Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 5057) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 5058) @@ -12,9 +12,11 @@ // array to store candidate PSF stars + int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS"); + if (!status) NSTARS = sources->n; stars = psArrayAlloc (sources->n); stars->n = 0; // select the candidate PSF stars (pointers to original sources) - for (int i = 0; i < sources->n; i++) { + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { pmSource *source = sources->data[i]; if (source->type != PM_SOURCE_PSFSTAR) continue; @@ -55,8 +57,8 @@ try = models->data[0]; int bestN = 0; - float bestM = try->dApResid; + float bestM = try->psf->dApResid; for (int i = 1; i < models->n; i++) { try = models->data[i]; - float M = try->dApResid; + float M = try->psf->dApResid; if (M < bestM) { bestM = M; @@ -67,10 +69,10 @@ // keep only the selected model: try = models->data[bestN]; - modelName = pmModelGetType (try->modelType); - psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, try->ApResid, try->dApResid); + pmPSF *psf = psMemCopy(try->psf); + psFree (models); // keep only the pmPSF resulting from this analysis - // keep only the pmPSF resulting from this analysis - pmPSF *psf = psMemCopy(try->psf); - psFree (models); + modelName = pmModelGetType (psf->type); + psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + return (psf); } Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 5057) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 5058) @@ -81,5 +81,5 @@ if (!status) { // if the fit fails, we need to change the classification - psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); + psLogMsg ("psphot", 5, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); source->type = PM_SOURCE_FAIL_FIT_GAL; // better choice? source->modelFLT = model; @@ -92,5 +92,5 @@ if (!goodfit) { // if the fit fails, we need to change the classification - psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); + psLogMsg ("psphot", 5, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); source->type = PM_SOURCE_POOR_FIT_GAL; // better choice? source->modelFLT = model; Index: /trunk/psphot/src/psphotImageStats.c =================================================================== --- /trunk/psphot/src/psphotImageStats.c (revision 5057) +++ /trunk/psphot/src/psphotImageStats.c (revision 5058) @@ -49,5 +49,5 @@ sky->sampleStdev = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE)); - psLogMsg ("psphot", 3, "stats: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", 4, "stats: %f sec\n", psTimerMark ("psphot")); psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev); Index: /trunk/psphot/src/psphotMarkPSF.c =================================================================== --- /trunk/psphot/src/psphotMarkPSF.c (revision 5057) +++ /trunk/psphot/src/psphotMarkPSF.c (revision 5058) @@ -36,5 +36,5 @@ if (source->modelPSF->params->data.F32[1] >= SATURATE) { if (source->type == PM_SOURCE_PSFSTAR) { - psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n"); + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); } source->type = PM_SOURCE_SATSTAR; @@ -42,5 +42,5 @@ } if (source->type == PM_SOURCE_SATSTAR) { - psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); + psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); source->type = PM_SOURCE_GOODSTAR; } @@ -70,5 +70,5 @@ if (source->type == PM_SOURCE_PSFSTAR) { - psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n"); + psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n"); } Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5057) +++ /trunk/psphot/src/psphotOutput.c (revision 5058) @@ -2,5 +2,5 @@ // output functions: we have several fixed modes (sx, obj, cmp) -void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources) { +void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { bool status; @@ -21,25 +21,25 @@ } if (!strcasecmp (outputMode, "TEXT")) { - pmSourcesWriteText (imdata, outputFile, sources); + pmSourcesWriteText (imdata, outputFile, sources, psf, sky); return; } if (!strcasecmp (outputMode, "OBJ")) { - pmSourcesWriteOBJ (imdata, outputFile, sources); + pmSourcesWriteOBJ (imdata, outputFile, sources, psf, sky); return; } if (!strcasecmp (outputMode, "SX")) { - pmSourcesWriteSX (imdata, outputFile, sources); + pmSourcesWriteSX (imdata, outputFile, sources, psf, sky); return; } if (!strcasecmp (outputMode, "CMP")) { - pmSourcesWriteCMP (imdata, outputFile, sources); + pmSourcesWriteCMP (imdata, outputFile, sources, psf, sky); return; } if (!strcasecmp (outputMode, "CMF")) { - pmSourcesWriteCMF (imdata, outputFile, sources); + pmSourcesWriteCMF (imdata, outputFile, sources, psf, sky); return; } @@ -47,5 +47,5 @@ } -bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources) { +bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky) { char *name = (char *) psAlloc (strlen(filename) + 10); @@ -68,10 +68,10 @@ // dophot-style output list with fixed line width -bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources) { +bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { int i, type, status; pmModel *model; psF32 *PAR, *dPAR; - float sky, dmag, apMag, fitMag; + float skyFlux, dmag, apMag, fitMag; psLine *line = psLineAlloc (104); // 104 is dophot-defined line length @@ -118,5 +118,5 @@ // save local sky for later - sky = model->params->data.F32[0]; + skyFlux = model->params->data.F32[0]; // replace source flux @@ -124,8 +124,9 @@ model->params->data.F32[0] = 0; pmSourceAddModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; + model->params->data.F32[0] = skyFlux; // measure object photometry status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + fitMag += psf->ApResid; // subtract object, leave local sky @@ -133,5 +134,5 @@ model->params->data.F32[0] = 0; pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; + model->params->data.F32[0] = skyFlux; if (status == FALSE) continue; @@ -159,5 +160,5 @@ // elixir/sextractor-style output list with fixed line width -bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources) { +bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { int i, status, flags; @@ -224,4 +225,5 @@ // measure object photometry status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + fitMag += psf->ApResid; // subtract object, leave local sky @@ -257,5 +259,5 @@ // elixir-style pseudo FITS table (header + ascii list) -bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources) { +bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { int i, type, status; @@ -328,4 +330,5 @@ // measure object photometry status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + fitMag += psf->ApResid; // subtract object, leave local sky @@ -360,5 +363,5 @@ // elixir-style FITS table output (header + table in 1st extension) -bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources) { +bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { // write imdata->header to file Index: /trunk/psphot/src/psphotSetup.c =================================================================== --- /trunk/psphot/src/psphotSetup.c (revision 5057) +++ /trunk/psphot/src/psphotSetup.c (revision 5058) @@ -28,6 +28,8 @@ psFits *file = psFitsAlloc (input); header = psFitsReadHeader (header, file); - image = psFitsReadImage (NULL, file, region, 0); + psImage *tmpimage = psFitsReadImage (NULL, file, region, 0); + image = psImageCopy (NULL, tmpimage, PS_TYPE_F32); psFree (file); + // psFree (input); // XXX EAM - does the current API expect to get a copy or a view? @@ -90,5 +92,5 @@ } - psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", 4, "load data: %f sec\n", psTimerMark ("psphot")); // save the image data & return it Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 5057) +++ /trunk/psphot/src/psphotSourceStats.c (revision 5058) @@ -47,7 +47,5 @@ } - psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot")); - psLogMsg ("psphot", 3, "%d sources\n", sources->n); - + psLogMsg ("psphot", 3, "%d moments: %f sec\n", sources->n, psTimerMark ("psphot")); return (sources); }