Index: /trunk/psphot/doc/notes.txt =================================================================== --- /trunk/psphot/doc/notes.txt (revision 4115) +++ /trunk/psphot/doc/notes.txt (revision 4115) @@ -0,0 +1,166 @@ + +Notes on psphot + +2005.06.04 + + progress on psphot has moved along well. At this point, the process + loads the image, finds the peaks, determines a best PSF model, + identifies objects which are consistent with that model, and + attempts to fit a galaxy model to the objects which are not + consistent with the model. + + a variety of minor issues remain, as well as some major issues. + + - use a proper noise image to keep the fits honest after other + objects have been subtracted. + + - define the masks in terms of a single mask image. this image + could be provided initially by the user. + + - add noise enhancement (couple with mask marking subtraction) + + - + +2005.04.12 + +- psPeak, psSource, etc: should be pmPeak, pmSource, etc +- pmCullPeaks: should be pmCullImagePeaks + +- does pmSourceMoments need an image argument? + +- psLibInit / p_psTimeInit + + we should not have to know about astronomy time concepts (ie, have a + timeConfig file) in order to use the time functions for basic work. + This means having a PS_TIME_UNIX which just works with the UNIX + clock and avoids all initialization issues. + +- psImageStats is much too slow for basic median + - this function does a complete sort on all pixel values. we need + to use the histogram method (spec a different type of median) + +- psImageStats, ROBUST_MEDIAN fails +- there are errors in psVectorStats for ROBUST MEAN, STDEV (& median?) + * this is ill-defined. revise the ADD: too sensitive to the + binning, smoothing scales + +- pmFindImagePeaks is not finding any peaks: + typo in interior row section: see bug 359 + after this fix: PASS + +* p_psGetRowVectorFromImage : PASS +* psGetRowVectorFromImage : PASS +- both have about the same time (~15 us +/- 5us) for 2048 pix on alala + +* pmFindVectorPeaks : PASS (returned row data in data.U32) + +- # include : needed in psMetadata.c + +* MyListAddPeak / pmPeakAlloc are inconsistent wrt col,row vs x,y + +- why does psMetadataItemAllocV allocate a string of length + MAX_STRING_LENGTH for the comment, then use strncpy? + + should be: + + // set metadata item comment + if (comment == NULL) { + // Per SDRS, null isn't allowed, must use "" instead + metadataItem->comment = psStringCopy (""); + } else { + metadataItem->comment = psStringCopy (comment); + } + +psMetadataAdd -> psMetadataAddItem is all confused: + - we end up with two psMetadataItemAlloc calls, when only one should + be used + - psMetadataItemAlloc was doing the wrong thing if the incoming data + was NULL + - I think the logic in this block is wrong as well. + + * I fixed this in psMetadata.c + +2005.04.14 : + + - psImageStats is still much too slow + + - pmObjects : changed to using psArray to carry peaks / sources, + rather than psList + + - added function pmPeaksSubset to replace pmCullPeaks + + - pmSourceLocalSky: added min/max fncs to force subimage to stay in + image + + - Sx, Sy can go negative: have forced limit to 0,0 + + - pmSourceRoughClass: implemented sigmaX, sigmaY search for stellar + clump (uses pmComparePeakDescend) + + - the search ignores the 0,0 pixel + + - pmSourceRoughClass uses all detected peaks. allow an exclusion + for peak > max, peak < min? + + - need to push clump stats on the metadata + + - the source classes are assigned so they are mutually exclusive + + - change PS_SOURCE_OTHER to STAR? + + - pmSourceSetPixelCircle: + + - name should be pmSourceSetPixelsCircle (geometry could be an + option...) + + - I cleaned the defined region to match the convensions of + pmSourceLocalSky + + - truncate and force subimage to lie on image + + - CheckRadius2 seems like an inefficient coding + + - valid pixels have mask value of 1 (inverse of other funcs) + + - pmSourceFitModel: + + - num iterations is much too large + + - + +mrq: + + - allocate alpha, beta arrays + - mrq2dcof (pars, alpha, beta) + - make guess pars', alpha', beta' + - mrq2dcof (pars', alpha', beta') + - if (chisq < old chisq), keep new guess + + +2005.04.19: + +- psMinimizeLMChi2 has some serious problems. I re-wrote it starting + from the mrqmin example in ohana. after a variety little bugs, it + seems to work quite well. the program 'modeltest.c' runs two + different gaussians through psMinimizeLMChi2, and gets the right + answer quite quickly. The first pass at the implementation had a + variety of problems. to avoid any worry about errors in the + psMatrix code, I implemented psGaussJordan, basing the code on the + elixir gaussj.c code. this worked fine, and let me find the errors + elsewhere in the code. the LUD version now works just the same as + the GaussJ version. + + +2005.04.21 + +- psMinimizeLMChi2 does an excellent job now in most cases. I am able + to run pmSourceFitModel on the PSF stars quite well. I added in the + (nearly) correct errors (actually, just using sqrt(N), limiting to + N>=1). I see one problem object, which is quite faint, and does not + converge well: I get very large values for x,y, Sx, Sy, etc. I need + to think about constraining these fits from running way off scale. + + fit: 0.123879 sec for 277 stars (gauss) + fit: 0.155940 sec for 277 stars (pgauss) (!) + fit: 0.153292 sec for 277 stars + (0.5 msec per star) Index: /trunk/psphot/doc/psphot.txt =================================================================== --- /trunk/psphot/doc/psphot.txt (revision 4114) +++ /trunk/psphot/doc/psphot.txt (revision 4115) @@ -1,144 +1,5 @@ -2005.04.12 : Notes on psphot +psphot is a stand-alone program which implements the Pan-STARRS IPP +object detection and analysis functions for using with individual +images. -- psPeak, psSource, etc: should be pmPeak, pmSource, etc -- pmCullPeaks: should be pmCullImagePeaks - -- does pmSourceMoments need an image argument? - -- psLibInit / p_psTimeInit - - we should not have to know about astronomy time concepts (ie, have a - timeConfig file) in order to use the time functions for basic work. - This means having a PS_TIME_UNIX which just works with the UNIX - clock and avoids all initialization issues. - -- psImageStats is much too slow for basic median - - this function does a complete sort on all pixel values. we need - to use the histogram method (spec a different type of median) - -- psImageStats, ROBUST_MEDIAN fails -- there are errors in psVectorStats for ROBUST MEAN, STDEV (& median?) - * this is ill-defined. revise the ADD: too sensitive to the - binning, smoothing scales - -- pmFindImagePeaks is not finding any peaks: - typo in interior row section: see bug 359 - after this fix: PASS - -* p_psGetRowVectorFromImage : PASS -* psGetRowVectorFromImage : PASS -- both have about the same time (~15 us +/- 5us) for 2048 pix on alala - -* pmFindVectorPeaks : PASS (returned row data in data.U32) - -- # include : needed in psMetadata.c - -* MyListAddPeak / pmPeakAlloc are inconsistent wrt col,row vs x,y - -- why does psMetadataItemAllocV allocate a string of length - MAX_STRING_LENGTH for the comment, then use strncpy? - - should be: - - // set metadata item comment - if (comment == NULL) { - // Per SDRS, null isn't allowed, must use "" instead - metadataItem->comment = psStringCopy (""); - } else { - metadataItem->comment = psStringCopy (comment); - } - -psMetadataAdd -> psMetadataAddItem is all confused: - - we end up with two psMetadataItemAlloc calls, when only one should - be used - - psMetadataItemAlloc was doing the wrong thing if the incoming data - was NULL - - I think the logic in this block is wrong as well. - - * I fixed this in psMetadata.c - -2005.04.14 : - - - psImageStats is still much too slow - - - pmObjects : changed to using psArray to carry peaks / sources, - rather than psList - - - added function pmPeaksSubset to replace pmCullPeaks - - - pmSourceLocalSky: added min/max fncs to force subimage to stay in - image - - - Sx, Sy can go negative: have forced limit to 0,0 - - - pmSourceRoughClass: implemented sigmaX, sigmaY search for stellar - clump (uses pmComparePeakDescend) - - - the search ignores the 0,0 pixel - - - pmSourceRoughClass uses all detected peaks. allow an exclusion - for peak > max, peak < min? - - - need to push clump stats on the metadata - - - the source classes are assigned so they are mutually exclusive - - - change PS_SOURCE_OTHER to STAR? - - - pmSourceSetPixelCircle: - - - name should be pmSourceSetPixelsCircle (geometry could be an - option...) - - - I cleaned the defined region to match the convensions of - pmSourceLocalSky - - - truncate and force subimage to lie on image - - - CheckRadius2 seems like an inefficient coding - - - valid pixels have mask value of 1 (inverse of other funcs) - - - pmSourceFitModel: - - - num iterations is much too large - - - - -mrq: - - - allocate alpha, beta arrays - - mrq2dcof (pars, alpha, beta) - - make guess pars', alpha', beta' - - mrq2dcof (pars', alpha', beta') - - if (chisq < old chisq), keep new guess - - -2005.04.19: - -- psMinimizeLMChi2 has some serious problems. I re-wrote it starting - from the mrqmin example in ohana. after a variety little bugs, it - seems to work quite well. the program 'modeltest.c' runs two - different gaussians through psMinimizeLMChi2, and gets the right - answer quite quickly. The first pass at the implementation had a - variety of problems. to avoid any worry about errors in the - psMatrix code, I implemented psGaussJordan, basing the code on the - elixir gaussj.c code. this worked fine, and let me find the errors - elsewhere in the code. the LUD version now works just the same as - the GaussJ version. - - -2005.04.21 - -- psMinimizeLMChi2 does an excellent job now in most cases. I am able - to run pmSourceFitModel on the PSF stars quite well. I added in the - (nearly) correct errors (actually, just using sqrt(N), limiting to - N>=1). I see one problem object, which is quite faint, and does not - converge well: I get very large values for x,y, Sx, Sy, etc. I need - to think about constraining these fits from running way off scale. - - fit: 0.123879 sec for 277 stars (gauss) - fit: 0.155940 sec for 277 stars (pgauss) (!) - fit: 0.153292 sec for 277 stars - (0.5 msec per star) Index: /trunk/psphot/src/apply_psf_model.c =================================================================== --- /trunk/psphot/src/apply_psf_model.c (revision 4114) +++ /trunk/psphot/src/apply_psf_model.c (revision 4115) @@ -22,5 +22,8 @@ keep = psRegionForImage (keep, image, keep); - // set the object surfact-brightness limit for fitted pixels + float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + float snFaint = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM"); + + // set the object surface-brightness limit for fitted pixels float FLUX_LIMIT = FIT_NSIGMA * sky->sampleStdev; psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); @@ -31,4 +34,7 @@ for (int i = 0; i < sources->n; i++) { psSource *source = sources->data[i]; + + // skip the existing PSF stars + // XXX drop this -- refit? if (source->modelPSF != NULL) continue; @@ -55,9 +61,12 @@ pmSourceMaskSaturated (source, SATURATE); - // fit as PSF, not FLT (drop poor fits) + // fit as PSF, not FLT (skip poor fits) if (!pmSourceFitModel (source, model, true)) continue; source->modelPSF = model; Niter += model[0].nIter; Nfit ++; + + mark_psf_source (source); + subtract_psf_source (source); } psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); Index: /trunk/psphot/src/by_SN.c =================================================================== --- /trunk/psphot/src/by_SN.c (revision 4115) +++ /trunk/psphot/src/by_SN.c (revision 4115) @@ -0,0 +1,16 @@ +# include "psphot.h" + +int by_SN (const void **a, const void **b) +{ + psSource *A = *a; + psSource *B = *b; + + psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN; + psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN; + + psF32 diff = fA - fB; + if (diff < FLT_EPSILON) return (-1); + if (diff > FLT_EPSILON) return (+1); + return (0); +} + Index: /trunk/psphot/src/choose_psf_model.c =================================================================== --- /trunk/psphot/src/choose_psf_model.c (revision 4114) +++ /trunk/psphot/src/choose_psf_model.c (revision 4115) @@ -68,5 +68,11 @@ for (int i = 0; i < test->sources->n; i++) { psSource *source = test->sources->data[i]; - source->modelPSF = test->modelPSF->data[i]; + // drop masked sources from PSFSTAR list + if (test->mask->data.U8[i]) { + source->type = PS_SOURCE_OTHER; + // XXX is this the right type to go to? + } else { + source->modelPSF = test->modelPSF->data[i]; + } } return (test->psf); Index: /trunk/psphot/src/find_peaks.c =================================================================== --- /trunk/psphot/src/find_peaks.c (revision 4114) +++ /trunk/psphot/src/find_peaks.c (revision 4115) @@ -9,5 +9,5 @@ // smooth the image - // should we also subtract a super-binned image? (as an option) + // should we also subtract a super-binned image? (as an option?) psTimerStart ("psphot"); @@ -24,5 +24,4 @@ psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); - // find the peaks in the smoothed image Index: /trunk/psphot/src/image_stats.c =================================================================== --- /trunk/psphot/src/image_stats.c (revision 4114) +++ /trunk/psphot/src/image_stats.c (revision 4115) @@ -13,22 +13,20 @@ // pass on the stats { - bool status = false; + bool status = false; + int Npixels = image->numRows*image->numCols; + int Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "NSUBSET")); + psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32); psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); - int Nsubset = psMetadataLookupF32 (&status, config, "NSUBSET"); - double fSubset = (double) Nsubset / (image->numRows*image->numCols); - psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32); - subset->n = 0; - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - double frnd = psRandomUniform (rnd); - if ((fSubset < 1.0) && (fSubset < frnd)) continue; - subset->data.F32[subset->n] = image->data.F32[j][i]; - subset->n ++; - if (subset->n == Nsubset) goto got_subset; - } + // choose Nsubset points between 0 and Nx*Ny, convert to coords + Npixels = image->numRows*image->numCols; + for (int i = 0; i < Nsubset; i++) { + double frnd = psRandomUniform (rnd); + int pixel = Npixels * frnd; + int ix = pixel / image->numCols; + int iy = pixel % image->numCols; + subset->data.F32[i] = image->data.F32[iy][ix]; } - got_subset: // this should use ROBUST not SAMPLE median // robust median is broken in pslib (0.5) Index: /trunk/psphot/src/mark_psf_source.c =================================================================== --- /trunk/psphot/src/mark_psf_source.c (revision 4115) +++ /trunk/psphot/src/mark_psf_source.c (revision 4115) @@ -0,0 +1,53 @@ +# include "psphot.h" + +// identify objects consistent with PSF shape/magnitude distribution +// we expect dparams[4],dparams[5] to have a scatter of: +// sigma_x / (S/N) * sqrt(2) +// 1 / (params[4],params[5])*(S/N) + +// any objects which is consistent with the PSF should have +// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) +// this should include a minimum dsx buffer for the brighter objects +// saturated stars should fall outside (but are already IDed) +// galaxies should be larger, cosmic rays smaller, but need to test? + +// any object which meets the criterion is marked as +// PS_SOURCE_BRIGHTSTAR +bool mark_psf_source (psArray *sources) +{ + bool status = false; + float dSX, dSY, SX, SY, SN; + float nSx, nSy; + + if (source->modelPSF == NULL) return (false); + + if (source->type == PS_SOURCE_PSFSTAR) return (true); + if (source->type == PS_SOURCE_SATSTAR) return (true); + + SN = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]; + SX = source->modelPSF->params->data.F32[4]; + SY = source->modelPSF->params->data.F32[5]; + dSX = source->modelPSF->dparams->data.F32[4]; + dSY = source->modelPSF->dparams->data.F32[5]; + + nSx = dSX * SX * SN; + nSy = dSY * SY * SN; + + // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[] + if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) { + source->type = PS_SOURCE_BRIGHTSTAR; + return (true); + } + + if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { + source->type = PS_SOURCE_GALAXY; + return (false); + } + // replace DEFECT with COSMIC? + if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { + source->type = PS_SOURCE_DEFECT; + return (false); + } + psTrace (".psphot.mark_psf_source", 2, "unexpected result: unmarked object\n"); + return (false); +} Index: /trunk/psphot/src/psPolynomials.c =================================================================== --- /trunk/psphot/src/psPolynomials.c (revision 4114) +++ /trunk/psphot/src/psPolynomials.c (revision 4115) @@ -191,7 +191,8 @@ // XXX EAM : test version of 1d fitting psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, - const psVector* x, - const psVector* y, - const psVector* yErr) + psVector *mask, + const psVector *x, + const psVector *y, + const psVector *yErr) { // I think this is 1 dimension down @@ -232,4 +233,5 @@ // Build the B and A data structs. for (int k = 0; k < x->n; k++) { + if ((mask != NULL) && mask->U8[k]) continue; xSums = BuildSums1D(xSums, x->data.F64[k], nTerm); @@ -291,8 +293,9 @@ // XXX EAM : test version of 2d fitting psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, - const psVector* x, - const psVector* y, - const psVector* z, - const psVector* zErr) + psVector* mask, + const psVector* x, + const psVector* y, + const psVector* z, + const psVector* zErr) { // I think this is 1 dimension down @@ -322,4 +325,5 @@ // Build the B and A data structs. for (int k = 0; k < x->n; k++) { + if ((mask != NULL) && mask->data.U8[i]) continue; Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm); @@ -391,5 +395,5 @@ } -psPolynomial2D* RobustFit2D(psPolynomial2D* poly, +psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, @@ -413,5 +417,5 @@ for (int N = 0; N < 3; N++) { // XXX EAM : this would be better defined with an element mask - poly = VectorFitPolynomial2DOrd_EAM (poly, X, Y, Z, dZ); + poly = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ); zFit = Polynomial2DEvalVectorD (poly, x, y); zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit); @@ -424,5 +428,7 @@ int n = 0; for (int i = 0; i < zResid->n; i++) { - if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) continue; + if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) { + continue; + } X->data.F64[n] = x->data.F64[i]; Y->data.F64[n] = y->data.F64[i]; @@ -435,4 +441,41 @@ Z->n = n; dZ->n = n; + } + return (poly); +} + +// XXX EAM : be careful here with F32 vs F64 vectors +psPolynomial2D* RobustFit2D(psPolynomial2D* poly, + psVector* mask, + const psVector* x, + const psVector* y, + const psVector* z, + const psVector* dz) +{ + psVector *zFit = NULL; + psVector *zResid = NULL; + psStats *stats = NULL; + + for (int N = 0; N < 3; N++) { + poly = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz); + zFit = Polynomial2DEvalVectorD (poly, x, y); + zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit); + + stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + stats = psVectorStats (stats, zResid, NULL, mask, 1); + psTrace (".psphot.RobustFit", 4, "residual stats for robust fit: %g +/- %g (%d pts)\n", + stats->clippedMean, stats->clippedStdev, stats->clippedNvalues); + + // set mask if pts are not valid + int n = 0; + for (int i = 0; i < zResid->n; i++) { + if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) { + mask->data.U8[i] = 1; + continue; + } else { + mask->data.U8[i] = 0; + } + n++; + } } return (poly); Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 4114) +++ /trunk/psphot/src/psphot.c (revision 4115) @@ -24,5 +24,9 @@ peaks = find_peaks (image, config, sky); + // construct sources and measure basic stats sources = source_moments (image, config, peaks); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, by_SN); // use stellar objects SN > PSF_SN_LIM @@ -35,12 +39,4 @@ // XXX it has some problems // test_psf_scatter (sources); - - // identify PSF-like objects (< PSF_SHAPE_NSIGMA from dSx,dSy model) - // SN <= FAINT_SN_LIM : FAINTSTAR - // SN > FAINT_SN_LIM : BRIGHTSTAR - mark_psf_sources (sources, config); - - // subtract BRIGHTSTAR, PSFSTAR, SATSTAR - subtract_psf_sources (sources); // fit_galaxies (image, config, sources); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 4114) +++ /trunk/psphot/src/psphot.h (revision 4115) @@ -13,4 +13,5 @@ psArray *modelFLT; // model fits, floating parameters psArray *modelPSF; // model fits, PSF parameters + psVector *mask; psVector *metric; psStats *metricStats; Index: /trunk/psphot/src/pspsf.c =================================================================== --- /trunk/psphot/src/pspsf.c (revision 4114) +++ /trunk/psphot/src/pspsf.c (revision 4115) @@ -36,8 +36,10 @@ test->metricStats = NULL; test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); + test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); for (int i = 0; i < test->modelFLT->n; i++) { test->modelFLT->data[i] = NULL; test->modelPSF->data[i] = NULL; + test->mask->data.U8[i] = 0; } return (test); @@ -57,6 +59,9 @@ psModel *model = pmSourceModelGuess (source, test->modelType); - // fit model as FLT, not PSF (skip poor fits) - if (!pmSourceFitModel (source, model, false)) continue; + // fit model as FLT, not PSF (mask & skip poor fits) + if (!pmSourceFitModel (source, model, false)) { + test->mask->data.U8[i] = 1; + continue; + } test->modelFLT->data[i] = model; Nflt ++; @@ -66,5 +71,5 @@ // stage 2: construct a psf (pmPSF) from this collection of model fits - pmPSFFromModels (test->psf, test->modelFLT); + pmPSFFromModels (test->psf, test->modelFLT, test->mask); // stage 3: refit with fixed shape parameters @@ -73,9 +78,18 @@ psSource *source = test->sources->data[i]; psModel *modelFLT = test->modelFLT->data[i]; - if (modelFLT == NULL) continue; + + // masked for: bad model fit, outlier in parameters + if (test->mask->data.U8[i]) continue; + if (modelFLT == NULL) { + psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n"); + continue; + } psModel *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model // fit model as PSF, not FLT (skip poor fits) - if (!pmSourceFitModel (source, modelPSF, true)) continue; + if (!pmSourceFitModel (source, modelPSF, true)) { + test->mask->data.U8[i] = 1; + continue; + } test->modelPSF->data[i] = modelPSF; Npsf ++; @@ -92,15 +106,21 @@ float fitSum = 0; + // is this metricMask redundant with test->mask ? metricMask->data.U8[i] = 1; test->metric->data.F64[i] = 0; + if (test->mask->data.U8[i]) continue; + psModel *model = test->modelPSF->data[i]; - if (model == NULL) continue; + if (model == NULL) { + psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n"); + continue; + } psImage *image = ((psSource *)test->sources->data[i])->pixels; psImage *mask = ((psSource *)test->sources->data[i])->mask; + // this metric is Ap-Fit float sky = model->params->data.F32[0]; - for (int ix = 0; ix < image->numCols; ix++) { for (int iy = 0; iy < image->numRows; iy++) { @@ -110,10 +130,12 @@ } } - // this metric is Ap-Fit // fprintf (stderr, "%d %f %f %f\n", i, obsSum, fitSum, test->metric->data.F64[i]); + // keep the good metrics if ((fitSum > 0) && (obsSum > 0)) { metricMask->data.U8[i] = 0; test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum); - } + } else { + test->mask->data.U8[i] = 1; + } } @@ -128,13 +150,15 @@ // input: an array of psModels, pre-allocated psf // some of the array entries may be NULL, ignore them -bool pmPSFFromModels (pmPSF *psf, psArray *models) { +bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) { int n; - // construct the x and y vectors from the collection of objects - // count the number of valid model measurements + // construct the fit vectors from the collection of objects + // use the mask to ignore missing fits psVector *x = psVectorAlloc (models->n, PS_TYPE_F64); psVector *y = psVectorAlloc (models->n, PS_TYPE_F64); - n = 0; + psVector *z = psVectorAlloc (models->n, PS_TYPE_F64); + psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64); + for (int i = 0; i < models->n; i++) { psModel *model = models->data[i]; @@ -142,30 +166,24 @@ // XXX EAM : this is fragile: x and y must be stored consistently at 2,3 - x->data.F64[n] = model->params->data.F32[2]; - y->data.F64[n] = model->params->data.F32[3]; - n++; - } - x->n = y->n = n; + x->data.F64[i] = model->params->data.F32[2]; + y->data.F64[i] = model->params->data.F32[3]; + } // we are doing a robust fit. after each pass, we drop points which are - // more deviant than three sigma. - psVector *z = psVectorAlloc (x->n, PS_TYPE_F64); - psVector *dz = psVectorAlloc (x->n, PS_TYPE_F64); + // more deviant than three sigma. + // mask is currently updated for each pass. this is inconsistent + for (int i = 0; i < psf->params->n; i++) { - n = 0; for (int j = 0; j < models->n; j++) { psModel *model = models->data[j]; if (model == NULL) continue; - z->data.F64[n] = model->params->data.F32[i + 4]; - dz->data.F64[n] = 1; - n++; + z->data.F64[j] = model->params->data.F32[i + 4]; + dz->data.F64[j] = 1; // XXX EAM : need to use actual errors? // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4 } - if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n"); - z->n = n; - dz->n = n; - - psf->params->data[i] = RobustFit2D (psf->params->data[i], x, y, z, dz); + // if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n"); + + psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz); // psPolynomial2DDump (psf->params->data[i]); } Index: /trunk/psphot/src/subtract_psf_source.c =================================================================== --- /trunk/psphot/src/subtract_psf_source.c (revision 4115) +++ /trunk/psphot/src/subtract_psf_source.c (revision 4115) @@ -0,0 +1,28 @@ +# include "psphot.h" + +// subtract model fits from image (only PSFSTAR && BRIGHTSTAR) +// need to control application of this with some thresholding + +bool subtract_psf_sources (psSource *source) +{ + float sky; + psSource *source = sources->data[i]; + + if (source->type == PS_SOURCE_GALAXY) return (false); + if (source->type == PS_SOURCE_DEFECT) return (false); + + psModel *model = source->modelPSF; + if (model == NULL) { + psLogMsg ("psphot.subtract_psf_source", "missing model for %f, %f\n", source->moments->x, source->moments->y); + return (false); + } + + psImage *pixels = source->pixels; + + // subtract object, leave local sky + sky = model->params->data.F32[0]; + model->params->data.F32[0] = 0; + pmSourceSubModel (pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + return (true); +} Index: /trunk/psphot/src/subtract_psf_sources.c =================================================================== --- /trunk/psphot/src/subtract_psf_sources.c (revision 4114) +++ /trunk/psphot/src/subtract_psf_sources.c (revision 4115) @@ -6,4 +6,7 @@ bool subtract_psf_sources (psArray *sources) { + + sources = psArraySort (sources, by_SN); + for (int i = 0; i < sources->n; i++) { float sky; Index: /trunk/psphot/src/test_psf_scatter.c =================================================================== --- /trunk/psphot/src/test_psf_scatter.c (revision 4114) +++ /trunk/psphot/src/test_psf_scatter.c (revision 4115) @@ -71,7 +71,7 @@ psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD); - psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL); + psVectorFitPolynomial1D (dsxLine, NULL, snBin, dsxBin, NULL); psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD); - psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL); + psVectorFitPolynomial1D (dsyLine, NULL, snBin, dsyBin, NULL); // these should have a slope of 1.0 psPolynomial1DDump (dsxLine);