Index: /tags/eam_tags/psphot_dev_06/psphot/.cvsignore =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/.cvsignore (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/.cvsignore (revision 22339) @@ -0,0 +1,1 @@ +bin Index: /tags/eam_tags/psphot_dev_06/psphot/Makefile =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/Makefile (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/Makefile (revision 22339) @@ -0,0 +1,136 @@ +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 -g +LFLAGS = $(LIBS) + +PSPHOT = \ +$(SRC)/psphot.$(ARCH).o \ +$(SRC)/psphotArguments.$(ARCH).o \ +$(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psphotImageStats.$(ARCH).o \ +$(SRC)/psphotSourceStats.$(ARCH).o \ +$(SRC)/psphotChoosePSF.$(ARCH).o \ +$(SRC)/psphotApplyPSF.$(ARCH).o \ +$(SRC)/psphotFixedPSF.$(ARCH).o \ +$(SRC)/psphotEnsemblePSF.$(ARCH).o \ +$(SRC)/psphotReapplyPSF.$(ARCH).o \ +$(SRC)/psphotFitGalaxies.$(ARCH).o \ +$(SRC)/psphotOutput.$(ARCH).o \ +$(SRC)/psphotMarkPSF.$(ARCH).o \ +$(SRC)/psphotSubtractPSF.$(ARCH).o \ +$(SRC)/psphotSortBySN.$(ARCH).o \ +$(SRC)/psphotDefinePixels.$(ARCH).o\ +$(SRC)/psphotMagnitudes.$(ARCH).o \ +$(SRC)/psphotImageBackground.$(ARCH).o \ +$(SRC)/psphotBasicDeblend.$(ARCH).o \ +$(SRC)/psphotModelGroupInit.$(ARCH).o \ +$(SRC)/pmSourceContour.$(ARCH).o \ +$(SRC)/psLine.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ +$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \ +$(SRC)/pmSourceFitFixed.$(ARCH).o \ +$(SRC)/psSparse.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o + +PSMODULES = \ +$(SRC)/psEllipse.$(ARCH).o \ +$(SRC)/pmPSF.$(ARCH).o \ +$(SRC)/pmPSFtry.$(ARCH).o \ +$(SRC)/pmModelGroup.$(ARCH).o \ +$(SRC)/pmObjects_EAM.$(ARCH).o + +MODELS = \ +$(SRC)/models/pmModel_GAUSS.c \ +$(SRC)/models/pmModel_PGAUSS.c \ +$(SRC)/models/pmModel_QGAUSS.c \ +$(SRC)/models/pmModel_SGAUSS.c \ +$(SRC)/models/pmModel_TGAUSS.c + +$(SRC)/pmModelGroup.$(ARCH).o: $(MODELS) + +MODELTEST = \ +$(SRC)/modelTestFitSource.$(ARCH).o \ +$(SRC)/modelTestArguments.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o \ +$(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psphotDefinePixels.$(ARCH).o \ +$(SRC)/psphotModelTest.$(ARCH).o + +TEST = \ +$(SRC)/psphotTest.$(ARCH).o \ +$(SRC)/psphotTestArguments.$(ARCH).o \ +$(SRC)/pmSourceContour.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o \ +$(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ +$(SRC)/psSparse.$(ARCH).o + +psphot: $(BIN)/psphot.$(ARCH) +$(BIN)/psphot.$(ARCH) : $(PSPHOT) +$(PSPHOT) $(PSMODULES): $(SRC)/psphot.h + +modeltest.install: psphotModelTest.install +modeltest: $(BIN)/psphotModelTest.$(ARCH) +$(BIN)/psphotModelTest.$(ARCH) : $(MODELTEST) +$(MODELTEST): $(SRC)/psphot.h + +test.install: psphotTest.install +test: $(BIN)/psphotTest.$(ARCH) +$(BIN)/psphotTest.$(ARCH) : $(TEST) +$(TEST): $(SRC)/psphot.h + +INSTALL = psphot psphotTest psphotModelTest + +# 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); fi + $(CC) $^ -o $@ $(LFLAGS) + @echo "done with $@" + +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$* + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* + +$(INSTALL): % : $(BIN)/%.$(ARCH) + +%.clean : + rm -f $(BIN)/$*.$(ARCH) + rm -f $(SRC)/*.$(ARCH).o + +%.install: + make $(DESTBIN)/$* + +# utilities ################################################# + +install: + for i in $(INSTALL); do make $$i.install; done + +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` Index: /tags/eam_tags/psphot_dev_06/psphot/doc/config.txt =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/doc/config.txt (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/doc/config.txt (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/doc/notes.txt =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/doc/notes.txt (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/doc/notes.txt (revision 22339) @@ -0,0 +1,293 @@ + +Notes on psphot + +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/eam_tags/psphot_dev_06/psphot/doc/psphot.txt =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/doc/psphot.txt (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/doc/psphot.txt (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/doc/versions.txt =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/doc/versions.txt (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/doc/versions.txt (revision 22339) @@ -0,0 +1,25 @@ + +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/eam_tags/psphot_dev_06/psphot/src/modelTestArguments.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/modelTestArguments.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/modelTestArguments.c (revision 22339) @@ -0,0 +1,73 @@ +# include "psphot.h" +static int usage (); + +psMetadata *modelTestArguments (int *argc, char **argv) { + + int N; + unsigned int Nfail; + int mode = PS_DATA_STRING | PS_META_REPLACE; + psF32 UseCoords_X, UseCoords_Y; + + // basic pslib options + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + // identify options in args list - these go on config + bool UseCoords = false; + if ((N = psArgumentGet (*argc, argv, "-coords"))) { + UseCoords = true; + psArgumentRemove (N, argc, argv); + UseCoords_X = atof (argv[N]); + psArgumentRemove (N, argc, argv); + UseCoords_Y = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + char *model = NULL; + if ((N = psArgumentGet (*argc, argv, "-model"))) { + psArgumentRemove (N, argc, argv); + model = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + char *fitmode = NULL; + if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { + psArgumentRemove (N, argc, argv); + fitmode = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + if (*argc != 3) usage (); + + // load config information + psMetadata *config = psMetadataAlloc (); + psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); + config = psMetadataConfigParse (config, &Nfail, argv[2], FALSE); + + // identify input image + psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); + + if (model != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODEL", mode, "", model); + } + if (fitmode != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODE", mode, "", fitmode); + } + if (UseCoords) { + mode = PS_DATA_F32 | PS_META_REPLACE; + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", mode, "", UseCoords_X); + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", mode, "", UseCoords_Y); + } + + return (config); +} + +static int usage () { + + fprintf (stderr, "USAGE: psphotModelTest (image.fits) (config)\n"); + fprintf (stderr, "options: \n"); + fprintf (stderr, " -coords x y : specify object center\n"); + fprintf (stderr, " -model (name) : specify model of interest\n"); + fprintf (stderr, " -fitmode (flt|psf) : specify how to fit model\n"); + exit (2); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/modelTestFitSource.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/modelTestFitSource.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/modelTestFitSource.c (revision 22339) @@ -0,0 +1,112 @@ +# include "psphot.h" + +int DumpImage (psImage *image, char *filename); + +bool modelTestFitSource (eamReadout *imdata, psMetadata *config) { + + bool status; + int modelType; + float obsMag, fitMag, value; + char name[64]; + + psMetadataItem *item = NULL; + + // 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); + DumpImage (source->pixels, "object.fits"); + + // 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; + + // 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++) { + sprintf (name, "TEST_FIT_PAR%d", i); + value = psMetadataLookupF32 (&status, config, name); + if (status) { + params[i] = value; + } + } + + // what fitting mode to use? + char *fitModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE"); + if (!status) { + fitModeWord = psStringCopy ("FLT"); + } + bool fitMode = !strcasecmp (fitModeWord, "PSF"); + + // define the pixels used for the fit + psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, model, fitMode); + + // 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); + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + // write out + DumpImage (source->pixels, "resid.fits"); + DumpImage (source->mask, "mask.fits"); + return true; +} + +int DumpImage (psImage *image, char *filename) { + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0); + psFitsClose (fits); + return (TRUE); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_GAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_GAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_GAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_PGAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_PGAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_PGAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_QGAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_QGAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_QGAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_RGAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_RGAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_RGAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_SGAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_SGAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_SGAUSS.c (revision 22339) @@ -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.001; + + 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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_TAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_TAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_TAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_TGAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_TGAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_TGAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_WAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_WAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_WAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/models/pmModel_ZGAUSS.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_ZGAUSS.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/models/pmModel_ZGAUSS.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/pmPeaksSigmaLimit.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/pmPeaksSigmaLimit.c (revision 22339) @@ -0,0 +1,48 @@ +# include "psphot.h" + +psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky) +{ + bool status = false; + float NSIGMA; + float SIGMA; + float threshold; + + // smooth the image + + psTimerStart ("psphot"); + + SIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_SIGMA"); + NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA"); + + psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); + psImageSmooth (smooth, SIGMA, NSIGMA); + psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); + + psTimerStart ("psphot"); + + // set peak threshold + NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT"); + + // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + threshold = NSIGMA*sky->sampleStdev; + psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold); + + // find the peaks in the smoothed image + psArray *peaks = pmFindImagePeaks (smooth, threshold); + if (peaks == NULL) psAbort ("find peaks", "no peaks found"); + psFree (smooth); + + // optional dump of all peak data + char *output = psMetadataLookupPtr (&status, config, "PEAKS_OUTPUT_FILE"); + if (status && (output != NULL) && (output[0])) { + pmPeaksWriteText (peaks, output); + psFree (output); + } + + psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); + return (peaks); +} + +// In this function, we smooth the image, then search for the peaks +// Should we also subtract a super-binned image? (as an option?) +// We need to gracefully handle no source detections Index: /tags/eam_tags/psphot_dev_06/psphot/src/pmSourceContour.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/pmSourceContour.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/pmSourceContour.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/pmSourceFitFixed.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/pmSourceFitFixed.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/pmSourceFitFixed.c (revision 22339) @@ -0,0 +1,90 @@ +# include "psphot.h" +# define USECLIP 0 + +bool pmSourceFitFixed (pmSource *source, pmModel *model) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->moments, 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); + + psF32 Ro, ymodel; + + psVector *params = model->params; + psVector *dparams = model->dparams; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + psF32 Sky = params->data.F32[0]; + + // construct the coordinate and value entries + psVector *x = psVectorAlloc(100, PS_TYPE_F32); + psVector *y = psVectorAlloc(100, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(100, PS_TYPE_F32); + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + x->n = y->n = yErr->n = 0; + psS32 count = 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; + } + + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + ymodel = modelFunc (NULL, params, coord); + + x->data.F32[count] = source->pixels->data.F32[i][j]; + y->data.F32[count] = ymodel; + + // this test modifies the weight based on deviation from the model flux + Ro = 1.0 + fabs (x->data.F32[count] - ymodel) / sqrt(PS_SQR(ymodel - Sky) + PS_SQR(Sky)); + yErr->data.F32[count] = sqrt (source->weight->data.F32[i][j] * Ro); + + psVectorExtend (x, 100, 1); + psVectorExtend (y, 100, 1); + psVectorExtend (yErr, 100, 1); + count++; + } + } + + psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); + +# if (USECLIP) + psVector *mask = psVectorAlloc (count, PS_TYPE_U8); + psVectorInit (mask, 0); + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + stats->clipIter = 3; + poly = psVectorClipFitPolynomial1D (poly, stats, mask, 0, y, yErr, x); +# else + poly = psVectorFitPolynomial1D (poly, NULL, 0, y, yErr, x); +# endif + + // fitted values are variations from current values + params->data.F32[0] += poly->coeff[0]; + params->data.F32[1] *= poly->coeff[1]; + + // XXX need to get errors back from psVectorClipFit... + dparams->data.F32[0] = 0; + dparams->data.F32[1] = 0; + + // XXX need to get chisq back from psVectorClipFit... + model->chisq = 0; + model->nIter = 0; + model->nDOF = 0; + + psFree(x); + psFree(y); + psFree(yErr); + +# if (USECLIP) + psFree(stats); + psFree(mask); +# endif + + return true; +} + Index: /tags/eam_tags/psphot_dev_06/psphot/src/psImageData.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psImageData.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psImageData.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psLibUtils.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psLibUtils.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psLibUtils.c (revision 22339) @@ -0,0 +1,65 @@ +# include // for strncasecmp +# include // for write +# include +# include "psLibUtils.h" + +// XXX EAM : These two functions (psLogArguments and psTraceArguments) have been +// rolled together into a single function in psLib +// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) +// 2 = default, -v = 3, -vv = 4, -vvv = 5 +psS32 psLogArguments (int *argc, char **argv) { + + int N, level; + + // default log level is 2 + level = 2; + + // set in order, so that -vvv overrides -vv overrides -v + if ((N = psArgumentGet (*argc, argv, "-v"))) { + psArgumentRemove (N, argc, argv); + level = 3; + } + if ((N = psArgumentGet (*argc, argv, "-vv"))) { + psArgumentRemove (N, argc, argv); + level = 4; + } + if ((N = psArgumentGet (*argc, argv, "-vvv"))) { + psArgumentRemove (N, argc, argv); + level = 5; + } + + if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { + if (*argc < N + 2) { + psAbort ("psLogArguments", "USAGE: -logfmt (format)"); + } + psArgumentRemove (N, argc, argv); + psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid + psArgumentRemove (N, argc, argv); + } + + // set the level, return the level + psLogSetLevel (level); + return (level); +} + +// set trace levels by facility +psS32 psTraceArguments (int *argc, char **argv) { + + int N; + + // argument format is: -trace (facil) (level) + while ((N = psArgumentGet (*argc, argv, "-trace"))) { + if (*argc < N + 3) { + psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); + } + psArgumentRemove (N, argc, argv); + psTraceSetLevel (argv[N], atoi(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { + psTracePrintLevels (); + exit (2); + } + return (TRUE); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psLibUtils.h =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psLibUtils.h (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psLibUtils.h (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psLine.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psLine.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psLine.c (revision 22339) @@ -0,0 +1,47 @@ +# 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); + + if (Nchar >= nMax) return (false); + return (true); +} + Index: /tags/eam_tags/psphot_dev_06/psphot/src/psModulesUtils.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psModulesUtils.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psModulesUtils.c (revision 22339) @@ -0,0 +1,102 @@ +# 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]; + + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + + 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; + if (fitSum <= 0) return false; + + *fitMag = -2.5*log10(fitSum); + *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; +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psModulesUtils.h =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psModulesUtils.h (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psModulesUtils.h (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psSparse.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psSparse.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psSparse.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psSparse.h =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psSparse.h (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psSparse.h (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psphot.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphot.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphot.c (revision 22339) @@ -0,0 +1,97 @@ +# include "psphot.h" + +int main (int argc, char **argv) { + + psMetadata *config = NULL; + eamReadout *imdata = NULL; + psArray *sources = NULL; + psArray *peaks = NULL; + pmPSF *psf = NULL; + psStats *sky = NULL; + pmPSFClump psfClump; + bool status; + + psphotModelGroupInit (); + + config = psphotArguments (&argc, argv); + + // load input data (image and config) + // create or load mask and weight images + // we have memory leaks here -- may be from psMetadata + imdata = psphotSetup (config); + char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT"); + if (!status) { + breakPt = psStringCopy ("NONE"); + } + + // measure image stats for initial guess + sky = psphotImageStats (imdata, config); + + // psPolynomial2D *skyModel = psphotImageBackground (imdata, config, sky); + psphotImageBackground (imdata, config, sky); + + // find the peaks in the image + peaks = pmPeaksSigmaLimit (imdata, config, sky); + if (!strcasecmp (breakPt, "PEAKS")) exit (0); + + // construct sources and measure basic stats + sources = psphotSourceStats (imdata, config, peaks); + if (!strcasecmp (breakPt, "STATS")) exit (0); + + // classify sources based on moments, brightness + // XXX - this has 1 leak + psfClump = pmSourcePSFClump (sources, config); + + // group into STAR, COSMIC, GALAXY, SATURATED + pmSourceRoughClass (sources, config, psfClump); + if (!strcasecmp (breakPt, "CLASS")) exit (0); + + psphotBasicDeblend (sources, config, sky); + if (!strcasecmp (breakPt, "DEBLEND")) exit (0); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // use bright stellar objects to measure PSF + psf = psphotChoosePSF (config, sources, sky); + if (!strcasecmp (breakPt, "PSFFIT")) exit (0); + + int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE"); + if (!status) FITMODE = 2; + + switch (FITMODE) { + case -2: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + psphotReapplyPSF (imdata, config, sources, psf, sky); + break; + + case -1: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + break; + + case 0: + psphotFixedPSF (imdata, config, sources, psf, sky); + break; + + case 1: + // test PSF on all except SATURATE and DEFECT (artifacts) + psphotApplyPSF (imdata, config, sources, psf, sky); + break; + + case 2: + // fit extended objects with galaxy models + psphotApplyPSF (imdata, config, sources, psf, sky); + psphotFitGalaxies (imdata, config, sources, sky); + break; + + case 3: + // fit extended objects with galaxy models + psphotFixedPSF (imdata, config, sources, psf, sky); + psphotFitGalaxies (imdata, config, sources, sky); + break; + } + + // write out data in appropriate format + psphotOutput (imdata, config, sources, psf, sky); + exit (0); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphot.h =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphot.h (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphot.h (revision 22339) @@ -0,0 +1,75 @@ +# include +# include // for strcasecmp +# include // for unlink +# 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); +psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); +psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); +bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); +void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); + +bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); +bool psphotSubtractPSF (pmSource *source); +int psphotSortBySN (const void **a, const void **b); +int psphotSortByY (const void **a, const void **b); +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); +bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); +void psphotModelGroupInit (void); + +bool pmSourceFitFixed (pmSource *source, pmModel *model); +psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky); +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); +pmModel *pmSourceSelectModel (pmSource *source); + +// eamReadout functions +eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); + +// 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 pmPeaksWriteText (psArray *sources, char *filename); +bool pmMomentsWriteText (psArray *sources, char *filename); +bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); +bool pmModelWriteFLTs (psArray *sources, char *filename); +bool pmModelWriteNULLs (psArray *sources, char *filename); + +int pmSourcesDophotType (pmSource *source); + +// psphotModelTest functions +psMetadata *modelTestArguments (int *argc, char **argv); +bool modelTestFitSource (eamReadout *imdata, psMetadata *config); +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +float psphotCrossProduct (pmSource *Mi, pmSource *Mj); +psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky); +bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); + +pmModel *pmModelCopy (pmModel *model); +psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); +psMetadata *psphotTestArguments (int *argc, char **argv); +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky); Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotApplyPSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotApplyPSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotApplyPSF.c (revision 22339) @@ -0,0 +1,98 @@ +# include "psphot.h" + +// fit psf model to all objects +// PSFSTAR objects will be refitted +// run this function to a specific flux limit? + +bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool status; + float x; + float y; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + // we may set this differently here from the value used to mark likely saturated stars + float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + + float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); + float PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); + float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); + float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); + float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + + // set the object surface-brightness limit for fitted pixels + float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; + psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + + // this function specifies the radius at this the model hits the given flux + pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); + + 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->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // use the source moments, etc to guess basic model parameters + pmModel *modelFLT = pmSourceModelGuess (source, psf->type); + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelFLT, psf); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + psFree (modelFLT); + + // 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 = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) { + psAbort ("apply_psf_model", "error in radius"); + } + + // check if we need to redefine the pixels + // XXX EAM : a better test would examine the source pixels + if (model->radius > OUTER_RADIUS) { + // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) + psphotDefinePixels (source, imdata, x, y, model->radius); + } + + // if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5); + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, model, true); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); + + if (!status || (model->params->data.F32[1] < 0)) { + psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + psFree (model); + continue; + } + // XXX EAM : this was an attempt to correct for fit-sky biases + // pmModelSkyOffset (model, x, y, model->radius); + + // model succeeded : tentatively keep it + source->modelPSF = model; + Niter += model[0].nIter; + Nfit ++; + + // is it a good model? + psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION); + if (psphotSubtractPSF (source)) { + 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); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotArguments.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotArguments.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotArguments.c (revision 22339) @@ -0,0 +1,86 @@ +# include "psphot.h" + +static void usage (void) ; + +psMetadata *psphotArguments (int *argc, char **argv) { + + int N; + unsigned int Nfail; + int mode = PS_DATA_STRING | PS_META_REPLACE; + + // basic pslib options + fprintf (stderr, "starting... %s\n", psLibVersion()); + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + // optional mask image - add to config + char *mask = NULL; + if ((N = psArgumentGet (*argc, argv, "-mask"))) { + psArgumentRemove (N, argc, argv); + mask = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + // optional weight image - add to config + char *weight = NULL; + if ((N = psArgumentGet (*argc, argv, "-weight"))) { + psArgumentRemove (N, argc, argv); + weight = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + // optional output residual image - add to config + char *resid = NULL; + if ((N = psArgumentGet (*argc, argv, "-resid"))) { + psArgumentRemove (N, argc, argv); + resid = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + // optional output residual image - add to config + char *photcode = NULL; + if ((N = psArgumentGet (*argc, argv, "-photcode"))) { + psArgumentRemove (N, argc, argv); + photcode = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + if (*argc != 4) usage (); + + // load config information + psMetadata *config = psMetadataAlloc (); + psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); + config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE); + + // identify input image & optional weight & mask images + // command-line entries override config-file entries + psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); + psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]); + + if (mask != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask); + } + if (weight != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight); + } + if (resid != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid); + } + if (photcode != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode); + } else { + psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE"); + } + 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/eam_tags/psphot_dev_06/psphot/src/psphotBasicDeblend.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotBasicDeblend.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotBasicDeblend.c (revision 22339) @@ -0,0 +1,134 @@ +# include "psphot.h" + +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky) { + + int N; + bool status; + float threshold; + pmSource *source, *testSource; + + FILE *f = fopen ("deblend.dat", "w"); + if (f == NULL) psAbort ("psphot", "can't open deblend.dat output file"); + int Nblend = 0; + + psTimerStart ("psphot"); + + // this should be added to the PM_SOURCE flags: + int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + + float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION"); + float NSIGMA = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA"); + float minThreshold = NSIGMA*sky->sampleStdev; + fprintf (stderr, "min threshold: %f\n", minThreshold); + + // 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); + + // XXX make sure this results in decreasing, not increasing, SN + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; + + if (source->type == 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) { + // XXX EAM : make the 1/4 user-defined. + // XXX keep threshold from dropping too low (N*sky.sigma) + threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; + threshold = PS_MAX (threshold, minThreshold); + + // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now + psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold); + if (contour == NULL) { + fprintf (stderr, "below threshold? invalid peak?\n"); + 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; + + int xp0 = source->moments->x - source->pixels->col0; + int xp1 = source->peak->x - source->pixels->col0; + int xp2 = testSource->moments->x - source->pixels->col0; + int xp3 = testSource->peak->x - source->pixels->col0; + + int yp0 = source->moments->y - source->pixels->row0; + int yp1 = source->peak->y - source->pixels->row0; + int yp2 = testSource->moments->y - testSource->pixels->row0; + int yp3 = testSource->peak->y - testSource->pixels->row0; + + fprintf (f, "%d %d (%f, %f) : %d %d (%f, %f) vs %f\n", + source->peak->x, source->peak->y, + source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1], + testSource->peak->x, testSource->peak->y, + testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold + ); + + testSource->type = PM_SOURCE_BLEND; + Nblend ++; + j = xv->n; + } + } + } + } + fclose (f); + psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f)\n", Nblend, psTimerMark ("psphot")); + return true; +} + + + Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotChoosePSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotChoosePSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotChoosePSF.c (revision 22339) @@ -0,0 +1,78 @@ +# 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; + + // array to store candidate PSF stars + int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS"); + if (!status) NSTARS = sources->n; + stars = psArrayAlloc (sources->n); + stars->n = 0; + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { + pmSource *source = sources->data[i]; + if (source->type != PM_SOURCE_PSFSTAR) continue; + psArrayAdd (stars, 200, source); + } + psTrace (".psphot.pspsf", 3, "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 EAM - 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; + } + } + + // keep only the selected model: + try = models->data[bestN]; + pmPSF *psf = psMemCopy(try->psf); + psFree (models); // keep only the pmPSF resulting from this analysis + + modelName = pmModelGetType (psf->type); + psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + + return (psf); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotDefinePixels.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotDefinePixels.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotDefinePixels.c (revision 22339) @@ -0,0 +1,21 @@ +# 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); + + return(mySource); +} + Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotEnsemblePSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotEnsemblePSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotEnsemblePSF.c (revision 22339) @@ -0,0 +1,191 @@ +# include "psphot.h" + +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { + + bool status; + float x; + float y; + float Sky; + + psTimerStart ("psphot"); + + // source analysis is done in spatial order + sources = psArraySort (sources, psphotSortByY); + + // this should be added to the PM_SOURCE flags: + int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_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; + + 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->type == PM_SOURCE_BLEND) continue; + if (inSource->type == PM_SOURCE_DEFECT) continue; + if (inSource->type == PM_SOURCE_SATURATED) continue; + + pmSource *otSource = pmSourceAlloc (); + + // XXX EAM : add option to use FLT or PSF form + // use the source moments, etc to guess basic model parameters + pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelFLT, psf); + psFree (modelFLT); + + // save the original coords + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // get function which specifies the radius at this 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 + if (model->radius > OUTER_RADIUS) { + psphotDefinePixels (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; + + // XXX EAM : set model to unit peak, zero sky, maybe use peak (x,y) + // model->params->data.F32[2] = inSource->peak->x; + // model->params->data.F32[3] = inSource->peak->y; + 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 = psphotCrossProduct (Mi, Mi); + psSparseMatrixElement (sparse, i, i, f); + + // find the image x model value + f = psphotCrossProduct (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 = psphotCrossProduct (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; + Fi->modelPSF->params->data.F32[1] = norm->data.F32[i]; + + // subtract object + pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); + + Fi->modelPSF->params->data.F32[0] = Sky; + // need to set this! + } + + // XXX EAM : need to free up many things here + + psLogMsg ("psphot.emsemble", 4, "apply models: %f\n", psTimerMark ("psphot")); + return true; +} + + +float psphotCrossProduct (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/eam_tags/psphot_dev_06/psphot/src/psphotFitGalaxies.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotFitGalaxies.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotFitGalaxies.c (revision 22339) @@ -0,0 +1,109 @@ +# include "psphot.h" + +bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) +{ + bool status, goodfit; + float x; + float y; + int Nfit = 0; + int Nfail = 0; + int Niter = 0; + + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + + float GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); + float GAL_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA"); + float GAL_FIT_PADDING = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING"); + float GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); + + float FLUX_LIMIT = GAL_FIT_NSIGMA * skyStats->sampleStdev; + + char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); + pmModelType modelType = pmModelSetType (modelName); + pmModelRadius modelRadius = pmModelRadius_GetFunction (modelType); + + psTimerStart ("psphot"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // sources which should not be fitted + // skip all valid stars + if (source->type == PM_SOURCE_PSFSTAR) continue; + if (source->type == PM_SOURCE_SATSTAR) continue; + if (source->type == PM_SOURCE_GOODSTAR) continue; + + // skip all likely defects + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // skip poorly fitted stars + if (source->type == PM_SOURCE_FAINTSTAR) continue; + if (source->type == PM_SOURCE_POOR_FIT_PSF) continue; + + // XXX EAM when do we pick these up again? + if (source->moments->SN < GAL_MIN_SN) { + source->type = PM_SOURCE_FAINT_GALAXY; // better choice? + continue; + } + + // recalculate the source moments using the larger galaxy moments radius + status = pmSourceMoments (source, GAL_MOMENTS_RAD); + if (!status) { + source->type = PM_SOURCE_DROP_GALAXY; // better choice? + continue; + } + + // use the source moments, etc to guess basic model parameters + pmModel *model = pmSourceModelGuess (source, modelType); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // 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 = modelRadius (model->params, FLUX_LIMIT) + GAL_FIT_PADDING; + if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); + + // check if we need to redefine the pixels + // XXX EAM : a better test would examine the source pixels + if (model->radius > OUTER_RADIUS) { + // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) + psphotDefinePixels (source, imdata, x, y, model->radius); + } + + // fit FLT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, model, false); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); + if (!status) { + // if the fit fails, we need to change the classification + psLogMsg ("psphot", 5, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); + source->type = PM_SOURCE_FAIL_FIT_GAL; // better choice? + source->modelFLT = model; + Nfail ++; + continue; + } + + // check if model fit is acceptable + goodfit = pmModelFitStatus (model); + if (!goodfit) { + // if the fit fails, we need to change the classification + psLogMsg ("psphot", 5, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); + source->type = PM_SOURCE_POOR_FIT_GAL; // better choice? + source->modelFLT = model; + Nfail ++; + continue; + } + + // accept the model + source->type = PM_SOURCE_GALAXY; + source->modelFLT = model; + Niter += model[0].nIter; + Nfit++; + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, model, false, false); + } + psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter); + return (true); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotFixedPSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotFixedPSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotFixedPSF.c (revision 22339) @@ -0,0 +1,83 @@ +# include "psphot.h" + +// fit psf model to all objects with centroid fixed + +bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool status; + float x; + float y; + int Nfit = 0; + + psTimerStart ("psphot"); + + // we may set this differently here from the value used to mark likely saturated stars + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_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.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + + // this function specifies the radius at this the model hits the given flux + pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); + + 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->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // use the source moments, etc to guess basic model parameters + pmModel *modelFLT = pmSourceModelGuess (source, psf->type); + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelFLT, psf); + // model->params->data.F32[2] = source->peak->x; + // model->params->data.F32[3] = source->peak->y; + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + psFree (modelFLT); + + // 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 = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) { + psAbort ("apply_psf_model", "error in radius"); + } + + // check if we need to redefine the pixels + // XXX EAM : a better test would examine the source pixels + if (model->radius > OUTER_RADIUS) { + // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) + psphotDefinePixels (source, imdata, x, y, model->radius); + } + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitFixed (source, model); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); + + if (!status || (model->params->data.F32[1] < 0)) { + psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (peak %f\n", x, y, model->params->data.F32[1]); + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + psFree (model); + continue; + } + + // model succeeded : tentatively keep it + source->modelPSF = model; + Nfit ++; + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + } + + psLogMsg ("psphot", 3, "fit fixed PSF models: %f sec for %d objects\n", psTimerMark ("psphot"), Nfit); + return (true); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotImageBackground.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotImageBackground.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotImageBackground.c (revision 22339) @@ -0,0 +1,72 @@ +# 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; + psLogMsg ("psphot", 5, "back: %f sec (select %d points)\n", psTimerMark ("psphot"), 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); + psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot")); + fprintf (stderr, "model: %f %f %f %f\n", + skyModel->coeff[0][0], skyModel->coeff[1][0], + skyModel->coeff[0][1], skyModel->coeff[1][1]); + + // 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; + } + } + } + + psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot")); + + psphotSaveImage (imdata->header, imdata->image, "backsub.fits"); + + return (skyModel); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotImageStats.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotImageStats.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotImageStats.c (revision 22339) @@ -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, "stats: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev); + + psFree (stats); + return (sky); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotMagnitudes.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotMagnitudes.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotMagnitudes.c (revision 22339) @@ -0,0 +1,73 @@ +# include "psphot.h" + +// XXX EAM : this aperture correction business is invalid (& wrong) for galaxies +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) { + + int status; + float x, y; + float sky, rflux, apMag, fitMag; + + // use the correct model (PSF vs FLT) + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) return NULL; + + x = source->peak->x; + y = source->peak->y; + + // we have already (psphotApplyPSF) defined pixels at least OUTER_RADIUS from source + // we need to mask pixels to measure the aperture magnitude + // set aperture mask circle of PSF_FIT_RADIUS + psImageKeepCircle (source->mask, x, y, apRadius, "OR", PSPHOT_MASK_MARKED); + + // save local sky for later + sky = model->params->data.F32[0]; + + // replace source flux + pmSourceAddModel (source->pixels, source->mask, model, false, false); + + // measure object photometry + status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + rflux = pow (10.0, 0.4*fitMag); + source->apMag = apMag - rflux * psf->skyBias * (M_PI * PS_SQR(apRadius)); + source->fitMag = fitMag + psf->ApResid; + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + // unmask aperture + psImageKeepCircle (source->mask, x, y, apRadius, "AND", ~PSPHOT_MASK_MARKED); + return model; +} + +pmModel *pmSourceSelectModel (pmSource *source) { + + pmModel *model; + + switch (source->type) { + + // use PSF model with stars + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_SATSTAR: + case PM_SOURCE_GOODSTAR: + model = source->modelPSF; + break; + + // use FLT model with galaxies + case PM_SOURCE_GALAXY: + model = source->modelFLT; + break; + + // skip defects and poorly fitted stars or galaxies + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + case PM_SOURCE_FAINTSTAR: + case PM_SOURCE_POOR_FIT_PSF: + case PM_SOURCE_POOR_FIT_GAL: + case PM_SOURCE_FAIL_FIT_GAL: + default: + model = NULL; + } + return (model); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotMarkPSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotMarkPSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotMarkPSF.c (revision 22339) @@ -0,0 +1,126 @@ +# include "psphot.h" + +// identify objects consistent with PSF shape/magnitude distribution +// we expect dparams[4],dparams[5] to have a scatter of: +// sigma_x / (S/N) * sqrt(2) +// 1 / (params[4],params[5])*(S/N) + +// 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 (but are already IDed) +// galaxies should be larger, cosmic rays smaller, but need to test? +// we also reject objects with S/N too low or ChiSquare to high + +// any object which meets the criterion is marked as PM_SOURCE_BRIGHTSTAR + +// floor for DS value +// XXX EAM : add to configuration? +# define MIN_DS 0.01 + +bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE) +{ + int keep; + float dSX, dSY, SX, SY, SN; + float nSx, nSy, Chi; + + pmModel *model = source->modelPSF; + + // do we actually have a valid PSF model? + if (model == NULL) return (false); + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_BADARGS: + case PM_MODEL_UNTRIED: + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + return false; + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", model->params->data.F32[2], model->params->data.F32[3], model->nIter); + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + return false; + default: + return false; + } + + // 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) + if (model->params->data.F32[1] >= SATURATE) { + if (source->type == PM_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->type = 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->type = PM_SOURCE_FAIL_FIT_PSF; + return (false); + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->type == PM_SOURCE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); + source->type = PM_SOURCE_GOODSTAR; + } + + 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) < shapeNsigma); + keep &= (fabs(nSy) < shapeNsigma); + keep &= (SN > minSN); + keep &= (Chi < maxChi); + if (keep) { + if (source->type == PM_SOURCE_PSFSTAR) return (true); + source->type = PM_SOURCE_GOODSTAR; + return (true); + } + + if (source->type == PM_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n"); + } + + // object appears to be small, suspected defect + if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { + source->type = PM_SOURCE_DEFECT; + return (false); + } + + // object appears to be large, suspected galaxy + if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { + source->type = PM_SOURCE_GALAXY; + return (false); + } + + // object appears to be extremely faint: what is this? + if (SN < minSN) { + source->type = PM_SOURCE_FAINTSTAR; + return (false); + } + + // these are pooly fitted, probable stars near other stars? + source->type = PM_SOURCE_POOR_FIT_PSF; + return (false); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotModelGroupInit.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotModelGroupInit.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotModelGroupInit.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psphotModelTest.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotModelTest.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotModelTest.c (revision 22339) @@ -0,0 +1,11 @@ +# include "psphot.h" + +int main (int argc, char **argv) { + + psMetadata *config = modelTestArguments (&argc, argv); + eamReadout *imdata = psphotSetup (config); + + modelTestFitSource (imdata, config); + + exit (0); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotOutput.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotOutput.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotOutput.c (revision 22339) @@ -0,0 +1,655 @@ +# 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; + + char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE"); + char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE"); + char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE"); + + if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage); + + 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, "TEXT")) { + pmSourcesWriteText (imdata, config, outputFile, sources, psf, sky); + 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; + } + psAbort ("psphot", "unknown output mode %s", outputMode); +} + +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.flt.dat", filename); + pmModelWriteFLTs (sources, name); + + sprintf (name, "%s.nul.dat", filename); + pmModelWriteNULLs (sources, name); + + sprintf (name, "%s.mnt.dat", filename); + pmMomentsWriteText (sources, name); + + psFree (name); + return true; +} + +// 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, "PSF_FIT_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", 32.0); + 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, "PSF_FIT_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", 32.0); // 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, apResid, lsky; + bool status; + + // find config information for output header + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); + char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE"); + + // write necessary information to output header + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT", PS_DATA_F32 | PS_META_REPLACE, "zero point", ZERO_POINT); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code", PHOTCODE); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + + // 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); + apResid = source->apMag - source->fitMag; + 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 (32.0, source->fitMag + ZERO_POINT)); + psLineAdd (line, "%03d ", (int)(1000*dmag)); + psLineAdd (line, "%2d ", type); + psLineAdd (line, "%3.1f ", lsky); + psLineAdd (line, "%6.3f ", 32.0); // should be 'Mgal + psLineAdd (line, "%6.3f ", PS_MIN (32.0, 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, "PSF_FIT_RADIUS"); + float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); + char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE"); + + // write necessary information to output header + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT", PS_DATA_F32 | PS_META_REPLACE, "zero point", ZERO_POINT); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code", PHOTCODE); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + + // set NAXIS to 0 + fprintf (stderr, "setting naxis\n"); + mdi = psMetadataLookup (imdata->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + table = psArrayAlloc (sources->n); + sources->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_HEAD, "X_PIX", PS_DATA_F32, "", PAR[2]); + psMetadataAdd (row, PS_LIST_HEAD, "Y_PIX", PS_DATA_F32, "", PAR[3]); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_GAL", PS_DATA_F32, "", type); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_AP", PS_DATA_F32, "", lsky); + psMetadataAdd (row, PS_LIST_HEAD, "LOG_SKY", PS_DATA_F32, "", 32.0); // should be 'Mgal' + psMetadataAdd (row, PS_LIST_HEAD, "FWHM_X", PS_DATA_F32, "", source->apMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_HEAD, "FWHM_Y", PS_DATA_F32, "", PAR[4]); // should be 'FHWM x' + psMetadataAdd (row, PS_LIST_HEAD, "THETA", PS_DATA_F32, "", PAR[5]); // should be 'FHWM y' + psMetadataAdd (row, PS_LIST_HEAD, "DOPHOT", PS_TYPE_S8, "", 0); // should be theta + // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY", PS_DATA_STRING, "", NULL); + + psArrayAdd (table, 100, row); + psFree (row); + } + + // psAbort ("psphotOutput", "FITS table output not finished"); + + // psImage *image = psImageAlloc (1, 1, PS_TYPE_F32); + // XXX missing: theader = psFitsHeaderFromTable (NULL, table, "SMPFILE"); + // psFitsWriteImage (fits, imdata->header, image, 0); + + theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteHeader (imdata->header, fits); + psFitsWriteTable (fits, theader, table); + psFitsClose (fits); + + return true; +} + +/***** Text Output Methods *****/ + +// 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; +} + +// dump the sources to an output file +bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf) { + + double dP, dmag; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + bool status; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename); + return false; + } + + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // valid source types for this function + if (source->type == PM_SOURCE_SATSTAR) goto valid; + if (source->type == PM_SOURCE_PSFSTAR) goto valid; + if (source->type == PM_SOURCE_GOODSTAR) goto valid; + continue; + + valid: + model = pmSourceMagnitudes (source, psf, RADIUS); + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[2]); + dP += PS_SQR(dPAR[3]); + dP = sqrt (dP); + + dmag = dPAR[1] / PAR[1]; + + fprintf (f, "%7.1f %7.1f %5.1f %7.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP); + + for (j = 0; j < model->params->n - 4; j++) { + fprintf (f, "%9.6f ", PAR[j+4]); + } + fprintf (f, " : "); + for (j = 0; j < model->params->n - 4; j++) { + fprintf (f, "%9.6f ", dPAR[j+4]); + } + fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", + source[0].type, + 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 pmModelWriteFLTs (psArray *sources, char *filename) { + + double dP; + int i, j; + FILE *f; + psVector *params; + psVector *dparams; + pmModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWriteFLTs", 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]; + model = (pmModel *) source->modelFLT; + if (model == NULL) continue; + if (source->type == PM_SOURCE_GALAXY) goto valid; + continue; + + valid: + params = model->params; + dparams = model->dparams; + + // XXX these are hardwired for SGAUSS : this should be pushed into the + // model functions as an abstract function + dP = 0; + dP += PS_SQR(dparams[0].data.F32[4] / params[0].data.F32[4]); + dP += PS_SQR(dparams[0].data.F32[5] / params[0].data.F32[5]); + dP += PS_SQR(dparams[0].data.F32[7] / params[0].data.F32[7]); + dP = sqrt (dP); + + fprintf (f, "%7.1f %7.1f %5.1f %7.3f %7.4f %7.4f ", + params[0].data.F32[2], + params[0].data.F32[3], + params[0].data.F32[0], + -2.5*log10(params[0].data.F32[1]), + (dparams[0].data.F32[1]/params[0].data.F32[1]), + dP); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", params[0].data.F32[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dparams[0].data.F32[j]); + } + fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", + source[0].type, 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; + + 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++) { + pmSource *source = (pmSource *) sources->data[i]; + + // skip these sources (in PSF or FLT) + if (source->type == PM_SOURCE_GALAXY) continue; + if (source->type == PM_SOURCE_PSFSTAR) continue; + if (source->type == PM_SOURCE_SATSTAR) continue; + if (source->type == PM_SOURCE_GOODSTAR) continue; + + if (source->moments == NULL) { + 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, + empty->x, empty->y, + empty->Sx, empty->Sy, + empty->Sum, empty->Peak, + empty->Sky, empty->SN, + empty->nPixels, source->type); + } else { + 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; + + 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 = (pmSource *) sources->data[i]; + if (source->moments == NULL) + continue; + 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; +} + +// 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_SATSTAR: + return (10); + + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_GOODSTAR: + return (1); + + case PM_SOURCE_POOR_FIT_PSF: + return (7); + + case PM_SOURCE_FAIL_FIT_PSF: + case PM_SOURCE_FAINTSTAR: + return (4); + + case PM_SOURCE_GALAXY: + case PM_SOURCE_FAINT_GALAXY: + case PM_SOURCE_DROP_GALAXY: + case PM_SOURCE_FAIL_FIT_GAL: + case PM_SOURCE_POOR_FIT_GAL: + return (2); + + case PM_SOURCE_OTHER: + 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 *modelFLT, pmPSF *psf, int x, int y, int dx, int dy) { + + psRegion region = {x - dx, x + dx, y - dy, y + dy}; + psImage *sample = psImageSubset (image, region); + psImageInit (sample, 0); + modelFLT->params->data.F32[2] = x; + modelFLT->params->data.F32[3] = y; + pmModel *modelPSF = pmModelFromPSF (modelFLT, psf); + pmSourceAddModel (sample, NULL, modelPSF, false, false); + psFree (modelPSF); + return (sample); +} + + +bool psphotSamplePSFs (pmPSF *psf, psImage *image) { + + // make sample PSFs for 4 corners and the center + psImage *sample; + + pmModel *modelFLT = pmModelAlloc (psf->type); + modelFLT->params->data.F32[0] = 0; + modelFLT->params->data.F32[1] = 1; + + psFits *fits = psFitsOpen ("sample.psf.fits", "w"); + + sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, image->numRows - 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelFLT, psf, 25, image->numRows - 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols / 2, image->numRows / 2, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + + psFitsClose (fits); + return (TRUE); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotReapplyPSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotReapplyPSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotReapplyPSF.c (revision 22339) @@ -0,0 +1,86 @@ +# include "psphot.h" + +bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool status; + float x, y, r; + 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); + + // this should be added to the PM_SOURCE flags: + int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + + // we may set this differently here from the value used to mark likely saturated stars + float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); + float PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); + float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + + 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->type == PM_SOURCE_BLEND) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + pmModel *model = source->modelPSF; + pmModel *save = pmModelCopy (model); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + r = model->radius; + + // mask outside radius to limit the following actions to the source + psImageKeepCircle (source->mask, x, y, r, "OR", PSPHOT_MASK_MARKED); + + // replace object in image + pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false); + + // fit PSF model (set/unset the pixel mask) + pmSourceFitModel (source, model, true); + + // model succeeded : tentatively keep it + Niter += model[0].nIter; + Nfit ++; + + // is it a good PSF model? + if (psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION)) { + psFree (save); + Nsub ++; + } else { + psFree (model); + source->modelPSF = save; + } + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + + // clear the masked region + psImageKeepCircle (source->mask, x, y, r, "AND", ~PSPHOT_MASK_MARKED); + } + + 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); +} + +pmModel *pmModelCopy (pmModel *model) { + + pmModel *new = pmModelAlloc (model->type); + + new->chisq = model->chisq; + new->nIter = model->nIter; + + 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); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotSetup.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotSetup.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotSetup.c (revision 22339) @@ -0,0 +1,109 @@ +# 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"); + + // 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] = image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN); + } + } + } + + // 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", 4, "load data: %f sec\n", psTimerMark ("psphot")); + + // save the image data & return it + eamReadout *imdata = eamReadoutAlloc(image, weight, mask, header); + return (imdata); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotSortBySN.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotSortBySN.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotSortBySN.c (revision 22339) @@ -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/eam_tags/psphot_dev_06/psphot/src/psphotSourceStats.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotSourceStats.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotSourceStats.c (revision 22339) @@ -0,0 +1,53 @@ +# include "psphot.h" + +psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *peaks) +{ + bool status = false; + psArray *sources = NULL; + + 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 + // XXX EAM : choose between these two versions + status = pmSourceMoments (source, RADIUS); + if (!status) { + psFree (source); + continue; + } + + // add to the source array + psArrayAdd (sources, 100, source); + psFree (source); + } + + psLogMsg ("psphot", 3, "%d moments: %f sec\n", sources->n, psTimerMark ("psphot")); + return (sources); +} + +// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky) Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotSparseMatrix.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotSparseMatrix.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotSparseMatrix.c (revision 22339) @@ -0,0 +1,116 @@ +# include "psphot.h" + +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { + + bool status; + float x; + float y; + float Sky; + int Nfit = 0; + + psTimerStart ("psphot"); + + // source analysis is done in spatial order + sources = psArraySort (sources, psphotSortByY); + + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_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.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + + // vector for model input coordinates + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + + // this function specifies the radius at this the model hits the given flux + pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); + + // this is the function to get the model function + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + // pre-calculate all model pixels + psArray *pixlist = psArrayAlloc (sources->n); + for (int k = 0; k < sources->n; k++) { + source = sources->data[k]; + + // use the source moments, etc to guess basic model parameters + pmModel *modelFLT = pmSourceModelGuess (source, psf->type); + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelFLT, psf); + // model->params->data.F32[2] = source->peak->x; + // model->params->data.F32[3] = source->peak->y; + model->params->data.F32[1] = 1.0; + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + psFree (modelFLT); + + // 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 = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) { + psAbort ("apply_psf_model", "error in radius"); + } + + // build the model image + nx = src->pixels->numCols; + ny = src->pixels->numRows; + im = psImageAlloc (nx, ny, PS_DATA_F32); + mk = psImageAlloc (nx, ny, PS_DATA_U8); + + // fill in the model pixel values + psImageKeepCircle (mk, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + psVector *params = model->params; + for (int i = 0; i < nx; i++) { + for (int j = 0; j < ny; j++) { + if (mk->data.U8[j][i]) continue; + x->data.F32[0] = (float) i + pix->col0; + x->data.F32[1] = (float) j + pix->row0; + im->data.F32[j][i] = modelFunc (NULL, params, x); + } + } + pixlist->data[k] = im; + } + + // fill out the sparse matrix + psSparse *sparse = psSparseAlloc (sources->n, 100); + psVector *Aij = sparse->Aij; + psVector *Bfj = sparse->Bfj; + psVector *Si = sparse->Si; + psVector *Sj = sparse->Sj; + + int Ri, Rj; + + Nelem = 0; + for (i = 0; i < pixlist->n; i++) { + pf = ((psSource *)sources->data[k])->pixels; + pi = pixlist->data[i]; + + // diagonal element (auto-cross-product) + f = psphotCrossProduct (pi, pi); + psSparseMatrixElement (sparse, i, i, f); + + // find the image x model value + f = psphotCrossProduct (pi, pf); + psSparseVectorElement (sparse, i, f); + + // loop over all other stars following this one + for (j = i + 1; j < sources->n; j++) { + pj = ((pmSource *)sources->data[j])->pixels; + + if (pi->col0 + pi->numCols < pj->col0) break; + if (pj->col0 + pj->numCols < pi->col0) continue; + if (pi->row0 + pi->numRows < pj->row0) continue; + if (pj->row0 + pj->numRows < pi->row0) continue; + + // got an overlap; calculate cross-product and add to output array + f = psphotCrossProduct (pi, pj); + psSparseMatrixElement (sparse, j, i, f); + } + } + + psSparseResort (sparse); + return (match); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotSubtractPSF.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotSubtractPSF.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotSubtractPSF.c (revision 22339) @@ -0,0 +1,29 @@ +# include "psphot.h" + +// subtract model fits from image (only PSFSTAR && BRIGHTSTAR) +// need to control application of this with some thresholding + +bool psphotSubtractPSF (pmSource *source) +{ + pmModel *model; + + // only subtract successful fits + if (source->type == PM_SOURCE_SATSTAR) goto valid; + if (source->type == PM_SOURCE_PSFSTAR) goto valid; + if (source->type == PM_SOURCE_GOODSTAR) goto valid; + return (false); + +valid: + model = source->modelPSF; + if (model == NULL) { + psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y); + return (false); + } + + // subtract object, leave local sky + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + // XXX EAM : amplify the weight matrix a la dophot? + + return (true); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotTest.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotTest.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotTest.c (revision 22339) @@ -0,0 +1,34 @@ +# include "psphot.h" + +int main (int argc, char **argv) { + + bool status; + + psMetadata *config = psphotTestArguments (&argc, argv); + eamReadout *imdata = psphotSetup (config); + + int x = psMetadataLookupF32 (&status, config, "TEST_FIT_X"); + int y = psMetadataLookupF32 (&status, config, "TEST_FIT_Y"); + float fraction = psMetadataLookupF32 (&status, config, "TEST_FRACTION"); + float threshold = psMetadataLookupF32 (&status, config, "TEST_THRESHOLD"); + + if (threshold == -1) { + threshold = fraction * imdata->image->data.F32[y][x]; + } + + psArray *contour = pmSourceContour_EAM (imdata->image, x, y, threshold); + if (contour == NULL) { + fprintf (stderr, "error: invalid contour\n"); + exit (1); + } + + psVector *xV = contour->data[0]; + psVector *yV = contour->data[1]; + + for (int i = 0; i < xV->n; i++) { + fprintf (stdout, "%d %f %f\n", i, xV->data.F32[i], yV->data.F32[i]); + } + + // psSparseMatrixTest (); + exit (0); +} Index: /tags/eam_tags/psphot_dev_06/psphot/src/psphotTestArguments.c =================================================================== --- /tags/eam_tags/psphot_dev_06/psphot/src/psphotTestArguments.c (revision 22339) +++ /tags/eam_tags/psphot_dev_06/psphot/src/psphotTestArguments.c (revision 22339) @@ -0,0 +1,70 @@ +# include "psphot.h" +static int usage (); + +psMetadata *psphotTestArguments (int *argc, char **argv) { + + int N; + unsigned int Nfail; + int mode = PS_DATA_STRING | PS_META_REPLACE; + psF32 UseCoords_X, UseCoords_Y; + + // basic pslib options + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + // identify options in args list - these go on config + bool UseCoords = false; + if ((N = psArgumentGet (*argc, argv, "-coords"))) { + UseCoords = true; + psArgumentRemove (N, argc, argv); + UseCoords_X = atof (argv[N]); + psArgumentRemove (N, argc, argv); + UseCoords_Y = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + float threshold = -1; + if ((N = psArgumentGet (*argc, argv, "-threshold"))) { + psArgumentRemove (N, argc, argv); + threshold = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + float fraction = 0.5; + if ((N = psArgumentGet (*argc, argv, "-frac"))) { + psArgumentRemove (N, argc, argv); + fraction = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + if (*argc != 3) usage (); + + // load config information + psMetadata *config = psMetadataAlloc (); + psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); + config = psMetadataConfigParse (config, &Nfail, argv[2], FALSE); + + // identify input image + psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); + + mode = PS_DATA_F32 | PS_META_REPLACE; + psMetadataAdd (config, PS_LIST_HEAD, "TEST_THRESHOLD", mode, "", threshold); + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FRACTION", mode, "", fraction); + + if (UseCoords) { + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", mode, "", UseCoords_X); + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", mode, "", UseCoords_Y); + } + + return (config); +} + +static int usage () { + + fprintf (stderr, "USAGE: psphotModelTest (image.fits) (config)\n"); + fprintf (stderr, "options: \n"); + fprintf (stderr, " -coords x y : specify object center\n"); + fprintf (stderr, " -threshold (threshold) : (defaults to FRAC of peak\n"); + fprintf (stderr, " -frac (FRAC) : specify how to determine threshold\n"); + exit (2); +}