Index: /tags/psphot_dev_07/psphot/.cvsignore =================================================================== --- /tags/psphot_dev_07/psphot/.cvsignore (revision 5994) +++ /tags/psphot_dev_07/psphot/.cvsignore (revision 5994) @@ -0,0 +1,1 @@ +bin Index: /tags/psphot_dev_07/psphot/Makefile =================================================================== --- /tags/psphot_dev_07/psphot/Makefile (revision 5994) +++ /tags/psphot_dev_07/psphot/Makefile (revision 5994) @@ -0,0 +1,117 @@ +default: psphot +help: + @echo "USAGE: make psphot" + +CC = gcc +SRC = src +BIN = bin + +DESTBIN = $(HOME)/src/bin/$(ARCH) + +LPSLIB := $(shell pslib-config --libs) +IPSLIB := $(shell pslib-config --cflags) + +LPSMOD := $(shell psmodule-config --libs) +IPSMOD := $(shell psmodule-config --cflags) + +INCS = $(IPSLIB) $(IPSMOD) +LIBS = $(LPSLIB) $(LPSMOD) +CFLAGS = $(INCS) -std=c99 -Wall -Werror -O2 +LFLAGS = $(LIBS) + +PSPHOT = \ +$(SRC)/psphot.$(ARCH).o \ +$(SRC)/psphotArguments.$(ARCH).o \ +$(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psphotModelTest.$(ARCH).o \ +$(SRC)/psphotImageStats.$(ARCH).o \ +$(SRC)/psphotImageBackground.$(ARCH).o \ +$(SRC)/psphotFindPeaks.$(ARCH).o \ +$(SRC)/psphotSourceStats.$(ARCH).o \ +$(SRC)/psphotRoughClass.$(ARCH).o \ +$(SRC)/psphotBasicDeblend.$(ARCH).o \ +$(SRC)/psphotChoosePSF.$(ARCH).o \ +$(SRC)/psphotEnsemblePSF.$(ARCH).o \ +$(SRC)/psphotFullFit.$(ARCH).o \ +$(SRC)/psphotBlendFit.$(ARCH).o \ +$(SRC)/psphotApplyPSF.$(ARCH).o \ +$(SRC)/psphotFitGalaxies.$(ARCH).o \ +$(SRC)/psphotReplaceUnfit.$(ARCH).o \ +$(SRC)/psphotApResid.$(ARCH).o \ +$(SRC)/psphotOutput.$(ARCH).o + +SUPPORT = \ +$(SRC)/psphotModelGroupInit.$(ARCH).o \ +$(SRC)/psphotGrowthCurve.$(ARCH).o \ +$(SRC)/psphotEvalPSF.$(ARCH).o \ +$(SRC)/psphotEvalFLT.$(ARCH).o \ +$(SRC)/psphotSourceFits.$(ARCH).o \ +$(SRC)/psphotSortBySN.$(ARCH).o \ +$(SRC)/psphotDefinePixels.$(ARCH).o \ +$(SRC)/psphotMagnitudes.$(ARCH).o \ +$(SRC)/psphotRadiusChecks.$(ARCH).o + +PSMODULES = \ +$(SRC)/psLine.$(ARCH).o \ +$(SRC)/psSparse.$(ARCH).o \ +$(SRC)/psPolynomialUtils.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ +$(SRC)/pmSourceContour.$(ARCH).o \ +$(SRC)/pmModelFitSet.$(ARCH).o \ +$(SRC)/pmSourceFitSet.$(ARCH).o + +TEST = \ +$(SRC)/psphotTest.$(ARCH).o \ +$(SRC)/psBicube.$(ARCH).o \ +$(SRC)/psphotTestArguments.$(ARCH).o + +psphot: $(BIN)/psphot.$(ARCH) +$(BIN)/psphot.$(ARCH) : $(PSPHOT) $(PSMODULES) $(SUPPORT) +$(PSPHOT) : $(SRC)/psphot.h + +test.install: psphotTest.install +test: $(BIN)/psphotTest.$(ARCH) +$(DESTBIN)/test : $(DESTBIN)/psphotTest + +$(BIN)/psphotTest.$(ARCH) : $(TEST) +$(TEST): $(SRC)/psphot.h + +INSTALL = psphot + +# dependancy rules for binary code ######################### +.PRECIOUS: %.$(ARCH).o +.PRECIOUS: $(BIN)/%.$(ARCH) + +%.$(ARCH).o : %.c + $(CC) $(CFLAGS) -c $< -o $@ + +$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o + @if [ ! -d $(BIN) ]; then mkdir -p $(BIN) || exit; fi + $(CC) $^ -o $@ $(LFLAGS) || exit + @echo "done with $@" + +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN) || exit; fi + rm -f $(DESTBIN)/$* || exit + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* || exit + +$(INSTALL): % : $(BIN)/%.$(ARCH) + +%.clean : + rm -f $(BIN)/$*.$(ARCH) + rm -f $(SRC)/*.$(ARCH).o + +%.install: + make $(DESTBIN)/$* || exit + +# utilities ################################################# + +install: + for i in $(INSTALL); do make $$i.install || exit; done + +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` Index: /tags/psphot_dev_07/psphot/doc/config.txt =================================================================== --- /tags/psphot_dev_07/psphot/doc/config.txt (revision 5994) +++ /tags/psphot_dev_07/psphot/doc/config.txt (revision 5994) @@ -0,0 +1,56 @@ + +# input data options +# IMAGE STR file.fits +# MASK_IMAGE STR mask.fits +# WEIGHT_IMAGE STR weight.fits + +# output data options +OUTPUT_FILE STR output.smp # output object file +OUTPUT_MODE STR TEXT # output mode: TEXT, OBJ, SX, CMP, CMF +# RESID_IMAGE STR resid.fits # output residual image + +# image noise parameters +RDNOISE STR HD:RDNOISE # read-noise in electrons +GAIN STR HD:GAIN # electrons / DN + +# masking parameters +XMIN F32 32 # minimum valid x-coord +XMAX F32 0 # maximum valid x-coord +YMIN F32 1 # minimum valid y-coord +YMAX F32 0 # maximum valid y-coord +SATURATION F32 50000 # saturation level on this chip + +# image statistics parameters +IMSTATS_NPIX S32 100000 # number of pixels to use for sky estimate: + +# peak finding +PEAKS_SMOOTH_SIGMA F32 1.0 # smoothing kernel sigma in pixels +PEAKS_SMOOTH_NSIGMA F32 3.0 # smoothing kernel width in sigmas +PEAKS_NSIGMA_LIMIT F32 10.0 # peak significance threshold + +# basic object statistics +SKY_INNER_RADIUS F32 15 # square annulus for local sky measurement +SKY_OUTER_RADIUS F32 25 # square annulus for local sky measurement +PSF_MOMENTS_RADIUS F32 5 +PSF_SN_LIM F32 100 + +# PSF model parameters : choose the PSF model +# list as many PSF_MODEL options as desired +PSF_MODEL STR PS_MODEL_QGAUSS +#PSF_MODEL STR PS_MODEL_PGAUSS +#PSF_MODEL STR PS_MODEL_GAUSS +PSF_FIT_RADIUS F32 25 # fitting radius for test PSF model + +# PSF model parameters : apply the PSF model +PSF_FIT_NSIGMA F32 1 # significance for pixel included in fit +PSF_FIT_PADDING F32 5 # extra annulus to use for fit +PSF_SHAPE_NSIGMA F32 3.0 # max significance for shape variation +PSF_MIN_SN F32 2.0 # reject objects below this significance +PSF_MAX_CHI F32 10.0 # reject objects worse that this + +# Galaxy model parameters +GAL_MODEL STR PS_MODEL_SGAUSS +GAL_MIN_SN F32 3 +GAL_FIT_NSIGMA F32 1 # significance for pixel included in fit +GAL_FIT_PADDING F32 5 # extra annulus to use for fit +GAL_MOMENTS_RADIUS F32 9 Index: /tags/psphot_dev_07/psphot/doc/notes.txt =================================================================== --- /tags/psphot_dev_07/psphot/doc/notes.txt (revision 5994) +++ /tags/psphot_dev_07/psphot/doc/notes.txt (revision 5994) @@ -0,0 +1,324 @@ + +2005.12.23 + + Improvements in psphot as of psphot_dev_07 + + - updated to work with psModule release 9.0 + - moved model test functions into psphot as optional mode + - added pmSourceFitSet : simultaneous model group fitting + - added bicubic interpolation to tweak up the peak coordinates + - added aperture residual analysis after complete object model & subtraction + - added 3D (x,y,rflux) aperture residual fitting + - added blend vs blob fitting (BlendFit) + - added ensemble PSF fitting, with sparse matrix inversion + - improved the flagging / marking of model results + - better abstraction for setting fit parameters + - various low-level fixes + - added PSF output option + - added break points + - added more optional output formats + - fixed FITS table output + + Things still to be added to psphot: + + - better initial background model (NxM median grid?) + - include local background model in EnsemblePSF Fit (use sparse matrix) + - 2nd pass (after object subtraction) to smooth and detect faint objects + - re-add background to object-subtracted image + - cleanup function names + - push more functions into psLib / psModule + - fix extended / galaxy models to improve stability + - add mode with input PSF + - add mode with input fixed sources (PS_SOURCE_FIXED) + - cleanup memory leaks + +2005.11.25 + + I've updated psphot to work with the current psLib v8, though a + number of psLib fixes were needed. These are pushed under + psLib:eam_rel8_b2. The psphot code which works with that version of + psLib is tagged psphot_dev_04. + + I'm working on converting psphot to work with the current release of + psModules, which basically includes all of the object functions. + + - change comments are relative to psModules (vs psphot versions) + + psEllipse.c : OK + psEllipse.h : OK + pmModelGroup.c : OK (missing TGAUSS entry) + pmModelGroup.h : OK (adds many comments) + + pmPSFtry.c : a bunch of formatting weirdness (try-> got line breaks + added pmPSFtryMetric_Alt + fixed up usage of psVectorClipFitPolynomial1D + fixed up usage of psPolynomial..Alloc + + pmPSFtry.h : OK (Added extensive comments) + + pmPSF.c : fixed up usage of psPolynomial..Alloc + + pmPSF.h : OK (Added extensive comments) + + pmObjects.c : cleaned up some formatting, + fixed usage of psImageCountPixelMask + + pmObjects.h : added apMag and fitMag to pmSource + dropped PS_MODEL_name defines + + psLibUtils.[ch] should be dropped from psModules + psModulesUtils.[ch] should be dropped from psModules + + Makefile.am : dropped psModulesUtils.[ch], psLibUtils.[ch] + + - I have successfully tested psphot with the psModules.v8 code + (eam_rel8_b1 tag). I only needed to make a few modes in the + psModules code. + +2005.11.16 + + I have made some fixes to make psphot work with pslib rel8_0. I + need to merge psphot with the rel8_0 version of the psMinimize and + pmObjects code. + +2005.09.06 + + I have built a working version of PSPhot using my own copy of the + pmObjects.[ch] files. I need to minimize the difference between the + v0.7.0, v0.5.0, and my own version of pmObjects.c + + additions to psLib based on the new version: + - extended the minimization to include parameter limits + + changes v0.5.0 / v0.7.0 / EAM + - changed all ps* to pm* + - fixed naming for local static functions (eg modelFree not p_psModelFree) + - dropped all hard-coded model names (vs v0.7.0) + - using current psMemSetDellocator name + - added noise image to pmSource + - converted old asserts to new asserts (PS_.._CHECK..) + - add return from error in pmFindImagePeaks + - allow isItInThisRegion to have NULL region (true) + - fixed subImages to correspond to region definition (end is + exclusive) in LocalSky and SetPixelCircle + - dropped old, redundant code in pmSourceLocalSky + - fixed pmSourceMoments use of mask (0 is valid, not 1) + - added pmSourceMoments validity tests (numPixels, Sum, centroid shift) + - added correction for Sxy (XY/Sum - x*y) + - added pmSourcePSFClump function (and pmPSFClump structure) + - allow missing sources in pmSourcePSFClump (not all peaks yield + sources) + - moved pmSourceRoughClass metadata lookups out of inner loop + - added saturated pixel test using pmImageCountPixelMask + - various redefinitions of initial classes + - dropped old, redundant code in pmSourceSetPixelsCircle + - mask image set to type psU8 (was psF32) + - abstract object model functions for pmSourceModelGuess + - converted evalModel to a public function psModelEval, modified params + - pmSourceFitModel operates on a specified model + - pmSourceFitModel had option to apply/ignore PSF information + - pmSourceFitModel calculates yErr + - pmSourceFitModel calculates the covariance matrix + - added paramMask to pmSourceFitModel + - pmSourceFitModel_EAM uses covar to carry in beta and param + limits... + - pmSourceFitModel requires solution to stay within image + - pmSourceFitModel calculates and saves chisq, nIter, nDOF + - pmSourceFitModel uses GaussNewtonDelta for frozen params + - added mask to pmSourceAdd, pmSourceSub + - model abstractions in pmSourceAdd, pmSourceSub + - fixed 'center' option in pmSourceAdd/Sub + + - pmPSF params should be F64 + + lingering concerns: + - Polynomials and related functions are STILL messed up. this is a + major priority!!! + - pmSourceFitModel calculates sqrt(var), psMinimize calculates sq(sqrt(var)) + - does psImageAlloc zero the image? if so, drop init routine + from pmSourcePSFClump + - pmSourceRoughClass uses RDNOISE, GAIN to calculate SN: should use + the source->noise image? part of pmSourceMoments? + - pmSourceMoments_EAM uses macro for checkRegion + - pmObjects_EAM.c uses some hardcoded mask values (LocalSky) + - why does LocalSky_EAM not need subImageMask->col0 = subImage->col0?? + - psGetRowVectorFromImage should be in psLib (in psLibUtils) + - pmSourceSetPixelsCircle and pmSourceLocalSky probably could use + the psImageMaskRegion, etc functions. + - pmSourceSetPixelsCircle uses a hard-coded mask value + - define a mask registration process?? + - pmSourceContour is using the model flux, not the image flux + +2005.09.05 + + I have a working version of PSPhot which handles PSF and galaxy + models. I am tagging the module with an alpha version number. this + version compiles with psLib-0.5.0, at least my version of it. My + current goal is to make it work with the current psLib head (nearly + 0.7.0) with a mimimum of _EAM versions of functions. + + the tag for the working version is dev_01 + +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: /tags/psphot_dev_07/psphot/doc/psphot.txt =================================================================== --- /tags/psphot_dev_07/psphot/doc/psphot.txt (revision 5994) +++ /tags/psphot_dev_07/psphot/doc/psphot.txt (revision 5994) @@ -0,0 +1,160 @@ + +TBD priorities: + +- test psImageSmooth / remove psImageSmooth_EAM +- unify pmSourceLocalSky / pmSourceLocalSky_EAM +- unify pmSourceMoments_EAM / pmSourceMoments_EAM +- merge with image data hierarchy stuff +- better name for psPSF_Test + +Defined APIs: + +top-level psphot functions : these do not need to be in psModules + psMetadata *psphotArguments (); + psImageData *psphotSetup (); + psStats *psphotImageStats (); + psArray *pmPeaksSigmaLimit (); + psArray *psphotSourceStats (); + pmPSF *psphotChoosePSF (); + bool psphotApplyPSF (); + bool psphotFitGalaxies (); + void psphotOutput (); + +psphot-specific functions + bool psphotMarkPSF (); + bool psphotSubtractPSF (); + int psphotSortBySN (); + int psphotSaveImage (); + +** add the following entries to the psLibSDRS + +psLib extra utilities + bool psTimerStart (); - already in SDRS + void psTimerFree (); - already in SDRS + bool psTimerClear (); - add to psLibSDRS + psF64 psTimerMark (); - already in SDRS + psS32 psLogArguments (); - now: psArgumentVerbosity in SDRS + psS32 psTraceArguments (); - now: psArgumentVerbosity in SDRS + int psArgumentGet (); - already in SDRS + int psArgumentRemove (); - already in SDRS + psVector *psVectorCreate (); - add to psLibSDRS + int psImageCountPixelMask ();- add to psLibSDRS + psVector *psGetRowVectorFromImage(); - add to psLibSDRS? + bool p_psVectorPrintRow (); - add to psLibSDRS? + +// basic image functions + bool psImageInit (); - already in SDRS + void psImageSmooth_EAM (); - test psLib version and drop + +// psLine functions + psLine *psLineAlloc (); - add to psLibSDRS + bool psLineInit (); - add to psLibSDRS + bool psLineAdd (); - add to psLibSDRS + +** compare with current implementation / update psLib + we need to define a mechanism to carry the beta values and the + parameter limits + +// minimize + psBool p_psMinLM_GuessABP_EAM (); + psBool psMinimizeLMChi2_EAM(); + psF64 p_psMinLM_dLinear (); -- not included in pslib.h files + +** already in psLib : some need to be cleaned up by MHPCC + +// polynomial functions + psF64 Polynomial2DEval(); + psImage *psBuildSums2D(); + psPolynomial2D *VectorFitPolynomial2DOrd_EAM(); + psPolynomial2D *Polynomial2DAlloc(); + void psPolynomial2DDump (); + psPolynomial2D *RobustFit2D_nomask(); + psPolynomial2D *RobustFit2D(); + psVector *Polynomial2DEvalVector(); + + psVector *psBuildSums1D(); + void psPolynomial1DDump (); + psF64 Polynomial1DEval_EAM(); + psVector *Polynomial1DEvalVector_EAM(); + psPolynomial1D *Polynomial1DAlloc(); + psPolynomial1D *VectorFitPolynomial1DOrd_EAM(); + +** make sure these are all in the psModules SDRS: + +pmObjects.h: + pmMoments *pmMomentsAlloc(); + pmModel *pmModelAlloc(); + pmSource *pmSourceAlloc(); + psVector *pmFindVectorPeaks() + psArray *pmFindImagePeaks() + psList *pmCullPeaks() + pmSource *pmSourceLocalSky() + bool pmSourceMoments() + pmPSFClump pmSourcePSFClump() + bool pmSourceRoughClass() + bool pmSourceSetPixelsCircle() + pmModel *pmSourceModelGuess() + psArray *pmSourceContour() + bool pmSourceFitModel() + bool pmSourceAddModel() + bool pmSourceSubModel() + int pmModelParameterCount (); + char *pmModelGetType (); + pmModelType pmModelSetType (); + pmModelFunc pmModelFunc_GetFunction (); + pmModelFlux pmModelFlux_GetFunction (); + pmModelRadius pmModelRadius_GetFunction (); + pmModelLimits pmModelLimits_GetFunction (); + pmModelGuessFunc pmModelGuessFunc_GetFunction (); + pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (); + pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (); + +model function abstractions: + typedef psMinimizeLMChi2Func pmModelFunc; + typedef psF64 (*pmModelFlux)(const psVector *params); + typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); + typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); + typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf); + typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +conversions between elliptical shape representations + EllipseAxes EllipseMomentsToAxes (); + EllipseShape EllipseAxesToShape (); + EllipseAxes EllipseShapeToAxes (); + +output functions + bool pmSourcesWriteText (); -- be careful about psImageData... + bool pmSourcesWriteOBJ (); -- be careful about psImageData... + bool pmSourcesWriteCMP (); -- be careful about psImageData... + bool pmSourcesWriteCMF (); -- be careful about psImageData... + bool pmSourcesWriteSX (); -- be careful about psImageData... + int pmSourcesDophotType (); + bool pmPeaksWriteText (); + bool pmMomentsWriteText (); + bool pmModelWritePSFs (); + bool pmModelWriteFLTs (); + bool pmModelWriteNULLs (); + +// psf utilities + pmPSF *pmPSFAlloc (); + pmPSF_Test *pmPSF_TestAlloc (); + pmPSF_Test *pmPSF_TestModel (); + bool pmPSFFromModels (); + pmModel *pmModelFromPSF (); + bool pmSourcePhotometry (); + bool pmPSFMetricModel (); + +// psModule extra utilities + bool pmSourceFitModel_EAM(); -- requires updated version of psMinimize... + bool pmSourceLocalSky_EAM (); -- compare with pmSourceLocalSky & choose + bool pmSourceMoments_EAM(); -- compare with pmSourceMoments & choose + bool pmSourceDefinePixels(); -- be careful about psImageData... + bool pmModelFitStatus (); + int pmSourceDophotType (); + psF32 pmConfigLookupF32 (); -- merge with Paul's image container hierarchy + +** this needs to be confronted with the phase2 / image container hierarchy + +// psImageData functions + psImageData *psImageDataAlloc (); Index: /tags/psphot_dev_07/psphot/doc/speed.txt =================================================================== --- /tags/psphot_dev_07/psphot/doc/speed.txt (revision 5994) +++ /tags/psphot_dev_07/psphot/doc/speed.txt (revision 5994) @@ -0,0 +1,129 @@ + +psphot -g +psModule -g +psLib -O2 +on alala: + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 0.983568 sec + System random seed value used seed = 25904F7C7454AE25 hex + background: 1229.280029 +/- 27.700180 + image stats: 0.063085 sec + System random seed value used seed = 411E5A0600A49AF3 hex + fit background model: 0.580152 sec + smooth: 1.842415 sec + threshold: 83.100540 DN + 4356 peaks: 0.571587 sec + 4230 moments: 3.965099 sec + SN range: 3.355829 - 1043.401001 + identified 107 blended objects (0.121737 sec) + selected candidate 300 PSF objects + fit flt: 3.682850 sec for 300 sources + fit psf: 1.704141 sec for 300 sources + fit stats: 0.000884 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717999 + select psf model: 5.402742 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.991345 (4230 objects) + built matrix: 1.472869 (27998 elements) + measure ensemble of PSFs: 4.334556 + fit PSF models: 12.205738 sec for 4230 objects + replace unfitted models: 0.021473 sec (4230 objects) + measure aperture residuals : 0.451886 sec + measure full-frame aperture residual: 0.454032 sec + aperture residual: -0.002842 +/- 0.002729 : 0.166072 bias + wrote output: 4.648821 sec + complete psphot run: 35.334080 sec + + +psphot -g +psModule -O2 +psLib -O2 +on alala: + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 1.001502 sec + System random seed value used seed = 5F43A63F532D4B3F hex + background: 1229.280029 +/- 27.700180 + image stats: 0.063169 sec + System random seed value used seed = F5F5FF99A6578F61 hex + fit background model: 0.530384 sec + smooth: 1.856414 sec + threshold: 83.100540 DN + 4365 peaks: 0.237292 sec + 4239 moments: 4.087611 sec + SN range: 3.356469 - 1043.400879 + identified 107 blended objects (0.123372 sec) + selected candidate 300 PSF objects + fit flt: 3.741444 sec for 300 sources + fit psf: 1.676478 sec for 300 sources + fit stats: 0.000883 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717985 + select psf model: 5.458156 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.791613 (4239 objects) + built matrix: 1.267909 (28101 elements) + measure ensemble of PSFs: 3.551364 + fit PSF models: 11.240197 sec for 4239 objects + replace unfitted models: 0.018051 sec (4239 objects) + measure aperture residuals : 0.411856 sec + measure full-frame aperture residual: 0.413843 sec + aperture residual: -0.002831 +/- 0.002717 : 0.164772 bias + wrote output: 4.355429 sec + complete psphot run: 33.126283 sec + +after reducing indirection in pmSourceMoments: + 4268 moments: 3.762452 sec + +psphot -O2 +psModule -O2 +psLib -O2 + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 0.484324 sec + System random seed value used seed = 980050732A392EC7 hex + background: 1229.280029 +/- 27.700180 + image stats: 0.058759 sec + System random seed value used seed = 44E1A6EE77B0E9CE hex + fit background model: 0.395125 sec + smooth: 1.861776 sec + threshold: 83.100540 DN + 4351 peaks: 0.234874 sec + 4264 moments: 3.719435 sec + SN range: 2.023345 - 1043.401123 + identified 124 blended objects (0.115230 sec) + selected candidate 300 PSF objects + fit flt: 3.525501 sec for 300 sources + fit psf: 1.595031 sec for 300 sources + fit stats: 0.000884 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717976 + select psf model: 5.160442 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.760805 (4264 objects) + built matrix: 0.908674 (23703 elements) + measure ensemble of PSFs: 3.115425 + fit PSF models: 10.162397 sec for 4264 objects + replace unfitted models: 0.017491 sec (4264 objects) + measure aperture residuals : 0.408242 sec + measure full-frame aperture residual: 0.410207 sec + aperture residual: -0.002908 +/- 0.002734 : 0.170692 bias + wrote output: 4.225560 sec + complete psphot run: 30.149281 sec + + + fit PSF models: 10.162397 sec for 4264 objects + select psf model: 5.160442 sec + wrote output: 4.225560 sec + 4264 moments: 3.719435 sec + measure ensemble of PSFs: 3.115425 sec + + other : 3.5 sec + + fixed per objects : 4 msec / object (17 of 30 sec) + + +** still can re-compile with -DPS_NO_TRACE +** ASSERTS cannot be compiled out at the moment! Index: /tags/psphot_dev_07/psphot/doc/versions.txt =================================================================== --- /tags/psphot_dev_07/psphot/doc/versions.txt (revision 5994) +++ /tags/psphot_dev_07/psphot/doc/versions.txt (revision 5994) @@ -0,0 +1,33 @@ + +2005.12.23 + psphot_dev_07 should be compiled against psLib:eam_rel9_b1 and + psModule:eam_rel9_b1. This version now handles the ensemble PSF + fitting, detection of blended sources, fitting of 'blends' (groups + of PSFs with distinct be 'blended' peaks) and 'blobs' (unresolved + extended objects, with a test for double PSF vs extended source + model. + +2005.11.25 + + psphot_dev_05 should be compiled against psLib tag eam_rel8_b2 and + psModules tag eam_rel8_b1. this version removes all _EAM versions + of psLib and psModules code. It now uses the psModules version of + the object code. + + psphot_dev_04 should be compiled against psLib tag eam_rel8_b2 this + version removes all _EAM versions of psLib code, so that it is now + consistent with psLib.v8 (fixes are pushed into the eam_rel8_b2 + branch). + +2005.11.16 + + psphot_dev_03 should be compiled against psLib tag eam_rel8_b1, a + branch starting from psLib rel8_0. I needed to make some minor + fixes in psLib to work successfully with psLib rel8_0. + + psphot_dev_02 was getting segfaults from a bug added into psLib in + rel8_0 in p_psVectorClippedStats. I also fixed some errors arising + from changing types (argument to psMetadataParseConfig). + + + Index: /tags/psphot_dev_07/psphot/src/models/pmModel_GAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_GAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_GAUSS.c (revision 5994) @@ -0,0 +1,152 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF32 pmModelFunc_GAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 r = exp(-z); + psF32 q = params->data.F32[1]*r; + psF32 f = q + params->data.F32[0]; + + if (deriv != NULL) { + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +psF64 pmModelFlux_GAUSS(const psVector *params) +{ + psF64 A1 = PS_SQR(params->data.F32[4]); + psF64 A2 = PS_SQR(params->data.F32[5]); + psF64 A3 = PS_SQR(params->data.F32[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// return the radius which yields the requested flux +psF64 pmModelRadius_GAUSS (const psVector *params, psF64 flux) +{ + + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + if (flux >= params->data.F32[1]) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + return (radius); +} + +// define the parameter limits +bool pmModelLimits_GAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (7, PS_TYPE_F32); + *params_min = psVectorAlloc (7, PS_TYPE_F32); + *params_max = psVectorAlloc (7, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + + return (TRUE); +} + +// make an initial guess for parameters +bool pmModelGuess_GAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = moments->x; + params[3] = moments->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + return(true); +} + +// construct the PSF model from the FLT model and the psf +bool pmModelFromPSF_GAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 7; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = psPolynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +// check the status of the fitted model +bool pmModelFitStatus_GAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_PGAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_PGAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_PGAUSS.c (revision 5994) @@ -0,0 +1,167 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF32 pmModelFunc_PGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + psF32 t = 1 + z + z*z/2.0; + psF32 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = PAR[1]*r + psF32 q = PAR[1]*r*r*t; + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +bool pmModelLimits_PGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (7, PS_TYPE_F32); + *params_min = psVectorAlloc (7, PS_TYPE_F32); + *params_max = psVectorAlloc (7, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + + return (TRUE); +} + +psF64 pmModelFlux_PGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + z + z*z/2 + z*z*z/6); + norm += f; + } + norm *= 0.01; + + psF64 Flux = params->data.F32[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_PGAUSS (const psVector *params, psF64 flux) +{ + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + if (flux >= params->data.F32[1]) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = moments->x; + params[3] = moments->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + return(true); +} + +bool pmModelFromPSF_PGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 7; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = psPolynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_PGAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_QGAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_QGAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_QGAUSS.c (revision 5994) @@ -0,0 +1,197 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^M + z^N) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = + params->data.F32[8] = +*****************************************************************************/ + +/* XXX EAM : we need a way to have user-set values for fix parameters */ +# define QG_S1 2.5 +# define dQG_S1 1.5 +/* dQG_S1 is QG_S1 - 1.0 */ + +psF32 pmModelFunc_QGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + psF32 r = 1.0 / (1 + PAR[7]*z + pow(z, QG_S1)); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(PAR[7] + QG_S1*pow(z, dQG_S1)); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -t*z; + } + return(f); +} + +bool pmModelLimits_QGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + beta_lim[0][0].data.F32[7] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + params_min[0][0].data.F32[7] = 0.1; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + params_max[0][0].data.F32[7] = 10.0; + + return (TRUE); +} + +bool pmModelGuess_QGAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + params[7] = 1.0; + return(true); +} + +psF64 pmModelFlux_QGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + PAR[7]*z + pow(z, QG_S1)); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_QGAUSS (const psVector *params, psF64 flux) +{ + psF64 z, f; + psF32 *PAR = params->data.F32; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // if Sx == Sy, sigma = Sx == Sy + psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); + psF64 dz = 1.0 / (2.0 * sigma*sigma); + psF64 limit = flux / PAR[1]; + + // we can do this much better with intelligent choices here + for (z = 0.0; z < 20.0; z += dz) { + f = 1.0 / (1 + PAR[7]*z + pow(z, QG_S1)); + if (f < limit) break; + } + psF64 radius = sigma * sqrt (2.0 * z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelFromPSF_QGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = psPolynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_QGAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_RGAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_RGAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_RGAUSS.c (revision 5994) @@ -0,0 +1,151 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z + z^Npow) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = 1 / SigmaX; + params->data.F32[5] = 1 / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = Npow +*****************************************************************************/ + +psF64 psModelFunc_RGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + // psF32 FACT = 1 + 5*exp(-5*PAR[7]); + + psF32 p = pow(z, PAR[7] - 1.0); + psF32 r = 1.0 / (1 + z + z*p); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(1 + PAR[7]*p); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -q*px*X*2; + deriv->data.F32[5] = -q*py*Y*2; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -5.0*t*log(z)*p*z; + + // deriv->data.F32[4] = -1.8*q*px*X*2; + // deriv->data.F32[5] = -1.8*q*py*Y*2; + // deriv->data.F32[6] = -1.5*q*X*Y; + // deriv->data.F32[7] = -5.0*t*log(z)*p*z; + } + return(f); +} + +psF64 psModelFlux_RGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + z + pow(z, PAR[7])); + norm += f; + } + norm *= 0.01; + + psF64 Flux = params->data.F32[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 psModelRadius_RGAUSS (const psVector *params, psF64 flux) +{ + psF64 z, f, p; + psF32 *PAR = params->data.F32; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // if Sx == Sy, sigma = Sx == Sy + psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); + psF64 dz = 1.0 / (2.0 * sigma*sigma); + psF64 limit = flux / PAR[1]; + + // we can do this much better with intelligent choices here + for (z = 0.0; z < 20.0; z += dz) { + p = pow(z, PAR[7]); + f = 1.0 / (1 + z + p); + if (f < limit) break; + } + psF64 radius = sigma * sqrt (2.0 * z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool psModelGuess_RGAUSS (psModel *model, psSource *source) { + + psVector *params = model->params; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(source->moments->Sx); + moments.y2 = PS_SQR(source->moments->Sy); + moments.xy = source->moments->Sxy; + + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = 1.0 / shape.sx; + params->data.F32[5] = 1.0 / shape.sy; + params->data.F32[6] = shape.sxy; + params->data.F32[7] = 2.0; + return(true); +} + +bool psModelFromPSF_RGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_SGAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_SGAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_SGAUSS.c (revision 5994) @@ -0,0 +1,342 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^Npow + St*z^Ntide) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = 1 / SigmaX; + params->data.F32[5] = 1 / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = Npow + params->data.F32[8] = St +*****************************************************************************/ + +# define SQ(A)((A)*(A)) +psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image); +psF64 p_psImageGetElementF64(psImage *a, int i, int j); + +psF32 pmModelFunc_SGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y), 1e-8); + // note that if z -> 0, dPAR[7] -> -inf + // also z^(PAR[7]-1) -> Inf + + psF32 pr = z*PAR[8]; + psF32 pr3 = pr*pr*pr; + psF32 p = pow(z, PAR[7] - 1.0); + psF32 r = 1.0 / (1 + z*p + pr*pr3); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(PAR[7]*p + 4*PAR[8]*pr3); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -2*t*log(z)*z*p; + deriv->data.F32[8] = -2*t*4*z*pr3; + } + return(f); +} + +bool pmModelLimits_SGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (9, PS_TYPE_F32); + *params_min = psVectorAlloc (9, PS_TYPE_F32); + *params_max = psVectorAlloc (9, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + beta_lim[0][0].data.F32[7] = 0.5; + beta_lim[0][0].data.F32[8] = 0.05; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + params_min[0][0].data.F32[7] = 0.5; + params_min[0][0].data.F32[8] = 0.05; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +3.0; + params_max[0][0].data.F32[7] = 5.0; + params_max[0][0].data.F32[8] = 0.5; + + return (TRUE); +} + +// measure the flux for the elliptical contour +psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image) { + + double t, dt, ct, st, xo, yo, value; + int N, Nt, x, y; + + // choose dt to uniformly divide contour, with ~1 pix spacing at most + dt = asin (1 / axes.minor); + Nt = (int)(2*M_PI / dt) + 1; + dt = 2*M_PI / Nt; + + ct = cos(axes.theta); + st = sin(axes.theta); + xo = xc - image->col0; + yo = yc - image->row0; + + psVector *contour = psVectorAlloc (Nt, PS_TYPE_F32); + for (t = 0, N = 0; (t < 2*M_PI) && (N < Nt); t += dt) { + x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo; + y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo; + value = p_psImageGetElementF64(image, x, y); + if (isfinite(value)) { + contour->data.F32[N] = value; + N++; + } + } + contour->n = N; + // accept every pixel: double counting is not so problematic here... + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, contour, NULL, NULL, 0); + value = stats->sampleMedian; + + psFree (stats); + psFree (contour); + + return (value); +} + +bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source) { + + pmMoments *sMoments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(sMoments->Sx); + moments.y2 = PS_SQR(sMoments->Sy); + moments.xy = sMoments->Sxy; + + // solve the math to go from Moments To Shape + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params[0] = sMoments->Sky; + params[1] = sMoments->Peak - sMoments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.0 / shape.sx; + params[5] = 1.0 / shape.sy; + params[6] = shape.sxy; + +# if (0) + f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + axes.major *= 2.0; + axes.minor *= 2.0; + f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + + if (f1 > f2) { + params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); + } else { + params[7] = 0.5; + } +# endif + + params[7] = 1.8; + params[8] = 0.1; + + + return(true); +} + +// XXX EAM : test version using flux contours to guess slope +bool pmModelGuess_SGAUSS_HARD (pmModel *model, pmSource *source) { + + pmMoments *sMoments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + float f1, f2; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(sMoments->Sx); + moments.y2 = PS_SQR(sMoments->Sy); + moments.xy = sMoments->Sxy; + + // solve the math to go from Moments To Shape + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params[0] = sMoments->Sky; + params[1] = sMoments->Peak - sMoments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.0 / shape.sx; + params[5] = 1.0 / shape.sy; + params[6] = shape.sxy; + + f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + axes.major *= 2.0; + axes.minor *= 2.0; + f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + + if (f1 > f2) { + params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); + } else { + params[7] = 0.5; + } + params[8] = 0.1; + + return(true); +} + +psF64 pmModelFlux_SGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + psF32 pr = PAR[8]*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// XXX need to define the radius along the major axis +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_SGAUSS (const psVector *params, psF64 flux) +{ + psF64 r, z, pr, f; + psF32 *PAR = params->data.F32; + + EllipseAxes axes; + EllipseShape shape; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // convert Sx,Sy,Sxy to major/minor axes + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = EllipseShapeToAxes (shape); + psF64 dr = 1.0 / axes.major; + psF64 limit = flux / PAR[1]; + + // XXX : we can do this faster with an intelligent starting choice + for (r = 0.0; r < 20.0; r += dr) { + z = SQ(r); + pr = PAR[8]*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + if (f < limit) break; + } + psF64 radius = 2.0 * axes.major * sqrt (z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelFromPSF_SGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 9; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = psPolynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_SGAUSS (pmModel *model) { + + psF32 dP; + bool status; + EllipseAxes axes; + EllipseShape shape; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = EllipseShapeToAxes (shape); + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP += PS_SQR(dPAR[7] / PAR[7]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + status &= (fabs(PAR[8]) < 0.5); + status &= (dPAR[8] < 0.1); + status &= (axes.major > 1.41); + status &= (axes.minor > 1.41); + status &= ((axes.major / axes.minor) < 5.0); + status &= (PAR[7] > 0.5); + + if (status) return true; + return false; +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_TAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_TAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_TAUSS.c (revision 5994) @@ -0,0 +1,164 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF32 pmModelFunc_TAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 r = exp(-z); + psF32 q = params->data.F32[1]*r; + psF32 f = q + params->data.F32[0]; + + if (deriv != NULL) { + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +psF64 pmModelFlux_TAUSS(const psVector *params) +{ + psF64 A1 = PS_SQR(params->data.F32[4]); + psF64 A2 = PS_SQR(params->data.F32[5]); + psF64 A3 = PS_SQR(params->data.F32[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// return the radius which yields the requested flux +psF64 pmModelRadius_TAUSS (const psVector *params, psF64 flux) +{ + + if (flux <= 0) + return (1.0); + if (params->data.F32[1] <= 0) + return (1.0); + if (flux >= params->data.F32[1]) + return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + return (radius); +} + +// define the parameter limits +bool pmModelLimits_TAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) +{ + + *beta_lim = psVectorAlloc (7, PS_TYPE_F32); + *params_min = psVectorAlloc (7, PS_TYPE_F32); + *params_max = psVectorAlloc (7, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + + return (TRUE); +} + +// make an initial guess for parameters +bool pmModelGuess_TAUSS (pmModel *model, pmSource *source) +{ + + pmMoments *moments = source->moments; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = moments->x; + params[3] = moments->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + return(true); +} + +// construct the PSF model from the FLT model and the psf +bool pmModelFromPSF_TAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +{ + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 7; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + // XXX: be wary of a bug here. EAM had his own version of psPolynomial2DEval(). + // I think the reason was that hewas using nX=xOrder. Not sure. I changed this + // without verifying. + // out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval(poly, out[2], out[3]); + } + return(true); +} + +// check the status of the fitted model +bool pmModelFitStatus_TAUSS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) + return true; + return false; +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_TGAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_TGAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_TGAUSS.c (revision 5994) @@ -0,0 +1,200 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^M + z^N) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = + params->data.F32[8] = +*****************************************************************************/ + +/* XXX EAM : we need a way to have user-set values for fix parameters */ +# define TG_S 2.60 +# define dTG_S 1.60 +# define TRF 0.50 +/* dTG_S is TG_S - 1.0 */ + +psF32 pmModelFunc_TGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + psF32 er = exp(TRF*z); + + psF32 r = 1.0 / (er + PAR[7]*z + pow(z, TG_S)); // (1/R) + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + psF32 t = PAR[1]*r*r; // df/dR + psF32 q = t*(TRF*er + PAR[7] + TG_S*pow(z, dTG_S)); // (df/dR)(dR/dz) + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -t*z; + } + return(f); +} + +bool pmModelLimits_TGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + beta_lim[0][0].data.F32[7] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + params_min[0][0].data.F32[7] = 0.1; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + params_max[0][0].data.F32[7] = 10.0; + + return (TRUE); +} + +bool pmModelGuess_TGAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + params[7] = 1.0; + return(true); +} + +psF64 pmModelFlux_TGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + psF32 er = exp(TRF*z); + f = 1.0 / (er + PAR[7]*z + pow(z, TG_S)); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_TGAUSS (const psVector *params, psF64 flux) +{ + psF64 z, f; + psF32 *PAR = params->data.F32; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // if Sx == Sy, sigma = Sx == Sy + psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); + psF64 dz = 1.0 / (2.0 * sigma*sigma); + psF64 limit = flux / PAR[1]; + + // we can do this much better with intelligent choices here + for (z = 0.0; z < 20.0; z += dz) { + psF32 er = exp(TRF*z); + f = 1.0 / (er + PAR[7]*z + pow(z, TG_S)); + if (f < limit) break; + } + psF64 radius = sigma * sqrt (2.0 * z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelFromPSF_TGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = psPolynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_TGAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_WAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_WAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_WAUSS.c (revision 5994) @@ -0,0 +1,99 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF64 psModelFunc_WAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[2]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); + psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ + psF32 f = params->data.F32[1]*r + params->data.F32[0]; + + if (deriv != NULL) { + // note difference from gaussian: q = params->data.F32[1]*r + psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; + deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; + // The values of 100 dampen the swing of params->data.F32[7,8] */ + } + return(f); +} + +// this is probably wrong since it uses the gauss integral 2 pi sigma^2 +psF64 psModelFlux_WAUSS(const psVector *params) +{ + psF64 A1 = 1 / PS_SQR(params->data.F32[4]); + psF64 A2 = 1 / PS_SQR(params->data.F32[5]); + psF64 A3 = params->data.F32[6]; + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// return the radius which yields the requested flux +psF64 psModelRadius_WAUSS (const psVector *params, psF64 flux) +{ + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + return (radius); +} + +bool psModelGuess_WAUSS (psModel *model, pmSource *source) { + + psVector *params = model->params; + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = sqrt(2.0) / source->moments->Sx; + params->data.F32[5] = sqrt(2.0) / source->moments->Sy; + params->data.F32[6] = source->moments->Sxy; + // XXX: What are these? + // params->data.F32[7] = B2; + // params->data.F32[8] = B3; + return(true); +} + +bool psModelFromPSF_WAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 9; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: /tags/psphot_dev_07/psphot/src/models/pmModel_ZGAUSS.c =================================================================== --- /tags/psphot_dev_07/psphot/src/models/pmModel_ZGAUSS.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/models/pmModel_ZGAUSS.c (revision 5994) @@ -0,0 +1,158 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^Npow + PAR8*z^4) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = 1 / SigmaX; + params->data.F32[5] = 1 / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = Npow +*****************************************************************************/ + +# define SQ(A)((A)*(A)) +# define PAR8 0.1 + +psF64 psModelFunc_ZGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + psF32 pr = PAR8*z; + psF32 p = pow(z, PAR[7] - 1.0); + psF32 r = 1.0 / (1 + z*p + SQ(SQ(pr))); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(PAR[7]*p + 4*PAR8*pr*pr*pr); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -q*px*X; + deriv->data.F32[5] = -q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -t*log(z)*z*p; + } + return(f); +} + +psF64 psModelFlux_ZGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + psF32 pr = PAR8*z; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// XXX need to define the radius along the major axis +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 psModelRadius_ZGAUSS (const psVector *params, psF64 flux) +{ + psF64 r, z, pr, f; + psF32 *PAR = params->data.F32; + + EllipseAxes axes; + EllipseShape shape; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // convert Sx,Sy,Sxy to major/minor axes + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = EllipseShapeToAxes (shape); + psF64 dr = 1.0 / axes.major; + psF64 limit = flux / PAR[1]; + + // XXX : we can do this faster with an intelligent starting choice + for (r = 0.0; r < 20.0; r += dr) { + z = SQ(r); + pr = PAR8*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + if (f < limit) break; + } + psF64 radius = 2.0 * axes.major * sqrt (z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool psModelGuess_ZGAUSS (psModel *model, psSource *source) { + + psVector *params = model->params; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(source->moments->Sx); + moments.y2 = PS_SQR(source->moments->Sy); + moments.xy = source->moments->Sxy; + + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = 1.0 / shape.sx; + params->data.F32[5] = 1.0 / shape.sy; + params->data.F32[6] = shape.sxy; + params->data.F32[7] = 1.9; + return(true); +} + +bool psModelFromPSF_ZGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: /tags/psphot_dev_07/psphot/src/pmModelFitSet.c =================================================================== --- /tags/psphot_dev_07/psphot/src/pmModelFitSet.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/pmModelFitSet.c (revision 5994) @@ -0,0 +1,56 @@ +# include "psphot.h" + +// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, +// nPar = nSrc*(nOnePar - 1) + 1 + +static pmModelFunc mFunc; +static int nPar; +static psVector *onePar; +static psVector *oneDeriv; + +bool pmModelFitSetInit (pmModelType type) { + + mFunc = pmModelFunc_GetFunction (type); + nPar = pmModelParameterCount (type); + + onePar = psVectorAlloc (nPar, PS_DATA_F32); + oneDeriv = psVectorAlloc (nPar, PS_DATA_F32); + + return true; +} + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x) +{ + + psF32 value; + psF32 model; + + psF32 *PAR = onePar->data.F32; + psF32 *dPAR = oneDeriv->data.F32; + + int nSrc = (params->n - 1) / (nPar - 1); + + model = params->data.F32[0]; + + PAR[0] = 0; + for (int i = 0; i < nSrc; i++) { + for (int n = 1; n < nPar; n++) { + PAR[n] = params->data.F32[i*nPar - i + n]; + } + value = mFunc (oneDeriv, onePar, x); + model += value; + for (int n = 1; n < nPar; n++) { + deriv->data.F32[i*nPar - i + n] = dPAR[n]; + } + } + return (model); +} + +/* +i: 0 1 2 +n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 +i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 +*/ + Index: /tags/psphot_dev_07/psphot/src/pmSourceContour.c =================================================================== --- /tags/psphot_dev_07/psphot/src/pmSourceContour.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/pmSourceContour.c (revision 5994) @@ -0,0 +1,213 @@ +# include "psphot.h" + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +# define LEFT false +# define RIGHT true + +// return the first coordinate at or below the threshold in the requested direction +static int findContourNeg( + psImage *image, + float threshold, + int x, + int y, + bool right) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = image->numCols - 1; + } else { + incr = -1; + lastColumn = 0; + } + + subCol = x; + + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value <= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +// return the last coordinate at or below the threshold in the requested direction +static int findContourPos( + psImage *image, + float threshold, + int x, + int y, + bool right, + int xEnd) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = PS_MIN (image->numCols - 1, xEnd); + } else { + incr = -1; + lastColumn = PS_MAX (0, xEnd); + } + + subCol = x; + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value >= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + if (subCol == x) { + return (subCol); + } + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +/****************************************************************************** +new implementation of source contour function +*****************************************************************************/ +psArray *pmSourceContour_EAM (psImage *image, int xc, int yc, float threshold) { + + int xR, yR, x0, x1, x0s, x1s; + + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + + int x = xc - image->col0; + int y = yc - image->row0; + + // the requested point must be within the contour + if (image->data.F32[y][x] < threshold) return NULL; + + // Ensure that the starting column is allowable. + if (x < 0) return NULL; + if (y < 0) return NULL; + if (x >= image->numCols) return NULL; + if (y >= image->numRows) return NULL; + + // Allocate data for x/y pairs. + psVector *xVec = psVectorAlloc(100, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(100, PS_TYPE_F32); + + // First row: find the left and right end-points + int Npt = 0; + + x0 = findContourNeg (image, threshold, x, y, LEFT); + x1 = findContourNeg (image, threshold, x, y, RIGHT); + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + y; + yVec->data.F32[Npt + 1] = image->row0 + y; + Npt += 2; + + x0s = x0; + x1s = x1; + + // look for contour outline above row + xR = x0s; + yR = y + 1; + while (yR < image->numRows) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt1; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR ++; + xR = x0; + } + +pt1: + // look for contour outline below row + xR = x0s; + x1 = x1s; + yR = y - 1; + while (yR >= 0) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt2; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + // fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR --; + } +pt2: + xVec->n = Npt; + yVec->n = Npt; + + // fprintf (stderr, "done\n"); + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) xVec; + tmpArray->data[1] = (psPtr *) yVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); +} Index: /tags/psphot_dev_07/psphot/src/pmSourceFitSet.c =================================================================== --- /tags/psphot_dev_07/psphot/src/pmSourceFitSet.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/pmSourceFitSet.c (revision 5994) @@ -0,0 +1,224 @@ +# include "psphot.h" + +bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, bool PSF) { + + double x, y, Io; + + FILE *f = fopen (fitset, "r"); + if (f == NULL) return false; + + psArray *modelSet = psArrayAlloc (16); + modelSet->n = 0; + + while (fscanf (f, "%lf %lf %lf", &x, &y, &Io) == 3) { + pmModel *model = pmModelAlloc (oneModel->type); + + for (psS32 i = 0; i < model->params->n; i++) { + model->params->data.F32[i] = oneModel->params->data.F32[i]; + model->dparams->data.F32[i] = oneModel->dparams->data.F32[i]; + } + model->params->data.F32[1] = Io; + model->params->data.F32[2] = x; + model->params->data.F32[3] = y; + psArrayAdd (modelSet, 16, model); + } + + pmSourceFitSet (source, modelSet, PSF); + + // write out positive object + psphotSaveImage (NULL, source->pixels, "object.fits"); + + // subtract object, leave local sky + for (int i = 0; i < modelSet->n; i++) { + pmModel *model = modelSet->data[i]; + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + fprintf (stderr, "output parameters (obj %d):\n", i); + for (int n = 0; n < model->params->n; n++) { + fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]); + } + } + + // write out + psphotSaveImage (NULL, source->pixels, "resid.fits"); + psphotSaveImage (NULL, source->mask, "mask.fits"); + return true; +} + +# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 +# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 + +bool pmSourceFitSet (pmSource *source, + psArray *modelSet, + const bool PSF) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // base values on first model + pmModel *model = modelSet->data[0]; + + // set the static variables + pmModelFitSetInit (model->type); + + int nSrc = modelSet->n; + int nPar = model->params->n - 1; // number of object parameters (excluding sky) + + psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + psVector *paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL; + + // all but the sky are allowed to vary independently (subject to PSF) + params->data.F32[0] = model->params->data.F32[0]; + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + for (int n = 1; n < nPar + 1; n++) { + params->data.F32[i*nPar + n] = model->params->data.F32[n]; + dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n]; + if (PSF) { + paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1; + } + } + } + + // PSF model only fits first 4 parameters, EXT model fits all + int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1; + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j]) { + continue; + } + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (source->weight->data.F32[i][j] == 0) { + continue; + } + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + nPix++; + } + } + if (nPix < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + return(false); + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // Set the parameter range checks + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + psVector *beta_lim = NULL; + psVector *params_min = NULL; + psVector *params_max = NULL; + + // XXX EAM : in this implementation, I pass in the limits with the covar matrix. + // in the SDRS, I define a new psMinimization which will take these in + psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); + modelLimits (&beta_lim, ¶ms_min, ¶ms_max); + for (int i = 0; i < nSrc; i++) { + for (int n = 1; n < nPar; n++) { + covar->data.F64[0][i*nPar + n] = beta_lim->data.F32[n]; + covar->data.F64[1][i*nPar + n] = params_min->data.F32[n]; + covar->data.F64[2][i*nPar + n] = params_max->data.F32[n]; + } + } + covar->data.F64[0][0] = beta_lim->data.F32[0]; + covar->data.F64[1][0] = params_min->data.F32[0]; + covar->data.F64[2][0] = params_max->data.F32[0]; + + psTrace (".pmObjects.pmSourceFitSet", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, pmModelFitSet); + + // parameter errors from the covariance matrix + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + } + + // assign back the parameters to the models + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + model->params->data.F32[0] = params->data.F32[0]; + for (int n = 1; n < nPar + 1; n++) { + model->params->data.F32[n] = params->data.F32[i*nPar + n]; + model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; + } + // save the resulting chisq, nDOF, nIter + // these are not unique for any one source + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // set the model success or failure status + model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE; + + // models can go insane: reject these + onPic &= (model->params->data.F32[2] >= source->pixels->col0); + onPic &= (model->params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (model->params->data.F32[3] >= source->pixels->row0); + onPic &= (model->params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + } + + source->mode |= PM_SOURCE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(paramMask); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + Index: /tags/psphot_dev_07/psphot/src/psImageData.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psImageData.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psImageData.c (revision 5994) @@ -0,0 +1,24 @@ +# include "psphot.h" + +static void eamReadoutFree (eamReadout *imdata) { + + if (imdata == NULL) return; + + psFree (imdata->image); + psFree (imdata->header); + psFree (imdata->weight); + psFree (imdata->mask); + return; +} + +eamReadout *eamReadoutAlloc (psImage *image, psImage *weight, psImage *mask, psMetadata *header) { + + eamReadout *imdata = psAlloc(sizeof(eamReadout)); + imdata->image = image; + imdata->header = header; + imdata->weight = weight; + imdata->mask = mask; + + psMemSetDeallocator(imdata, (psFreeFunc) eamReadoutFree); + return (imdata); +} Index: /tags/psphot_dev_07/psphot/src/psLibUtils.h =================================================================== --- /tags/psphot_dev_07/psphot/src/psLibUtils.h (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psLibUtils.h (revision 5994) @@ -0,0 +1,17 @@ + +# ifndef PS_LIB_UTILS +# define PS_LIB_UTILS + +// structure to carry a dynamic string +typedef struct { + int NLINE; + int Nline; + char *line; +} psLine; + +// psLine functions -- keep out for now? +psLine *psLineAlloc (int Nline); +bool psLineInit (psLine *line); +bool psLineAdd (psLine *line, char *format, ...); + +# endif Index: /tags/psphot_dev_07/psphot/src/psLine.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psLine.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psLine.c (revision 5994) @@ -0,0 +1,48 @@ +# include +# include "psLibUtils.h" + +static void psLineFree (psLine *line) { + + if (line == NULL) return; + + psFree (line->line); + return; +} + +// allocate a psLine structrue +psLine *psLineAlloc (int Nline) { + + psLine *line; + line = psAlloc (sizeof(psLine)); + psMemSetDeallocator(line, (psFreeFunc) psLineFree); + + line->Nline = 0; + line->NLINE = Nline; + line->line = psAlloc (Nline); + return (line); +} + +bool psLineInit (psLine *line) { + if (line == NULL) return (false); + line->Nline = 0; + return (true); +} + +bool psLineAdd (psLine *line, char *format, ...) { + + int Nchar; + va_list ap; + + if (line == NULL) return (false); + + int nMax = line->NLINE - line->Nline; + + va_start (ap, format); + Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap); + line->Nline += PS_MIN (nMax - 1, Nchar); + va_end (ap); + + if (Nchar >= nMax) return (false); + return (true); +} + Index: /tags/psphot_dev_07/psphot/src/psModulesUtils.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psModulesUtils.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psModulesUtils.c (revision 5994) @@ -0,0 +1,163 @@ +# include "psModulesUtils.h" + +// extract config informatin from config data or from image header, if specified +// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme +psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) { + + char *source; + char *keyword; + psF32 value; + psMetadataItem *item; + + // look for the entry in the config collection + item = psMetadataLookup (config, name); + if (item == NULL) { + psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name); + value = psMetadataLookupF32 (status, header, name); + if (*status == false) { + psAbort ("pmConfigLookupF32", "no key %s in header", name); + } + *status = true; + return (value); + } + + // I'm either expecting a string, with the name "HD:keyword"... + if (item->type == PS_DATA_STRING) { + source = item->data.V; + if (!strncasecmp (source, "HD:", 3)) { + keyword = &source[3]; + value = psMetadataLookupF32 (status, header, keyword); + if (*status == false) { + psAbort ("pmConfigLookupF32", "no key %s in config", name); + } + *status = true; + // psFree (item); + return (value); + } + } + + // ... or a value (F32?) + if (item->type == PS_DATA_F32) { + value = item->data.F32; + // psFree (item); + return (value); + } + + *status = false; + psError(PS_ERR_UNKNOWN, true, "invalid item"); + return (0); +} + +bool pmModelFitStatus (pmModel *model) { + + bool status; + + pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); + status = statusFunc (model); + + return (status); +} + +bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) { + + float obsSum = 0; + float fitSum = 0; + float sky = model->params->data.F32[0]; + + *fitMag = 99.0; + *obsMag = 99.0; + + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + if (fitSum <= 0) return false; + if (!isfinite(fitSum)) return false; + *fitMag = -2.5*log10(fitSum); + + for (int ix = 0; ix < image->numCols; ix++) { + for (int iy = 0; iy < image->numRows; iy++) { + if (mask->data.U8[iy][ix]) continue; + obsSum += image->data.F32[iy][ix] - sky; + } + } + if (obsSum <= 0) return false; + *obsMag = -2.5*log10(obsSum); + + return (true); +} + +// force the fitted sky to meet the source flux at outer radius +bool pmModelSkyOffset (pmModel *model, float x, float y, float radius) { + + float flux; + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + psVector *params = model->params; + + coord->data.F32[0] = x + radius; + coord->data.F32[1] = y; + flux = modelFunc (NULL, params, coord); + params->data.F32[0] = flux; + + psFree (coord); + + return true; +} + +pmModel *pmModelCopy (pmModel *model) { + + pmModel *new = pmModelAlloc (model->type); + + new->chisq = model->chisq; + new->nDOF = model->nDOF; + new->nIter = model->nIter; + new->status = model->status; + new->radius = model->radius; + + for (int i = 0; i < new->params->n; i++) { + new->params->data.F32[i] = model->params->data.F32[i]; + new->dparams->data.F32[i] = model->dparams->data.F32[i]; + } + + return (new); +} + +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) { + + int Xs, Xe, Ys, Ye; + int xi, xj, yi, yj; + int xIs, xJs, yIs, yJs; + int xIe, yIe; + float flux; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Ti = Mi->mask; + psImage *Tj = Mj->mask; + + Xs = PS_MAX (Pi->col0, Pj->col0); + Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); + + Ys = PS_MAX (Pi->row0, Pj->row0); + Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); + + xIs = Xs - Pi->col0; + xJs = Xs - Pj->col0; + yIs = Ys - Pi->row0; + yJs = Ys - Pj->row0; + + xIe = Xe - Pi->col0; + yIe = Ye - Pi->row0; + + flux = 0; + for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { + for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { + if (Ti->data.U8[yi][xi]) continue; + if (Tj->data.U8[yj][xj]) continue; + flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]; + } + } + return (flux); +} Index: /tags/psphot_dev_07/psphot/src/psModulesUtils.h =================================================================== --- /tags/psphot_dev_07/psphot/src/psModulesUtils.h (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psModulesUtils.h (revision 5994) @@ -0,0 +1,21 @@ + +# ifndef PS_MODULE_UTILS +# define PS_MODULE_UTILS + +# include // for strcasecmp +# include +# include +# include +# include "psLibUtils.h" + +// psModule extra utilities +// bool pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF); +bool pmModelFitStatus (pmModel *model); +int pmSourceDophotType (pmSource *source); +bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask); +bool pmModelSkyOffset (pmModel *model, float x, float y, float radius); + +// unify with paul's image/header/metadata functions +psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); + +# endif Index: /tags/psphot_dev_07/psphot/src/psPolynomialUtils.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psPolynomialUtils.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psPolynomialUtils.c (revision 5994) @@ -0,0 +1,168 @@ +# include +# include // for strcasecmp +# include // for unlink +# include + +psPolynomial4D *psVectorChiClipFitPolynomial4D( + psPolynomial4D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_PTR_NON_NULL(stats, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL); + PS_ASSERT_VECTOR_NON_NULL(fErr, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + + // clipping range defined by min and max and/or clipSigma + float minClipSigma; + float maxClipSigma; + if (isfinite(stats->max)) { + maxClipSigma = +fabs(stats->max); + } else { + maxClipSigma = +fabs(stats->clipSigma); + } + if (isfinite(stats->min)) { + minClipSigma = -fabs(stats->min); + } else { + minClipSigma = -fabs(stats->clipSigma); + } + psVector *fit = NULL; + psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); + + // eventual expansion: user supplies one of various stats option pairs, + // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to + // evaluate the clipping sigma + // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used + stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + for (int N = 0; N < stats->clipIter; N++) { + int Nkeep = 0; + + poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); + fit = psPolynomial4DEvalVector (poly, x, y, z, t); + resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); + + stats = psVectorStats (stats, resid, NULL, mask, maskValue); + psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + + // set mask if pts are not valid + // we are masking out any point which is out of range + // recovery is not allowed with this scheme + for (int i = 0; i < resid->n; i++) { + if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { + continue; + } + float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev); + if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + Nkeep ++; + } + + psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", + Nkeep, x->n); + + stats->clippedNvalues = Nkeep; + psFree (fit); + } + // Free local temporary variables + psFree (resid); + + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); + return(NULL); + } + return(poly); +} + +psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) { + + int ix = x - image->col0; + int iy = y - image->row0; + + psF32 *Fm = &image->data.F32[iy - 1][ix]; + psF32 *Fo = &image->data.F32[iy + 0][ix]; + psF32 *Fp = &image->data.F32[iy + 1][ix]; + + double Fxm = Fm[-1] + Fo[-1] + Fp[-1]; + double Fxp = Fm[+1] + Fo[+1] + Fp[+1]; + double Fym = Fm[-1] + Fm[+0] + Fm[+1]; + double Fyp = Fp[-1] + Fp[+0] + Fp[+1]; + double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1]; + + psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD); + poly->mask[2][2] = 1; + poly->mask[1][2] = 1; + poly->mask[2][1] = 1; + + poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ; + + poly->coeff[1][0] = (Fxp - Fxm)/6.0; + poly->coeff[0][1] = (Fyp - Fym)/6.0; + + poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0; + poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0; + + poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0; + + return (poly); +} + +psPlane psImageBicubeMin (psPolynomial2D *poly) { + + psPlane min; + + min.xErr = min.yErr = 0; + + double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); + + min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det; + min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det; + return (min); +} + +// may need to extend this to change the keyname in the copy +bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key) { + + psMetadataItem *item = psMetadataLookup (in, key); + if (item == NULL) return false; + + psMetadataAddItem (out, item, PS_LIST_TAIL, PS_META_REPLACE); + return true; +} Index: /tags/psphot_dev_07/psphot/src/psSparse.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psSparse.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psSparse.c (revision 5994) @@ -0,0 +1,221 @@ +# include "psphot.h" + +void psSparseMatrixTest () { + + // build a sparse matrix + psSparse *sparse = psSparseAlloc (3, 9); + + psSparseMatrixElement (sparse, 0, 0, 3.0); + psSparseMatrixElement (sparse, 1, 1, 2.0); + psSparseMatrixElement (sparse, 2, 2, 1.0); + + psSparseMatrixElement (sparse, 1, 0, 0.1); + psSparseMatrixElement (sparse, 2, 0, -0.1); + + psSparseResort (sparse); + for (int i = 0; i < sparse->Nelem; i++) { + fprintf (stderr, "%d %d %f\n", + sparse->Si->data.S32[i], + sparse->Sj->data.S32[i], + sparse->Aij->data.F32[i]); + } + + psVector *x = psVectorAlloc (3, PS_DATA_F32); + x->data.F32[0] = 3; + x->data.F32[1] = 5; + x->data.F32[2] = 7; + fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); + + psVector *B = psSparseMatrixTimesVector (NULL, sparse, x); + fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]); + + sparse->Bfj->data.F32[0] = B->data.F32[0]; + sparse->Bfj->data.F32[1] = B->data.F32[1]; + sparse->Bfj->data.F32[2] = B->data.F32[2]; + + x = psSparseSolve (x, sparse, 0); + fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); + + x = psSparseSolve (x, sparse, 1); + fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); + + x = psSparseSolve (x, sparse, 2); + fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); + + x = psSparseSolve (x, sparse, 3); + fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); + + x = psSparseSolve (x, sparse, 4); + fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); + return; +} + +void psSparseResort (psSparse *sparse) { + + int Nelem = sparse->Nelem; + + psVector *index = psVectorSortIndex (NULL, sparse->Sj); + psVector *Aij = sparse->Aij; + psVector *Si = sparse->Si; + psVector *Sj = sparse->Sj; + + // allocate new temporary vectors + psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32); + psVector *tSi = psVectorAlloc (Nelem, PS_DATA_S32); + psVector *tSj = psVectorAlloc (Nelem, PS_DATA_S32); + for (int i = 0; i < Nelem; i++) { + int j = index->data.U32[i]; + tAij->data.F32[i] = Aij->data.F32[j]; + tSi->data.S32[i] = Si->data.S32[j]; + tSj->data.S32[i] = Sj->data.S32[j]; + } + psFree (Aij); + psFree (Si); + psFree (Sj); + + sparse->Aij = tAij; + sparse->Si = tSi; + sparse->Sj = tSj; + return; +} + +// user should only add elements above the diagonal, but we don't check this +void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) { + + int k; + + if (i < j) { + fprintf (stderr, "*** error: subdiagonal element ***\n"); + return; + } + + if (i == j) { + // add to the diagonal + sparse->Qii->data.F32[i] = value; + + // check vectors lengths and extend if needed + if (sparse->Nelem >= sparse->Aij->nalloc) { + psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100); + psVectorRealloc (sparse->Si, sparse->Si->nalloc + 100); + psVectorRealloc (sparse->Sj, sparse->Sj->nalloc + 100); + } + + k = sparse->Nelem; + sparse->Aij->data.F32[k] = value; + sparse->Si->data.S32[k] = i; + sparse->Sj->data.S32[k] = j; + + sparse->Nelem ++; + sparse->Aij->n ++; + sparse->Si->n ++; + sparse->Sj->n ++; + } else { + // check vectors lengths and extend if needed + if (sparse->Nelem >= sparse->Aij->nalloc - 1) { + psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100); + psVectorRealloc (sparse->Si, sparse->Si->nalloc + 100); + psVectorRealloc (sparse->Sj, sparse->Sj->nalloc + 100); + } + + k = sparse->Nelem; + sparse->Aij->data.F32[k] = value; + sparse->Si->data.S32[k] = i; + sparse->Sj->data.S32[k] = j; + k++; + + sparse->Aij->data.F32[k] = value; + sparse->Si->data.S32[k] = j; + sparse->Sj->data.S32[k] = i; + + sparse->Nelem += 2; + sparse->Aij->n += 2; + sparse->Si->n += 2; + sparse->Sj->n += 2; + } + return; +} + +void psSparseVectorElement (psSparse *sparse, int i, float value) { + + sparse->Bfj->data.F32[i] = value; + return; +} + +// multiple A * x +psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) { + + int i, Nelem; + float F; + + if (output == NULL) { + output = psVectorAlloc (vector->n, PS_DATA_F32); + } + + Nelem = 0; + for (int j = 0; j < vector->n; j++) { + F = 0; + while (matrix->Sj->data.S32[Nelem] == j) { + i = matrix->Si->data.S32[Nelem]; + F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem]; + Nelem++; + } + output->data.F32[j] = F; + } + return (output); +} + +psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) { + + psF32 dG; + + psVector *Qii = sparse->Qii; + psVector *Bfj = sparse->Bfj; + + guess = psVectorCopy (guess, Bfj, PS_DATA_F32); + + // temporary storage for intermediate results + psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32); + + for (int j = 0; j < Niter; j++) { + dQ = psSparseMatrixTimesVector (dQ, sparse, guess); + for (int i = 0; i < dQ->n; i++) { + dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i]; + guess->data.F32[i] -= dG; + } + } + psFree (dQ); + return (guess); +} + +static void psSparseFree (psSparse *sparse) { + if (sparse == NULL) return; + psFree (sparse->Aij); + psFree (sparse->Bfj); + psFree (sparse->Qii); + psFree (sparse->Si); + psFree (sparse->Sj); + return; +} + +// allocate a sparse matrix container for Nrows, with Nelem slots allocated +psSparse *psSparseAlloc (int Nrows, int Nelem) { + + psSparse *sparse = (psSparse *) psAlloc (sizeof(psSparse)); + sparse->Aij = psVectorAlloc (Nelem, PS_DATA_F32); + sparse->Si = psVectorAlloc (Nelem, PS_DATA_S32); + sparse->Sj = psVectorAlloc (Nelem, PS_DATA_S32); + + sparse->Aij->n = 0; + sparse->Si->n = 0; + sparse->Sj->n = 0; + sparse->Nelem = 0; + + sparse->Bfj = psVectorAlloc (Nrows, PS_DATA_F32); + sparse->Qii = psVectorAlloc (Nrows, PS_DATA_F32); + + sparse->Nrows = Nrows; + psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree); + return (sparse); +} + + Index: /tags/psphot_dev_07/psphot/src/psSparse.h =================================================================== --- /tags/psphot_dev_07/psphot/src/psSparse.h (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psSparse.h (revision 5994) @@ -0,0 +1,18 @@ + +typedef struct { + psVector *Aij; + psVector *Bfj; + psVector *Qii; + psVector *Si; + psVector *Sj; + int Nelem; + int Nrows; +} psSparse; + +void psSparseMatrixTest (); +void psSparseResort (psSparse *sparse); +void psSparseMatrixElement (psSparse *sparse, int i, int j, float value); +void psSparseVectorElement (psSparse *sparse, int i, float value); +psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector); +psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter); +psSparse *psSparseAlloc (int Nrows, int Nelem); Index: /tags/psphot_dev_07/psphot/src/psphot.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphot.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphot.c (revision 5994) @@ -0,0 +1,83 @@ +# include "psphot.h" + +// XXX need a better structure for handling optional sequence +int main (int argc, char **argv) { + + psMetadata *config = NULL; + eamReadout *imdata = NULL; + psArray *sources = NULL; + psArray *peaks = NULL; + pmPSF *psf = NULL; + psStats *sky = NULL; + bool status; + + psTimerStart ("complete"); + + // load command-line arguments and options + config = psphotArguments (&argc, argv); + + // load input data (config and images (signal, noise, mask) + imdata = psphotSetup (config); + + // run a single-model test if desired + psphotModelTest (imdata, config); + + // measure image stats for initial guess + sky = psphotImageStats (imdata, config); + + // generate a background model (currently, 2D polynomial) + // XXX this should be available to be re-added to the original image + psphotImageBackground (imdata, config, sky); + + // find the peaks in the image + peaks = psphotFindPeaks (imdata, config, sky); + + // construct sources and measure basic stats + sources = psphotSourceStats (imdata, config, peaks); + + // classify sources based on moments, brightness + psphotRoughClass (sources, config); + + // mark blended peaks PS_SOURCE_BLEND + psphotBasicDeblend (sources, config, sky); + + // use bright stellar objects to measure PSF + psf = psphotChoosePSF (config, sources, sky); + + // XXX change FITMODE to a string + int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE"); + switch (FITMODE) { + case 0: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + break; + + case 1: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + psphotFullFit (imdata, config, sources, psf, sky); + psphotReplaceUnfit (sources); + psphotApResid (imdata, sources, config, psf); + break; + + case 2: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + psphotBlendFit (imdata, config, sources, psf, sky); + psphotReplaceUnfit (sources); + psphotApResid (imdata, sources, config, psf); + break; + + case 3: + psphotApplyPSF (imdata, config, sources, psf, sky); + break; + + case 4: + psphotApplyPSF (imdata, config, sources, psf, sky); + psphotFitExtended (imdata, config, sources, sky); + break; + } + + // write out data in appropriate format + psphotOutput (imdata, config, sources, psf, sky); + psLogMsg ("psphot", 3, "wrote output: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); + exit (0); +} Index: /tags/psphot_dev_07/psphot/src/psphot.h =================================================================== --- /tags/psphot_dev_07/psphot/src/psphot.h (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphot.h (revision 5994) @@ -0,0 +1,122 @@ +# include +# include // for strcasecmp +# include // for unlink +# include +# include +# include +# include +# include +# include +# include "psLibUtils.h" +# include "psModulesUtils.h" +# include "psSparse.h" + +# define psMemCopy(A)(psMemIncrRefCounter((A))) + +typedef struct { + psImage *image; + psImage *mask; + psImage *weight; + psMetadata *header; +} eamReadout; + +// top-level psphot functions +psMetadata *psphotArguments (int *argc, char **argv); +eamReadout *psphotSetup (psMetadata *config); +bool psphotModelTest (eamReadout *imdata, psMetadata *config); +psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); +psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky); +psArray *psphotFindPeaks (eamReadout *imdata, psMetadata *config, psStats *sky); +psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); +bool psphotRoughClass (psArray *sources, psMetadata *config); +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky); +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); +void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); + +// optional object analysis steps +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotReplaceUnfit (psArray *sources); +bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); +bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf); + +// basic support functions +pmModel *pmModelCopy (pmModel *model); +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); +psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); +void psphotModelGroupInit (void); +int psphotSortBySN (const void **a, const void **b); +int psphotSortByY (const void **a, const void **b); +bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf); +void psphotTestArguments (int *argc, char **argv); + +// functions to set the correct source pixels +bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type); +bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model); +bool psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type); +bool psphotCheckRadiusEXT (eamReadout *imdata, pmSource *source, pmModel *model); +bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); +bool psphotRedefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); + +// output functions +bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); + +bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); +bool pmModelWriteEXTs (psArray *sources, char *filename); +bool pmModelWriteNULLs (psArray *sources, char *filename); +bool pmPeaksWriteText (psArray *sources, char *filename); +bool pmMomentsWriteText (psArray *sources, char *filename); + +bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image); +bool psphotDumpMoments (psMetadata *config, psArray *sources); +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); +bool psphotUpdateHeader (psMetadata *header, psMetadata *config); +int pmSourcesDophotType (pmSource *source); +bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key); + +// PSF / DBL / EXT evaluation functions +bool psphotEvalPSF (pmSource *source, pmModel *model); +bool psphotEvalDBL (pmSource *source, pmModel *model); +bool psphotEvalEXT (pmSource *source, pmModel *model); + +// functions to support simultaneous multi-source fitting +bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF); +bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF); +psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x); +bool pmModelFitSetInit (pmModelType type); + +// functions to support the source fitting process +bool psphotInitLimitsPSF (psMetadata *config); +bool psphotInitLimitsEXT (psMetadata *config, psStats *sky); +bool psphotFitBlend (eamReadout *imdata, pmSource *source); +bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources); +bool psphotFitPSF (eamReadout *imdata, pmSource *source); +pmModel *psphotFitEXT (eamReadout *imdata, pmSource *source); +psArray *psphotFitDBL (eamReadout *imdata, pmSource *source); + +// eamReadout functions +eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); + +// bicubic interpolation +psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); +psPlane psImageBicubeMin (psPolynomial2D *poly); + +// optional mode for clip fit? +psPolynomial4D *psVectorChiClipFitPolynomial4D( + psPolynomial4D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t); Index: /tags/psphot_dev_07/psphot/src/psphotApResid.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotApResid.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotApResid.c (revision 5994) @@ -0,0 +1,271 @@ +# include "psphot.h" + +// measure the aperture residual statistics +bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf) { + + int Npsf; + bool status; + pmModel *model; + pmSource *source; + + psTimerStart ("psphot"); + + // measure the aperture loss as a function of radius for PSF + float REF_RADIUS = psMetadataLookupF32 (&status, config, "PSF_REF_RADIUS"); + psf->growth = pmGrowthCurveAlloc (3.0, REF_RADIUS, 0.1); + psphotGrowthCurve (imdata, psf); + + psVector *mask = psVectorAlloc (300, PS_TYPE_U8); + psVector *xPos = psVectorAlloc (300, PS_TYPE_F64); + psVector *yPos = psVectorAlloc (300, PS_TYPE_F64); + psVector *flux = psVectorAlloc (300, PS_TYPE_F64); + psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64); + psVector *apResid = psVectorAlloc (300, PS_TYPE_F64); + psVector *dMag = psVectorAlloc (300, PS_TYPE_F64); + mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0; + Npsf = 0; + + // select all good PM_SOURCE_STAR entries + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + if (source->type != PM_SOURCE_STAR) continue; + if (source->mode & PM_SOURCE_SATSTAR) continue; + if (source->mode & PM_SOURCE_BLEND) continue; + if (source->mode & PM_SOURCE_FAIL) continue; + if (source->mode & PM_SOURCE_POOR) continue; + + // get uncorrected magnitudes in scaled apertures + model = pmSourceMagnitudes (source, NULL, 0); + if (model == NULL) continue; + + mask->data.U8[Npsf] = 0; + xPos->data.F64[Npsf] = model->params->data.F32[2]; + yPos->data.F64[Npsf] = model->params->data.F32[3]; + + flux->data.F64[Npsf] = pow(10.0, -0.4*source->fitMag); + r2rflux->data.F64[Npsf] = PS_SQR(model->radius) / flux->data.F64[Npsf]; + + apResid->data.F64[Npsf] = source->apMag + pmGrowthCurveCorrect (psf->growth, model->radius) - source->fitMag ; + + // XXX sanity clip? + // XXX need to see if all data were tossed? + // XXX need to subtract median? + if (fabs(apResid->data.F64[Npsf]) > 0.2) continue; + + dMag->data.F64[Npsf] = model->dparams->data.F32[1] / model->params->data.F32[1]; + + psVectorExtend (mask, 100, 1); + psVectorExtend (xPos, 100, 1); + psVectorExtend (yPos, 100, 1); + psVectorExtend (flux, 100, 1); + psVectorExtend (r2rflux, 100, 1); + psVectorExtend (dMag, 100, 1); + psVectorExtend (apResid, 100, 1); + Npsf ++; + } + psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects\n", psTimerMark ("psphot"), Npsf); + + // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL + char *ApTrendOption = psMetadataLookupPtr (&status, config, "APTREND"); + if (!status) ApTrendOption = psStringCopy ("SKYBIAS"); + + // 3hi/1lo sigma clipping on the rflux vs metric fit + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + stats->min = 3.0; + stats->max = 3.0; + + // constant only + if (!strcasecmp (ApTrendOption, "CONSTANT")) { + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and skybias only + if (!strcasecmp (ApTrendOption, "SKYBIAS")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + if (!strcasecmp (ApTrendOption, "SKYSAT")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and linear X,Y only + if (!strcasecmp (ApTrendOption, "XY_LIN")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_XY_LIN); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and quadratic X,Y only + if (!strcasecmp (ApTrendOption, "XY_QUAD")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_XY_QUAD); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and sky, linear X,Y only + if (!strcasecmp (ApTrendOption, "SKY_XY_LIN")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKY_XY_LIN); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and sky, linear X,Y only + if (!strcasecmp (ApTrendOption, "SKYSAT_XY_LIN")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT_XY_LIN); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + if (!strcasecmp (ApTrendOption, "ALL")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // fit just SkyBias and clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // finally, fit x, y, SkyBias and clip out objects which are too far from the median + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_ALL); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + +# if (1) + psPolynomial4D *poly = psf->ApTrend; + for (int nt = 0; nt <= poly->nT; nt++) { + for (int nz = 0; nz <= poly->nZ; nz++) { + for (int ny = 0; ny <= poly->nY; ny++) { + for (int nx = 0; nx <= poly->nX; nx++) { + if (poly->mask[nx][ny][nz][nt]) continue; + fprintf (stderr, "%d %d %d %d : %22.15g\n", nx, ny, nz, nt, poly->coeff[nx][ny][nz][nt]); + } + } + } + } +# endif + + // construct the fitted values and the residuals + psVector *fit = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux); + psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit); + +# if (0) + FILE *fout = fopen ("resid.dat", "w"); + for (int i = 0; i < resid->n; i++) { + fprintf (fout, "%d %f %f %f %f %f %f %f %f %d\n", + i, + (float) psVectorGet(xPos, i), (float) psVectorGet(yPos, i), (float) psVectorGet(r2rflux, i), (float) psVectorGet(flux, i), + (float) psVectorGet(apResid, i), (float) psVectorGet(fit, i), (float) psVectorGet(resid, i), (float) psVectorGet(dMag, i), + mask->data.U8[i]); + } + fclose (fout); +# endif + + // measure scatter for sources with dMag < 0.01 (S/N = 100) + psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); + for (int i = 0; i < dMag->n; i++) { + if (dMag->data.F64[i] > 0.01) { + mask->data.U8[i] |= 0x02; + } + } + residStats = psVectorStats (residStats, resid, NULL, mask, 0x03); + + // apply ApTrend results + psf->skyBias = psf->ApTrend->coeff[0][0][1][0]; + psf->skySat = psf->ApTrend->coeff[0][0][0][1]; + psf->ApResid = psf->ApTrend->coeff[0][0][0][0]; + psf->dApResid = residStats->sampleStdev; + psf->ApTrend->coeff[0][0][1][0] = 0; + psf->ApTrend->coeff[0][0][0][1] = 0; + psf->nApResid = residStats->clippedNvalues; + + // save results for later output + psMetadataAdd (config, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (config, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skySat); + psMetadataAdd (config, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (config, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (config, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->growth->apLoss); + psMetadataAdd (config, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid); + + psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias); + psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", + 1e3*psf->ApTrend->coeff[1][0][0][0], + 1e6*psf->ApTrend->coeff[2][0][0][0], + 1e6*psf->ApTrend->coeff[1][1][0][0], + 1e3*psf->ApTrend->coeff[0][1][0][0], + 1e6*psf->ApTrend->coeff[0][2][0][0]); + + // psFree (stats); + // psFree (mask); + // psFree (rflux); + // psFree (apResid); + + return true; +} + +/* + (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) + (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) + (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) +*/ + Index: /tags/psphot_dev_07/psphot/src/psphotApplyPSF.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotApplyPSF.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotApplyPSF.c (revision 5994) @@ -0,0 +1,61 @@ +# include "psphot.h" + +// fit psf model to all objects +bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + float x; + float y; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + psphotInitLimitsPSF (config); + psphotInitRadiusPSF (config, sky, psf->type); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->mode & PM_SOURCE_BLEND) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // use the source moments, etc to guess basic model parameters + pmModel *modelEXT = pmSourceModelGuess (source, psf->type); + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelEXT, psf); + psFree (modelEXT); + + source->modelPSF = model; + + // sets the model radius (via source->model) and adjusts pixels as needed + psphotCheckRadiusPSF (imdata, source, model); + + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, model, true); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); + + Niter += model[0].nIter; + Nfit ++; + + // check if model fit is acceptable + if (psphotEvalPSF (source, source->modelPSF)) { + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + Nsub ++; + } + } + + psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d PSF models\n", Nsub); + return (true); +} Index: /tags/psphot_dev_07/psphot/src/psphotArguments.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotArguments.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotArguments.c (revision 5994) @@ -0,0 +1,153 @@ +# include "psphot.h" + +static void usage (void) ; + +psMetadata *psphotArguments (int *argc, char **argv) { + + int N; + unsigned int Nfail; + + psMetadata *options = psMetadataAlloc (); + + // basic pslib options + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + // command-line options -- place on options MD + // mask image + if ((N = psArgumentGet (*argc, argv, "-mask"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "MASK_IMAGE", 0, "", psStringCopy(argv[N])); + psArgumentRemove (N, argc, argv); + } + + // weight image + if ((N = psArgumentGet (*argc, argv, "-weight"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "WEIGHT_IMAGE", 0, "", psStringCopy(argv[N])); + psArgumentRemove (N, argc, argv); + } + + // output residual image + if ((N = psArgumentGet (*argc, argv, "-resid"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "RESID_IMAGE", 0, "", psStringCopy(argv[N])); + psArgumentRemove (N, argc, argv); + } + + // analysis region + if ((N = psArgumentGet (*argc, argv, "-region"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", psStringCopy(argv[N])); + psArgumentRemove (N, argc, argv); + } + + // output residual image + psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", 0, "", psStringCopy("NONE")); + if ((N = psArgumentGet (*argc, argv, "-photcode"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", psStringCopy(argv[N])); + psArgumentRemove (N, argc, argv); + } + + // input PSF file + if ((N = psArgumentGet (*argc, argv, "-psf"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "PSF_INPUT_FILE", 0, "", psStringCopy (argv[N])); + psArgumentRemove (N, argc, argv); + } + + // run the 0ltest model (requires X,Y coordinate) + if ((N = psArgumentGet (*argc, argv, "-modeltest"))) { + psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT", 0, "", true); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); + + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + + // specify the modeltest model + if ((N = psArgumentGet (*argc, argv, "-model"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", psStringCopy (argv[N])); + psArgumentRemove (N, argc, argv); + } + + // specify the test fit mode + if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", psStringCopy (argv[N])); + psArgumentRemove (N, argc, argv); + } + if ((N = psArgumentGet (*argc, argv, "-fitset"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", psStringCopy (argv[N])); + psArgumentRemove (N, argc, argv); + } + } + + // other arbitrary options: -D key value (all added as string) + while ((N = psArgumentGet (*argc, argv, "-D"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + // other arbitrary options: -Df key value (all added as float) + while ((N = psArgumentGet (*argc, argv, "-Df"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + // other arbitrary options: -Di key value (all added as int) + while ((N = psArgumentGet (*argc, argv, "-Di"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + if (*argc != 4) usage (); + + // load config information + 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 file values override defaults set here + config = psMetadataConfigParse (config, &Nfail, argv[3], TRUE); + + psMetadataItem *item = NULL; + psMetadataIterator *iter = psMetadataIteratorAlloc (options, PS_LIST_HEAD, NULL); + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + psMetadataAddItem (config, item, PS_LIST_TAIL, PS_META_REPLACE); + psFree (item); + } + psFree (iter); + + // identify input image & optional weight & mask images + // command-line entries override config-file entries + psMetadataAddStr (config, PS_LIST_TAIL, "IMAGE", PS_META_REPLACE, "", argv[1]); + psMetadataAddStr (config, PS_LIST_TAIL, "OUTPUT_FILE", PS_META_REPLACE, "", argv[2]); + + return (config); +} + +static void usage (void) { + + fprintf (stderr, "USAGE: psphot (image) (output) (config)\n"); + fprintf (stderr, "options: \n"); + fprintf (stderr, " -mask (filename)\n"); + fprintf (stderr, " -weight (filename)\n"); + fprintf (stderr, " -resid (filename)\n"); + fprintf (stderr, " -photcode (photcode)\n"); + exit (2); +} Index: /tags/psphot_dev_07/psphot/src/psphotBasicDeblend.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotBasicDeblend.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotBasicDeblend.c (revision 5994) @@ -0,0 +1,119 @@ +# include "psphot.h" + +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky) { + + int N; + bool status; + float threshold; + pmSource *source, *testSource; + + int Nblend = 0; + + psTimerStart ("psphot"); + + float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION"); + if (!status) FRACTION = 0.25; + + float NSIGMA = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA"); + if (!status) NSIGMA = 5.0; + + float minThreshold = NSIGMA*sky->sampleStdev; + + // we need sources spatially-sorted to find overlaps + sources = psArraySort (sources, psphotSortByY); + + // source analysis is done in peak order (brightest first) + // we use an index for this so the spatial sorting is kept + psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32); + for (int i = 0; i < SN->n; i++) { + source = sources->data[i]; + SN->data.F32[i] = source->peak->counts; + } + psVector *index = psVectorSortIndex (NULL, SN); + // this results in an index of increasing SN + + // temporary array for overlapping objects we find + psArray *overlap = psArrayAlloc (100); + + // examine sources in decreasing SN order + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; + + if (source->mode & PM_SOURCE_BLEND) continue; + + overlap->n = 0; + + // search backwards for overlapping sources + for (int j = N - 1; j >= 0; j--) { + testSource = sources->data[j]; + if (testSource->peak->x < source->pixels->col0) continue; + if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; + if (testSource->peak->y < source->pixels->row0) break; + if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) { + fprintf (stderr, "warning: invalid condition\n"); + continue; + } + psArrayAdd (overlap, 100, testSource); + } + + // search forwards for overlapping sources + for (int j = N + 1; j < sources->n; j++) { + testSource = sources->data[j]; + if (testSource->peak->x < source->pixels->col0) continue; + if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; + if (testSource->peak->y < source->pixels->row0) { + fprintf (stderr, "warning: invalid condition\n"); + continue; + } + if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break; + psArrayAdd (overlap, 100, testSource); + } + + // generate source contour (1/4 peak counts) + if (overlap->n > 0) { + // set the threshold based on user inputs + threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; + threshold = PS_MAX (threshold, minThreshold); + + // generate a basic contour (set of x,y coordinates at-or-below flux level) + psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold); + if (contour == NULL) continue; + + // the source contour consists of two vectors, xv and yv. the contour is + // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]). both + // coordinate pairs have the same yv[] value, with xv[i] corresponding to + // the left boundary and xv[i+1] corresponding to the right boundary + + psVector *xv = contour->data[0]; + psVector *yv = contour->data[1]; + for (int k = 0; k < overlap->n; k++) { + testSource = overlap->data[k]; + if (testSource->peak->counts > source->peak->counts) continue; + for (int j = 0; j < xv->n; j+=2) { + if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; + if (xv->data.F32[j+0] > testSource->peak->x) break; + if (xv->data.F32[j+1] < testSource->peak->x) break; + + testSource->mode |= PM_SOURCE_BLEND; + + // add this to the list of source->blends + if (source->blends == NULL) { + source->blends = psArrayAlloc (16); + source->blends->n = 0; + } + psArrayAdd (source->blends, 16, testSource); + + Nblend ++; + j = xv->n; + } + } + } + } + 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: /tags/psphot_dev_07/psphot/src/psphotBlendFit.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotBlendFit.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotBlendFit.c (revision 5994) @@ -0,0 +1,45 @@ +# include "psphot.h" + +// XXX I don't like this name +bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { + + psTimerStart ("psphot"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + psphotInitLimitsPSF (config); + psphotInitLimitsEXT (config, sky); + psphotInitRadiusPSF (config, sky, psf->type); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->mode & PM_SOURCE_BLEND) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // if model is NULL, we don't have a starting guess + if (source->modelPSF == NULL) continue; + + // replace object in image + pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false); + + psTrace ("psphot.blend", 5, "trying source at %f, %f\n", source->moments->x, source->moments->y); + + if (psphotFitBlend (imdata, source)) continue; + if (psphotFitBlob (imdata, source, sources)) continue; + + psTrace ("psphot.blend", 5, "failed fits\n"); + + // re-subtract PSF for object, leave local sky + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode |= PM_SOURCE_TEMPSUB; + } + + psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects\n", psTimerMark ("psphot"), sources->n); + return (true); +} Index: /tags/psphot_dev_07/psphot/src/psphotChoosePSF.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotChoosePSF.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotChoosePSF.c (revision 5994) @@ -0,0 +1,98 @@ +# include "psphot.h" + +// try PSF models and select best option +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) { + + bool status; + char *modelName; + pmPSFtry *try = NULL; + psArray *stars = NULL; + psMetadataItem *item = NULL; + + psTimerStart ("psphot"); + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // array to store candidate PSF stars + int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS"); + if (!status) NSTARS = PS_MIN (sources->n, 200); + + stars = psArrayAlloc (sources->n); + stars->n = 0; + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { + pmSource *source = sources->data[i]; + if (source->mode & PM_SOURCE_PSFSTAR) psArrayAdd (stars, 200, source); + } + psLogMsg ("psphot.pspsf", 4, "selected candidate %d PSF objects\n", stars->n); + + // get the fixed PSF fit radius + // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + + // get the list pointers for the PSF_MODEL entries + psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + + psList *list = (psList *) mdi->data.list; + psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); + + // set up an array to store the results + psArray *models = psArrayAlloc (list->n); + + // try each model option listed in config + for (int i = 0; i < models->n; i++) { + + item = psListGetAndIncrement (iter); + modelName = item->data.V; + + models->data[i] = pmPSFtryModel (stars, modelName, RADIUS); + psFree (modelName); + psFree (item); + } + psFree (iter); + // psFree (list); XXX is list freed with iter? + psFree (stars); + + // select the best of the models + // here we are using the clippedStdev on the metric as the indicator + try = models->data[0]; + int bestN = 0; + float bestM = try->psf->dApResid; + for (int i = 1; i < models->n; i++) { + try = models->data[i]; + float M = try->psf->dApResid; + if (M < bestM) { + bestM = M; + bestN = i; + } + } + + // use the best model: + try = models->data[bestN]; + + // unset the PSFSTAR flag for stars not used for PSF model + for (int i = 0; i < try->sources->n; i++) { + pmSource *source = try->sources->data[i]; + if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) { + source->mode &= ~PM_SOURCE_PSFSTAR; + } + } + + // save only the best model; + pmPSF *psf = psMemCopy(try->psf); + psFree (models); + + psMetadataAdd (config, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "number of stars used to make PSF", psf->nPSFstars); + + modelName = pmModelGetType (psf->type); + psLogMsg ("psphot.pspsf", 3, "select psf model: %f sec\n", psTimerMark ("psphot")); + 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: /tags/psphot_dev_07/psphot/src/psphotDefinePixels.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotDefinePixels.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotDefinePixels.c (revision 5994) @@ -0,0 +1,62 @@ +# include "psphot.h" + +bool psphotDefinePixels(pmSource *mySource, + const eamReadout *imdata, + psF32 x, + psF32 y, + psF32 Radius) +{ + psRegion srcRegion; + + // Grab a subimage of the original image of size (2 * outerRadius). + srcRegion = psRegionForSquare (x, y, Radius); + srcRegion = psRegionForImage (imdata->image, srcRegion); + + mySource->pixels = psImageSubset(imdata->image, srcRegion); + mySource->weight = psImageSubset(imdata->weight, srcRegion); + mySource->mask = psImageSubset(imdata->mask, srcRegion); + mySource->region = srcRegion; + + return true; +} + +bool psphotRedefinePixels(pmSource *mySource, + const eamReadout *imdata, + psF32 x, + psF32 y, + psF32 Radius) +{ + bool extend; + psRegion newRegion; + + if (Radius == 0) return false; + + // check to see if new region is completely contained within old region + newRegion = psRegionForSquare (x, y, Radius); + newRegion = psRegionForImage (imdata->image, newRegion); + + extend = false; + extend |= (int)(newRegion.x0) < (int)(mySource->region.x0); + extend |= (int)(newRegion.x1) > (int)(mySource->region.x1); + extend |= (int)(newRegion.y0) < (int)(mySource->region.y0); + extend |= (int)(newRegion.y1) > (int)(mySource->region.y1); + + extend |= (mySource->pixels == NULL); + extend |= (mySource->weight == NULL); + extend |= (mySource->mask == NULL); + + // extend = true; + if (extend) { + // Grab a new subimage + psFree (mySource->pixels); + psFree (mySource->weight); + psFree (mySource->mask); + + mySource->pixels = psImageSubset(imdata->image, newRegion); + mySource->weight = psImageSubset(imdata->weight, newRegion); + mySource->mask = psImageSubset(imdata->mask, newRegion); + mySource->region = newRegion; + } + + return extend; +} Index: /tags/psphot_dev_07/psphot/src/psphotEnsemblePSF.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotEnsemblePSF.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotEnsemblePSF.c (revision 5994) @@ -0,0 +1,190 @@ +# include "psphot.h" + +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { + + bool status; + float x; + float y; + + psTimerStart ("psphot"); + + // source analysis is done in spatial order + sources = psArraySort (sources, psphotSortByY); + + float INNER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS"); + float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); + float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); + + // set the object surface-brightness limit for fitted pixels + float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; + psLogMsg ("psphot.ensemble", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + + // pre-calculate all model pixels + psArray *models = psArrayAlloc (sources->n); + psVector *index = psVectorAlloc (sources->n, PS_TYPE_U32); + models->n = 0; + index->n = 0; + + bool UseAnalysisRegion = false; + psRegion AnalysisRegion = {0, 0, 0, 0}; + char *region = psMetadataLookupPtr (&status, config, "ANALYSIS_REGION"); + if (status) { + UseAnalysisRegion = true; + AnalysisRegion = psRegionFromString (region); + psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", + AnalysisRegion.x0, AnalysisRegion.y0, + AnalysisRegion.x1, AnalysisRegion.y1); + } + + for (int i = 0; i < sources->n; i++) { + pmSource *inSource = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (inSource->mode & PM_SOURCE_BLEND) continue; + if (inSource->type == PM_SOURCE_DEFECT) continue; + if (inSource->type == PM_SOURCE_SATURATED) continue; + + if (UseAnalysisRegion) { + if (inSource->moments->x < AnalysisRegion.x0) continue; + if (inSource->moments->y < AnalysisRegion.y0) continue; + if (inSource->moments->x > AnalysisRegion.x1) continue; + if (inSource->moments->y > AnalysisRegion.y1) continue; + } + + pmSource *otSource = pmSourceAlloc (); + + // really saturated stars should be re-measured for a better centroid + // XXX EAM : move this to a 'clear satstar function' + if (inSource->mode & PM_SOURCE_SATSTAR) { + status = pmSourceMoments (inSource, INNER_RADIUS); + } + + // XXX EAM : add option to use EXT or PSF form + // use the source moments, etc to guess basic model parameters + pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); + if (inSource->mode & PM_SOURCE_SATSTAR) { + modelEXT->params->data.F32[2] = inSource->moments->x; + modelEXT->params->data.F32[3] = inSource->moments->y; + } else { + // peak-up on peak (for non-sat objects) + + int ix = inSource->peak->x; + int iy = inSource->peak->y; + + psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy); + psPlane min = psImageBicubeMin (bicube); + + psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", + modelEXT->params->data.F32[2], modelEXT->params->data.F32[3], min.x + ix, min.y + iy); + + // if min point is too deviant, keep the old value + if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { + modelEXT->params->data.F32[2] = min.x + ix; + modelEXT->params->data.F32[3] = min.y + iy; + } + psFree (bicube); + } + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelEXT, psf); + psFree (modelEXT); + + // save the original coords + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // get function which specifies the radius at which the model hits the given flux + pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); + + // set the fit radius based on the object flux limit and the model + model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) psAbort ("psphotEnsemblePSF", "error in radius"); + + // if needed, ask for more object pixels + psphotRedefinePixels (inSource, imdata, x, y, model->radius); + + // make temporary copies of the image pixels and mask + otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); + otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); + + // build the model image + psImage *flux = otSource->pixels; + psImage *mask = otSource->mask; + + // set model to unit peak, zero sky (we assume sky is constant) + model->params->data.F32[0] = 0.0; + model->params->data.F32[1] = 1.0; + + // fill in the model pixel values + psImageInit (flux, 0.0); + psImageKeepCircle (mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceAddModel (flux, mask, model, false, false); + + // save source in list + otSource->modelPSF = model; + index->data.U32[models->n] = i; + psArrayAdd (models, 100, otSource); + } + psLogMsg ("psphot.emsemble", 4, "built models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n); + + float f; + + // fill out the sparse matrix + psSparse *sparse = psSparseAlloc (models->n, 100); + for (int i = 0; i < models->n; i++) { + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; + pmSource *Mi = models->data[i]; + + // diagonal element (auto-cross-product) + f = pmSourceCrossProduct (Mi, Mi); + psSparseMatrixElement (sparse, i, i, f); + + // find the image x model value + f = pmSourceCrossProduct (Fi, Mi); + psSparseVectorElement (sparse, i, f); + + // loop over all other stars following this one + for (int j = i + 1; j < models->n; j++) { + pmSource *Mj = models->data[j]; + + // skip over disjoint source images, break after last possible overlap + if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; + if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; + if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; + if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; + + // got an overlap; calculate cross-product and add to output array + f = pmSourceCrossProduct (Mi, Mj); + psSparseMatrixElement (sparse, j, i, f); + } + } + psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); + + // solve for normalization terms (need include local sky) + psSparseResort (sparse); + psVector *norm = psSparseSolve (NULL, sparse, 3); + + // adjust models, set sources and subtract + for (int i = 0; i < models->n; i++) { + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; + pmSource *Mi = models->data[i]; + + Fi->modelPSF = Mi->modelPSF; + + // assign linearly-fitted normalization + Fi->modelPSF->params->data.F32[1] = norm->data.F32[i]; + + // subtract object + pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); + Fi->mode |= PM_SOURCE_SUBTRACTED; + Fi->mode |= PM_SOURCE_TEMPSUB; + } + + // XXX EAM : need to free up many things here + + psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot")); + return true; +} Index: /tags/psphot_dev_07/psphot/src/psphotEvalFLT.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotEvalFLT.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotEvalFLT.c (revision 5994) @@ -0,0 +1,48 @@ +# include "psphot.h" + +bool psphotEvalEXT (pmSource *source, pmModel *model) +{ + int keep; + + // do we actually have a valid EXT model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + psLogMsg ("psphot", 5, "EXT fail fit\n"); + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // unless we prove otherwise, this object is extended + source->type = PM_SOURCE_EXTENDED; + + // the following source->mode information pertains to modelEXT: + source->mode |= PM_SOURCE_EXTMODEL; + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] <= 0) { + source->mode |= PM_SOURCE_FAIL; + return false; + } + + keep = pmModelFitStatus (model); + if (keep) return true; + + // poor-quality fit; only keep if nothing else works... + psLogMsg ("psphot", 5, "EXT poor fit\n"); + source->mode |= PM_SOURCE_POOR; + return false; +} Index: /tags/psphot_dev_07/psphot/src/psphotEvalPSF.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotEvalPSF.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotEvalPSF.c (revision 5994) @@ -0,0 +1,207 @@ +# include "psphot.h" + +// given a pmSource which has been fitted using modelPSF, evaluate the +// resulting fit: did the fit succeed? is this object PSF-like? is this object +// extended? return status is TRUE for a valid PSF, FALSE otherwise. + +// 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) + +// sigma_x : 1 / SX +// dsx : 1 / (SX * SN) +// dsx_o = hypot (1/(SX*SN), MIN_DSX) + +// any objects which is consistent with the PSF should have +// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) +// this includes a minimum buffer (DS) for the brighter objects + +// saturated stars should fall outside (larger), but have peaks above SATURATION +// extended sources should be larger, cosmic rays smaller +// we also reject objects with S/N too low or ChiSquare to high + +// floor for DS value +// XXX EAM : add to configuration? +# define MIN_DS 0.01 + +static float SATURATION; +static float PSF_MIN_SN; +static float PSF_MAX_CHI; +static float PSF_SHAPE_NSIGMA; + +bool psphotInitLimitsPSF (psMetadata *config) { + + bool status; + + // we may set this differently here from the value used to mark likely saturated stars + SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); + PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); + PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + + return true; +} + +// examine the model->status, fit parameters, etc and decide if the model succeeded +// set the source->type and source->mode appropriately +bool psphotEvalPSF (pmSource *source, pmModel *model) { + + int keep; + float dSX, dSY, SX, SY, SN; + float nSx, nSy, Chi; + + // do we actually have a valid PSF model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_STAR; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_PSFMODEL; + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests + // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... + if (model->params->data.F32[1] >= SATURATION) { + if (source->mode & PM_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->mode |= PM_SOURCE_SATSTAR; + return true; + } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] <= 0) { + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->mode & PM_SOURCE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_SATSTAR; + } + + SN = model->params->data.F32[1]/model->dparams->data.F32[1]; + SX = model->params->data.F32[4]; + SY = model->params->data.F32[5]; + dSX = model->dparams->data.F32[4]; + dSY = model->dparams->data.F32[5]; + Chi = model->chisq / model->nDOF; + + // swing of sigma_x,y in sigmas + nSx = dSX / hypot (MIN_DS, 1 / (SX * SN)); + nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); + + // assign PM_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] + keep = TRUE; + keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); + keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); + keep &= (SN > PSF_MIN_SN); + keep &= (Chi < PSF_MAX_CHI); + if (keep) return true; + + // this source is not a star, warn if it was a PSFSTAR + if (source->mode & PM_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality (%f, %f : %f %f %f %f)\n", + model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi); + } else { + psLogMsg ("psphot", 5, "fails PSF fit (%f, %f : %f %f %f %f)\n", + model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi); + } + + // object appears to be small, suspected defect + if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_DEFECT; + return false; + } + + // object appears to be large, suspected extended source + if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_EXTENDED; + return false; + } + + // poor-quality fit; only keep if nothing else works... + source->mode |= PM_SOURCE_POOR; + return false; +} + +// examine the model->status, fit parameters, etc and decide if the model succeeded +// set the source->type and source->mode appropriately +bool psphotEvalDBL (pmSource *source, pmModel *model) { + + // do we actually have a valid PSF model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_STAR; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_PSFMODEL; + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests + // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... + if (model->params->data.F32[1] >= SATURATION) { + if (source->mode & PM_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->mode |= PM_SOURCE_SATSTAR; + return true; + } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] <= 0) { + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->mode & PM_SOURCE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_SATSTAR; + } + return true; +} Index: /tags/psphot_dev_07/psphot/src/psphotFindPeaks.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotFindPeaks.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotFindPeaks.c (revision 5994) @@ -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: /tags/psphot_dev_07/psphot/src/psphotFitGalaxies.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotFitGalaxies.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotFitGalaxies.c (revision 5994) @@ -0,0 +1,66 @@ +# include "psphot.h" + +bool psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) +{ + bool status; + float x; + float y; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + float EXT_MIN_SN = psMetadataLookupF32 (&status, config, "EXT_MIN_SN"); + float EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS"); + char *modelName = psMetadataLookupPtr (&status, config, "EXT_MODEL"); + pmModelType modelType = pmModelSetType (modelName); + + psphotInitRadiusEXT (config, skyStats, modelType); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // only fit suspected extended sources + if (source->type != PM_SOURCE_EXTENDED) continue; + if (source->mode & PM_SOURCE_BLEND) continue; + + // skip possible extended sources below fit threshold + if (source->moments->SN < EXT_MIN_SN) continue; + + // recalculate the source moments using the larger extended-source moments radius + status = pmSourceMoments (source, EXT_MOMENTS_RAD); + if (!status) { + source->mode |= PM_SOURCE_FAIL; // better choice? + continue; + } + + // use the source moments, etc to guess basic model parameters + pmModel *model = pmSourceModelGuess (source, modelType); + source->modelEXT = model; + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // sets the model radius (via source->model) and adjusts pixels as needed + psphotCheckRadiusEXT (imdata, source, model); + + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, model, false); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); + + Niter += model[0].nIter; + Nfit++; + + // check if model fit is acceptable + if (pmModelFitStatus (model)) { + pmSourceSubModel (source->pixels, source->mask, model, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + Nsub ++; + } + } + psLogMsg ("psphot", 3, "fit EXT models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d EXT objects\n", Nsub); + return (true); +} Index: /tags/psphot_dev_07/psphot/src/psphotFullFit.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotFullFit.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotFullFit.c (revision 5994) @@ -0,0 +1,137 @@ +# include "psphot.h" + +bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool status; + float x, y; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // extended source model parameters + float EXT_MIN_SN = psMetadataLookupF32 (&status, config, "EXT_MIN_SN"); + float EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS"); + + // extended source model descriptions + char *modelNameEXT = psMetadataLookupPtr (&status, config, "EXT_MODEL"); + pmModelType modelTypeEXT = pmModelSetType (modelNameEXT); + + psphotInitLimitsPSF (config); + psphotInitRadiusPSF (config, sky, psf->type); + psphotInitRadiusEXT (config, sky, modelTypeEXT); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (source->mode & PM_SOURCE_BLEND) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // save the PSF model from the Ensemble fit + pmModel *LIN = pmModelCopy (source->modelPSF); + + // attempt to fit the PSF model + pmModel *modelPSF = source->modelPSF; + + // replace object in image + pmSourceAddModel (source->pixels, source->mask, modelPSF, false, false); + source->mode &= ~PM_SOURCE_SUBTRACTED; + + psphotCheckRadiusPSF (imdata, source, modelPSF); + + x = modelPSF->params->data.F32[2]; + y = modelPSF->params->data.F32[3]; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, modelPSF->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, modelPSF, true); + psImageKeepCircle (source->mask, x, y, modelPSF->radius, "AND", ~PSPHOT_MASK_MARKED); + + // track model evaluations + Niter += modelPSF[0].nIter; + Nfit ++; + + // does the PSF model succeed? + if (psphotEvalPSF (source, source->modelPSF)) { + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + psFree (LIN); + Nsub ++; + continue; + } + + // use the moments to get the extended source angle + // displace on either side of the main peak by 1sigma or so + // cut the peak in half? use the main peak? + // also try fit set on (unsaturated?) objects with a blend within 2-3 sigma + // allow all possible nearby blends to be fit? + // EllipseMoments moments; + // moments.x2 = source->moments->Sx; + // moments.y2 = source->moments->Sy; + // moments.xy = source->moments->Sxy; + // EllipseAxes axes = EllipseMomentsToAxes (moments); + + // skip the source if we don't think it is extended + if (source->type != PM_SOURCE_EXTENDED) goto subLINEAR; + if (source->moments->SN < EXT_MIN_SN) goto subLINEAR; + + // add the double-PSF fit mode + // how do we compare the results? chisq? + + // recalculate the source moments using the larger extended-source moments radius + if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) goto subLINEAR; + + // use the source moments, etc to guess basic model parameters + source->modelEXT = pmSourceModelGuess (source, modelTypeEXT); + pmModel *modelEXT = source->modelEXT; + + psphotCheckRadiusEXT (imdata, source, modelEXT); + + x = modelEXT->params->data.F32[2]; + y = modelEXT->params->data.F32[3]; + + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, modelEXT->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, modelEXT, false); + psImageKeepCircle (source->mask, x, y, modelEXT->radius, "AND", ~PSPHOT_MASK_MARKED); + + // track model evaluations + Niter += modelEXT[0].nIter; + Nfit++; + + // does the EXT model succeed? + if (psphotEvalEXT (source, source->modelEXT)) { + pmSourceSubModel (source->pixels, source->mask, source->modelEXT, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + psFree (LIN); + Nsub ++; + continue; + } + + // re-subtract PSF for object, leave local sky + subLINEAR: + psFree (source->modelPSF); + source->modelPSF = LIN; + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode |= PM_SOURCE_TEMPSUB; + Nsub ++; + } + + psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub); + return (true); +} + + +// XXX EAM : can we distinguish the outcomes enough to re-add poor fits, etc? Index: /tags/psphot_dev_07/psphot/src/psphotGrowthCurve.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotGrowthCurve.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotGrowthCurve.c (revision 5994) @@ -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: /tags/psphot_dev_07/psphot/src/psphotImageBackground.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotImageBackground.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotImageBackground.c (revision 5994) @@ -0,0 +1,65 @@ +# include "psphot.h" + +psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky) +{ + + // fit background to subset of image points within 3 sigma of sky + psTimerStart ("psphot"); + + psImage *image = imdata->image; + psImage *mask = imdata->mask; + + bool status = false; + int Npixels = image->numRows*image->numCols; + int Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "IMSTATS_NPIX")); + psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); + + psVector *x = psVectorAlloc (Nsubset, PS_TYPE_F32); + psVector *y = psVectorAlloc (Nsubset, PS_TYPE_F32); + psVector *z = psVectorAlloc (Nsubset, PS_TYPE_F32); + + float Sky = sky->sampleMean; + float dSky = 3*sky->sampleStdev; + + // choose Nsubset points between 0 and Nx*Ny, convert to coords + int Nout = 0; + Npixels = image->numRows*image->numCols; + for (int i = 0; (i < 3*Nsubset) && (Nout < Nsubset); i++) { + double frnd = psRandomUniform (rnd); + int pixel = Npixels * frnd; + int ix = pixel / image->numCols; + int iy = pixel % image->numCols; + float value = image->data.F32[iy][ix]; + if (fabs(value - Sky) > dSky) continue; + z->data.F32[Nout] = value; + x->data.F32[Nout] = ix; + y->data.F32[Nout] = iy; + Nout ++; + } + x->n = y->n = z->n = Nout; + + psPolynomial2D *skyModel = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); + skyModel->mask[1][1] = 1; + + // psVector *mask = psVectorAlloc (Nout, PS_TYPE_U8); + // psVectorInit (mask, 0); + // psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + // skyModel = psVectorFitPolynomial2D (skyModel, mask, 1, z, NULL, x, y); + + skyModel = psVectorFitPolynomial2D (skyModel, NULL, 0, z, NULL, x, y); + + // this is a very inefficient way to evaluate the function.. + // lame, temporary step to zero masked pixels + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + if (!mask->data.U8[j][i]) { + Sky = psPolynomial2DEval (skyModel, i, j); + image->data.F32[j][i] -= Sky; + } + } + } + sky->sampleMean = 0.0; + + psLogMsg ("psphot", 3, "fit background model: %f sec\n", psTimerMark ("psphot")); + return (skyModel); +} Index: /tags/psphot_dev_07/psphot/src/psphotImageStats.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotImageStats.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotImageStats.c (revision 5994) @@ -0,0 +1,56 @@ +# include "psphot.h" + +psStats *psphotImageStats (eamReadout *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) + psTimerStart ("psphot"); + + // this section fixes a bug in the psImageStats definition: + // it should limit the number of elements used for the initial + // pass on the stats + { + bool status = false; + int Npixels = image->numRows*image->numCols; + int Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "IMSTATS_NPIX")); + psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32); + psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); + + // 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]; + } + + // this should use ROBUST not SAMPLE median + // robust median is broken in pslib (0.5) + stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + stats = psVectorStats (stats, subset, NULL, NULL, 0); // use a mask? + psFree (subset); + psFree (rnd); + } + + // we store these values in mean,stdev + bool status = false; + float NOISE = psMetadataLookupF32 (&status, config, "RDNOISE"); + float GAIN = psMetadataLookupF32 (&status, config, "GAIN"); + + // convert instrumental background to poisson stats + sky = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + sky->sampleMean = stats->sampleMedian; + sky->sampleStdev = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE)); + + psLogMsg ("psphot", 4, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev); + psLogMsg ("psphot", 3, "image stats: %f sec\n", psTimerMark ("psphot")); + + psFree (stats); + return (sky); +} Index: /tags/psphot_dev_07/psphot/src/psphotMagnitudes.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotMagnitudes.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotMagnitudes.c (revision 5994) @@ -0,0 +1,72 @@ +# include "psphot.h" + +// XXX EAM : the apMag should only be calculated for the brighter sources? +// XXX EAM : SN limit set by user? +// XXX EAM : masked region should be (optionally) elliptical +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) { + + int status; + bool isPSF; + float x, y; + float rflux; + float radius; + pmModel *model; + + switch (source->type) { + case PM_SOURCE_STAR: + model = source->modelPSF; + if (model == NULL) return NULL; + radius = (apRadius > 0) ? apRadius : model->radius; + isPSF = true; + break; + + case PM_SOURCE_EXTENDED: + model = source->modelEXT; + if (model == NULL) return NULL; + radius = model->radius; + isPSF = false; + break; + + default: + return NULL; + } + + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // replace source flux + pmSourceAddModel (source->pixels, source->mask, model, false, false); + + // set aperture mask circle of PSF_FIT_RADIUS + psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED); + + // measure object photometry + status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask); + + // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + if (isPSF && (psf != NULL)) { + if (psf->growth != NULL) { + source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius); + } + + rflux = pow (10.0, 0.4*source->fitMag); + source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux; + source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); + } + + // unmask aperture + psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + if (!status) return NULL; + return model; +} + +/* +aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) + +*/ Index: /tags/psphot_dev_07/psphot/src/psphotModelGroupInit.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotModelGroupInit.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotModelGroupInit.c (revision 5994) @@ -0,0 +1,20 @@ +# include "psphot.h" + +# include "models/pmModel_TAUSS.c" + +static pmModelGroup userModels[] = { + {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS, pmModelFlux_TAUSS, pmModelRadius_TAUSS, pmModelLimits_TAUSS, pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS}, +}; + + +void psphotModelGroupInit (void) +{ + + pmModelGroupInit (); + + int Nmodels = sizeof (userModels) / sizeof (pmModelGroup); + for (int i = 0; i < Nmodels; i++) { + pmModelGroupAdd (&userModels[i]); + } + return; +} Index: /tags/psphot_dev_07/psphot/src/psphotModelTest.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotModelTest.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotModelTest.c (revision 5994) @@ -0,0 +1,175 @@ +# include "psphot.h" +# include "psEllipse.h" + +bool psphotModelTest (eamReadout *imdata, psMetadata *config) { + + bool status; + int modelType; + unsigned int Nfail; + float obsMag, fitMag, value; + char name[64]; + pmPSF *psf = NULL; + + psMetadataItem *item = NULL; + + // run model fitting tests on a single source + if (!psMetadataLookupBool (&status, config, "TEST_FIT")) return false; + + // what fitting mode to use? + char *psfModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE"); + if (!status) { + psfModeWord = psStringCopy ("EXT"); + } + bool psfMode = !strcasecmp (psfModeWord, "PSF"); + + // in psfMode, psf sets the model type + if (psfMode) { + char *psfFile = psMetadataLookupPtr (&status, config, "PSF_INPUT_FILE"); + if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied"); + psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE); + psf = pmPSFfromMD (psfData); + modelType = psf->type; + } else { + // find the model: supplied by user or first in the PSF_MODEL list + char *modelName = psMetadataLookupPtr (&status, config, "TEST_FIT_MODEL"); + if (modelName == NULL) { + // get the list pointers for the PSF_MODEL entries + psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + + // take the first list element + psList *list = (psList *) mdi->data.list; + item = psListGet (list, PS_LIST_HEAD); + modelName = item->data.V; + } + modelType = pmModelSetType (modelName); + if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName); + } + + // find the fitting parameters (try test values first) + float INNER = psMetadataLookupF32 (&status, config, "TEST_FIT_INNER_RADIUS"); + if (!status) { + INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS"); + } + + float OUTER = psMetadataLookupF32 (&status, config, "TEST_FIT_OUTER_RADIUS"); + if (!status) { + OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + } + + float RADIUS = psMetadataLookupF32 (&status, config, "TEST_FIT_RADIUS"); + if (!status) { + RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + } + + float mRADIUS = psMetadataLookupF32 (&status, config, "TEST_MOMENTS_RADIUS"); + if (!status) { + mRADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS"); + } + + // define the source of interest + float xObj = psMetadataLookupF32 (&status, config, "TEST_FIT_X"); + float yObj = psMetadataLookupF32 (&status, config, "TEST_FIT_Y"); + + // construct the source structures + pmSource *source = pmSourceAlloc(); + source->peak = pmPeakAlloc (xObj, yObj, 0, 0); + psphotDefinePixels (source, imdata, xObj, yObj, OUTER); + + // find the local sky + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER); + if (!status) psAbort ("fitsource", "pmSourceLocalSky error"); + + // get the source moments + status = pmSourceMoments (source, mRADIUS); + if (!status) psAbort ("fitsource", "pmSourceLocalSky error"); + source->peak->counts = source->moments->Peak; + + fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y); + + EllipseMoments moments; + moments.x2 = source->moments->Sx; + moments.y2 = source->moments->Sy; + moments.xy = source->moments->Sxy; + EllipseAxes axes = EllipseMomentsToAxes (moments); + + fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + + // get the initial model parameter guess + pmModel *model = pmSourceModelGuess (source, modelType); + // if any parameters are defined, use those values + int nParams = pmModelParameterCount (modelType); + psF32 *params = model->params->data.F32; + for (int i = 0; i < nParams; i++) { + if (i == 2) { + params[i] = xObj; + continue; + } + if (i == 3) { + params[i] = yObj; + continue; + } + sprintf (name, "TEST_FIT_PAR%d", i); + value = psMetadataLookupF32 (&status, config, name); + if (status) { + params[i] = value; + } + } + + float area = params[4]*params[5]; + fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y); + + if (psfMode) { + pmModel *modelPSF = pmModelFromPSF (model, psf); + psFree (model); + model = modelPSF; + params = model->params->data.F32; + } + + // list model input shape + EllipseShape shape; + shape.sx = 1.4 / model->params->data.F32[4]; + shape.sy = 1.4 / model->params->data.F32[5]; + shape.sxy = model->params->data.F32[6]; + axes = EllipseShapeToAxes (shape); + + fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + + + fprintf (stderr, "input parameters: \n"); + for (int i = 0; i < nParams; i++) { + fprintf (stderr, "%d : %f\n", i, params[i]); + } + + // define the pixels used for the fit + psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED); + + char *fitset = psMetadataLookupPtr (&status, config, "TEST_FIT_SET"); + if (status) { + status = psphotFitSet (source, model, fitset, psfMode); + return status; + } + + status = pmSourceFitModel (source, model, psfMode); + + // measure the source mags + pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask); + fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag); + + // write out positive object + psphotSaveImage (NULL, source->pixels, "object.fits"); + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + fprintf (stderr, "output parameters: \n"); + for (int i = 0; i < nParams; i++) { + fprintf (stderr, "%d : %f\n", i, params[i]); + } + + // write out + psphotSaveImage (NULL, source->pixels, "resid.fits"); + psphotSaveImage (NULL, source->mask, "mask.fits"); + + exit (0); +} Index: /tags/psphot_dev_07/psphot/src/psphotOutput.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotOutput.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotOutput.c (revision 5994) @@ -0,0 +1,658 @@ +# include "psphot.h" + +// output functions: we have several fixed modes (sx, obj, cmp) +void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { + + bool status; + + psTimerStart ("psphot"); + + char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE"); + char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE"); + char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE"); + char *psfFile = psMetadataLookupPtr (&status, config, "PSF_OUTPUT_FILE"); + + if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage); + psphotSamplePSFs (config, psf, imdata->image); + + if (psfFile != NULL) { + psMetadata *psfData = pmPSFtoMD (NULL, psf); + psMetadataConfigWrite (psfData, psfFile); + psFree (psfData); + } + if (outputFile == NULL) { + psLogMsg ("output", 3, "no data output file selected"); + return; + } + if (outputMode == NULL) { + psLogMsg ("output", 3, "no data output mode selected"); + return; + } + if (!strcasecmp (outputMode, "OBJ")) { + pmSourcesWriteOBJ (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "SX")) { + pmSourcesWriteSX (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "CMP")) { + pmSourcesWriteCMP (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "CMF")) { + pmSourcesWriteCMF (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "TEXT")) { + pmSourcesWriteText (imdata, config, outputFile, sources, psf, sky); + return; + } + + psAbort ("psphot", "unknown output mode %s", outputMode); +} + +// dophot-style output list with fixed line width +bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { + + int type; + pmModel *model; + psF32 *PAR, *dPAR; + float dmag, apResid; + bool status; + + psLine *line = psLineAlloc (104); // 104 is dophot-defined line length + + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); + + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + model = pmSourceMagnitudes (source, psf, RADIUS); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + apResid = source->apMag - source->fitMag; + + psLineInit (line); + psLineAdd (line, "%3d", type); + psLineAdd (line, "%8.2f", PAR[2]); + psLineAdd (line, "%8.2f", PAR[3]); + psLineAdd (line, "%8.3f", source->fitMag); + psLineAdd (line, "%6.3f", dmag); + psLineAdd (line, "%9.2f", PAR[0]); + psLineAdd (line, "%9.3f", PAR[4]); + psLineAdd (line, "%9.3f", PAR[5]); + psLineAdd (line, "%7.2f", PAR[6]); + psLineAdd (line, "%8.3f", 99.999); + psLineAdd (line, "%8.3f", source->apMag); + psLineAdd (line, "%8.2f\n", apResid); + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +// elixir-mode / sextractor-style output list with fixed line width +bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { + + pmModel *model; + psF32 *PAR, *dPAR; + float dmag; + bool status; + + psLine *line = psLineAlloc (110); // 110 is sextractor line length + + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("WriteSourceSX", 3, "can't open output file for output %s\n", filename); + return false; + } + + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); + + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + model = pmSourceMagnitudes (source, psf, RADIUS); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // pmSourceSextractType (source, &type, &flags); XXX EAM : implement me... + + dmag = dPAR[1] / PAR[1]; + + psLineInit (line); + psLineAdd (line, "%5.2f", 0.0); // should be type + psLineAdd (line, "%11.3f", PAR[2]); + psLineAdd (line, "%11.3f", PAR[3]); + psLineAdd (line, "%9.4f", source->fitMag); + psLineAdd (line, "%9.4f", dmag); + psLineAdd (line, "%13.4f", PAR[0]); + psLineAdd (line, "%9.2f", 0.0); // should be FWHMx + psLineAdd (line, "%9.2f", 0.0); // should be FWHMy + psLineAdd (line, "%6.1f", 0.0); // should be Theta + psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO + psLineAdd (line, "%9.4f", source->apMag); + psLineAdd (line, "%4d\n", 0); // should be flags + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +// elixir-style pseudo FITS table (header + ascii list) +bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { + + int i, type; + pmModel *model; + psMetadataItem *mdi; + psF32 *PAR, *dPAR; + float dmag, lsky; + bool status; + + // find config information for output header + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); + float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); + if (!status) ZERO_POINT = 25.0; + + // write necessary information to output header + psphotUpdateHeader (imdata->header, config); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); + + // create file, write-out header + psFits *fits = psFitsOpen (filename, "w"); + + // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 + mdi = psMetadataLookup (imdata->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // psMetadataAdd (imdata->header, PS_LIST_HEAD, "NAXIS", PS_DATA_S32 | PS_META_REPLACE, "", 0); + psFitsWriteHeader (imdata->header, fits); + psFitsClose (fits); + + // re-open, add data to end of file + FILE *f = fopen (filename, "a+"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + fseek (f, 0, SEEK_END); + + psLine *line = psLineAlloc (67); // 66 is imclean-defined line length + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + model = pmSourceMagnitudes (source, psf, RADIUS); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + psLineInit (line); + psLineAdd (line, "%6.1f ", PAR[2]); + psLineAdd (line, "%6.1f ", PAR[3]); + psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT)); + psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); + psLineAdd (line, "%2d ", type); + psLineAdd (line, "%3.1f ", lsky); + psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal + psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); + psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' + psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' + psLineAdd (line, "%5.1f\n", 0); // should be theta + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +// elixir-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data +bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { + + psArray *table; + psMetadataItem *mdi; + psMetadata *row; + psMetadata *theader; + int i, type; + pmModel *model; + psF32 *PAR, *dPAR; + float dmag, lsky; + bool status; + + // find config information for output header + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); + float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); + + // write necessary information to output header + psphotUpdateHeader (imdata->header, config); + + table = psArrayAlloc (sources->n); + table->n = 0; + + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + model = pmSourceMagnitudes (source, psf, RADIUS); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", type); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", lsky); + psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", source->apMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_STRING, "", "0"); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "123"); + + psArrayAdd (table, 100, row); + psFree (row); + } + + // write the correct number of sources in the header + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true); + + // set NAXIS to 0 (we don't write out the data array) + mdi = psMetadataLookup (imdata->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create the basic table header + theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); + + // write out the header and table + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteHeader (imdata->header, fits); + psFitsWriteTable (fits, theader, table); + psFitsClose (fits); + + return true; +} + +/***** Text Output Methods *****/ + +bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky) { + + char *name = (char *) psAlloc (strlen(filename) + 10); + + sprintf (name, "%s.psf.dat", filename); + pmModelWritePSFs (sources, config, name, psf); + + sprintf (name, "%s.ext.dat", filename); + pmModelWriteEXTs (sources, name); + + sprintf (name, "%s.nul.dat", filename); + pmModelWriteNULLs (sources, name); + + sprintf (name, "%s.mnt.dat", filename); + pmMomentsWriteText (sources, name); + + psFree (name); + return true; +} + +// write the PSF sources to an output file +bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf) { + + double dPos, dMag; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + bool status; + + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename); + return false; + } + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + if (source->type != PM_SOURCE_STAR) continue; + model = pmSourceMagnitudes (source, psf, RADIUS); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // dPos is positional error, dMag is mag error + dPos = hypot (dPAR[2], dPAR[3]); + dMag = dPAR[1] / PAR[1]; + + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", PAR[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dPAR[j]); + } + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + source[0].apMag, source[0].type, source[0].mode, + log10(model[0].chisq/model[0].nDOF), + source[0].moments->SN, + model[0].radius, + model[0].nDOF, + model[0].nIter); + } + fclose (f); + return true; +} + +// dump the sources to an output file +bool pmModelWriteEXTs (psArray *sources, char *filename) { + + double dPos, dMag; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename); + return false; + } + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + if (source->type != PM_SOURCE_EXTENDED) continue; + model = pmSourceMagnitudes (source, NULL, 0.0); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // dPos is shape error + // XXX these are hardwired for SGAUSS + dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5])); + dMag = dPAR[1] / PAR[1]; + + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", PAR[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dPAR[j]); + } + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + source->apMag, + source[0].type, source[0].mode, + log10(model[0].chisq/model[0].nDOF), + source[0].moments->SN, + model[0].radius, + model[0].nDOF, + model[0].nIter); + } + fclose (f); + return true; +} + +// dump the sources to an output file +bool pmModelWriteNULLs (psArray *sources, char *filename) +{ + + int i; + FILE *f; + pmMoments *moment = NULL; + pmSource *source = NULL; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); + return false; + } + + pmMoments *empty = pmMomentsAlloc (); + + // write sources with models first + for (i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // skip these sources (in PSF or EXT) + if (source->type == PM_SOURCE_STAR) continue; + if (source->type == PM_SOURCE_EXTENDED) continue; + + if (source->moments == NULL) { + moment = empty; + } else { + moment = source->moments; + } + + fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", + source->peak->x, source->peak->y, source->peak->counts, + source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, + source->moments->Sum, source->moments->Peak, + source->moments->Sky, source->moments->SN, + source->moments->nPixels, source->type); + } + fclose (f); + return true; +} + +// write the moments to an output file +bool pmMomentsWriteText (psArray *sources, char *filename) +{ + + int i; + FILE *f; + pmSource *source = NULL; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + source = sources->data[i]; + if (source->moments == NULL) + continue; + fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", + source->peak->x, source->peak->y, source->peak->counts, + source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, + source->moments->Sum, source->moments->Peak, + source->moments->Sky, source->moments->SN, + source->moments->nPixels, source->type, source->mode); + } + fclose (f); + return true; +} + +// write the peaks to an output file +bool pmPeaksWriteText (psArray *peaks, char *filename) { + + int i; + FILE *f; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); + return false; + } + + for (i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + if (peak == NULL) continue; + fprintf (f, "%5d %5d %7.1f\n", + peak->x, peak->y, peak->counts); + } + fclose (f); + return true; +} + +// translations between psphot object types and dophot object types +int pmSourceDophotType (pmSource *source) { + + switch (source->type) { + + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + return (8); + + case PM_SOURCE_STAR: + if (source->mode & PM_SOURCE_SATSTAR) return (10); + if (source->mode & PM_SOURCE_POOR) return (7); + if (source->mode & PM_SOURCE_FAIL) return (4); + return (1); + + case PM_SOURCE_EXTENDED: + return (2); + + default: + return (0); + } + return (0); +} + +int psphotSaveImage (psMetadata *header, psImage *image, char *filename) { + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0); + psFitsClose (fits); + return (TRUE); +} + +psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) { + + psRegion region = {x - dx, x + dx, y - dy, y + dy}; + psImage *view = psImageSubset (image, region); + psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32); + psImageInit (sample, 0); + modelEXT->params->data.F32[2] = x; + modelEXT->params->data.F32[3] = y; + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); + pmSourceAddModel (sample, NULL, modelPSF, false, false); + psFree (modelPSF); + return (sample); +} + + +bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image) { + + bool status; + + // make sample PSFs for 4 corners and the center + psImage *sample; + + // optional dump of all rough source data + char *output = psMetadataLookupPtr (&status, config, "PSF_SAMPLE_FILE"); + if (!status) return false; + if (output == NULL) return false; + if (output[0] == 0) return false; + + pmModel *modelEXT = pmModelAlloc (psf->type); + modelEXT->params->data.F32[0] = 0; + modelEXT->params->data.F32[1] = 1; + + psFits *fits = psFitsOpen (output, "w"); + + sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + + psFitsClose (fits); + return (TRUE); +} + +bool psphotDumpMoments (psMetadata *config, psArray *sources) { + + bool status; + + // optional dump of all rough source data + char *output = psMetadataLookupPtr (&status, config, "MOMENTS_OUTPUT_FILE"); + if (!status) return false; + if (output == NULL) return false; + if (output[0] == 0) return false; + + pmMomentsWriteText (sources, output); + psFree (output); + return true; +} + +bool psphotUpdateHeader (psMetadata *header, psMetadata *config) { + + // write necessary information to output header + psMetadataItemTransfer (header, config, "NSTARS"); + psMetadataItemTransfer (header, config, "ZERO_PT"); + psMetadataItemTransfer (header, config, "APMIFIT"); + psMetadataItemTransfer (header, config, "dAPMIFIT"); + psMetadataItemTransfer (header, config, "SKYBIAS"); + psMetadataItemTransfer (header, config, "PHOTCODE"); + + // write necessary information to output header + psMetadataItemTransfer (header, config, "NPSFSTAR"); + psMetadataItemTransfer (header, config, "SKYBIAS"); + psMetadataItemTransfer (header, config, "SKYSAT"); + psMetadataItemTransfer (header, config, "APMIFIT"); + psMetadataItemTransfer (header, config, "DAPMIFIT"); + psMetadataItemTransfer (header, config, "NAPMIFIT"); + psMetadataItemTransfer (header, config, "APLOSS"); + + // XXX these need to be defined from elsewhere + psMetadataAdd (header, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + + return true; +} Index: /tags/psphot_dev_07/psphot/src/psphotRadiusChecks.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotRadiusChecks.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotRadiusChecks.c (revision 5994) @@ -0,0 +1,65 @@ +# include "psphot.h" + +static float PSF_OUTER_RADIUS; +static float PSF_FIT_PADDING; +static float PSF_FLUX_LIMIT; +static pmModelRadius modelRadiusPSF; + +bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type) { + + bool status; + + float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); + PSF_OUTER_RADIUS = psMetadataLookupF32 (&status, config, "PSF_OUTER_RADIUS"); + PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); + PSF_FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; + + // this function specifies the radius at this the model hits the given flux + modelRadiusPSF = pmModelRadius_GetFunction (type); + return true; +} + +bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model) { + + // set the fit radius based on the object flux limit and the model + // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) + model->radius = modelRadiusPSF (model->params, PSF_FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius"); + + if (source->mode & PM_SOURCE_SATSTAR) { + model->radius *= 2; + } + + bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); + return status; +} + +static float EXT_OUTER_RADIUS; +static float EXT_FIT_PADDING; +static float EXT_FLUX_LIMIT; +static pmModelRadius modelRadiusEXT; + +bool psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type) { + + bool status; + + float EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "EXT_FIT_NSIGMA"); + EXT_OUTER_RADIUS = psMetadataLookupF32 (&status, config, "EXT_OUTER_RADIUS"); + EXT_FIT_PADDING = psMetadataLookupF32 (&status, config, "EXT_FIT_PADDING"); + EXT_FLUX_LIMIT = EXT_FIT_NSIGMA * sky->sampleStdev; + + // this function specifies the radius at this the model hits the given flux + modelRadiusEXT = pmModelRadius_GetFunction (type); + return true; +} + +bool psphotCheckRadiusEXT (eamReadout *imdata, pmSource *source, pmModel *model) { + + // set the fit radius based on the object flux limit and the model + model->radius = modelRadiusEXT (model->params, EXT_FLUX_LIMIT) + EXT_FIT_PADDING; + if (isnan(model->radius)) psAbort (__func__, "error in radius"); + + // redefine the pixels if needed + bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); + return status; +} Index: /tags/psphot_dev_07/psphot/src/psphotReplaceUnfit.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotReplaceUnfit.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotReplaceUnfit.c (revision 5994) @@ -0,0 +1,21 @@ +# include "psphot.h" + +bool psphotReplaceUnfit (psArray *sources) { + + pmSource *source; + + psTimerStart ("psphot"); + + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_TEMPSUB)) continue; + if (source->modelPSF == NULL) continue; + + pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode &= ~PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + } + psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%d objects)\n", psTimerMark ("psphot"), sources->n); + return true; +} Index: /tags/psphot_dev_07/psphot/src/psphotRoughClass.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotRoughClass.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotRoughClass.c (revision 5994) @@ -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, EXTENDED, 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: /tags/psphot_dev_07/psphot/src/psphotSetup.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotSetup.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotSetup.c (revision 5994) @@ -0,0 +1,112 @@ +# include "psphot.h" + +// load the image +// load or construct the mask +// apply X/Y border info +// apply the SATURATE mask? +// load or construct the weight + +eamReadout *psphotSetup (psMetadata *config) +{ + psMetadata *header = NULL; + psImage *image = NULL; + psImage *weight = NULL; + psImage *mask = NULL; + psRegion region = {0,0,0,0}; // a region representing the entire array + bool status; + + psTimerStart ("psphot"); + + // load implementation-specific models + psphotModelGroupInit (); + + // setup header template, specify COMMENT and HISTORY as MULTI + // XXX EAM - shouldn't this be part of the basic FITS concept? + header = psMetadataAlloc (); + psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_DATA_METADATA_MULTI, "folder for comment", NULL); + psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_DATA_METADATA_MULTI, "folder for history", NULL); + + // read header and image data from INPUT + char *input = psMetadataLookupPtr (&status, config, "IMAGE"); + if (!status) psAbort ("psphot", "input image not specified"); + psFits *file = psFitsOpen (input, "r"); + header = psFitsReadHeader (header, file); + psImage *tmpimage = psFitsReadImage (NULL, file, region, 0); + image = psImageCopy (NULL, tmpimage, PS_TYPE_F32); + psFitsClose (file); + + // psFree (input); + // XXX EAM - does the current API expect to get a copy or a view? + + // grab these values from the approrpiate location (image header if necessary) + float RDNOISE = pmConfigLookupF32 (&status, config, header, "RDNOISE"); + float GAIN = pmConfigLookupF32 (&status, config, header, "GAIN"); + psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_DATA_F32 | PS_META_REPLACE, "read noise value used", RDNOISE); + psMetadataAdd (config, PS_LIST_TAIL, "GAIN", PS_DATA_F32 | PS_META_REPLACE, "gain value used", GAIN); + + // load the weight image if it is supplied, otherwise build from input + // by 'weight' we really mean variance image + char *weightName = psMetadataLookupPtr (&status, config, "WEIGHT"); + if (status == true) { + file = psFitsOpen (weightName, "r"); + weight = psFitsReadImage (NULL, file, region, 0); + psFitsClose (file); + // psFree (weightName); XXX - see psFree (input) + } else { + // build the weight image from the input image, RDNOISE, and GAIN + weight = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32); + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + weight->data.F32[iy][ix] = PS_MAX (image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN), 0.0); + } + } + } + + // load the mask if specified, otherwise construct + char *maskName = psMetadataLookupPtr (&status, config, "MASK"); + if (status == true) { + // XXX EAM require / convert mask to psU8? + file = psFitsOpen (maskName, "r"); + mask = psFitsReadImage (NULL, file, region, 0); + psFitsClose (file); + // psFree (maskName); XXX - see psFree (input) + } else { + mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); + } + + // mask the excluded outer pixels + float XMIN = psMetadataLookupF32 (&status, config, "XMIN"); + float XMAX = psMetadataLookupF32 (&status, config, "XMAX"); + float YMIN = psMetadataLookupF32 (&status, config, "YMIN"); + float YMAX = psMetadataLookupF32 (&status, config, "YMAX"); + psRegion keep = psRegionSet (XMIN, XMAX, YMIN, YMAX); + keep = psRegionForImage (image, keep); + psImageKeepRegion (mask, keep, "OR", PSPHOT_MASK_INVALID); + + // mask the saturated pixels + // XXX EAM does the mask need to grow? + float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + for (int i = 0; status && (i < image->numRows); i++) { + for (int j = 0; j < image->numCols; j++) { + if (image->data.F32[i][j] >= SATURATION) { + mask->data.U8[i][j] |= PSPHOT_MASK_SATURATED; + } + } + } + + // mask the pixels below min threshold + float MIN_VALID = psMetadataLookupF32 (&status, config, "MIN_VALID_PIXEL"); + for (int i = 0; status && (i < image->numRows); i++) { + for (int j = 0; j < image->numCols; j++) { + if (image->data.F32[i][j] < MIN_VALID) { + mask->data.U8[i][j] |= PSPHOT_MASK_INVALID; + } + } + } + + psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); + + // return image data + eamReadout *imdata = eamReadoutAlloc(image, weight, mask, header); + return (imdata); +} Index: /tags/psphot_dev_07/psphot/src/psphotSortBySN.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotSortBySN.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotSortBySN.c (revision 5994) @@ -0,0 +1,34 @@ +# include "psphot.h" + +// sort by SN (descending) +int psphotSortBySN (const void **a, const void **b) +{ + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; + + psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN; + psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN; + if (isnan (fA)) fA = 0; + if (isnan (fB)) fB = 0; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} + +// sort by Y (ascending) +int psphotSortByY (const void **a, const void **b) +{ + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; + + psF32 fA = (A->moments == NULL) ? 0 : A->moments->y; + psF32 fB = (B->moments == NULL) ? 0 : B->moments->y; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); +} + Index: /tags/psphot_dev_07/psphot/src/psphotSourceFits.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotSourceFits.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotSourceFits.c (revision 5994) @@ -0,0 +1,287 @@ +# include "psphot.h" +# include "psEllipse.h" + +// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful + +bool psphotFitPSF (eamReadout *imdata, pmSource *source) { + + float x, y; + + // save the PSF model from the Ensemble fit + pmModel *PSF = pmModelCopy (source->modelPSF); + + // extend source radius as needed + psphotCheckRadiusPSF (imdata, source, PSF); + + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, PSF, true); + psImageKeepCircle (source->mask, x, y, PSF->radius, "AND", ~PSPHOT_MASK_MARKED); + + // does the PSF model succeed? + if (!psphotEvalPSF (source, PSF)) { + psFree (PSF); + return false; + } + + psTrace ("psphot.blend", 5, "fitted as PSF\n"); + pmSourceSubModel (source->pixels, source->mask, PSF, false, false); + + // free old model, save new model + psFree (source->modelPSF); + source->modelPSF = PSF; + + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + return true; +} + + +static float EXT_MIN_SN; +static float EXT_MOMENTS_RAD; +static pmModelType modelTypeEXT; + +bool psphotInitLimitsEXT (psMetadata *config, psStats *sky) { + + bool status; + + // extended source model parameters + EXT_MIN_SN = psMetadataLookupF32 (&status, config, "EXT_MIN_SN"); + EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS"); + + // extended source model descriptions + char *modelNameEXT = psMetadataLookupPtr (&status, config, "EXT_MODEL"); + modelTypeEXT = pmModelSetType (modelNameEXT); + psphotInitRadiusEXT (config, sky, modelTypeEXT); + + return true; +} + +bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources) { + + bool okEXT, okDBL; + float chiEXT, chiDBL; + + // skip the source if we don't think it is extended + if (source->type == PM_SOURCE_UNKNOWN) return false; + if (source->type == PM_SOURCE_DEFECT) return false; + if (source->type == PM_SOURCE_SATURATED) return false; + if (source->moments->SN < EXT_MIN_SN) return false; + + // recalculate the source moments using the larger extended-source moments radius + if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false; + + psTrace ("psphot.blend", 5, "trying blob...\n"); + + pmSource *tmpSrc = pmSourceAlloc (); + + pmModel *EXT = psphotFitEXT (imdata, source); + okEXT = psphotEvalEXT (tmpSrc, EXT); + chiEXT = EXT->chisq / EXT->nDOF; + + psArray *DBL = psphotFitDBL (imdata, source); + okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); + okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); + pmModel *ONE = DBL->data[0]; + chiDBL = ONE->chisq / ONE->nDOF; + + if (okEXT && okDBL) { + psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL); + // XXX EAM : a bogus bias: need to examine this better + if (3*chiEXT > chiDBL) goto keepDBL; + goto keepEXT; + } + + if (okEXT && !okDBL) goto keepEXT; + if (!okEXT && okDBL) goto keepDBL; + return false; + +keepEXT: + // sub EXT + pmSourceSubModel (source->pixels, source->mask, EXT, false, false); + psTrace ("psphot.blend", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[2], EXT->params->data.F32[3]); + + // save new model + source->modelEXT = EXT; + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + return true; + +keepDBL: + // sub DLB + pmSourceSubModel (source->pixels, source->mask, (pmModel *) DBL->data[0], false, false); + pmSourceSubModel (source->pixels, source->mask, (pmModel *) DBL->data[1], false, false); + psTrace ("psphot.blend", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]); + + source->modelPSF = (pmModel *) DBL->data[0]; + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + + // save new model + // tmpSrc->modelPSF = (pmModel *) modelSet->data[1]; + // psArrayAdd (sources, 100, tmpSrc); + + return true; +} + +// fit a double PSF source to an extended blob +psArray *psphotFitDBL (eamReadout *imdata, pmSource *source) { + + float x, y, dx, dy; + pmModel *DBL; + pmModel *PSF; + EllipseAxes axes; + EllipseMoments moments; + psArray *modelSet; + + // make a guess at the position of the two sources + moments.x2 = source->moments->Sx; + moments.y2 = source->moments->Sy; + moments.xy = source->moments->Sxy; + axes = EllipseMomentsToAxes (moments); + + dx = 2 * cos (axes.theta); + dy = 2 * sin (axes.theta); + + // save the PSF model from the Ensemble fit + PSF = source->modelPSF; + psphotCheckRadiusPSF (imdata, source, PSF); + + modelSet = psArrayAlloc (2); + + DBL = pmModelCopy (PSF); + DBL->params->data.F32[1] *= 0.5; + DBL->params->data.F32[2] += dx; + DBL->params->data.F32[3] += dy; + modelSet->data[0] = DBL; + + DBL = pmModelCopy (PSF); + DBL->params->data.F32[1] *= 0.5; + DBL->params->data.F32[2] -= dx; + DBL->params->data.F32[3] -= dy; + modelSet->data[1] = DBL; + + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; + + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitSet (source, modelSet, true); + psImageKeepCircle (source->mask, x, y, PSF->radius, "AND", ~PSPHOT_MASK_MARKED); + + return (modelSet); +} + +pmModel *psphotFitEXT (eamReadout *imdata, pmSource *source) { + + float x, y; + + // use the source moments, etc to guess basic model parameters + pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); + + psphotCheckRadiusEXT (imdata, source, EXT); + + x = EXT->params->data.F32[2]; + y = EXT->params->data.F32[3]; + + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, EXT->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, EXT, false); + psImageKeepCircle (source->mask, x, y, EXT->radius, "AND", ~PSPHOT_MASK_MARKED); + + return (EXT); +} + +bool psphotFitBlend (eamReadout *imdata, pmSource *source) { + + float x, y, dR; + + // if this source is not a possible blend, just fit as PSF + if ((source->blends == NULL) || (source->mode & PM_SOURCE_SATSTAR)) { + bool status = psphotFitPSF (imdata, source); + return status; + } + + psTrace ("psphot.blend", 5, "trying blend...\n"); + + // save the PSF model from the Ensemble fit + pmModel *PSF = pmModelCopy (source->modelPSF); + + // extend source radius as needed + psphotCheckRadiusPSF (imdata, source, PSF); + + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; + + psArray *modelSet = psArrayAlloc (source->blends->n + 1); + modelSet->n = 0; + psArrayAdd (modelSet, 16, PSF); + + psArray *sourceSet = psArrayAlloc (source->blends->n + 1); + sourceSet->n = 0; + psArrayAdd (sourceSet, 16, source); + + for (int i = 0; i < source->blends->n; i++) { + pmSource *blend = source->blends->data[i]; + + // is this object close enough to include in the fit + // XXX this test is bogus: need to think about this + dR = hypot (blend->peak->x - x, blend->peak->y - y); + if (dR > PSF->radius - 3) continue; + + // create the model and guess parameters for this blend + blend->modelPSF = pmModelAlloc (PSF->type); + pmModel *model = blend->modelPSF; + + for (psS32 i = 0; i < model->params->n; i++) { + model->params->data.F32[i] = PSF->params->data.F32[i]; + model->dparams->data.F32[i] = PSF->dparams->data.F32[i]; + } + + model->params->data.F32[1] = blend->peak->counts; + model->params->data.F32[2] = blend->peak->x; + model->params->data.F32[3] = blend->peak->y; + + // add this blend to the list + psArrayAdd (modelSet, 16, model); + psArrayAdd (sourceSet, 16, blend); + } + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitSet (source, modelSet, true); + psImageKeepCircle (source->mask, x, y, PSF->radius, "AND", ~PSPHOT_MASK_MARKED); + + // evaluate the blend objects, subtract if good, free otherwise + for (int i = 1; i < modelSet->n; i++) { + pmSource *src = sourceSet->data[i]; + pmModel *model = modelSet->data[i]; + + if (!psphotEvalPSF (src, model)) { + psFree (model); + src->modelPSF = NULL; + continue; + } + + psTrace ("psphot.blend", 5, "fitted blend as PSF\n"); + pmSourceSubModel (source->pixels, source->mask, model, false, false); + src->mode |= PM_SOURCE_SUBTRACTED; + src->mode &= ~PM_SOURCE_TEMPSUB; + } + + // evaluate the primary object + if (!psphotEvalPSF (source, PSF)) { + psFree (PSF); + return false; + } + + psTrace ("psphot.blend", 5, "fitted primary as PSF\n"); + pmSourceSubModel (source->pixels, source->mask, PSF, false, false); + source->modelPSF = PSF; + source->mode |= PM_SOURCE_SUBTRACTED; + source->mode &= ~PM_SOURCE_TEMPSUB; + return true; +} Index: /tags/psphot_dev_07/psphot/src/psphotSourceStats.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotSourceStats.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotSourceStats.c (revision 5994) @@ -0,0 +1,70 @@ +# include "psphot.h" + +psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *peaks) +{ + bool status = false; + psArray *sources = NULL; + float BIG_RADIUS; + + psTimerStart ("psphot"); + + // determine properties (sky, moments) of initial sources + float INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS"); + float OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS"); + + sources = psArrayAlloc (peaks->n); + sources->n = 0; + + for (int i = 0; i < peaks->n; i++) { + + // create a new source, add peak + pmSource *source = pmSourceAlloc(); + source->peak = (pmPeak *)psMemCopy(peaks->data[i]); + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + // XXX EAM : this has 12 leaks (v.5) + psphotDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER); + + // measure a local sky value + // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER); + if (!status) { + psFree (source); + continue; + } + + // measure basic source moments + status = pmSourceMoments (source, RADIUS); + if (status) { + // add to the source array + psArrayAdd (sources, 100, source); + psFree (source); + continue; + } + + // if no valid pixels, or massive swing, likely saturated source, + // try a much larger box + BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS); + if (status) { + // add to the source array + psArrayAdd (sources, 100, source); + psFree (source); + continue; + } + + psFree (source); + continue; + } + + 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); +} + +// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky) Index: /tags/psphot_dev_07/psphot/src/psphotTest.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotTest.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotTest.c (revision 5994) @@ -0,0 +1,66 @@ +# include "psphot.h" + +void psExit (int status, char *process, char *format, ...) { + + va_list ap; + + va_start (ap, format); + fprintf (stderr, "exiting %s\n", process); + vfprintf (stderr, format, ap); + va_end (ap); + + exit (status); +} + +int main (int argc, char **argv) { + + psMetadata *row; + psArray *table; + + psMetadataItem *mdi; + psRegion region = {0,0,0,0}; // a region representing the entire array + + psphotTestArguments (&argc, argv); + + psFits *file = psFitsOpen (argv[1], "r"); + psMetadata *header = psFitsReadHeader (NULL, file); + psImage *image = psFitsReadImage (NULL, file, region, 0); + psFitsClose (file); + + psMetadataConfigWrite (header, argv[2]); + + // attempt to write image with NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create a test image + // psImage *tmpimage = psImageAlloc (10, 10, PS_DATA_F32); + + // create a test table + table = psArrayAlloc (10); + table->n = 0; + + for (int i = 0; i < 10; i++) { + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "ROW", PS_DATA_S32, "", i); + psMetadataAdd (row, PS_LIST_TAIL, "FROW", PS_TYPE_F32, "", 0.1*i); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "test line"); + + table->data[i] = row; + } + table->n = 10; + + psMetadata *theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); + + psFits *fits = psFitsOpen (argv[3], "w"); + // psFitsWriteImage (fits, header, tmpimage, 0); + psFitsWriteHeader (header, fits); + psFitsWriteTable (fits, theader, table); + psFitsClose (fits); + + psFree (header); + psFree (image); + exit (0); +} Index: /tags/psphot_dev_07/psphot/src/psphotTestArguments.c =================================================================== --- /tags/psphot_dev_07/psphot/src/psphotTestArguments.c (revision 5994) +++ /tags/psphot_dev_07/psphot/src/psphotTestArguments.c (revision 5994) @@ -0,0 +1,19 @@ +# include "psphot.h" +static int usage (); + +void psphotTestArguments (int *argc, char **argv) { + + // basic pslib options + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + if (*argc != 4) usage (); + + return; +} + +static int usage () { + + fprintf (stderr, "USAGE: psphotTest (input.fits) (output.hdr) (out.fits)\n"); + exit (2); +}