Index: /tags/RHL_v1_2branch_0/psphot/.cvsignore =================================================================== --- /tags/RHL_v1_2branch_0/psphot/.cvsignore (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/.cvsignore (revision 10106) @@ -0,0 +1,18 @@ +bin +Makefile +Makefile.in +aclocal.m4 +autom4te.cache +compile +config.log +config.status +configure +depcomp +install-sh +missing +config.guess +libtool +ltmain.sh +stamp-h1 +config.sub +psphot.pc Index: /tags/RHL_v1_2branch_0/psphot/Makefile.am =================================================================== --- /tags/RHL_v1_2branch_0/psphot/Makefile.am (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/Makefile.am (revision 10106) @@ -0,0 +1,10 @@ +SUBDIRS = src + +CLEANFILES = *.pyc *~ core core.* + +pkgconfigdir = $(libdir)/pkgconfig +pkgconfig_DATA= psphot.pc + +EXTRA_DIST = \ + psphot.pc.in \ + autogen.sh Index: /tags/RHL_v1_2branch_0/psphot/autogen.sh =================================================================== --- /tags/RHL_v1_2branch_0/psphot/autogen.sh (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/autogen.sh (revision 10106) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=psphot +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL=aclocal +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOlIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /tags/RHL_v1_2branch_0/psphot/configure.ac =================================================================== --- /tags/RHL_v1_2branch_0/psphot/configure.ac (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/configure.ac (revision 10106) @@ -0,0 +1,58 @@ +dnl Process this file with autoconf to produce a configure script. +AC_PREREQ(2.59) + +AC_INIT([psphot], [0.0.1], [http://pan-starrs.ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC +AC_PROG_INSTALL +AC_PROG_LIBTOOL + +dnl enable largefile supports +AC_SYS_LARGEFILE + +dnl handle debug building +AC_ARG_ENABLE(optimize, + [AS_HELP_STRING(--enable-optimize,enable compiler optimization)], + [AC_MSG_RESULT(compile optimization enabled) + CFLAGS="-O2"], + [AC_MSG_RESULT([compile optimization disabled]) + CFLAGS="-O0 -g"] +) + +dnl handle profiler building +AC_ARG_ENABLE(profile, + [AS_HELP_STRING(--enable-profile,enable compiler profiler information inclusion)], + [AC_MSG_RESULT(compile optimization enabled) + CFLAGS="${CFLAGS=} -pg"] +) + +dnl ------------------------------------------------------------ + +AC_PATH_PROG([ERRORCODES], [psParseErrorCodes], [missing]) +if test "$ERRORCODES" = "missing" ; then + AC_MSG_ERROR([psParseErrorCodes is required]) +fi + +PKG_CHECK_MODULES(PSLIB, pslib >= 0.12.99) +PKG_CHECK_MODULES(PSMODULE, psmodules >= 0.12.99) + +dnl Set CFLAGS for build +CFLAGS="${CFLAGS} -Wall -Werror -std=c99" + +AC_SUBST([PSPHOT_CFLAGS]) +AC_SUBST([PSPHOT_LIBS]) + +AC_CONFIG_FILES([ + Makefile + src/Makefile + psphot.pc +]) + +AC_OUTPUT Index: /tags/RHL_v1_2branch_0/psphot/doc/config.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/config.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/config.txt (revision 10106) @@ -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/RHL_v1_2branch_0/psphot/doc/notes.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/notes.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/notes.txt (revision 10106) @@ -0,0 +1,593 @@ + +ensemble: + + * first block: + * select sources of interest / skip + * generate a guess model + * re-measure moments for saturated stars (MOVE) + * update peak position (MOVE to peaks) + * set the object radius + * build temporary source copies (including pixels, mask, weight) + + * second block: + * measure the cross-products and weights + * add to the sparse matrix + + * build full-image mask + * measure coord weights + * build order elements + + * solve for source amplitudes + * update models + +2006.11.16 + + * create psSparseBorder to solve matrix equations which have a large + sparse region, with a completely filled border of a smaller number + of rows. + + * trying to modify psphotEnsemble to (1) be more general (not just + PSFs) and (2) fit for the sky offset. + + - how is the radius for each source set? + + +2006.11.11 + + I am trying to define a complete set of tests to validate the + measurement of the photometry compared with the simulated images. + Part of my confusion is controlling exactly which corrections are + applied to the output magnitudes. Here is the complete list: + + psfMag_out = psfMag_meas + ApTrend(x,y,0.0,0.0) + apMag_out = apMag_meas + growth(radius,refRadius) + f(skyBias,skySat) + + these corrections are modified by the following config settings: + + IGNORE_GROWTH (sets all growth magnitudes to the fitMag value, + effectively turning off the growth correction) + + APTREND : set the level of correction + NONE : no correction is fitted + CONSTANT : just a single constant (applied to get psfMag_out) + SKYBIAS : the skybias term is fitted (constant applied to psfMag, + slope applied apMag) + +2006.11.08 + + I have been testing the pmSourceMagnitudes interpolation of the + source position when measuing the aperture magnitudes. I built a + test suite for pmSourceMagnitudes and pmGrowthCurve. even with + bicubic interpolation, there are noticable errors for small + apertures (ie, small compared to the PSF sigma). I used a gaussian + PSF to make the tests. it seems that, if the radius is less than 2 + sigma or so, the interpolated aperture mag is in error by as much as + 12mmag. The size of this effect must depend on the smoothing + algorithm, the shape of the PSF, the size of the PSF relative to the + pixels, and the size of the subpixel interpolation. I believe the + effect is a result of the smoothing introduced by the interpolation: + the smoothing pushes the flux in the inner portions of the PSF down + and enhances the flux in the outer portions. thus, there is + actually a turn over in the error: at very large values, the error + is very small because the aperture contains all of the flux, and the + interpolation has little effect. as the radius shrinks, the error + grows (in the sense of M_raw - M_interpolated; ie, the flux in the + interpolated source is too large). around 2-3sigma, the error turns + over and becomes negative (too little flux in the interpolated + source). + + The sigmas I used for these tests (1.0, 1.5, 2.0 pix) are + appropriate for the range we expect to see from PS1, or any + well-sampled detector. Thus, even though the effect is probably + less important for poorer seeing, we have to be careful to avoid it + in our Pan-STARRS analysis. + + Currently, the minimum aperture is defined only by the config + variable PSF_FIT_PAD, which really implies a padding. My proposal + is to require the minimum aperture to be a number of sigma (an + additional parameter). Then, the formula for the radius would be + something like: + + MAX (modelRadiusPSF (PSF_FIT_NSIGMA), PSF_FIT_RADIUS_NSIG) + + PSF_FIT_PAD + + with PSF_FIT_RADIUS_NSIG (choose a better name!) set to something + like 2 or 3. + +2006.11.06 + + I have added the aperture interpolation as a component of the + pmSourceMagnitudes function. + + I am trying to make things a bit more organized, with more + consistent APIs. + + * changed pmPSFtry to use the pmSourceMagnitudes API + * changed pmPSFtry to use new pmSource entries for each try (saving + * the pmModels generated) + + other notes: + * moved psphotGrowthCurve -> pmGrowthCurveGenerate + + things I know need to be added / fixed: + * measure bicube fit to peak from smoothed image + * use psphotEnsemble fit to re-measure normalization? + * extended object statistical measures + +2006.11.02 + + Looking into the growth curve problem. RHL sent his recipe for + correcting to true circular apertures. I've done a test fix of + interpolating the stellar image to the desired center. it looks + like bilinear interpolation is insufficient (about 8% error for 2 + pixel radius aperture). Unless I've got the formula wrong, the + effective smoothing of the gaussian is moving flux from the center + out and making small apertures in error that depends on the offset. + Bicubic could be enough to clean this up? + + I think the current code may need to be careful modified to include + the interpolation for every aperture measurement. there are paired + 'pmModelAdd' and pmModelSub steps to add and remove the model flux. + I need to be careful about how the mask affects these steps. + + I also need to be careful in interpreting the test results: since + the simulated images create objects by filling in a pixel and then + smoothing, all objects have centroids which are exactly 0.5 pixel + offsets. These all then have a fixed growth correction error, which + fortunately happens to be substantial (because of a 0.5 pixel offset + in the growth curve reference object!) + + I need to double check the 0.5 pixel error issues. There are + problems with the pslib version of psImageShift and + psImageInterpolate (or at least inconsistencies). the pmModelAdd + and Sub functions and the centroid functions are probably also + inconsistent. + +2006.11.01 + + Trying to understand low-level errors in psphot using simulated + data. the main issue I am concerned with is the bias RHL is seeing + in SDSS-psphot comparisons. But there are other issues that are + creeping in and probably point at logical errors in the code. + + * in pmPSFtry, I measure the first pass on the aperture residual (to + choose between PSF model options, if that is being tested). This + analysis currently uses a polynomial fit (1st order) for the offset + and the sky bias (using r^2/flux as an independent variable). + + Notes: + + - is the psf-fit error distribution reasonable? perhaps a little + on the high side: at 10000 DN, noise is 0.0138, should be + 0.0107. + + - scatter of ap-fit is surprisingly high. M_input - M_fit has a + median chisq of 2 (mean of 5), while M_input - M_ap has a median + chisq of 125 (mean of 380). This is using a radius of 15. It + cleans up hugely when I go to a radius of 7 and even more so + with radius = 4. (In this case, the PSF sigma is 1.5, so 4 is + well matched and 7 is only slightly large). + + * make the initial PSF fitting radius and aperture radius + different? + * make the initial PSF fitting radius a function of the measured + PSF sigma? + + - I have been fitting ap-fit mag without weights; the resulting + ap-fit errors are highly over-estimated since the ap-fit scatter + is too high. I have added code to carry the weights (error from + the fit only at the moment, ie, not the ap error as well). This + is good, but not as important if the radius is chosen well above. + + * a bias in the psphotEnsemble results depending on the sky level + but not in the PSF model fits? turning on/off CONST_PHOT_WT has + an equivalent effect, but is quite small. why is one sensitive to + the sky and the other is not? + + * the error measured by psphotEnsemble was wrong (and inconsistent + for a given weighting situation): fixed the def of the error. + + * still trying to understand the corrections being made. running a + few tests: + + keeping PSF_FIT_RADIUS at 4.0 + + t1 : no growth correction, constant weights + t2 : growth correction, constant weights + t3 : no growth correction, no constant weights + t4 : growth correction, no constant weights + t5 : growth correction, constant weights, constant ap (4) + t6 : growth correction, constant weights, integer radius + t7 : growth correction, constant weights, integer radius + (for t7, I've fixed the bug that the growth correction include the + min radius). + + * psphotGrowthCurve is messed up : the measurement does not correct + for the actual object center, and at small radii this is an + important source of noise. + +2006.06.28 + + I am adding a few minor features. First up is the ability to break + the processing at interesting stages. These will be: + + - PEAKS: after sources are generated, but before moments are + measured + - MOMENTS: after moments are measured, before the PSF is measured + - PSFMODEL: after the psf model is measured, before sources are + fitted in bulk + - ENSEMBLE: after the linear ensemble fitting + - DEFAULT: complete processing + +2006.04.22 + + I have updated psphot to work with the cycle 11 release of psLib and + psModules. At this point, the IfA group is now working off of the + CVS head rather than from a separate IfA branch of psLib and + psModules. The cycle 11 release of psLib needed a few updates to + fix minor problems. The psModules code needed significant work to + fix discrepancies. The new version of psphot has also been leak + checked, at least for basic configurations. This new version of + psphot should be used with the new camera configuration system, + which can be found in the ipp/config directory. Here are the + necessary tags: + + psphot: psphot_dev_11_0 + psModules: rel11_ifa_pre0 + psLib: rel11_ifa_0 + +2006.04.16 + + Examining some of the throughput issues. the flux measurement speed + is strongly dependent on the interation of the function step size: + + dz = 0.01 : 3.24s for 4400 objects + dz = 0.1 : 0.44s for 4400 objects + + there is an increased error in the result with the larger step size, + but it is smaller than the phot error. Still, we could do better + with a smarter integration function. + + Including the aperture photometry on all sources, but not the + pixWeight, the source magnitudes takes 4.8 sec for 4400 objects. + + Adding the pixWeight increases the time to 7.5 sec for 4400 objects + +2006.02.12 + + Last week, I made a new tag (psphot_dev_08) after finishing the work + to clean memory leaks and to modify psphot to work with Paul's + ppImage infrastructure. + + I have added a new background model function which generates a + median map based on Robert's suggestion of a) half-step windows and + b) linear interpolation. + + Now, I am working to remove the constant sky level concepts and make + every relevant test use the signal/noise ratio as the trigger. This + has impacts on the peak detections, the threshold for blends, and + the threshold for choosing pixels in the fitting routines. + +2005.12.23 + + Improvements in psphot as of psphot_dev_07 + + - updated to work with psModule release 9.0 + - moved model test functions into psphot as optional mode + - added pmSourceFitSet : simultaneous model group fitting + - added bicubic interpolation to tweak up the peak coordinates + - added aperture residual analysis after complete object model & subtraction + - added 3D (x,y,rflux) aperture residual fitting + - added blend vs blob fitting (BlendFit) + - added ensemble PSF fitting, with sparse matrix inversion + - improved the flagging / marking of model results + - better abstraction for setting fit parameters + - various low-level fixes + - added PSF output option + - added break points + - added more optional output formats + - fixed FITS table output + + Things still to be added to psphot: + + - better initial background model (NxM median grid?) + - include local background model in EnsemblePSF Fit (use sparse matrix) + - 2nd pass (after object subtraction) to smooth and detect faint objects + - re-add background to object-subtracted image + - cleanup function names + - push more functions into psLib / psModule + - fix extended / galaxy models to improve stability + - add mode with input PSF + - add mode with input fixed sources (PS_SOURCE_FIXED) + - cleanup memory leaks + +2005.11.25 + + I've updated psphot to work with the current psLib v8, though a + number of psLib fixes were needed. These are pushed under + psLib:eam_rel8_b2. The psphot code which works with that version of + psLib is tagged psphot_dev_04. + + I'm working on converting psphot to work with the current release of + psModules, which basically includes all of the object functions. + + - change comments are relative to psModules (vs psphot versions) + + psEllipse.c : OK + psEllipse.h : OK + pmModelGroup.c : OK (missing TGAUSS entry) + pmModelGroup.h : OK (adds many comments) + + pmPSFtry.c : a bunch of formatting weirdness (try-> got line breaks + added pmPSFtryMetric_Alt + fixed up usage of psVectorClipFitPolynomial1D + fixed up usage of psPolynomial..Alloc + + pmPSFtry.h : OK (Added extensive comments) + + pmPSF.c : fixed up usage of psPolynomial..Alloc + + pmPSF.h : OK (Added extensive comments) + + pmObjects.c : cleaned up some formatting, + fixed usage of psImageCountPixelMask + + pmObjects.h : added apMag and fitMag to pmSource + dropped PS_MODEL_name defines + + psLibUtils.[ch] should be dropped from psModules + psModulesUtils.[ch] should be dropped from psModules + + Makefile.am : dropped psModulesUtils.[ch], psLibUtils.[ch] + + - I have successfully tested psphot with the psModules.v8 code + (eam_rel8_b1 tag). I only needed to make a few modes in the + psModules code. + +2005.11.16 + + I have made some fixes to make psphot work with pslib rel8_0. I + need to merge psphot with the rel8_0 version of the psMinimize and + pmObjects code. + +2005.09.06 + + I have built a working version of PSPhot using my own copy of the + pmObjects.[ch] files. I need to minimize the difference between the + v0.7.0, v0.5.0, and my own version of pmObjects.c + + additions to psLib based on the new version: + - extended the minimization to include parameter limits + + changes v0.5.0 / v0.7.0 / EAM + - changed all ps* to pm* + - fixed naming for local static functions (eg modelFree not p_psModelFree) + - dropped all hard-coded model names (vs v0.7.0) + - using current psMemSetDellocator name + - added noise image to pmSource + - converted old asserts to new asserts (PS_.._CHECK..) + - add return from error in pmFindImagePeaks + - allow isItInThisRegion to have NULL region (true) + - fixed subImages to correspond to region definition (end is + exclusive) in LocalSky and SetPixelCircle + - dropped old, redundant code in pmSourceLocalSky + - fixed pmSourceMoments use of mask (0 is valid, not 1) + - added pmSourceMoments validity tests (numPixels, Sum, centroid shift) + - added correction for Sxy (XY/Sum - x*y) + - added pmSourcePSFClump function (and pmPSFClump structure) + - allow missing sources in pmSourcePSFClump (not all peaks yield + sources) + - moved pmSourceRoughClass metadata lookups out of inner loop + - added saturated pixel test using pmImageCountPixelMask + - various redefinitions of initial classes + - dropped old, redundant code in pmSourceSetPixelsCircle + - mask image set to type psU8 (was psF32) + - abstract object model functions for pmSourceModelGuess + - converted evalModel to a public function psModelEval, modified params + - pmSourceFitModel operates on a specified model + - pmSourceFitModel had option to apply/ignore PSF information + - pmSourceFitModel calculates yErr + - pmSourceFitModel calculates the covariance matrix + - added paramMask to pmSourceFitModel + - pmSourceFitModel_EAM uses covar to carry in beta and param + limits... + - pmSourceFitModel requires solution to stay within image + - pmSourceFitModel calculates and saves chisq, nIter, nDOF + - pmSourceFitModel uses GaussNewtonDelta for frozen params + - added mask to pmSourceAdd, pmSourceSub + - model abstractions in pmSourceAdd, pmSourceSub + - fixed 'center' option in pmSourceAdd/Sub + + - pmPSF params should be F64 + + lingering concerns: + - Polynomials and related functions are STILL messed up. this is a + major priority!!! + - pmSourceFitModel calculates sqrt(var), psMinimize calculates sq(sqrt(var)) + - does psImageAlloc zero the image? if so, drop init routine + from pmSourcePSFClump + - pmSourceRoughClass uses RDNOISE, GAIN to calculate SN: should use + the source->noise image? part of pmSourceMoments? + - pmSourceMoments_EAM uses macro for checkRegion + - pmObjects_EAM.c uses some hardcoded mask values (LocalSky) + - why does LocalSky_EAM not need subImageMask->col0 = subImage->col0?? + - psGetRowVectorFromImage should be in psLib (in psLibUtils) + - pmSourceSetPixelsCircle and pmSourceLocalSky probably could use + the psImageMaskRegion, etc functions. + - pmSourceSetPixelsCircle uses a hard-coded mask value + - define a mask registration process?? + - pmSourceContour is using the model flux, not the image flux + +2005.09.05 + + I have a working version of PSPhot which handles PSF and galaxy + models. I am tagging the module with an alpha version number. this + version compiles with psLib-0.5.0, at least my version of it. My + current goal is to make it work with the current psLib head (nearly + 0.7.0) with a mimimum of _EAM versions of functions. + + the tag for the working version is dev_01 + +2005.06.04 + + progress on psphot has moved along well. At this point, the process + loads the image, finds the peaks, determines a best PSF model, + identifies objects which are consistent with that model, and + attempts to fit a galaxy model to the objects which are not + consistent with the model. + + a variety of minor issues remain, as well as some major issues. + + - use a proper noise image to keep the fits honest after other + objects have been subtracted. + + - define the masks in terms of a single mask image. this image + could be provided initially by the user. + + - add noise enhancement (couple with mask marking subtraction) + + - + +2005.04.12 + +- psPeak, psSource, etc: should be pmPeak, pmSource, etc +- pmCullPeaks: should be pmCullImagePeaks + +- does pmSourceMoments need an image argument? + +- psLibInit / p_psTimeInit + + we should not have to know about astronomy time concepts (ie, have a + timeConfig file) in order to use the time functions for basic work. + This means having a PS_TIME_UNIX which just works with the UNIX + clock and avoids all initialization issues. + +- psImageStats is much too slow for basic median + - this function does a complete sort on all pixel values. we need + to use the histogram method (spec a different type of median) + +- psImageStats, ROBUST_MEDIAN fails +- there are errors in psVectorStats for ROBUST MEAN, STDEV (& median?) + * this is ill-defined. revise the ADD: too sensitive to the + binning, smoothing scales + +- pmFindImagePeaks is not finding any peaks: + typo in interior row section: see bug 359 + after this fix: PASS + +* p_psGetRowVectorFromImage : PASS +* psGetRowVectorFromImage : PASS +- both have about the same time (~15 us +/- 5us) for 2048 pix on alala + +* pmFindVectorPeaks : PASS (returned row data in data.U32) + +- # include : needed in psMetadata.c + +* MyListAddPeak / pmPeakAlloc are inconsistent wrt col,row vs x,y + +- why does psMetadataItemAllocV allocate a string of length + MAX_STRING_LENGTH for the comment, then use strncpy? + + should be: + + // set metadata item comment + if (comment == NULL) { + // Per SDRS, null isn't allowed, must use "" instead + metadataItem->comment = psStringCopy (""); + } else { + metadataItem->comment = psStringCopy (comment); + } + +psMetadataAdd -> psMetadataAddItem is all confused: + - we end up with two psMetadataItemAlloc calls, when only one should + be used + - psMetadataItemAlloc was doing the wrong thing if the incoming data + was NULL + - I think the logic in this block is wrong as well. + + * I fixed this in psMetadata.c + +2005.04.14 : + + - psImageStats is still much too slow + + - pmObjects : changed to using psArray to carry peaks / sources, + rather than psList + + - added function pmPeaksSubset to replace pmCullPeaks + + - pmSourceLocalSky: added min/max fncs to force subimage to stay in + image + + - Sx, Sy can go negative: have forced limit to 0,0 + + - pmSourceRoughClass: implemented sigmaX, sigmaY search for stellar + clump (uses pmComparePeakDescend) + + - the search ignores the 0,0 pixel + + - pmSourceRoughClass uses all detected peaks. allow an exclusion + for peak > max, peak < min? + + - need to push clump stats on the metadata + + - the source classes are assigned so they are mutually exclusive + + - change PS_SOURCE_OTHER to STAR? + + - pmSourceSetPixelCircle: + + - name should be pmSourceSetPixelsCircle (geometry could be an + option...) + + - I cleaned the defined region to match the convensions of + pmSourceLocalSky + + - truncate and force subimage to lie on image + + - CheckRadius2 seems like an inefficient coding + + - valid pixels have mask value of 1 (inverse of other funcs) + + - pmSourceFitModel: + + - num iterations is much too large + + - + +mrq: + + - allocate alpha, beta arrays + - mrq2dcof (pars, alpha, beta) + - make guess pars', alpha', beta' + - mrq2dcof (pars', alpha', beta') + - if (chisq < old chisq), keep new guess + + +2005.04.19: + +- psMinimizeLMChi2 has some serious problems. I re-wrote it starting + from the mrqmin example in ohana. after a variety little bugs, it + seems to work quite well. the program 'modeltest.c' runs two + different gaussians through psMinimizeLMChi2, and gets the right + answer quite quickly. The first pass at the implementation had a + variety of problems. to avoid any worry about errors in the + psMatrix code, I implemented psGaussJordan, basing the code on the + elixir gaussj.c code. this worked fine, and let me find the errors + elsewhere in the code. the LUD version now works just the same as + the GaussJ version. + + +2005.04.21 + +- psMinimizeLMChi2 does an excellent job now in most cases. I am able + to run pmSourceFitModel on the PSF stars quite well. I added in the + (nearly) correct errors (actually, just using sqrt(N), limiting to + N>=1). I see one problem object, which is quite faint, and does not + converge well: I get very large values for x,y, Sx, Sy, etc. I need + to think about constraining these fits from running way off scale. + + fit: 0.123879 sec for 277 stars (gauss) + fit: 0.155940 sec for 277 stars (pgauss) (!) + fit: 0.153292 sec for 277 stars + (0.5 msec per star) Index: /tags/RHL_v1_2branch_0/psphot/doc/outline.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/outline.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/outline.txt (revision 10106) @@ -0,0 +1,21 @@ + +Here is a summary outline of the steps taken by psphotReadout: + +* setup (choose recipe, choose readout, define break-point) +* create mask and weight images if needed, set mask for analysis region +* psphotImageMedian : construct a background model image and subtract it +* psphotFindPeaks : smooth and find the peak pixels from the smoothed image + - no pmSource defined yet, only pmPeak +* psphotSourceStats : create sources for each peak and measure their moments + - pmSource defined here, with pixels covering SKY_OUTER_RADIUS + - no pmModel is defined yet. +* psphotBasicDeblend : identify blended sources by proximity and valley depth +* psphotRoughClass : source classification guess based on moments + - set the source.type,mode. + - re-calculate moments for saturated stars +* psphotChoosePSF : define the psf model + - this generates model fits to the identified psf stars + (one for each tested model) + - these models are not kept (seems like a waste, but later fits must be consistent) +* psphotGuessModels : set the initial PSF model for each object (even EXTENDED) + - creates modelPSF entries Index: /tags/RHL_v1_2branch_0/psphot/doc/output.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/output.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/output.txt (revision 10106) @@ -0,0 +1,15 @@ + +we have several output formats: + - TEXT : separate ASCII tables for PSF, EXT, MNT, and NULL + - SX : sextractor-style ASCII file + - OBJ : dophot-style ASCII file + - CMP : old-style elixir format (header + text) + - CMF : fits-table + +we have several options for organization: + - single chip/readout : outroot.ext + - split chips : outroot.NN.ext + - mef chips : outroot.ext (extensions with names) + + * only CMF will work with MEF model + Index: /tags/RHL_v1_2branch_0/psphot/doc/psfmodel.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/psfmodel.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/psfmodel.txt (revision 10106) @@ -0,0 +1,22 @@ + +2006.10.27 + + I have been working to fix the PSF modeling in psphot. The PSF + model consists of a flux model for an individual object using an + analytical model with a number of parameters. For a collection of + PSF objects, the variation of the parameters as a function of + position are then themselves fitted with a model. The PSF model for + a single object consists of a radial profile with a functional form + f(z) where a given value of z defines an elliptical coutour of the + form z = \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + + \sigma_{xy}xy. + + + + The term \sigma_{xy} is difficult to model as a simple function of x + and y (eg, a low-order polynomial). + + A better model can be constructed for + \frac{\sigma_{xy}}{(\sigma_x^{-2} + \sigma_y^{-2})^2}, which varies + like a^2 \sin 2\theta + Index: /tags/RHL_v1_2branch_0/psphot/doc/psphot.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/psphot.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/psphot.txt (revision 10106) @@ -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/RHL_v1_2branch_0/psphot/doc/regions.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/regions.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/regions.txt (revision 10106) @@ -0,0 +1,159 @@ + +I am having some trouble with functions that refer to both subimages +and regions. there are inconsistencies between the coordinate +conventions. I am going to go through the list of image functions are +decide if we can choose either 1) external coordinates always refer to +parent coordinates or 2) something else! + +bool psImageSet(const psImage *image, int x, int y, double complex value); +double complex psImageGet(const psImage *image, int x, int y); + - sdrs does not specify + - code uses subimage + +psRegion psRegionForImage(psImage *image, psRegion in); + - sdrs says subimage + - code uses subimage + +psImage *psFitsReadImage(psImage *out, const psFits *fits, psRegion region, int z); + - region here specifies a subimage of the disk image + +bool psFitsUpdateImage(psFits *fits, const psImage *input, int x0, int y0, int z); + - x0,y0 here refer to the disk image coordinates + +bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage +*input, int depth, const char *extname); + - note that the output disk image is truncated to the subimage + +bool psFitsInsertImage(psFits *fits, psMetadata *header, const psImage *input, int depth, + const char *extname, bool after); + - same as above + +psImage *psImageSubset(psImage *image, psRegion region); + - code uses region as coordinates in parent image + - sdrs does not specify + - XXX I think the range tests for region fail if the input is a + subimage + +psImage *psImageCopy(psImage *output, const psImage *input, psElemType type); + - maintains the col0,row0 values + +psVector *psImageRow(psVector *out, const psImage *input, int row); +psVector *psImageCol(psVector *out, const psImage *input, int column); + - code uses subimage coordinates for row and column + +psVector *psImageSlice() +psVector* psImageCut() + - code uses subimage coordinates (no correction) + - sdrs does not specify + +psImage *psImageRebin() + - mask and input image must currently agree + +psImage *psImageTransform() + - sdrs is unclear: I think the region defines the output image such + that the pixel in the input image *parent* coordinates which + corresponds to region.x0,y0 in turn corresponds to the output image + 0,0 pixel; the output image always has col0,row0 = 0,0 + +long psImageCountPixelMask (psImage *mask, psRegion region, psMaskType value); + - sdrs says the region corresponds to subimage coordinates + +psPolynomial2D *psImageFitPolynomial(psPolynomial2D *coeffs, const psImage *input); +psImage *psImageEvalPolynomial(psImage *input, const psPolynomial2D *coeffs); + - sdrs does not specify + - polynomial coordinates should refer to the *parent* image + +double complex psImagePixelInterpolate(input, x, y, ...) + - sdrs does not specify + +int psImageOverlaySection(psImage *image, const psImage *overlay, ) + - sdrs does not specify + +psImage *psPixelsToMask() +psPixels *psPixelsFromMask() + - sdrs does not specify + +void psImageMaskRegion() +void psImageKeepRegion() + - sdrs does not specify + +void psImageMaskCircle() +void psImageKeepCircle() + - sdrs does not specify + +pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, psU8 maskVal, + int binFactor, psStats *stats, float clipSD); + - choice of polynomial coordinates + +stats *pmFringeStats (psArray *fringePoints, psImage *image, psMetadata *config); + - unspecified + +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row); + - sdrs says col,row are in *subimage* coords + +psArray *pmFindImagePeaks(const psImage *image, float threshold); + - does not specify coordinate system + - code uses subimage? + +bool pmSourceDefinePixels() +bool pmSourceRedefinePixels() + - sdrs does not specify very clearly + + + +--------- + +functions which are not affected +psFitsReadImage() +psFitsUpdateImage() +psFitsWriteImage() +psFitsInsertImage() +psImageCopy() -- maintains col0,row0 +psImageRebin() -- mask and image must match + +functions which currently use parent coordinates (CORRECT) +psImageSubset() -- check on range tests +psImageMaskRegion() +psImageKeepRegion() +psImageMaskCircle() +psImageKeepCircle() +pmSourceMoments -- (OK if input peaks are in parent coords) + +functions which currently use subimage coordinates (FIX) +psImageSet() +psImageGet() +psRegionForImage() *fixed* +psImageRow() +psImageCol() +psImageSlice() +psImageCut() +psImageCountPixelMask (); +pmModelEval() +pmFindImagePeaks() *fixed* +pmSourceDefinePixels() *fixed by psRegionForImage fix* +pmSourceRedefinePixels() *fixed by psRegionForImage fix* + +functions which are unclear (CHECK) +psImageTransform() +psImageFitPolynomial() +psImageEvalPolynomial() +psImagePixelInterpolate() +psPixelsToMask() +pmSubtractSky() +pmFringeStats() + + +***** one ambiguity: if we have a subimage, but we do not have the + parent image data (subimage does not carry dimensions of parent + array), then a 0 or negative upper-edge of a region cannot be + defined relative to the parent. we can specify it relative to + the subimage. so, for example, given a subimage [32,100:32,100] + of a much larger image, then psRegionForImage with 0,0,0,0 as + input would return a region with values 32,100:32,100... + + + + + +*** things I've modified + psRegionForImage expects parent (not subimage) coordinates Index: /tags/RHL_v1_2branch_0/psphot/doc/speed.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/speed.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/speed.txt (revision 10106) @@ -0,0 +1,129 @@ + +psphot -g +psModule -g +psLib -O2 +on alala: + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 0.983568 sec + System random seed value used seed = 25904F7C7454AE25 hex + background: 1229.280029 +/- 27.700180 + image stats: 0.063085 sec + System random seed value used seed = 411E5A0600A49AF3 hex + fit background model: 0.580152 sec + smooth: 1.842415 sec + threshold: 83.100540 DN + 4356 peaks: 0.571587 sec + 4230 moments: 3.965099 sec + SN range: 3.355829 - 1043.401001 + identified 107 blended objects (0.121737 sec) + selected candidate 300 PSF objects + fit flt: 3.682850 sec for 300 sources + fit psf: 1.704141 sec for 300 sources + fit stats: 0.000884 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717999 + select psf model: 5.402742 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.991345 (4230 objects) + built matrix: 1.472869 (27998 elements) + measure ensemble of PSFs: 4.334556 + fit PSF models: 12.205738 sec for 4230 objects + replace unfitted models: 0.021473 sec (4230 objects) + measure aperture residuals : 0.451886 sec + measure full-frame aperture residual: 0.454032 sec + aperture residual: -0.002842 +/- 0.002729 : 0.166072 bias + wrote output: 4.648821 sec + complete psphot run: 35.334080 sec + + +psphot -g +psModule -O2 +psLib -O2 +on alala: + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 1.001502 sec + System random seed value used seed = 5F43A63F532D4B3F hex + background: 1229.280029 +/- 27.700180 + image stats: 0.063169 sec + System random seed value used seed = F5F5FF99A6578F61 hex + fit background model: 0.530384 sec + smooth: 1.856414 sec + threshold: 83.100540 DN + 4365 peaks: 0.237292 sec + 4239 moments: 4.087611 sec + SN range: 3.356469 - 1043.400879 + identified 107 blended objects (0.123372 sec) + selected candidate 300 PSF objects + fit flt: 3.741444 sec for 300 sources + fit psf: 1.676478 sec for 300 sources + fit stats: 0.000883 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717985 + select psf model: 5.458156 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.791613 (4239 objects) + built matrix: 1.267909 (28101 elements) + measure ensemble of PSFs: 3.551364 + fit PSF models: 11.240197 sec for 4239 objects + replace unfitted models: 0.018051 sec (4239 objects) + measure aperture residuals : 0.411856 sec + measure full-frame aperture residual: 0.413843 sec + aperture residual: -0.002831 +/- 0.002717 : 0.164772 bias + wrote output: 4.355429 sec + complete psphot run: 33.126283 sec + +after reducing indirection in pmSourceMoments: + 4268 moments: 3.762452 sec + +psphot -O2 +psModule -O2 +psLib -O2 + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 0.484324 sec + System random seed value used seed = 980050732A392EC7 hex + background: 1229.280029 +/- 27.700180 + image stats: 0.058759 sec + System random seed value used seed = 44E1A6EE77B0E9CE hex + fit background model: 0.395125 sec + smooth: 1.861776 sec + threshold: 83.100540 DN + 4351 peaks: 0.234874 sec + 4264 moments: 3.719435 sec + SN range: 2.023345 - 1043.401123 + identified 124 blended objects (0.115230 sec) + selected candidate 300 PSF objects + fit flt: 3.525501 sec for 300 sources + fit psf: 1.595031 sec for 300 sources + fit stats: 0.000884 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717976 + select psf model: 5.160442 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.760805 (4264 objects) + built matrix: 0.908674 (23703 elements) + measure ensemble of PSFs: 3.115425 + fit PSF models: 10.162397 sec for 4264 objects + replace unfitted models: 0.017491 sec (4264 objects) + measure aperture residuals : 0.408242 sec + measure full-frame aperture residual: 0.410207 sec + aperture residual: -0.002908 +/- 0.002734 : 0.170692 bias + wrote output: 4.225560 sec + complete psphot run: 30.149281 sec + + + fit PSF models: 10.162397 sec for 4264 objects + select psf model: 5.160442 sec + wrote output: 4.225560 sec + 4264 moments: 3.719435 sec + measure ensemble of PSFs: 3.115425 sec + + other : 3.5 sec + + fixed per objects : 4 msec / object (17 of 30 sec) + + +** still can re-compile with -DPS_NO_TRACE +** ASSERTS cannot be compiled out at the moment! Index: /tags/RHL_v1_2branch_0/psphot/doc/versions.txt =================================================================== --- /tags/RHL_v1_2branch_0/psphot/doc/versions.txt (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/doc/versions.txt (revision 10106) @@ -0,0 +1,49 @@ + +2006.02.02 + I'm tagging psphot after I updated it to work with the readouts + using the structures being used for ppImage. In this release, I + have also moved to autoconf for psphot, and the psphotReadout + portion is being defined as a library against which the other + elements are linked. It should be possible to use this library with + ppImage, but I have not yet tested that. + + link psphot_dev_08 against psLib:eam_rel9_b2 and + psModule:eam_rel9_b2. + +2006.01.18 + the tag below (psphot_dev_07) was moved to a later release. It + should be used with psLib:eam_rel9_b1 from 2006/1/18 and + psModule:eam_rel9_p0. + +2005.12.23 + psphot_dev_07 should be compiled against psLib:eam_rel9_b1 and + psModule:eam_rel9_b1. This version now handles the ensemble PSF + fitting, detection of blended sources, fitting of 'blends' (groups + of PSFs with distinct be 'blended' peaks) and 'blobs' (unresolved + extended objects, with a test for double PSF vs extended source + model. + +2005.11.25 + + psphot_dev_05 should be compiled against psLib tag eam_rel8_b2 and + psModules tag eam_rel8_b1. this version removes all _EAM versions + of psLib and psModules code. It now uses the psModules version of + the object code. + + psphot_dev_04 should be compiled against psLib tag eam_rel8_b2 this + version removes all _EAM versions of psLib code, so that it is now + consistent with psLib.v8 (fixes are pushed into the eam_rel8_b2 + branch). + +2005.11.16 + + psphot_dev_03 should be compiled against psLib tag eam_rel8_b1, a + branch starting from psLib rel8_0. I needed to make some minor + fixes in psLib to work successfully with psLib rel8_0. + + psphot_dev_02 was getting segfaults from a bug added into psLib in + rel8_0 in p_psVectorClippedStats. I also fixed some errors arising + from changing types (argument to psMetadataParseConfig). + + + Index: /tags/RHL_v1_2branch_0/psphot/psphot.pc.in =================================================================== --- /tags/RHL_v1_2branch_0/psphot/psphot.pc.in (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/psphot.pc.in (revision 10106) @@ -0,0 +1,10 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ + +Name: libpsphot +Description: Pan-STARRS Photometry Library +Version: @VERSION@ +Libs: -L${libdir} -lpsphot @PSPHOT_LIBS@ +Cflags: -I${includedir} @PSPHOT_CFLAGS@ Index: /tags/RHL_v1_2branch_0/psphot/src/.cvsignore =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/.cvsignore (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/.cvsignore (revision 10106) @@ -0,0 +1,18 @@ +*.o +*.lo +.libs +.deps +Makefile +Makefile.in +libpsphot.la +psphot +psphotTest +psphot.loT +psphotErrorCodes.h +psphotErrorCodes.c +config.h +config.h.in +stamp-h1 + +polyfitTest +growthTest Index: /tags/RHL_v1_2branch_0/psphot/src/Makefile.am =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/Makefile.am (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/Makefile.am (revision 10106) @@ -0,0 +1,78 @@ + +lib_LTLIBRARIES = libpsphot.la +libpsphot_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) + +# bin_PROGRAMS = psphot polyfitTest growthTest +bin_PROGRAMS = psphot + +psphot_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) +psphot_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) +psphot_LDADD = libpsphot.la +psphot_SOURCES = psphot.c + +polyfitTest_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) +polyfitTest_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) +polyfitTest_SOURCES = polyfitTest.c + +growthTest_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) +growthTest_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) +growthTest_SOURCES = growthTest.c psphotGrowthCurve.c + +libpsphot_la_SOURCES = \ + psphotErrorCodes.c \ + psphotVersion.c \ + psphotModelGroupInit.c \ + psphotArguments.c \ + psphotParseCamera.c \ + psphotImageLoop.c \ + psphotMaskReadout.c \ + psphotReadout.c \ + psphotImageMedian.c \ + psphotFindPeaks.c \ + psphotSourceStats.c \ + psphotRoughClass.c \ + psphotBasicDeblend.c \ + psphotChoosePSF.c \ + psphotGuessModels.c \ + psphotEnsemblePSF.c \ + psphotBlendFit.c \ + psphotReplaceUnfit.c \ + psphotApResid.c \ + psphotMagnitudes.c \ + psphotSkyReplace.c \ + psphotEvalPSF.c \ + psphotEvalFLT.c \ + psphotSourceFits.c \ + psphotRadiusChecks.c \ + psphotSortBySN.c \ + psphotOutput.c \ + psphotGrowthCurve.c \ + psphotFakeSources.c \ + psphotModelTest.c \ + psphotFitSet.c \ + psphotWeightBias.c \ + psphotSourceFreePixels.c \ + psphotTestPSF.c \ + psphotCleanup.c + +include_HEADERS = \ + psphot.h \ + psphotErrorCodes.h + +clean-local: + -rm -f TAGS + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` +# +# Error codes. +# +BUILT_SOURCES = psphotErrorCodes.h psphotErrorCodes.c +CLEANFILES = psphotErrorCodes.h psphotErrorCodes.c + +psphotErrorCodes.h : psphotErrorCodes.dat psphotErrorCodes.h.in + $(ERRORCODES) --data=psphotErrorCodes.dat --outdir=. psphotErrorCodes.h + +psphotErrorCodes.c : psphotErrorCodes.dat psphotErrorCodes.c.in psphotErrorCodes.h + $(ERRORCODES) --data=psphotErrorCodes.dat --outdir=. psphotErrorCodes.c Index: /tags/RHL_v1_2branch_0/psphot/src/growthTest.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/growthTest.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/growthTest.c (revision 10106) @@ -0,0 +1,89 @@ +# include "psphot.h" + +int main (int argc, char **argv) { + + if (argc != 3) { + fprintf (stderr, "USAGE: growthTest\n"); + exit (2); + } + + // does the growth curve correctly fix aperture mags? + pmModelGroupInit (); + + // generate a simple readout + pmReadout *readout = pmReadoutAlloc (NULL); + readout->image = psImageAlloc (512, 512, PS_TYPE_F32); + readout->mask = psImageAlloc (512, 512, PS_TYPE_U8); + + // generate a simple psf + pmModelType type = pmModelSetType ("PS_MODEL_GAUSS"); + psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0); + pmPSF *psf = pmPSFAlloc (type, true, psfTrend); + psf->growth = pmGrowthCurveAlloc (2.0, 100.0, 15.0); + + psPolynomial2D *poly; + + poly = psf->params_NEW->data[PM_PAR_SXX]; + poly->coeff[0][0] = M_SQRT2*1.5; + + poly = psf->params_NEW->data[PM_PAR_SYY]; + poly->coeff[0][0] = M_SQRT2*1.5; + + poly = psf->params_NEW->data[PM_PAR_SXY]; + poly->coeff[0][0] = 0; + + // create template model + pmModel *modelRef = pmModelAlloc(psf->type); + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + + // generate the growth curve (element of psf) + psphotGrowthCurve (readout, psf, false); + + for (float radius = 2; radius < 8; radius += 1.0) { + for (float dx = -1; dx <= 1; dx += 0.1) { + for (float dy = -1; dy <= 1; dy += 0.1) { + + float fitMag, apMag, deltaMag; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_XPOS] = 200+dx; + modelRef->params->data.F32[PM_PAR_YPOS] = 200+dy; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + model->radiusFit = radius; + + // measure the fitMag for this model + pmSourcePhotometryModel (&fitMag, model); + + psRegion region = {180, 220, 180, 220}; + psImage *imageView = psImageSubset (readout->image, region); + psImage *maskView = psImageSubset (readout->mask, region); + + psImageInit (imageView, 0.0); + psImageInit (maskView, 0); + pmModelAdd (imageView, maskView, model, false, false); + + float dX = 0.5 - model->params->data.F32[PM_PAR_XPOS] + (int)model->params->data.F32[PM_PAR_XPOS]; + float dY = 0.5 - model->params->data.F32[PM_PAR_YPOS] + (int)model->params->data.F32[PM_PAR_YPOS]; + psImage *tmpImage = psImageShift (NULL, imageView, dX, dY, 0.0, PS_INTERPOLATE_BICUBE); + // psImage *tmpImage = tmp_psImageShift (imageView, dX, dY); + + psImageKeepCircle (maskView, 200+dx+dX, 200+dy+dY, radius, "OR", PM_MASK_MARK); + pmSourcePhotometryAper (&apMag, model, tmpImage, maskView); + psImageKeepCircle (maskView, 200+dx+dX, 200+dy+dY, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); + + deltaMag = pmGrowthCurveCorrect (psf->growth, radius); + + fprintf (stdout, "%f %f %f %f %f %f\n", radius, dx, dy, fitMag, apMag, deltaMag); + + psFree (tmpImage); + psFree (imageView); + psFree (maskView); + psFree (model); + } + } + } +} Index: /tags/RHL_v1_2branch_0/psphot/src/models/pmModel_TEST1.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/models/pmModel_TEST1.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/models/pmModel_TEST1.c (revision 10106) @@ -0,0 +1,236 @@ +/****************************************************************************** + * this file defines the TEST1 source shape model. Note that these model functions are loaded + * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own. The + * models use a psVector to represent the set of parameters, with the sequence used to specify + * the meaning of the parameter. The meaning of the parameters may thus vary depending on the + * specifics of the model. All models which are used a PSF representations share a few + * parameters, for which # define names are listed in pmModel.h: + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) * SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) * SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + *****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_TEST1 +# define PM_MODEL_FLUX pmModelFlux_TEST1 +# define PM_MODEL_GUESS pmModelGuess_TEST1 +# define PM_MODEL_LIMITS pmModelLimits_TEST1 +# define PM_MODEL_RADIUS pmModelRadius_TEST1 +# define PM_MODEL_FROM_PSF pmModelFromPSF_TEST1 +# define PM_MODEL_FIT_STATUS pmModelFitStatus_TEST1 + +// XXX consider changing to form PAR[SXY]*(1/PAR[SXX]^2 + 1/PAR[SYY]^2)*X*Y +// this would provide natural limits on PAR[SXY] of -0.25 : +0.25 + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +psF32 PM_MODEL_FUNC(psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + // XXX this is fitting sigma_x/sqrt(2), sigma_y/sqrt(2) + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*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[PM_PAR_I0]*r + PAR[PM_PAR_SKY]; + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + psF32 q = PAR[PM_PAR_I0]*r*r*t; + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; + } + return(f); +} + +// define the parameter limits +bool PM_MODEL_LIMITS (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[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; + beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; + + params_min[0][0].data.F32[PM_PAR_SKY] = -1000; + params_min[0][0].data.F32[PM_PAR_I0] = 0; + params_min[0][0].data.F32[PM_PAR_XPOS] = -100; + params_min[0][0].data.F32[PM_PAR_YPOS] = -100; + params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; + + params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; + params_max[0][0].data.F32[PM_PAR_I0] = 1e8; + params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; + params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; // XXX if we change the fitted term as above, this is naturally +/-1 + + return (TRUE); +} + +// make an initial guess for parameters +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) +{ + + pmMoments *moments = source->moments; + psF32 *PAR = model->params->data.F32; + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + PAR[PM_PAR_XPOS] = moments->x; + PAR[PM_PAR_YPOS] = moments->y; + PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy); + PAR[PM_PAR_SXY] = 0.0; // XXX we can get this right if we do the integral + + return(true); +} + +psF64 PM_MODEL_FLUX(const psVector *params) +{ + float norm, z; + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); + z += DZ; + f2 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psEllipseShape shape; + + if (flux <= 0) + return (1.0); + if (params->data.F32[PM_PAR_I0] <= 0) + return (1.0); + if (flux >= params->data.F32[PM_PAR_I0]) + return (1.0); + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + + // this estimates the radius assuming f(z) is roughly exp(-z) + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); + + if (isnan(radius)) psAbort ("psphot.model", "error in code: never return invalid radius"); + if (radius < 0) psAbort ("psphot.model", "error in code: never return invalid radius"); + + return (radius); +} + +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +{ + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + // we require these two parameters to exist + assert (psf->params_NEW->n > PM_PAR_YPOS); + assert (psf->params_NEW->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) { + out[i] = in[i]; + } else { + psPolynomial2D *poly = psf->params_NEW->data[i]; + out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here + out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + + return(true); +} + +// XXX double-check these definitions below +// this test is invalid if the parameters are derived +// from the PSF model +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); + dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + if (status) + return true; + return false; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /tags/RHL_v1_2branch_0/psphot/src/polyfitTest.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/polyfitTest.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/polyfitTest.c (revision 10106) @@ -0,0 +1,93 @@ +# include "psphot.h" + +int main (int argc, char **argv) { + + float X, Y, SX, SY, SXY; + + if (argc != 3) { + fprintf (stderr, "USAGE: polyfitTest (input) (output)\n"); + exit (2); + } + + FILE *f = fopen (argv[1], "r"); + if (f == NULL) exit (1); + + psVector *x = psVectorAlloc (100, PS_TYPE_F32); + psVector *y = psVectorAlloc (100, PS_TYPE_F32); + psVector *sx = psVectorAlloc (100, PS_TYPE_F32); + psVector *sy = psVectorAlloc (100, PS_TYPE_F32); + psVector *sxy = psVectorAlloc (100, PS_TYPE_F32); + + int i = 0; + while (fscanf (f, "%f %f %f %f %f %*f %*d", &X, &Y, &SX, &SY, &SXY) != EOF) { + x->data.F32[i] = X; + y->data.F32[i] = Y; + sx->data.F32[i] = SX; + sy->data.F32[i] = SY; + sxy->data.F32[i] = SXY; + + psVectorExtend (x, 100, 1); + psVectorExtend (y, 100, 1); + psVectorExtend (sx, 100, 1); + psVectorExtend (sy, 100, 1); + psVectorExtend (sxy, 100, 1); + + i++; + } + fclose (f); + + fprintf (stderr, "loaded %ld pts\n", x->n); + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + psVector *mask = psVectorAlloc (x->n, PS_TYPE_MASK); + mask->n = x->n; + psVectorInit (mask, 0); + + psPolynomial2D *polySX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); + polySX->mask[2][1] = 1; + polySX->mask[2][2] = 1; + polySX->mask[1][2] = 1; + + psVectorClipFitPolynomial2D (polySX, stats, mask, 0xff, sx, NULL, x, y); + fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); + fprintf (stderr, "SX: %g %g %g\n", polySX->coeff[0][0], polySX->coeff[1][0], polySX->coeff[0][1]); + fprintf (stderr, "SX: %g %g %g\n", polySX->coeff[0][2], polySX->coeff[1][1], polySX->coeff[2][0]); + + psPolynomial2D *polySY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); + polySY->mask[2][1] = 1; + polySY->mask[2][2] = 1; + polySY->mask[1][2] = 1; + + psVectorClipFitPolynomial2D (polySY, stats, mask, 0xff, sy, NULL, x, y); + fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); + fprintf (stderr, "SY: %g %g %g\n", polySY->coeff[0][0], polySY->coeff[1][0], polySY->coeff[0][1]); + fprintf (stderr, "SY: %g %g %g\n", polySY->coeff[0][2], polySY->coeff[1][1], polySY->coeff[2][0]); + + psPolynomial2D *polySXY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); + polySXY->mask[2][1] = 1; + polySXY->mask[2][2] = 1; + polySXY->mask[1][2] = 1; + + psVectorClipFitPolynomial2D (polySXY, stats, mask, 0xff, sxy, NULL, x, y); + fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); + fprintf (stderr, "SXY: %g %g %g\n", polySXY->coeff[0][0], polySXY->coeff[1][0], polySXY->coeff[0][1]); + fprintf (stderr, "SXY: %g %g %g\n", polySXY->coeff[0][2], polySXY->coeff[1][1], polySXY->coeff[2][0]); + + psVector *sxf = psPolynomial2DEvalVector (polySX, x, y); + psVector *syf = psPolynomial2DEvalVector (polySY, x, y); + psVector *sxyf = psPolynomial2DEvalVector (polySXY, x, y); + + f = fopen (argv[2], "w"); + for (i = 0; i < x->n; i++) { + fprintf (f, "%f %f %f %f %f %f %f %f\n", + x->data.F32[i], y->data.F32[i], + sx->data.F32[i], sy->data.F32[i], sxy->data.F32[i], + sxf->data.F32[i], syf->data.F32[i], sxyf->data.F32[i]); + } + fclose (f); + exit (0); +} + + // 0,0 1,0 2,0 + // 0,1 1,1 2,1 + // 0,2 1,2 2,2 Index: /tags/RHL_v1_2branch_0/psphot/src/psphot.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphot.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphot.c (revision 10106) @@ -0,0 +1,40 @@ +# include "psphot.h" + +static void usage (void) { + fprintf (stderr, "USAGE: ppImage [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); + exit (2); +} + +int main (int argc, char **argv) { + + psTimerStart ("complete"); + + psphotErrorRegister(); // register our error codes/messages + + // basic pslib options + psLogSetFormat ("M"); + + // load implementation-specific models + psphotModelGroupInit (); + + // load command-line arguments, options, and system config data + pmConfig *config = psphotArguments (argc, argv); + if (config == NULL) usage (); + + // load input data (config and images (signal, noise, mask) + if (!psphotParseCamera (config)) { + psErrorStackPrint(stderr, " "); + exit(1); + } + + // call psphot for each readout + if (!psphotImageLoop (config)) { + psErrorStackPrint(stderr, " "); + exit(1); + } + + psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); + + psphotCleanup (config); + exit (EXIT_SUCCESS); +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphot.h =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphot.h (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphot.h (revision 10106) @@ -0,0 +1,88 @@ +# ifdef HAVE_CONFIG_H +# include +# endif + +# include +# include // for strcasecmp +# include // for unlink +# include +# include + +# include "psphotErrorCodes.h" + +# define PSPHOT_RECIPE "PSPHOT" + +// top-level psphot functions +const char *psphotCVSName(void); + +// XXX test functions +bool psphotTestPSF (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool pmPSFtestModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask); + +pmConfig *psphotArguments (int argc, char **argv); +bool psphotParseCamera (pmConfig *config); +bool psphotImageLoop (pmConfig *config); + +bool psphotModelTest (pmReadout *readout, psMetadata *recipe); +bool psphotReadout (pmConfig *config, pmFPAview *view); +void psphotCleanup (pmConfig *config); + +psArray *psphotFakeSources (); + +// psphotReadout functions +bool psphotImageMedian (pmConfig *config, pmFPAview *view); +psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe); +psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *allpeaks); +bool psphotRoughClass (psArray *sources, psMetadata *recipe); +bool psphotBasicDeblend (psArray *sources, psMetadata *recipe); +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); +bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); +bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); +bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotReplaceUnfit (psArray *sources); +bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotSkyReplace (pmConfig *config, pmFPAview *view); +bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources); + +// basic support functions +void psphotModelGroupInit (void); +int psphotSortBySN (const void **a, const void **b); +int psphotSortByY (const void **a, const void **b); +bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore); +void psphotTestArguments (int *argc, char **argv); +bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe); +void psphotSourceFreePixels (psArray *sources); + +// functions to set the correct source pixels +bool psphotInitRadiusPSF (psMetadata *recipe, pmModelType type); +bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model); +bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR); +bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type); +bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model); + +// output functions +bool psphotDumpMoments (psMetadata *recipe, psArray *sources); +psMetadata *psphotDefineHeader (psMetadata *recipe); +bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); + +// PSF / DBL / EXT evaluation functions +bool psphotEvalPSF (pmSource *source, pmModel *model); +bool psphotEvalDBL (pmSource *source, pmModel *model); +bool psphotEvalEXT (pmSource *source, pmModel *model); + +// functions to support the source fitting process +bool psphotInitLimitsPSF (psMetadata *recipe); +bool psphotInitLimitsEXT (psMetadata *recipe); +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf); +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf); +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf); +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source); +psArray *psphotFitDBL (pmReadout *readout, pmSource *source); + +// functions to support simultaneous multi-source fitting +bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode); + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotApResid.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotApResid.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotApResid.c (revision 10106) @@ -0,0 +1,359 @@ +# include "psphot.h" +// XXXX this code fails if there are too few sources to measure the aperture residual +// the larger problem is that the rules for accepting more polynomial terms are weak. + +static pmPSFApTrendOptions DEFAULT_OPTION = PM_PSF_APTREND_SKYBIAS; + +// measure the aperture residual statistics +bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { + + int Nfail = 0; + int Nskip = 0; + int Npsf; + bool status; + pmModel *model; + pmSource *source; + + PS_ASSERT_PTR_NON_NULL(psf, false); + PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false); + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + + psTimerStart ("psphot"); + + // S/N limit to perform full non-linear fits + float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET"); + + // this is the smallest radius allowed: need to at least extend growth curve down to this... + float PSF_FIT_PAD = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING"); + bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); + bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); + + pmSourcePhotometryMode photMode = 0; + if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH; + if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP; + + // set limits on the aperture magnitudes + pmSourceMagnitudesInit (recipe); + + // measure the aperture loss as a function of radius for PSF + float REF_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_REF_RADIUS"); + psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS); + + pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH); + + psVector *mask = psVectorAllocEmpty (300, PS_TYPE_U8); + psVector *xPos = psVectorAllocEmpty (300, PS_TYPE_F64); + psVector *yPos = psVectorAllocEmpty (300, PS_TYPE_F64); + psVector *flux = psVectorAllocEmpty (300, PS_TYPE_F64); + psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64); + psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64); + psVector *dMag = psVectorAllocEmpty (300, PS_TYPE_F64); + Npsf = 0; + + // XXX test dump + FILE *f = fopen ("apresid.dat", "w"); + + // select all good PM_SOURCE_TYPE_STAR entries + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + model = source->modelPSF; + + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_FAIL) continue; + if (source->mode & PM_SOURCE_MODE_POOR) continue; + + // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures + // will fail if below S/N threshold or model is missing + if (!pmSourceMagnitudes (source, psf, photMode)) { + Nskip ++; + continue; + } + + apResid->data.F64[Npsf] = source->apMag - source->psfMag; + + xPos->data.F64[Npsf] = model->params->data.F32[PM_PAR_XPOS]; + yPos->data.F64[Npsf] = model->params->data.F32[PM_PAR_YPOS]; + + flux->data.F64[Npsf] = pow(10.0, -0.4*source->psfMag); + r2rflux->data.F64[Npsf] = PS_SQR(model->radiusFit) / flux->data.F64[Npsf]; + + mask->data.U8[Npsf] = 0; + + // XXX sanity clip? + // XXX need to see if all data were tossed? + // XXX need to subtract median? + // XXX need to put this in the config data... + // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue; + if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) { + Nfail ++; + continue; + } + + dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + + fprintf (f, "%d %f %f %f %f %f %f %f %f\n", i, + xPos->data.F64[Npsf], yPos->data.F64[Npsf], + flux->data.F64[Npsf], r2rflux->data.F64[Npsf], + source->apMag, source->psfMag, apResid->data.F64[Npsf],dMag->data.F64[Npsf]); + + psVectorExtend (mask, 100, 1); + psVectorExtend (xPos, 100, 1); + psVectorExtend (yPos, 100, 1); + psVectorExtend (flux, 100, 1); + psVectorExtend (r2rflux, 100, 1); + psVectorExtend (dMag, 100, 1); + psVectorExtend (apResid, 100, 1); + Npsf ++; + } + psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n", + psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail); + + fclose (f); + + // XXX choose a better value here? + if (Npsf < 5) { + psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up"); + return false; + } + + // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL + // APTREND options are used in the switch block below + pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION; + char *optionName = psMetadataLookupPtr (&status, recipe, "APTREND"); + if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName); + if (ApTrendOption == PM_PSF_APTREND_ERROR) { + psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName); + return false; + } + + // 3hi/1lo sigma clipping on the rflux vs metric fit + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + stats->min = 2.0; + stats->max = 3.0; + + // no correction + switch (ApTrendOption) { + case PM_PSF_APTREND_NONE: + // remove ApTrend fit from pmPSFtry + psf->ApTrend->coeff[0][0][0][0] = 0; + break; + case PM_PSF_APTREND_CONSTANT: + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction"); + return false; + } + break; + case PM_PSF_APTREND_SKYBIAS: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + return false; + } + break; + case PM_PSF_APTREND_SKYSAT: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat"); + return false; + } + break; + case PM_PSF_APTREND_XY_LIN: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN"); + return false; + } + break; + case PM_PSF_APTREND_XY_QUAD: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD"); + return false; + } + break; + case PM_PSF_APTREND_SKY_XY_LIN: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin"); + return false; + } + break; + case PM_PSF_APTREND_SKYSAT_XY_LIN: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin"); + return false; + } + break; + case PM_PSF_APTREND_ALL: + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend"); + return false; + } + // fit just SkyBias and clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias"); + return false; + } + // finally, fit x, y, SkyBias and clip out objects which are too far from the median + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting all"); + return false; + } + break; + default: + psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName); + return false; + } + + // construct the fitted values and the residuals + psVector *fit = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux); + psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit); + + // measure scatter for sources with dMag < 0.01 (S/N = 100) + int Nkeep = 0; + psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); + for (int i = 0; i < dMag->n; i++) { + if (dMag->data.F64[i] > 0.01) { + mask->data.U8[i] |= 0x02; + } + if (! mask->data.U8[i]) Nkeep ++; + } + residStats = psVectorStats (residStats, resid, NULL, mask, 0x03); + + // apply ApTrend results + psf->skyBias = psf->ApTrend->coeff[0][0][1][0]; + psf->skySat = psf->ApTrend->coeff[0][0][0][1]; + psf->ApResid = psf->ApTrend->coeff[0][0][0][0]; + psf->dApResid = residStats->sampleStdev; + psf->ApTrend->coeff[0][0][1][0] = 0; + psf->ApTrend->coeff[0][0][0][1] = 0; + psf->nApResid = Npsf; + + // save results for later output + psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skySat); + psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->growth->apLoss); + psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid); + + psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n", + psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf); + psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", + 1e3*psf->ApTrend->coeff[1][0][0][0], + 1e6*psf->ApTrend->coeff[2][0][0][0], + 1e6*psf->ApTrend->coeff[1][1][0][0], + 1e3*psf->ApTrend->coeff[0][1][0][0], + 1e6*psf->ApTrend->coeff[0][2][0][0]); + + psFree (mask); + psFree (xPos); + psFree (yPos); + psFree (flux); + psFree (r2rflux); + psFree (apResid); + psFree (dMag); + + psFree (fit); + psFree (resid); + psFree (stats); + psFree (residStats); + return true; +} + +/* + (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) + (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) + (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) +*/ + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotArguments.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotArguments.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotArguments.c (revision 10106) @@ -0,0 +1,136 @@ +# include "psphot.h" +# include + +pmConfig *psphotArguments(int argc, char **argv) { + + int N; + bool status; + + if (argc == 1) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc); + return NULL; + } + + // load config data from default locations + pmConfig *config = pmConfigRead(&argc, argv); + if (config == NULL) { + psError(PSPHOT_ERR_CONFIG, false, "Can't read site configuration"); + return NULL; + } + + // save the recipe options separately until we have loaded the recipe (psphotParseDetrend) + psMetadata *options = psMetadataAlloc (); + psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "PSPHOT.OPTIONS", PS_DATA_METADATA, "", options); + + // the following options override the PSPHOT recipe options + + // run the test model (requires X,Y coordinate) + if ((N = psArgumentGet (argc, argv, "-modeltest"))) { + psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT", 0, "", true); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); + + psArgumentRemove (N, &argc, argv); + psArgumentRemove (N, &argc, argv); + psArgumentRemove (N, &argc, argv); + + // specify the modeltest model + if ((N = psArgumentGet (argc, argv, "-model"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // specify the test fit mode + if ((N = psArgumentGet (argc, argv, "-fitmode"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-fitset"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + } + + // photcode : used in output to supplement header data (argument or recipe?) + if ((N = psArgumentGet (argc, argv, "-photcode"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // break : used from recipe throughout psphotReadout + if ((N = psArgumentGet (argc, argv, "-break"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // fitmode : used from recipe throughout psphotReadout + if ((N = psArgumentGet (argc, argv, "-fitmode"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF + if ((N = psArgumentGet (argc, argv, "-region"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // other arbitrary recipe values: -D key value (all added as string) + while ((N = psArgumentGet (argc, argv, "-D"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); + psArgumentRemove (N, &argc, argv); + psArgumentRemove (N, &argc, argv); + } + + // other arbitrary recipe values: -Df key value (all added as float) + while ((N = psArgumentGet (argc, argv, "-Df"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); + psArgumentRemove (N, &argc, argv); + psArgumentRemove (N, &argc, argv); + } + + // other arbitrary recipe values: -Di key value (all added as int) + while ((N = psArgumentGet (argc, argv, "-Di"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); + psArgumentRemove (N, &argc, argv); + psArgumentRemove (N, &argc, argv); + } + + // drop the local view on the options (saved on config->arguments) + psFree (options); + + // the input file is a required argument; if not found, we will exit + status = pmConfigFileSetsMD (config->arguments, config, "INPUT", "-file", "-list"); + if (!status) { + psError(PSPHOT_ERR_ARGUMENTS, false, "pmConfigFileSetsMD failed to parse arguments"); + return NULL; + } + + // chip selection is used to limit chips to be processed + if ((N = psArgumentGet (argc, argv, "-chip"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", psStringCopy(argv[N])); + psArgumentRemove (N, &argc, argv); + } + + if (argc != 2) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Expected to see one more argument; saw %d", argc - 1); + return NULL; + } + + // output position is fixed + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + psTrace("psphot", 1, "Done with psphotArguments...\n"); + return (config); +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotBasicDeblend.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotBasicDeblend.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotBasicDeblend.c (revision 10106) @@ -0,0 +1,127 @@ +# include "psphot.h" + +// 2006.02.07 : no leaks +bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) { + + int N; + bool status; + float threshold; + pmSource *source, *testSource; + + int Nblend = 0; + + psTimerStart ("psphot"); + + float FRACTION = psMetadataLookupF32 (&status, recipe, "DEBLEND_PEAK_FRACTION"); + if (!status) FRACTION = 0.25; + + float NSIGMA = psMetadataLookupF32 (&status, recipe, "DEBLEND_SKY_NSIGMA"); + if (!status) NSIGMA = 5.0; + + // 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 + + // examine sources in decreasing SN order + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; + + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + + // temporary array for overlapping objects we find + psArray *overlap = psArrayAllocEmpty (100); + + // 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); + } + + if (overlap->n == 0) { + psFree (overlap); + continue; + } + + // this source has overlapping neighbors, check for actual blends + // generate source contour (1/4 peak counts) + // set the threshold based on user inputs + + // threshold is fraction of the source peak flux + // image is background subtracted; source->moments->Sky should always be 0.0 + threshold = FRACTION * source->moments->Peak; + // threshold is no less than NSIGMA above the local median sigma? + threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky)); + + // generate a basic contour (set of x,y coordinates at-or-below flux level) + psArray *contour = pmSourceContour (source->pixels, source->peak->x, source->peak->y, threshold); + if (contour == NULL) { + psFree (overlap); + continue; + } + + // the source contour consists of two vectors, xv and yv. the contour is + // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]). both + // coordinate pairs have the same yv[] value, with xv[i] corresponding to + // the left boundary and xv[i+1] corresponding to the right boundary + + psVector *xv = contour->data[0]; + psVector *yv = contour->data[1]; + for (int k = 0; k < overlap->n; k++) { + testSource = overlap->data[k]; + if (testSource->peak->counts > source->peak->counts) continue; + for (int j = 0; j < xv->n; j+=2) { + if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; + if (xv->data.F32[j+0] > testSource->peak->x) break; + if (xv->data.F32[j+1] < testSource->peak->x) break; + + testSource->mode |= PM_SOURCE_MODE_BLEND; + + // add this to the list of source->blends + if (source->blends == NULL) { + source->blends = psArrayAllocEmpty (16); + } + psArrayAdd (source->blends, 16, testSource); + + Nblend ++; + j = xv->n; + } + } + psFree (overlap); + psFree (contour); + } + psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f sec)\n", Nblend, psTimerMark ("psphot")); + + psFree (SN); + psFree (index); + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotBlendFit.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotBlendFit.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotBlendFit.c (revision 10106) @@ -0,0 +1,90 @@ +# include "psphot.h" + +// XXX I don't like this name +bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { + + int Nfit = 0; + int Npsf = 0; + int Next = 0; + int Nfail = 0; + bool status; + + psTimerStart ("psphot"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // S/N limit to perform full non-linear fits + float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM"); + + psphotInitLimitsPSF (recipe); + psphotInitLimitsEXT (recipe); + psphotInitRadiusPSF (recipe, psf->type); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // skip DBL second sources (ie, added by psphotFitBlob) + if (source->mode & PM_SOURCE_MODE_PAIR) continue; + + // limit selection to some SN limit + if (source->moments == NULL) continue; + if (source->moments->SN < FIT_SN_LIM) continue; + + if (source->moments->x < AnalysisRegion.x0) continue; + if (source->moments->y < AnalysisRegion.y0) continue; + if (source->moments->x > AnalysisRegion.x1) continue; + if (source->moments->y > AnalysisRegion.y1) continue; + + // if model is NULL, we don't have a starting guess + if (source->modelPSF == NULL) continue; + + // skip sources which are insignificant flux? + if (source->modelPSF->params->data.F32[1] < 0.1) { + psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", + source->modelPSF->params->data.F32[1], + source->modelPSF->params->data.F32[2], + source->modelPSF->params->data.F32[3]); + continue; + } + + // replace object in image + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + } + Nfit ++; + + // try fitting PSFs, then try extended sources + if (psphotFitBlend (readout, source, psf)) { + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); + Npsf ++; + continue; + } + if (psphotFitBlob (readout, source, sources, psf)) { + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); + Next ++; + continue; + } + + psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y); + Nfail ++; + + // re-subtract PSF for object, leave local sky + pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->mode |= PM_SOURCE_MODE_TEMPSUB; + } + + psLogMsg ("psphot.psphotBlendFit", 3, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot"), Nfit, Npsf, Next, Nfail, sources->n - Nfit); + return (true); +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotChoosePSF.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotChoosePSF.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotChoosePSF.c (revision 10106) @@ -0,0 +1,277 @@ +# include "psphot.h" + +// try PSF models and select best option +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + char *modelName; + pmPSF *psf = NULL; + pmPSFtry *try = NULL; + psArray *stars = NULL; + + psTimerStart ("psphot"); + + // check if a PSF model is supplied by the user + psf = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.PSF"); + if (psf != NULL) return psf; + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // array to store candidate PSF stars + int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS"); + if (!status) { + NSTARS = PS_MIN (sources->n, 200); + psWarning("PSF_MAX_NSTARS is not set in the recipe --- defaulting to %d\n", NSTARS); + } + + // supply the measured sky variance for optional constant errors (non-poissonian) + float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV"); + if (!status) { + SKY_STDEV = 1.0; + psWarning("SKY_STDEV is not set --- defaulting to %f\n", SKY_STDEV); + } + // use poissonian errors or local-sky errors + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); + if (!status) { + POISSON_ERRORS = true; + psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n"); + } + pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_STDEV), POISSON_ERRORS); + + // how to model the PSF variations across the field + // XXX make a default value? or not? + psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK"); + psPolynomial2D *psfTrendMask; + if (!status || !md) { + psWarning("PSF.TREND.MASK is not set in the recipe --- defaulting to use zeroth order.\n"); + psfTrendMask = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0); + } else { + psfTrendMask = psPolynomial2DfromMetadata (md); + if (!psfTrendMask) { + psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); + return false; + } + } + + stars = psArrayAllocEmpty (sources->n); + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { + pmSource *source = sources->data[i]; + if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source); + } + psLogMsg ("psphot.pspsf", 4, "selected candidate %ld PSF objects\n", stars->n); + + if (stars->n == 0) { + psError(PSPHOT_ERR_PSF, true, "Failed to find any PSF candidates"); + return NULL; + } + + // get the fixed PSF fit radius + // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS + float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + if (!status) { + psWarning("PSF_FIT_RADIUS is not set in the recipe --- defaulting to 20.0\n"); + RADIUS = 20.0; + } + + // get the list pointers for the PSF_MODEL entries + psList *list = NULL; + psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + + if (mdi->type == PS_DATA_STRING) { + list = psListAlloc(NULL); + psListAdd (list, PS_LIST_HEAD, mdi); + } else { + if (mdi->type != PS_DATA_METADATA_MULTI) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + list = psMemIncrRefCounter(mdi->data.list); + } + + // set up an array to store the results + psArray *models = psArrayAlloc (list->n); + + // try each model option listed in config + psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); + for (int i = 0; i < models->n; i++) { + psMetadataItem *item = psListGetAndIncrement (iter); + modelName = item->data.V; + models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask); + } + + psFree (iter); + psFree (list); + psFree (stars); + psFree (psfTrendMask); + + // select the best of the models + // here we are using the clippedStdev on the metric as the indicator + int bestN = -1; + float bestM = 0.0; + for (int i = 0; i < models->n; i++) { + try = models->data[i]; + if (try == NULL) { + psError(PSPHOT_ERR_PSF, false, "PSF model %d is NULL", i); + continue; + } + float M = try->psf->dApResid; + if (bestN < 0 || M < bestM) { + bestM = M; + bestN = i; + } + } + + // use the best model: + if (bestN < 0) { + psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF"); + psFree (models); + return NULL; + } + + try = models->data[bestN]; + + // XXX test dump of psf star data and psf-subtracted image + if (psTraceGetLevel("psphot.psfstars") > 5) { + for (int i = 0; i < try->sources->n; i++) { + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = try->sources->data[i]; + float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + + // set the mask and subtract the PSF model + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK); + pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + } + + psphotSaveImage (NULL, readout->image, "psfstars.fits"); + pmSourcesWritePSFs (try->sources, "psfstars.dat"); + psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); + psMetadataConfigWrite (psfData, "psfmodel.dat"); + psFree (psfData); + psLogMsg ("psphot.choosePSF", 3, "wrote out psf-subtracted image, psf data, exiting\n"); + exit (0); + } + + // unset the PSFSTAR flag for stars not used for PSF model + for (int i = 0; i < try->sources->n; i++) { + pmSource *source = try->sources->data[i]; + if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) { + source->mode &= ~PM_SOURCE_MODE_PSFSTAR; + } + } + + // save only the best model; + // XXX we are not saving the fitted sources + // XXX do we want to keep them so we may optionally write them out? + psf = psMemIncrRefCounter(try->psf); + psFree (models); + + psMetadataAdd (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "number of stars used to make PSF", psf->nPSFstars); + + modelName = pmModelGetType (psf->type); + psLogMsg ("psphot.pspsf", 3, "select psf model: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + + return (psf); +} + +bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) { + + psEllipseShape shape; + psEllipseAxes axes; + + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + PS_ASSERT_PTR_NON_NULL(psf, false); + + psImage *image = readout->image; + PS_ASSERT_PTR_NON_NULL(image, false); + + pmModel *modelEXT = pmModelAlloc (psf->type); + PS_ASSERT_PTR_NON_NULL(modelEXT, false); + + // make a model with unit central intensity at the image center + modelEXT->params->data.F32[PM_PAR_SKY] = 0; + modelEXT->params->data.F32[PM_PAR_I0] = 1; + modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*image->numCols; + modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*image->numRows; + + // construct a PSF model at this coordinate + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); + PS_ASSERT_PTR_NON_NULL(modelPSF, false); + + // get the correct model-radius function + pmModelRadius modelRadiusFunc = pmModelRadius_GetFunction (psf->type); + + // get the model full-width at half-max + psF64 FWHM_X = 2*modelRadiusFunc (modelPSF->params, 0.5); + + // XXX make sure this is consistent with the re-definition of PM_PAR_SXX + shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; + shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; + shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape); + + psF64 FWHM_Y = FWHM_X * (axes.minor / axes.major); + + psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X); + psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); + psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF angle", axes.theta); + + psFree (modelEXT); + psFree (modelPSF); + + return true; +} + +// XXX fix this wrt re-definition of PM_PAR_SXX +bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) { + + // without the PSF model, we can only rely on the source->moments + // as a measure of the FWHM values + + int FWHM_N = 0; + double FWHM_X = 0.0; + double FWHM_Y = 0.0; + double FWHM_T = 0.0; + + psEllipseMoments moments; + psEllipseAxes axes; + + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) continue; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + // moments->Sx,Sy are roughly sigma_x,y + moments.x2 = PS_SQR(source->moments->Sx); + moments.y2 = PS_SQR(source->moments->Sy); + moments.xy = source->moments->Sxy; + + axes = psEllipseMomentsToAxes (moments); + + FWHM_X += axes.major * 2.35; + FWHM_Y += axes.minor * 2.35; + FWHM_T += axes.theta; + FWHM_N ++; + } + + FWHM_X /= FWHM_N; + FWHM_Y /= FWHM_N; + FWHM_T /= FWHM_N; + + psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X); + psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); + psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF angle", FWHM_T); + + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotCleanup.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotCleanup.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotCleanup.c (revision 10106) @@ -0,0 +1,16 @@ +# include "psphot.h" + +void psphotCleanup (pmConfig *config) { + + psFree (config); + + psTimerStop (); + psMemCheckCorruption (true); + pmModelGroupCleanup (); + psTimeFinalize (); + pmConceptsDone (); + pmConfigDone (); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); + // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot"); + return; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotEnsemblePSF.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotEnsemblePSF.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotEnsemblePSF.c (revision 10106) @@ -0,0 +1,306 @@ +# include "psphot.h" + +// fit flux (and optionally sky model) to all reasonable sources +// with the linear fitting process. sources must have an associated +// model with selected pixels, and the fit radius must be defined + +// given the set of sources, each of which points to the pixels in the +// science image, we construct a set of simulated sources with their own pixels. +// these are used to determine the simultaneous linear fit of fluxes. +// the analysis is performed wrt the simulated pixel values + +static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope, bool constant_weights); + +bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) { + + bool isPSF; + bool status; + float x; + float y; + float f; + float r; + + psTimerStart ("psphot"); + + // source analysis is done in spatial order + refSources = psArraySort (refSources, psphotSortByY); + + // storage arrays for fitSources and sequence index + psArray *fitSources = psArrayAllocEmpty (refSources->n); + psVector *index = psVectorAllocEmpty (refSources->n, PS_TYPE_U32); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionFromString (region); + AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion); + if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined"); + + const bool CONSTANT_PHOTOMETRIC_WEIGHTS = + psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + if (!status) { + psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); + } + + // construct source table with stand-alone fitSources + for (int i = 0; i < refSources->n; i++) { + pmSource *refSource = refSources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (refSource->type == PM_SOURCE_TYPE_DEFECT) continue; + if (refSource->type == PM_SOURCE_TYPE_SATURATED) continue; + if (final) { + if (refSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; + } else { + if (refSource->mode & PM_SOURCE_MODE_BLEND) continue; + } + + // which model to use? + pmModel *model = pmSourceGetModel (&isPSF, refSource); + if (model == NULL) continue; // model must be defined + + // save the original coords + x = model->params->data.F32[PM_PAR_XPOS]; + y = model->params->data.F32[PM_PAR_YPOS]; + + // is the source in the region of interest? + if (x < AnalysisRegion.x0) continue; + if (y < AnalysisRegion.y0) continue; + if (x > AnalysisRegion.x1) continue; + if (y > AnalysisRegion.y1) continue; + + pmSource *fitSource = pmSourceAlloc (); + + // make temporary copies of the image pixels and mask + // we need to have a copy which will not be modified by changes to its neighbor + fitSource->mask = psImageCopy (NULL, refSource->mask, PS_TYPE_U8); + fitSource->pixels = psImageCopy (NULL, refSource->pixels, PS_TYPE_F32); + fitSource->weight = psImageCopy (NULL, refSource->weight, PS_TYPE_F32); + + // set model to unit peak, zero sky (we assume sky is subtracted) + // XXX do this with the actual model or use keep the original guess? + model->params->data.F32[PM_PAR_SKY] = 0.0; + model->params->data.F32[PM_PAR_I0] = 1.0; + + // fill in the model pixel values + psImageInit (fitSource->pixels, 0.0); + psImageKeepCircle (fitSource->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); + pmModelAdd (fitSource->pixels, fitSource->mask, model, false, false); + + // save source in list + if (isPSF) { + fitSource->modelPSF = psMemIncrRefCounter (model); + } else { + fitSource->modelEXT = psMemIncrRefCounter (model); + } + index->data.U32[fitSources->n] = i; + psArrayAdd (fitSources, 100, fitSource); + psFree (fitSource); + } + psLogMsg (__func__, 4, "built fitSources: %f (%ld objects)\n", psTimerMark ("psphot"), refSources->n); + + // vectors to store stats for each object + psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32); + psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32); + + // create the border matrix (includes the sparse matrix) + // for just sky: 1 row; for x,y terms: 3 rows + psSparse *sparse = psSparseAlloc (fitSources->n, 100); + psSparseBorder *border = psSparseBorderAlloc (sparse, 1); + + // fill out the sparse matrix elements and border elements (B) + // Ri is the current refSource of interest (sci pixels) + // Fi is the current fitSource of interest (fit pixels) + for (int i = 0; i < fitSources->n; i++) { + int N = index->data.U32[i]; + pmSource *Ri = refSources->data[N]; + pmSource *Fi = fitSources->data[i]; + + // scale by diagonal element (auto-cross-product) + r = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS); + weight->data.F32[i] = r; + + // the formal error depends on the weighting scheme + if (CONSTANT_PHOTOMETRIC_WEIGHTS) { + float var = pmSourceCrossProduct (Fi, Fi, false); + errors->data.F32[i] = 1.0 / sqrt(var); + } else { + errors->data.F32[i] = 1.0 / sqrt(r); + } + + psSparseMatrixElement (sparse, i, i, 1.0); + + // find the image x model value + f = pmSourceCrossProduct (Ri, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseVectorElement (sparse, i, f / r); + + // add the per-source weights (border region) + float p = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); + // float px = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); + // float py = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 0, p); + // psSparseBorderElementB (border, i, 1, px); + // psSparseBorderElementB (border, i, 2, py); + + // loop over all other stars following this one + for (int j = i + 1; j < fitSources->n; j++) { + pmSource *Fj = fitSources->data[j]; + + // skip over disjoint source images, break after last possible overlap + if (Fi->pixels->row0 + Fi->pixels->numRows < Fj->pixels->row0) break; + if (Fj->pixels->row0 + Fj->pixels->numRows < Fi->pixels->row0) continue; + if (Fi->pixels->col0 + Fi->pixels->numCols < Fj->pixels->col0) continue; + if (Fj->pixels->col0 + Fj->pixels->numCols < Fi->pixels->col0) continue; + + // got an overlap; calculate cross-product and add to output array + f = pmSourceCrossProduct (Fi, Fj, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseMatrixElement (sparse, j, i, f / r); + } + } + psSparseResort (sparse); + psLogMsg (__func__, 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); + + // set the sky, sky_x, sky_y components of border matrix + // XXX ignore sky slope terms for now + SetBorderT (border, readout, refSources, fitSources, index, false, CONSTANT_PHOTOMETRIC_WEIGHTS); + + psSparseConstraint constraint; + constraint.paramMin = 0.0; + constraint.paramMax = 1e8; + constraint.paramDelta = 1e8; + + // solve for normalization terms (need include local sky?) + # if (1) + psVector *norm = NULL; + psVector *skyfit = NULL; + psSparseBorderSolve (&norm, &skyfit, constraint, border, 3); + fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); + # else + psVector *norm = psSparseSolve (NULL, constraint, sparse, 3); + psVector *skyfit = NULL; + # endif + + // adjust fitSources, set refSources and subtract + for (int i = 0; i < fitSources->n; i++) { + int N = index->data.U32[i]; + pmSource *Ri = refSources->data[N]; + pmModel *model = pmSourceGetModel (NULL, Ri); + + // assign linearly-fitted normalization + if (isnan(norm->data.F32[i])) { + psAbort ("psphot", "ensemble source is nan"); + } + model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; + // XXX note: is the value of 'errors' modified by the sky fit? + model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; + + // subtract object + pmModelSub (Ri->pixels, Ri->mask, model, false, false); + Ri->mode |= PM_SOURCE_MODE_SUBTRACTED; + if (!final) Ri->mode |= PM_SOURCE_MODE_TEMPSUB; + // XXX not sure about the use of TEMPSUB + } + + // measure chisq for each source + for (int i = 0; final && (i < fitSources->n); i++) { + int N = index->data.U32[i]; + pmSource *Ri = refSources->data[N]; + pmModel *model = pmSourceGetModel (NULL, Ri); + + x = model->params->data.F32[PM_PAR_XPOS]; + y = model->params->data.F32[PM_PAR_YPOS]; + + psImageKeepCircle (Ri->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); + pmSourceChisq (model, Ri->pixels, Ri->mask, Ri->weight); + psImageKeepCircle (Ri->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + } + + psFree (index); + psFree (sparse); + psFree (fitSources); + psFree (norm); + psFree (skyfit); + psFree (weight); + psFree (errors); + psFree (border); + + psLogMsg (__func__, 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot")); + return true; +} + +// calculate the weight terms for the sky fit component of the matrix +static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope, bool constant_weights) { + + // generate the image-wide weight terms + // turn on MARK for all image pixels + psRegion fullArray = psRegionSet (0, 0, 0, 0); + fullArray = psRegionForImage (readout->mask, fullArray); + psImageMaskRegion (readout->mask, fullArray, "OR", PM_MASK_MARK); + + // turn off MARK for all object pixels + // we must use the refSources here since their pixels point at the science image + for (int i = 0; i < fitSources->n; i++) { + int N = index->data.U32[i]; + pmSource *refSource = refSources->data[N]; + pmModel *model = pmSourceGetModel (NULL, refSource); + if (model == NULL) continue; + float x = model->params->data.F32[PM_PAR_XPOS]; + float y = model->params->data.F32[PM_PAR_YPOS]; + psImageMaskCircle (refSource->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + } + + // accumulate the image statistics from the masked regions + psF32 **image = readout->image->data.F32; + psF32 **weight = readout->weight->data.F32; + psU8 **mask = readout->mask->data.U8; + + double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; + w = x = y = x2 = xy = y2 = fo = fx = fy = 0; + + int col0 = readout->image->col0; + int row0 = readout->image->row0; + + for (int j = 0; j < readout->image->numRows; j++) { + for (int i = 0; i < readout->image->numCols; i++) { + if (mask[j][i]) continue; + if (constant_weights) { + wt = 1.0; + } else { + wt = weight[j][i]; + } + f = image[j][i]; + w += 1/wt; + fo += f/wt; + + if (!fitSlope) continue; + xc = i + col0; + yc = j + row0; + x += xc/wt; + y += yc/wt; + x2 += xc*xc/wt; + xy += xc*yc/wt; + y2 += yc*yc/wt; + fx += f*xc/wt; + fy += f*yc/wt; + } + } + + // turn off MARK for all image pixels + psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // set the Border T elements + psSparseBorderElementG (border, 0, fo); + psSparseBorderElementT (border, 0, 0, w); + if (fitSlope) { + psSparseBorderElementG (border, 0, fx); + psSparseBorderElementG (border, 0, fy); + psSparseBorderElementT (border, 1, 0, x); + psSparseBorderElementT (border, 2, 0, y); + psSparseBorderElementT (border, 0, 1, x); + psSparseBorderElementT (border, 1, 1, x2); + psSparseBorderElementT (border, 2, 1, xy); + psSparseBorderElementT (border, 0, 2, y); + psSparseBorderElementT (border, 1, 2, xy); + psSparseBorderElementT (border, 2, 2, y2); + } + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.c.in =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.c.in (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.c.in (revision 10106) @@ -0,0 +1,26 @@ +/* + * The line + { PSPHOT_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +#include "pslib.h" +#include "psphotErrorCodes.h" + +void psphotErrorRegister(void) +{ + static psErrorDescription errors[] = { + { PSPHOT_ERR_BASE, "First value we use; lower values belong to psLib" }, + { PSPHOT_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = PSPHOT_ERR_NERROR - PSPHOT_ERR_BASE; // number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + p_psMemSetPersistent(tmp, true); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.dat =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.dat (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.dat (revision 10106) @@ -0,0 +1,13 @@ +# +# This file is used to generate pmErrorClasses.h +# +BASE = 512 First value we use; lower values belong to psLib +UNKNOWN Unknown PM error code +NOT_IMPLEMENTED Desired feature is not yet implemented +ARGUMENTS Incorrect arguments +CONFIG Problem in configure files +FITS Problem in FITS I/O +FITS_WCS Error interpreting FITS WCS information +PHOTOM Problem in photometry +PSF Problem in PSF +APERTURE Problem with aperture photometry Index: /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.h.in =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.h.in (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotErrorCodes.h.in (revision 10106) @@ -0,0 +1,18 @@ +#if !defined(PSPHOT_ERROR_CODES_H) +#define PSPHOT_ERROR_CODES_H +/* + * The line + * PSPHOT_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + PSPHOT_ERR_BASE = 512, + PSPHOT_ERR_${ErrorCode}, + PSPHOT_ERR_NERROR +} psphotErrorCode; + +void psphotErrorRegister(void); + +#endif Index: /tags/RHL_v1_2branch_0/psphot/src/psphotEvalFLT.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotEvalFLT.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotEvalFLT.c (revision 10106) @@ -0,0 +1,48 @@ +# include "psphot.h" + +bool psphotEvalEXT (pmSource *source, pmModel *model) +{ + int keep; + + // do we actually have a valid EXT model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + psLogMsg ("psphot", 5, "EXT fail fit\n"); + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // unless we prove otherwise, this object is extended + source->type = PM_SOURCE_TYPE_EXTENDED; + + // the following source->mode information pertains to modelEXT: + source->mode |= PM_SOURCE_MODE_EXTMODEL; + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] <= 0) { + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + keep = pmModelFitStatus (model); + if (keep) return true; + + // poor-quality fit; only keep if nothing else works... + psLogMsg ("psphot", 5, "EXT poor fit\n"); + source->mode |= PM_SOURCE_MODE_POOR; + return false; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotEvalPSF.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotEvalPSF.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotEvalPSF.c (revision 10106) @@ -0,0 +1,204 @@ +# include "psphot.h" + +// given a pmSource which has been fitted using modelPSF, evaluate the +// resulting fit: did the fit succeed? is this object PSF-like? is this object +// extended? return status is TRUE for a valid PSF, FALSE otherwise. + +// identify objects consistent with PSF shape/magnitude distribution +// we expect dparams[PM_PAR_SXX],dparams[PM_PAR_SXX] to have a scatter of: +// sigma_x / (S/N) + +// any objects which is consistent with the PSF should have +// abs(dparams[PM_PAR_SXX]) < N * dsxLine(mag) & abs(dparams[PM_PAR_SYY]) < N * dsyLine(mag) +// this includes a minimum buffer (DS) for the brighter objects + +// EAM : 2006.10.20 : I have re-defined params[PM_PAR_SXX] : it is now sqrt(2)*sigma_x, not 1/sigma_x + +// saturated stars should fall outside (larger), but have peaks above SATURATION +// extended sources should be larger, cosmic rays smaller +// we also reject objects with S/N too low or ChiSquare to high + +// floor for DS value +// XXX EAM : add to configuration? +# define MIN_DS 0.01 + +static float SATURATION; +static float PSF_MIN_SN; +static float PSF_MAX_CHI; +static float PSF_SHAPE_NSIGMA; + +bool psphotInitLimitsPSF (psMetadata *recipe) { + + bool status; + + // we may set this differently here from the value used to mark likely saturated stars + SATURATION = psMetadataLookupF32 (&status, recipe, "SATURATION"); + PSF_MIN_SN = psMetadataLookupF32 (&status, recipe, "PSF_MIN_SN"); + PSF_MAX_CHI = psMetadataLookupF32 (&status, recipe, "PSF_MAX_CHI"); + PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_SHAPE_NSIGMA"); + + return true; +} + +// examine the model->status, fit parameters, etc and decide if the model succeeded +// set the source->type and source->mode appropriately +bool psphotEvalPSF (pmSource *source, pmModel *model) { + + int keep; + float dSX, dSY, SX, SY, SN; + float nSx, nSy, Chi; + + // do we actually have a valid PSF model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_TYPE_STAR; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_MODE_PSFMODEL; + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests + // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... + if (model->params->data.F32[PM_PAR_I0] >= SATURATION) { + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->mode |= PM_SOURCE_MODE_SATSTAR; + return true; + } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[PM_PAR_I0] <= 0) { + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_MODE_SATSTAR; + } + + SN = model->params->data.F32[PM_PAR_I0]/model->dparams->data.F32[1]; + SX = model->params->data.F32[PM_PAR_SXX]/M_SQRT2; + SY = model->params->data.F32[PM_PAR_SYY]/M_SQRT2; + dSX = model->dparams->data.F32[PM_PAR_SXX]; + dSY = model->dparams->data.F32[PM_PAR_SYY]; + Chi = model->chisqNorm / model->nDOF; + + // swing of sigma_x,y in sigmas + nSx = dSX / hypot (MIN_DS, (SX / SN)); + nSy = dSY / hypot (MIN_DS, (SY / SN)); + + // assign PM_SOURCE_MODE_GOODSTAR to bright objects within PSF region of dparams[] + keep = TRUE; + keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); + keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); + keep &= (SN > PSF_MIN_SN); + keep &= (Chi < PSF_MAX_CHI); + if (keep) return true; + + // this source is not a star, warn if it was a PSFSTAR + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality (%f, %f : %f %f %f %f)\n", + model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], nSx, nSy, SN, Chi); + } else { + psLogMsg ("psphot", 5, "fails PSF fit (%f, %f : %f %f %f %f)\n", + model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], nSx, nSy, SN, Chi); + } + + // object appears to be small, suspected defect + if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_TYPE_DEFECT; + return false; + } + + // object appears to be large, suspected extended source + if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + return false; + } + + // poor-quality fit; only keep if nothing else works... + source->mode |= PM_SOURCE_MODE_POOR; + return false; +} + +// examine the model->status, fit parameters, etc and decide if the model succeeded +// set the source->type and source->mode appropriately +bool psphotEvalDBL (pmSource *source, pmModel *model) { + + // do we actually have a valid PSF model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_TYPE_STAR; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_MODE_PSFMODEL; + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests + // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... + if (model->params->data.F32[1] >= SATURATION) { + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->mode |= PM_SOURCE_MODE_SATSTAR; + return true; + } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] <= 0) { + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_MODE_SATSTAR; + } + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotFakeSources.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotFakeSources.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotFakeSources.c (revision 10106) @@ -0,0 +1,30 @@ +# include "psphot.h" + +psArray *psphotFakeSources () { + + // psphotUpdateHeader (header, config); + + psArray *sources = psArrayAlloc (50); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = pmSourceAlloc (); + source->moments = pmMomentsAlloc (); + source->moments->x = 10; + source->moments->y = 10; + source->moments->Sx = 1; + source->moments->Sy = 1; + source->moments->Sxy = 0; + source->moments->Sum = 1000; + source->moments->Peak = 100; + source->moments->Sky = 10; + source->moments->nPixels = 10; + + source->peak = pmPeakAlloc (10, 10, 0, 0); + source->type = PM_SOURCE_TYPE_STAR; + + pmModelType modelType = pmModelSetType ("PS_MODEL_QGAUSS"); + source->modelPSF = pmSourceModelGuess (source, modelType); + sources->data[i] = source; + } + return sources; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotFindPeaks.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotFindPeaks.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotFindPeaks.c (revision 10106) @@ -0,0 +1,74 @@ +# include "psphot.h" + +// 2006.02.02 : no leaks +psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe) { + + bool status = false; + float NSIGMA; + float SIGMA; + float threshold; + float value; + + // smooth the image and weight map + psTimerStart ("psphot"); + + SIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); + NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + + // XXX we are somewhat doing the wrong thing here with the mask. + // it would be better to smooth the image using the mask information + psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); + psImageSmooth (smooth_im, SIGMA, NSIGMA); + psLogMsg ("psphot", 4, "smooth image: %f sec\n", psTimerMark ("psphot")); + + psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); + psImageSmooth (smooth_wt, SIGMA, NSIGMA); + psLogMsg ("psphot", 4, "smooth weight: %f sec\n", psTimerMark ("psphot")); + + psImage *mask = readout->mask; + + // optionally save example images under trace + if (psTraceGetLevel("psphot") > 5) { + psphotSaveImage (NULL, smooth_im, "imsmooth.fits"); + psphotSaveImage (NULL, smooth_wt, "wtsmooth.fits"); + } + + psTimerStart ("psphot"); + for (int j = 0; j < smooth_im->numRows; j++) { + for (int i = 0; i < smooth_im->numCols; i++) { + value = smooth_im->data.F32[j][i]; + smooth_im->data.F32[j][i] = 0.0; + if (mask->data.U8[j][i]) continue; + if (smooth_wt->data.F32[j][i] <= 0) continue; + smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; + } + } + psLogMsg ("psphot", 4, "built S/N map: %f sec\n", psTimerMark ("psphot")); + + psTimerStart ("psphot"); + + // set peak threshold + NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); + // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + threshold = PS_SQR(NSIGMA); + psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold); + + // find the peaks in the smoothed image + psArray *peaks = pmFindImagePeaks (smooth_im, threshold); + if (peaks == NULL) psAbort ("find peaks", "no peaks found"); + psFree (smooth_im); + psFree (smooth_wt); + + // optional dump of all peak data + char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE"); + if (status && (output != NULL) && (output[0])) { + pmPeaksWriteText (peaks, output); + } + psLogMsg ("psphot", 3, "%ld 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?) +// XXX : We need to gracefully handle no source detections Index: /tags/RHL_v1_2branch_0/psphot/src/psphotFitSet.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotFitSet.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotFitSet.c (revision 10106) @@ -0,0 +1,47 @@ +# include "psphot.h" + +// XXX this is not used in main psphot code +bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) { + + double x, y, Io; + + FILE *f = fopen (fitset, "r"); + if (f == NULL) return false; + + psArray *modelSet = psArrayAllocEmpty (16); + + while (fscanf (f, "%lf %lf %lf", &x, &y, &Io) == 3) { + pmModel *model = pmModelAlloc (oneModel->type); + + for (psS32 i = 0; i < model->params->n; i++) { + model->params->data.F32[i] = oneModel->params->data.F32[i]; + model->dparams->data.F32[i] = oneModel->dparams->data.F32[i]; + } + model->params->data.F32[1] = Io; + model->params->data.F32[2] = x; + model->params->data.F32[3] = y; + psArrayAdd (modelSet, 16, model); + } + + pmSourceFitSet (source, modelSet, mode); + + // write out positive object + psphotSaveImage (NULL, source->pixels, "object.fits"); + + // subtract object, leave local sky + for (int i = 0; i < modelSet->n; i++) { + pmModel *model = modelSet->data[i]; + pmModelSub (source->pixels, source->mask, model, false, false); + + fprintf (stderr, "output parameters (obj %d):\n", i); + for (int n = 0; n < model->params->n; n++) { + fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]); + } + } + + // write out + psphotSaveImage (NULL, source->pixels, "resid.fits"); + psphotSaveImage (NULL, source->mask, "mask.fits"); + return true; +} + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotGrowthCurve.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotGrowthCurve.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotGrowthCurve.c (revision 10106) @@ -0,0 +1,80 @@ +# include "psphot.h" + +// XXX we can probably move this into pmGrowthCurve.c +// XXX need to change the way we grab an image, or else use +// the 'center' option + +// XXX add a option to turn off the curve-of-growth (ie, make the apMag = fitMag everywhere); + +bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore) { + + // bool status; + float xc, yc, dx, dy; + float fitMag, apMag; + float radius; + + // create template model + pmModel *modelRef = pmModelAlloc(psf->type); + + // use the center of the center pixel of the image + xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; + yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; + dx = psf->growth->maxRadius + 1; + dy = psf->growth->maxRadius + 1; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + modelRef->params->data.F32[PM_PAR_XPOS] = xc; + modelRef->params->data.F32[PM_PAR_YPOS] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + + // measure the fitMag for this model + pmSourcePhotometryModel (&fitMag, model); + psf->growth->fitMag = fitMag; + + // generate working image for this source + psRegion region = {xc - dx, xc + dx, yc - dy, yc + dy}; + psImage *view = psImageSubset (readout->image, region); + psImage *image = psImageCopy (NULL, view, PS_TYPE_F32); + psImage *mask = psImageCopy (NULL, view, PS_TYPE_U8); + psImageInit (mask, 0); + + // loop over a range of source fluxes + // no need to interpolate since we have force the object center + // to 0.5, 0.5 above + for (int i = 0; i < psf->growth->radius->n; i++) { + + psImageInit (image, 0.0); + + radius = psf->growth->radius->data.F32[i]; + + psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK); + pmModelAdd (image, mask, model, false, false); + pmSourcePhotometryAper (&apMag, model, image, mask); + psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // the 'ignore' mode is for testing + if (ignore) { + psf->growth->apMag->data.F32[i] = fitMag; + } else { + psf->growth->apMag->data.F32[i] = apMag; + } + } + + psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius); + psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef; + + psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f\n", psf->growth->apLoss); + + psFree (view); + psFree (image); + psFree (mask); + psFree (model); + psFree (modelRef); + + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotGuessModels.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotGuessModels.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotGuessModels.c (revision 10106) @@ -0,0 +1,41 @@ +# include "psphot.h" + +// construct an initial PSF model for each object +bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { + + // setup the PSF fit radius details + psphotInitRadiusPSF (recipe, psf->type); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // use the source moments, etc to guess basic model parameters + pmModel *modelEXT = pmSourceModelGuess (source, psf->type); + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->x; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->y; + } else { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf; + } + + // set PSF parameters for this model (apply 2D shape model) + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); + psFree (modelEXT); + + // XXX need to define the guess flux? + // set the fit radius based on the object flux limit and the model + psphotCheckRadiusPSF (readout, source, modelPSF); + + // set the source PSF model + source->modelPSF = modelPSF; + } + psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n); + return true; +} + +// XXX do we always know which model is supposed to be used? Index: /tags/RHL_v1_2branch_0/psphot/src/psphotImageLoop.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotImageLoop.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotImageLoop.c (revision 10106) @@ -0,0 +1,98 @@ +# include "psphot.h" + +// XXX where do we load optional mask and weight input images? +bool psphotImageLoop (pmConfig *config) { + + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + if (!status) { + psError(PS_ERR_IO, false, "Can't find input data!"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for fpa in psphot\n"); + psFree(view); + return false; + } + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in psphot\n", view->chip); + psFree (view); + return false; + } + + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! cell->process || ! cell->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d, cell %d in psphot\n", view->chip, view->cell); + psFree (view); + return false; + } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d, cell %d, readout %d in psphot\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + + if (! readout->data_exists) { continue; } + + // run the actual photometry analysis + if (!psphotReadout (config, view)) { + psError(PS_ERR_UNKNOWN, false, "failed psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d, cell %d, readout %d in psphot\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d, cell %d in psphot\n", view->chip, view->cell); + psFree (view); + return false; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in psphot\n", view->chip); + psFree (view); + return false; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for fpa in psphot\n"); + psFree (view); + return false; + } + + psFree (view); + return true; +} + +// I/O files related to psphot: +// PSPHOT.INPUT : input image file(s) +// PSPHOT.RESID : residual image +// PSPHOT.OUTPUT : output object tables (object) + +// PSPHOT.BACKSUB : background subtracted image +// PSPHOT.BACKGND : background model (full-scale image?) +// PSPHOT.BACKMDL : background model (binned image?) +// PSPHOT.PSF : sample PSF images + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotImageMedian.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotImageMedian.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotImageMedian.c (revision 10106) @@ -0,0 +1,192 @@ +# include "psphot.h" + +// generate the median in NxN boxes, clipping heavily +// linear interpolation to generate full-scale model +bool psphotImageMedian (pmConfig *config, pmFPAview *view) +{ + bool status; + pmFPA *inFPA; + pmFPAfile *file; + psRegion region; + int MAX_SAMPLE_PIXELS; + + pmReadout *model = NULL; + pmReadout *readout = NULL; + pmReadout *background = NULL; + pmReadout *backSub = NULL; + + psTimerStart ("psphot"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); + + MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX"); + if (!status) { + MAX_SAMPLE_PIXELS = 1000; + } + unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED"); + if (!status) { + seed = 0; + } + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed); + + // subtract this amount extra from the sky + float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS"); + if (!status) { + SKY_BIAS = 0; + } + + // find the currently selected readout + file = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + inFPA = file->fpa; + readout = pmFPAviewThisReadout (view, inFPA); + + psImage *image = readout->image; + psImage *mask = readout->mask; + + // dimensions of input & output image + int Nx = image->numCols; + int Ny = image->numRows; + + // scaling factor + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + + // overhang : we will balance this evenly + int dx = (Nx % DX) / 2; + int dy = (Ny % DY) / 2; + + // dimensions of binned image + int nx = (Nx % DX) ? (int)(Nx / DX) + 1 : Nx / DX; + int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY; + + // select model pixels (from output background model file, or create internal file) + file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL"); + if (file == NULL) { + // we are not using PSPHOT.BACKMDL as an I/O file: define an internal version + model = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32); + } else { + // we are using PSPHOT.BACKMDL as an I/O file: select readout or create + model = pmFPAviewThisReadout (view, file->fpa); + if (model == NULL) { + // readout does not yet exist: create from input + pmFPAfileCopyStructureView (file->fpa, inFPA, DX, DY, view); + model = pmFPAviewThisReadout (view, file->fpa); + if (nx != model->image->numCols) psAbort ("psphot", "inconsistent size"); + if (ny != model->image->numRows) psAbort ("psphot", "inconsistent size"); + } + } + psF32 **modelData = model->image->data.F32; + assert(model->analysis != NULL); + psMetadataAdd(model->analysis, PS_LIST_TAIL, "XBIN", PS_DATA_S32 | PS_META_REPLACE, + "Background x-binsize", DX); + psMetadataAdd(model->analysis, PS_LIST_TAIL, "YBIN", PS_DATA_S32 | PS_META_REPLACE, + "Background x-binsize", DY); + psMetadataAdd(model->analysis, PS_LIST_TAIL, "x0", PS_DATA_S32 | PS_META_REPLACE, + "Background x-overhang", dx); + psMetadataAdd(model->analysis, PS_LIST_TAIL, "y0", PS_DATA_S32 | PS_META_REPLACE, + "Background y-overhang", dy); + + // measure clipped median for subimages + for (int iy = 0; iy < ny; iy++) { + for (int ix = 0; ix < nx; ix++) { + // sx, sy are in parent coords + int sx = ix*DX - dx + image->col0; + int sy = iy*DY - dy + image->row0; + region = psRegionSet (sx, sx + 2*DX, sy, sy + 2*DY); + region = psRegionForImage (image, region); + psImage *subset = psImageSubset (image, region); + if (!subset->numCols || !subset->numRows) { + psFree (subset); + continue; + } + + psImage *submask = psImageSubset (mask, region); + + // XXX the value of the upper and lower cuts probably should be studied... + psStats *stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS, rng); + modelData[iy][ix] = stats->robustMedian + SKY_BIAS; + + psFree (stats); + psFree (subset); + psFree (submask); + + // XXX psImageStats is still very inefficient and poorly coded... + // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + // stats = psImageStats (stats, subset, maskset, 0xff); + // backMdl->data.F32[iy][ix] = stats->clippedMean; + } + } + psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot")); + + // measure background stats and save for later output + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX); + stats = psImageStats (stats, model->image, NULL, 0); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_DATA_F32 | PS_META_REPLACE, "sky model mean", stats->sampleMean); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDEV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev", stats->sampleStdev); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MAX", PS_DATA_F32 | PS_META_REPLACE, "sky model maximum value", stats->max); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MIN", PS_DATA_F32 | PS_META_REPLACE, "sky model minimum value", stats->min); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_NX", PS_DATA_S32 | PS_META_REPLACE, "sky model size (x)", nx); + psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_NY", PS_DATA_S32 | PS_META_REPLACE, "sky model size (y)", ny); + psFree (stats); + + // select background pixels, from output background file, or create + file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND"); + if (file) { + // we are using PSPHOT.BACKMDL as an I/O file: select readout or create + background = pmFPAviewThisReadout (view, file->fpa); + if (background == NULL) { + // readout does not yet exist: create from input + pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view); + background = pmFPAviewThisReadout (view, file->fpa); + if (Nx != background->image->numCols) psAbort ("psphot", "inconsistent size"); + if (Ny != background->image->numRows) psAbort ("psphot", "inconsistent size"); + } + } else { + background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32); + } + psF32 **backData = background->image->data.F32; + + // linear interpolation to full-scale + psImageUnbin (background->image, model->image, DX, DY, dx, dy); + psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot")); + + // back-sub image pixels, from output background file (don't create if not requested) + file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKSUB"); + if (file) { + // we are using PSPHOT.BACKMDL as an I/O file: select readout or create + backSub = pmFPAviewThisReadout (view, file->fpa); + if (backSub == NULL) { + // readout does not yet exist: create from input + pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view); + backSub = pmFPAviewThisReadout (view, file->fpa); + } + } + +# ifdef TESTSAVE + psphotSaveImage (NULL, image, "image.fits"); + psphotSaveImage (NULL, background->image, "back.fits"); + psphotSaveImage (NULL, mask, "mask.fits"); + psphotSaveImage (NULL, model->image, "backmdl.fits"); +# endif + + // subtract the background model (save in backSub, if requested) + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + image->data.F32[j][i] -= backData[j][i]; + if (backSub) { + backSub->image->data.F32[j][i] = image->data.F32[j][i]; + } + } + } + +# ifdef TESTSAVE + psphotSaveImage (NULL, image, "backsub.fits"); +# endif + + psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot")); + psFree(rng); + + // the pmReadout selected in this function are all view on entries in config->files + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotMagnitudes.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotMagnitudes.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotMagnitudes.c (revision 10106) @@ -0,0 +1,30 @@ +# include "psphot.h" + +bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf) { + + bool status = false; + int Nap = 0; + + psTimerStart ("psphot"); + + pmSourceMagnitudesInit (recipe); + + bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); + bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); + + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT; + if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH; + if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + status = pmSourceMagnitudes (source, psf, photMode); + if (status) Nap ++; + } + + psLogMsg ("psphot.magnitudes", 4, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot"), sources->n, Nap); + return true; +} + +// XXX add in a measurement of the bright and faint completeness values +// XXX push these into the recipe Index: /tags/RHL_v1_2branch_0/psphot/src/psphotMaskReadout.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotMaskReadout.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotMaskReadout.c (revision 10106) @@ -0,0 +1,24 @@ +# include "psphot.h" + +bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe) { + + bool status; + + // mask the excluded outer pixels + // these coordinates refer to the parent image + // these bounds will saturate on the subimage + // negative upper bounds will subtract from the *subimage* + float XMIN = psMetadataLookupF32 (&status, recipe, "XMIN"); + float XMAX = psMetadataLookupF32 (&status, recipe, "XMAX"); + float YMIN = psMetadataLookupF32 (&status, recipe, "YMIN"); + float YMAX = psMetadataLookupF32 (&status, recipe, "YMAX"); + psRegion valid = psRegionSet (XMIN, XMAX, YMIN, YMAX); + + // restrict the supplied region above to the valid area on the image + psRegion keep = psRegionForImage (readout->image, valid); + + // psImageKeepRegion assumes the region refers to the parent coordinates + psImageKeepRegion (readout->mask, keep, "OR", PM_MASK_BAD); + + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotModelGroupInit.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotModelGroupInit.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotModelGroupInit.c (revision 10106) @@ -0,0 +1,23 @@ +# include "psphot.h" + +// Add locally-defined models here. As these mature, they can be moved to +// psModule/src/objects/models + +# include "models/pmModel_TEST1.c" + +static pmModelGroup userModels[] = { + {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1, pmModelFlux_TEST1, pmModelRadius_TEST1, pmModelLimits_TEST1, pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelFitStatus_TEST1}, +}; + +void psphotModelGroupInit (void) +{ + + // if pmModelGroupInit returns false, we have already init'ed + if (!pmModelGroupInit ()) return; + + int Nmodels = sizeof (userModels) / sizeof (pmModelGroup); + for (int i = 0; i < Nmodels; i++) { + pmModelGroupAdd (&userModels[i]); + } + return; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotModelTest.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotModelTest.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotModelTest.c (revision 10106) @@ -0,0 +1,191 @@ +# include "psphot.h" +static char DEFAULT_MODE[] = "EXT"; + +bool psphotModelTest (pmReadout *readout, psMetadata *recipe) { + + bool status; + int modelType; + unsigned int Nfail; + float obsMag, fitMag, value; + char name[64]; + pmPSF *psf = NULL; + pmSourceFitMode fitMode; + + psMetadataItem *item = NULL; + + // use poissonian errors or local-sky errors + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); + if (!status) POISSON_ERRORS = true; + pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS); + + // run model fitting tests on a single source + if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false; + + // what fitting mode to use? + char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE"); + if (!status) { + fitModeWord = DEFAULT_MODE; + } + fitMode = PM_SOURCE_FIT_EXT; + if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF; + + // in fitMode, psf sets the model type + if (fitMode == PM_SOURCE_FIT_PSF) { + char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE"); + if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied"); + psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, psfFile, FALSE); + psf = pmPSFfromMetadata (psfData); + modelType = psf->type; + } else { + // find the model: supplied by user or first in the PSF_MODEL list + char *modelName = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL"); + if (modelName == NULL) { + // get the list pointers for the PSF_MODEL entries + + psList *list = NULL; + psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotModelTest", "missing PSF_MODEL selection"); + if (mdi->type == PS_DATA_STRING) { + list = psListAlloc(NULL); + psListAdd (list, PS_LIST_HEAD, mdi); + } else { + if (mdi->type != PS_DATA_METADATA_MULTI) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + list = psMemIncrRefCounter(mdi->data.list); + } + + // take the first list element + 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, recipe, "TEST_FIT_INNER_RADIUS"); + if (!status) { + INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + } + + float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS"); + if (!status) { + OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + } + + float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS"); + if (!status) { + RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + } + + float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS"); + if (!status) { + mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + } + + // define the source of interest + float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); + float yObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y"); + + // construct the source structures + pmSource *source = pmSourceAlloc(); + source->peak = pmPeakAlloc (xObj, yObj, 0, 0); + pmSourceDefinePixels (source, readout, xObj, yObj, OUTER); + + // find the local sky + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER); + if (!status) psAbort ("fitsource", "pmSourceLocalSky error"); + + // get the source moments + status = pmSourceMoments (source, mRADIUS); + if (!status) psAbort ("fitsource", "pmSourceLocalSky error"); + source->peak->counts = source->moments->Peak; + + fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y); + + psEllipseMoments moments; + moments.x2 = source->moments->Sx; + moments.y2 = source->moments->Sy; + moments.xy = source->moments->Sxy; + psEllipseAxes axes = psEllipseMomentsToAxes (moments); + + fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + + // get the initial model parameter guess + pmModel *model = pmSourceModelGuess (source, modelType); + // if any parameters are defined, use those values + int nParams = pmModelParameterCount (modelType); + psF32 *params = model->params->data.F32; + for (int i = 0; i < nParams; i++) { + if (i == 2) { + params[i] = xObj; + continue; + } + if (i == 3) { + params[i] = yObj; + continue; + } + sprintf (name, "TEST_FIT_PAR%d", i); + value = psMetadataLookupF32 (&status, recipe, name); + if (status) { + params[i] = value; + } + } + + float area = params[4]*params[5]; + fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y); + + if (fitMode == PM_SOURCE_FIT_PSF) { + pmModel *modelPSF = pmModelFromPSF (model, psf); + psFree (model); + model = modelPSF; + params = model->params->data.F32; + } + + // list model input shape + psEllipseShape shape; + shape.sx = 1.4 / model->params->data.F32[4]; + shape.sy = 1.4 / model->params->data.F32[5]; + shape.sxy = model->params->data.F32[6]; + axes = psEllipseShapeToAxes (shape); + + fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + + + fprintf (stderr, "input parameters: \n"); + for (int i = 0; i < nParams; i++) { + fprintf (stderr, "%d : %f\n", i, params[i]); + } + + // define the pixels used for the fit + psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PM_MASK_MARK); + + char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET"); + if (status) { + status = psphotFitSet (source, model, fitset, fitMode); + exit (0); + } + + status = pmSourceFitModel (source, model, fitMode); + + // measure the source mags + pmSourcePhotometryModel (&fitMag, model); + pmSourcePhotometryAper (&obsMag, model, source->pixels, source->mask); + fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag); + + // write out positive object + psphotSaveImage (NULL, source->pixels, "object.fits"); + + // subtract object, leave local sky + pmModelSub (source->pixels, source->mask, model, false, false); + + fprintf (stderr, "output parameters: \n"); + for (int i = 0; i < nParams; i++) { + fprintf (stderr, "%d : %f\n", i, params[i]); + } + + // write out + psphotSaveImage (NULL, source->pixels, "resid.fits"); + psphotSaveImage (NULL, source->mask, "mask.fits"); + + exit (0); +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotOutput.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotOutput.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotOutput.c (revision 10106) @@ -0,0 +1,155 @@ +# include "psphot.h" + +int psphotSaveImage (psMetadata *header, psImage *image, char *filename) { + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + return (TRUE); +} + +bool psphotDumpMoments (psMetadata *recipe, psArray *sources) { + + bool status; + + // optional dump of all rough source data + char *output = psMetadataLookupStr (&status, recipe, "MOMENTS_OUTPUT_FILE"); + if (!status) return false; + if (output == NULL) return false; + if (output[0] == 0) return false; + + pmMomentsWriteText (sources, output); + return true; +} + +bool psphotDumpSource (pmSource *source, char *name) { + + FILE *f = fopen (name, "w"); + if (f == NULL) psAbort ("psphot", "can't open file"); + + for (int i = 0; i < source->pixels->numRows; i++) { + for (int j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + + fprintf (f, "%d %d %f %f %d\n", + (j + source->pixels->col0), + (i + source->pixels->row0), + source->pixels->data.F32[i][j], + 1.0 / source->weight->data.F32[i][j], + source->mask->data.U8[i][j]); + } + } + fclose (f); + return true; +} + +// these values are saved in an output header stub - they are added to either the +// PHU header (CMP) or the MEF table header (CMF) +psMetadata *psphotDefineHeader (psMetadata *recipe) { + + psMetadata *header = psMetadataAlloc (); + + // write necessary information to output header + psMetadataItemTransfer (header, recipe, "ZERO_PT"); + psMetadataItemTransfer (header, recipe, "PHOTCODE"); + + psMetadataItemTransfer (header, recipe, "APMIFIT"); + psMetadataItemTransfer (header, recipe, "DAPMIFIT"); + psMetadataItemTransfer (header, recipe, "NAPMIFIT"); + psMetadataItemTransfer (header, recipe, "SKYBIAS"); + psMetadataItemTransfer (header, recipe, "SKYSAT"); + + // PSF model parameters (shape values for image center) + psMetadataItemTransfer (header, recipe, "NPSFSTAR"); + psMetadataItemTransfer (header, recipe, "APLOSS"); + psMetadataItemTransfer (header, recipe, "FWHM_X"); + psMetadataItemTransfer (header, recipe, "FWHM_Y"); + psMetadataItemTransfer (header, recipe, "ANGLE"); + + // XXX these need to be defined from elsewhere + psMetadataAdd (header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + + // sky background model statistics + psMetadataItemTransfer (header, recipe, "SKY_MEAN"); + psMetadataItemTransfer (header, recipe, "SKY_STDEV"); + psMetadataItemTransfer (header, recipe, "SKY_MIN"); + psMetadataItemTransfer (header, recipe, "SKY_MAX"); + psMetadataItemTransfer (header, recipe, "SKY_NX"); + psMetadataItemTransfer (header, recipe, "SKY_NY"); + + return header; +} + +# if (0) +// output functions: we have several fixed modes (sx, obj, cmp) +void psphotOutput (pmReadout *readout, psMetadata *arguments) { + + bool status; + char *outputFile = NULL; + + psMetadata *header = pmReadoutGetHeader (readout); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + // sample PSF images?? + if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); + if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); + + if (psfFile != NULL) { + psMetadata *psfData = pmPSFtoMD (NULL, psf); + psMetadataConfigWrite (psfData, psfFile); + psFree (psfData); + } +} + +psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) { + + psRegion region = {x - dx, x + dx, y - dy, y + dy}; + psImage *view = psImageSubset (image, region); + psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32); + psImageInit (sample, 0); + modelEXT->params->data.F32[2] = x; + modelEXT->params->data.F32[3] = y; + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); + pmModelAdd (sample, NULL, modelPSF, false, false); + psFree (modelPSF); + return (sample); +} + +bool psphotSamplePSFs (pmPSF *psf, psImage *image, char *output) { + + // make sample PSFs for 4 corners and the center + psImage *sample; + + // optional dump of all rough source data + if (output[0] == 0) return false; + + pmModel *modelEXT = pmModelAlloc (psf->type); + modelEXT->params->data.F32[0] = 0; + modelEXT->params->data.F32[1] = 1; + + psFits *fits = psFitsOpen (output, "w"); + + // the centers are in parent coordinates; they do not need to correspond to valid pixels... + sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25); + psFitsWriteImage (fits, NULL, sample, 0); + + psFitsClose (fits); + return (TRUE); +} +# endif Index: /tags/RHL_v1_2branch_0/psphot/src/psphotParseCamera.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotParseCamera.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotParseCamera.c (revision 10106) @@ -0,0 +1,90 @@ +# include "psphot.h" + +// define the needed / desired I/O files +bool psphotParseCamera (pmConfig *config) { + + bool status = false; + + // the input image defines the camera + pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PSPHOT.INPUT", "INPUT"); + if (!status) { + psError(PSPHOT_ERR_CONFIG, false, "Failed to build FPA from PSPHOT.INPUT"); + return status; + } + + // select recipe options supplied on command line + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); + psMetadata *options = psMetadataLookupPtr (&status, config->arguments, "PSPHOT.OPTIONS"); + psMetadataCopy (recipe, options); + + // set default recipe values here + psMetadataAddStr (recipe, PS_LIST_TAIL, "FITMODE", PS_META_NO_REPLACE, "default fitting mode", "NONE"); + psMetadataAddStr (recipe, PS_LIST_TAIL, "PHOTCODE", PS_META_NO_REPLACE, "default photcode", "NONE"); + psMetadataAddStr (recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_NO_REPLACE, "default break point", "NONE"); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "ZERO_PT", PS_META_NO_REPLACE, "default zero point", 25.00); + psMetadataAddS32 (recipe, PS_LIST_TAIL, "BACKGROUND.XBIN", PS_META_NO_REPLACE, "default binning", 64); + psMetadataAddS32 (recipe, PS_LIST_TAIL, "BACKGROUND.YBIN", PS_META_NO_REPLACE, "default binning", 64); + + // these calls bind the I/O handle to the specified fpa + bool saveOutput = psMetadataLookupBool (&status, recipe, "SAVE.OUTPUT"); + if (!status || saveOutput) { + if (!pmFPAfileDefineOutput (config, input->fpa, "PSPHOT.OUTPUT")) { + psError(PS_ERR_IO, false, "Cannot find a rule for PSPHOT.OUTPUT"); + return false; + } + } + + // optionally save the residual image + if (psMetadataLookupBool(NULL, recipe, "SAVE.RESID")) { + pmFPAfileDefineOutput (config, input->fpa, "PSPHOT.RESID"); + } + + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + + // these calls construct a new fpa for the I/O handle + + // optionally save the background model (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { + pmFPAfileDefineFromFPA (config, input->fpa, DX, DY, "PSPHOT.BACKMDL"); + } + // optionally save the full background image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKGND")) { + pmFPAfileDefineFromFPA (config, input->fpa, 1, 1, "PSPHOT.BACKGND"); + } + // optionally save the background-subtracted image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKSUB")) { + pmFPAfileDefineFromFPA (config, input->fpa, 1, 1, "PSPHOT.BACKSUB"); + } + // optionally save the PSF Model + if (psMetadataLookupBool(NULL, recipe, "SAVE.PSF")) { + pmFPAfileDefineOutput (config, input->fpa, "PSPHOT.PSF.SAVE"); + } + // optionally load the PSF Model + if (psMetadataLookupBool(NULL, recipe, "LOAD.PSF")) { + pmFPAfileDefineInput (config, input->fpa, "PSPHOT.PSF.LOAD"); + } + + // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.PSF_SAMPLE"); + + // Chip selection: turn on only the chips specified (pass status to suppress missing-key log msg) + char *chipLine = psMetadataLookupStr(&status, config->arguments, "CHIP_SELECTIONS"); + psArray *chips = psStringSplitArray (chipLine, ",", false); + if (chips->n > 0) { + pmFPASelectChip (input->fpa, -1, true); // deselect all chips + for (int i = 0; i < chips->n; i++) { + int chipNum = atoi(chips->data[i]); + if (! pmFPASelectChip(input->fpa, chipNum, false)) { + psError(PS_ERR_IO, false, "Chip number %d doesn't exist in camera.\n", chipNum); + return false; + } + } + } + psFree (chips); + + psTrace("psphot", 1, "Done with psphotParseCamera...\n"); + + psErrorClear(); // some metadata lookup may have failed + + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotRadiusChecks.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotRadiusChecks.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotRadiusChecks.c (revision 10106) @@ -0,0 +1,84 @@ +# include "psphot.h" +# define RADIUS_TYPE int + +static float PSF_FIT_NSIGMA; +static float PSF_FIT_PADDING; +static pmModelRadius modelRadiusPSF; + +bool psphotInitRadiusPSF (psMetadata *recipe, pmModelType type) { + + bool status; + + PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_FIT_NSIGMA"); + PSF_FIT_PADDING = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING"); + + // this function specifies the radius at this the model hits the given flux + modelRadiusPSF = pmModelRadius_GetFunction (type); + return true; +} + +// call this function whenever you (re)-define the PSF model +bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model) { + + pmMoments *moments = source->moments; + if (moments == NULL) return false; + + // set the fit radius based on the object flux limit and the model + model->radiusFit = (RADIUS_TYPE) (modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + PSF_FIT_PADDING); + if (isnan(model->radiusFit)) psAbort ("apply_psf_model", "error in radius"); + + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + model->radiusFit *= 2; + } + + bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit); + return status; +} + +bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR) { + + pmMoments *moments = source->moments; + if (moments == NULL) return false; + + // set the fit radius based on the object flux limit and the model + model->radiusFit = (RADIUS_TYPE) (modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + dR + PSF_FIT_PADDING); + if (isnan(model->radiusFit)) psAbort ("apply_psf_model", "error in radius"); + + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + model->radiusFit *= 2; + } + + bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit); + return status; +} + +static float EXT_FIT_NSIGMA; +static float EXT_FIT_PADDING; +static pmModelRadius modelRadiusEXT; + +bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) { + + bool status; + + EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA"); + EXT_FIT_PADDING = psMetadataLookupF32 (&status, recipe, "EXT_FIT_PADDING"); + + // this function specifies the radius at this the model hits the given flux + modelRadiusEXT = pmModelRadius_GetFunction (type); + return true; +} + +// call this function whenever you (re)-define the EXT model +bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model) { + + pmMoments *moments = source->moments; + if (moments == NULL) return false; + + // set the fit radius based on the object flux limit and the model + model->radiusFit = (RADIUS_TYPE) (modelRadiusEXT (model->params, EXT_FIT_NSIGMA*moments->dSky) + EXT_FIT_PADDING); + if (isnan(model->radiusFit)) psAbort (__func__, "error in radius"); + + // redefine the pixels if needed + bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit); + return status; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotReadout.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotReadout.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotReadout.c (revision 10106) @@ -0,0 +1,159 @@ +# include "psphot.h" + +bool psphotReadout (pmConfig *config, pmFPAview *view) { + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT"); + + // find the currently selected readout + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); + + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + + // XXX does this need to invoke I/O? + // XXX move this input the psphotImageLoop level? + if (!pmReadoutGenerateMaskWeight (readout)) { + return false; + } + + // I have a valid mask, now mask in the analysis region of interest + psphotMaskReadout (readout, recipe); + + // run a single-model test if desired + psphotModelTest (readout, recipe); + + if (psTraceGetLevel("psphot") > 5) { + psphotSaveImage (NULL, readout->image, "image.fits"); + psphotSaveImage (NULL, readout->mask, "mask.fits"); + psphotSaveImage (NULL, readout->weight, "weight.fits"); + } + + // generate a background model (median, smoothed image) + psphotImageMedian (config, view); + + // find the peaks in the image + psArray *peaks = psphotFindPeaks (readout, recipe); + + // construct sources and measure basic stats + // limit moments analysis by S/N? + psArray *sources = psphotSourceStats (readout, recipe, peaks); + psFree (peaks); + + if (!strcasecmp (breakPt, "PEAKS")) { + if (!psphotReadoutCleanup (config, readout, recipe, NULL, sources)) return false; + return true; + } + + // mark blended peaks PS_SOURCE_BLEND + psphotBasicDeblend (sources, recipe); + + // classify sources based on moments, brightness + // faint sources not classified? + if (!psphotRoughClass (sources, recipe)) { + psFree (sources); + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + return false; + } + if (!strcasecmp (breakPt, "MOMENTS")) { + if (!psphotReadoutCleanup (config, readout, recipe, NULL, sources)) return false; + return true; + } + + // psphotTestPSF (readout, sources, recipe); + // exit (1); + + // use bright stellar objects to measure PSF + pmPSF *psf = psphotChoosePSF (readout, sources, recipe); + if (psf == NULL) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to construct a psf model"); + return false; + } + + if (!strcasecmp (breakPt, "PSFMODEL")) { + if (!psphotReadoutCleanup (config, readout, recipe, psf, sources)) return false; + return true; + } + + psphotGuessModels (readout, sources, recipe, psf); + + // linear PSF fit to peaks + psphotEnsemblePSF (readout, sources, recipe, psf, FALSE); + if (!strcasecmp (breakPt, "ENSEMBLE")) { + goto finish; + } + + // non-linear PSF and EXT fit to brighter sources + psphotBlendFit (readout, sources, recipe, psf); + + // XXX this was an attempt to measure bias due to the weighting. + // psphotWeightBias (readout, sources, recipe, psf); + + // replace fitted sources + psphotReplaceUnfit (sources); + + // XXX find remaining peaks after first source subtraction pass + // psphotFindPeaks (); + + // linear PSF fit to remaining peaks + // psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); + +finish: + + // measure aperture photometry corrections + if (!psphotApResid (readout, sources, recipe, psf)) { + psTrace ("psphot", 4, "failure on psphotApResid"); + psError(PSPHOT_ERR_PHOTOM, false, "Measure aperture photometry corrections"); + return false; + } + + // calculate source magnitudes + psphotMagnitudes (sources, recipe, psf); + + // replace background in residual image + psphotSkyReplace (config, view); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (sources); + // psphotSaveImage (NULL, readout->image, "resid.fits"); + + // create the exported-metadata and free local data + if (!psphotReadoutCleanup (config, readout, recipe, psf, sources)) return false; + return true; +} + +bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources) { + + // use the psf-model to measure FWHM stats + if (psf) { + if (!psphotPSFstats (readout, recipe, psf)) { + psError(PS_ERR_IO, true, "Failed to measure PSF shape parameters"); + return false; + } + } else if (sources) { + if (!psphotMomentsStats (readout, recipe, sources)) { + psError(PS_ERR_IO, true, "Failed to measure Moment shape parameters"); + return false; + } + } + + // create an output header with stats results + psMetadata *header = psphotDefineHeader (recipe); + + // save the results of the analysis + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER", PS_DATA_METADATA, "header stats", header); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); + + // remove internal pmFPAfiles, if created + pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + + psFree (psf); + psFree (header); + psFree (sources); + + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotReplaceUnfit.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotReplaceUnfit.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotReplaceUnfit.c (revision 10106) @@ -0,0 +1,24 @@ +# include "psphot.h" + +bool psphotReplaceUnfit (psArray *sources) { + + pmSource *source; + + psTimerStart ("psphot"); + + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_MODE_TEMPSUB)) continue; + if (source->modelPSF == NULL) continue; + + psTrace ("psphot", 3, "replacing object at %f,%f\n", + source->modelPSF->params->data.F32[2], source->modelPSF->params->data.F32[3]); + + pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; + source->mode &= ~PM_SOURCE_MODE_TEMPSUB; + } + psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot"), sources->n); + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotRoughClass.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotRoughClass.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotRoughClass.c (revision 10106) @@ -0,0 +1,26 @@ +# include "psphot.h" + +// 2006.02.02 : no leaks +bool psphotRoughClass (psArray *sources, psMetadata *recipe) { + + pmPSFClump psfClump; + + psfClump = pmSourcePSFClump (sources, recipe); + psLogMsg ("psphot", 3, "psf clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); + psLogMsg ("psphot", 3, "psf clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); + + // did we find a valid psfClump? + if (!psfClump.X || !psfClump.Y) { + psError(PSPHOT_ERR_PSF, true, "Failed to find a valid PSF clump"); + return false; + } + + // group into STAR, COSMIC, EXTENDED, SATURATED, etc. + pmSourceRoughClass (sources, recipe, psfClump); + + // optional printout of source moments only + psphotDumpMoments (recipe, sources); + + return true; +} + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotSkyReplace.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotSkyReplace.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotSkyReplace.c (revision 10106) @@ -0,0 +1,30 @@ +# include "psphot.h" + +// XXX make this an option? +// in order to successfully replace the sky, we must define a corresponding file... +bool psphotSkyReplace (pmConfig *config, pmFPAview *view) { + + // find the currently selected readout + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); + if (readout == NULL) psAbort ("psphot", "input not defined"); + + // select background pixels, from output background file, or create + pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND"); + if (background == NULL) psAbort ("psphot", "background not defined"); + + // select the corresponding images + psF32 **image = readout->image->data.F32; + psU8 **mask = readout->mask->data.U8; + psF32 **back = background->image->data.F32; + + // replace the background model + for (int j = 0; j < readout->image->numRows; j++) { + for (int i = 0; i < readout->image->numCols; i++) { + if (!mask[j][i]) { + image[j][i] += back[j][i]; + } + } + } + return true; +} + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotSortBySN.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotSortBySN.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotSortBySN.c (revision 10106) @@ -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/RHL_v1_2branch_0/psphot/src/psphotSourceFits.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotSourceFits.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotSourceFits.c (revision 10106) @@ -0,0 +1,350 @@ +# include "psphot.h" + +// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful + +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { + + float x, y, dR; + + // if this source is not a possible blend, just fit as PSF + if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { + bool status = psphotFitPSF (readout, source, psf); + return status; + } + psTrace ("psphot", 5, "trying blend...\n"); + + // save the PSF model from the Ensemble fit + pmModel *PSF = pmModelCopy (source->modelPSF); + if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in blend fit primary"); + + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; + + psArray *modelSet = psArrayAllocEmpty (source->blends->n + 1); + psArrayAdd (modelSet, 16, PSF); + + psArray *sourceSet = psArrayAllocEmpty (source->blends->n + 1); + psArrayAdd (sourceSet, 16, source); + + // we need to include all blends in the fit (unless primary is saturated?) + dR = 0; + for (int i = 0; i < source->blends->n; i++) { + pmSource *blend = source->blends->data[i]; + + // find the blend which is furthest from source + dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y)); + + // create the model and guess parameters for this blend + pmModel *model = pmModelAlloc (PSF->type); + for (int j = 0; j < model->params->n; j++) { + model->params->data.F32[j] = PSF->params->data.F32[j]; + model->dparams->data.F32[j] = PSF->dparams->data.F32[j]; + } + + // XXX assume local sky is 0.0? + model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky; + if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit"); + model->params->data.F32[2] = blend->peak->x; + model->params->data.F32[3] = blend->peak->y; + + // add this blend to the list + psArrayAdd (modelSet, 16, model); + psArrayAdd (sourceSet, 16, blend); + + // free to avoid double counting model + psFree (model); + } + + // extend source radius as needed + psphotCheckRadiusPSFBlend (readout, source, PSF, dR); + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); + pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // correct model chisq for flux trend + double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]); + PSF->chisqNorm = PSF->chisq / chiTrend; + + // evaluate the blend objects, subtract if good, free otherwise + for (int i = 1; i < modelSet->n; i++) { + pmSource *blend = sourceSet->data[i]; + pmModel *model = modelSet->data[i]; + + // correct model chisq for flux trend + chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]); + model->chisqNorm = model->chisq / chiTrend; + + // if this one failed, skip it + if (!psphotEvalPSF (blend, model)) { + psTrace ("psphot", 5, "failed on blend %d of %ld\n", i, modelSet->n); + continue; + } + + // otherwise, supply the resulting model to the corresponding blend + psFree(blend->modelPSF); + blend->modelPSF = psMemIncrRefCounter (model); + psTrace ("psphot", 5, "fitted blend as PSF\n"); + pmModelSub (source->pixels, source->mask, model, false, false); + blend->mode |= PM_SOURCE_MODE_SUBTRACTED; + blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; + } + psFree (modelSet); + psFree (sourceSet); + + // evaluate the primary object + if (!psphotEvalPSF (source, PSF)) { + psTrace ("psphot", 5, "failed on blend 0 of %ld\n", modelSet->n); + psFree (PSF); + return false; + } + + psTrace ("psphot", 5, "fitted primary as PSF\n"); + pmModelSub (source->pixels, source->mask, PSF, false, false); + psFree (source->modelPSF); + source->modelPSF = PSF; + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->mode &= ~PM_SOURCE_MODE_TEMPSUB; + return true; +} + +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { + + float x, y; + double chiTrend; + + // save the PSF model from the Ensemble fit + pmModel *PSF = pmModelCopy (source->modelPSF); + if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in psf fit"); + + // extend source radius as needed + psphotCheckRadiusPSF (readout, source, PSF); + + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); + pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF); + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // correct model chisq for flux trend + chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]); + PSF->chisqNorm = PSF->chisq / chiTrend; + + // does the PSF model succeed? + if (!psphotEvalPSF (source, PSF)) { + psFree (PSF); + return false; + } + + psTrace ("psphot", 5, "fitted as PSF\n"); + pmModelSub (source->pixels, source->mask, PSF, false, false); + + // free old model, save new model + psFree (source->modelPSF); + source->modelPSF = PSF; + + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->mode &= ~PM_SOURCE_MODE_TEMPSUB; + + return true; +} + +static float EXT_MIN_SN; +static float EXT_MOMENTS_RAD; +static pmModelType modelTypeEXT; + +bool psphotInitLimitsEXT (psMetadata *recipe) { + + bool status; + + // extended source model parameters + EXT_MIN_SN = psMetadataLookupF32 (&status, recipe, "EXT_MIN_SN"); + EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, recipe, "EXT_MOMENTS_RADIUS"); + + // extended source model descriptions + char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL"); + modelTypeEXT = pmModelSetType (modelNameEXT); + psphotInitRadiusEXT (recipe, modelTypeEXT); + + return true; +} + +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { + + bool okEXT, okDBL; + float chiEXT, chiDBL; + double chiTrend; + pmModel *ONE = NULL; + + // skip the source if we don't think it is extended + if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false; + if (source->type == PM_SOURCE_TYPE_DEFECT) return false; + if (source->type == PM_SOURCE_TYPE_SATURATED) return false; + if (source->moments->SN < EXT_MIN_SN) return false; + + // recalculate the source moments using the larger extended-source moments radius + if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false; + + psTrace ("psphot", 5, "trying blob...\n"); + + // this temporary source is used as a place-holder by the psphotEval functions below + pmSource *tmpSrc = pmSourceAlloc (); + + pmModel *EXT = psphotFitEXT (readout, source); + okEXT = psphotEvalEXT (tmpSrc, EXT); + chiEXT = EXT->chisq / EXT->nDOF; + + psArray *DBL = psphotFitDBL (readout, source); + okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); + okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); + // XXX should I keep / save the flags set in the eval functions? + + // correct first model chisqs for flux trend + ONE = DBL->data[0]; + chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); + ONE->chisqNorm = ONE->chisq / chiTrend; + + // save chisq for double-star/galaxy comparison + chiDBL = ONE->chisqNorm; + + // correct second model chisqs for flux trend + ONE = DBL->data[1]; + chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); + ONE->chisqNorm = ONE->chisq / chiTrend; + + psFree (tmpSrc); + + if (okEXT && okDBL) { + psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL); + // XXX EAM : a bogus bias: need to examine this better + if (3*chiEXT > chiDBL) goto keepDBL; + goto keepEXT; + } + + if (okEXT && !okDBL) goto keepEXT; + if (!okEXT && okDBL) goto keepDBL; + + // both models failed; reject them both + psFree (EXT); + psFree (DBL); + return false; + +keepEXT: + // sub EXT + psFree (DBL); + pmModelSub (source->pixels, source->mask, EXT, false, false); + psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[2], EXT->params->data.F32[3]); + + // save new model + source->modelEXT = EXT; + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->mode &= ~PM_SOURCE_MODE_TEMPSUB; + return true; + +keepDBL: + // sub DLB + psFree (EXT); + pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[0], false, false); + pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[1], false, false); + psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]); + + // drop old model, save new second model... + psFree (source->modelPSF); + source->modelPSF = psMemIncrRefCounter (DBL->data[0]); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->mode |= PM_SOURCE_MODE_PAIR; + source->mode &= ~PM_SOURCE_MODE_TEMPSUB; + + // copy most data from the primary source (modelEXT, blends stay NULL) + pmSource *newSrc = pmSourceAlloc (); + newSrc->peak = psMemIncrRefCounter (source->peak); + newSrc->pixels = psMemIncrRefCounter (source->pixels); + newSrc->weight = psMemIncrRefCounter (source->weight); + newSrc->mask = psMemIncrRefCounter (source->mask); + newSrc->moments = psMemIncrRefCounter (source->moments); + newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); + newSrc->type = source->type; + newSrc->mode = source->mode; + psArrayAdd (sources, 100, newSrc); + psFree (newSrc); + psFree (DBL); + return true; +} + +// fit a double PSF source to an extended blob +psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { + + float x, y, dx, dy; + pmModel *DBL; + pmModel *PSF; + psEllipseAxes axes; + psEllipseMoments moments; + psArray *modelSet; + + // make a guess at the position of the two sources + moments.x2 = source->moments->Sx; + moments.y2 = source->moments->Sy; + moments.xy = source->moments->Sxy; + axes = psEllipseMomentsToAxes (moments); + + dx = 2 * cos (axes.theta); + dy = 2 * sin (axes.theta); + + // save the PSF model from the Ensemble fit + PSF = source->modelPSF; + psphotCheckRadiusPSFBlend (readout, source, PSF, 8.0); + if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in dbl fit"); + + modelSet = psArrayAlloc (2); + + DBL = pmModelCopy (PSF); + DBL->params->data.F32[1] *= 0.5; + DBL->params->data.F32[2] = source->moments->x + dx; + DBL->params->data.F32[3] = source->moments->y + dy; + modelSet->data[0] = DBL; + + DBL = pmModelCopy (PSF); + DBL->params->data.F32[1] *= 0.5; + DBL->params->data.F32[2] = source->moments->x - dx; + DBL->params->data.F32[3] = source->moments->y - dy; + modelSet->data[1] = DBL; + + x = source->moments->x; + y = source->moments->y; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); + pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + + return (modelSet); +} + +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { + + float x, y; + + // use the source moments, etc to guess basic model parameters + pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); + // if (isnan(EXT->params->data.F32[1])) psAbort ("psphot", "nan in ext fit"); + + psphotCheckRadiusEXT (readout, source, EXT); + + x = EXT->params->data.F32[2]; + y = EXT->params->data.F32[3]; + + if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) { + psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy); + } + + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, EXT->radiusFit, "OR", PM_MASK_MARK); + pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT); + psImageKeepCircle (source->mask, x, y, EXT->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + + return (EXT); +} + Index: /tags/RHL_v1_2branch_0/psphot/src/psphotSourceFreePixels.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotSourceFreePixels.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotSourceFreePixels.c (revision 10106) @@ -0,0 +1,12 @@ +# include "psphot.h" + +void psphotSourceFreePixels (psArray *sources) { + + if (!sources) return; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmSourceFreePixels (source); + } + return; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotSourceStats.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotSourceStats.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotSourceStats.c (revision 10106) @@ -0,0 +1,81 @@ +# include "psphot.h" + +psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks) +{ + bool status = false; + psArray *sources = NULL; + float BIG_RADIUS; + + psTimerStart ("psphot"); + + // determine properties (sky, moments) of initial sources + float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); + + sources = psArrayAllocEmpty (peaks->n); + + for (int i = 0; i < peaks->n; i++) { + + // create a new source, add peak + pmSource *source = pmSourceAlloc(); + source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]); + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + if (!strcasecmp (breakPt, "PEAKS")) { + psArrayAdd (sources, 100, source); + psFree (source); + continue; + } + + // XXX skip faint sources? + + // 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 the local sky variance (needed if noise is not sqrt(signal)) + // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken + status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER); + if (!status) { + psFree (source); + continue; + } + + // measure basic source moments + status = pmSourceMoments (source, RADIUS); + if (status) { + // add to the source array + psArrayAdd (sources, 100, source); + psFree (source); + continue; + } + + // if no valid pixels, or massive swing, likely saturated source, + // try a much larger box + BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS); + if (status) { + // add to the source array + psArrayAdd (sources, 100, source); + psFree (source); + continue; + } + + psFree (source); + continue; + } + + psLogMsg ("psphot", 3, "%ld 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/RHL_v1_2branch_0/psphot/src/psphotTest.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotTest.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotTest.c (revision 10106) @@ -0,0 +1,83 @@ +# include "psphot.h" + +void psExit (int status, char *process, char *format, ...) { + + va_list ap; + + va_start (ap, format); + fprintf (stderr, "exiting %s\n", process); + vfprintf (stderr, format, ap); + va_end (ap); + + exit (status); +} + +int main (int argc, char **argv) { + + psRegion region = {0,0,0,0}; // a region representing the entire array + psphotTestArguments (&argc, argv); + + psFits *file = psFitsOpen (argv[1], "r"); + psMetadata *header = psFitsReadHeader (NULL, file); + psImage *image = psFitsReadImage (NULL, file, region, 0); + psFitsClose (file); + + psImageJpegColormap (argv[5]); + + // psImage *fimage = psImageCopy (NULL, image, PS_TYPE_F32); + + int binning = atof(argv[6]); + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN); + psImage *fimage = psImageRebin (NULL, image, NULL, 0, binning, stats); + + float min = atof(argv[3]); + float max = atof(argv[4]); + + psImageJpeg (fimage, argv[2], min, max); + + psFree (header); + psFree (image); + exit (0); +} + + +# if (0) + + psMetadata *row; + psArray *table; + + psMetadataItem *mdi; + + psMetadataConfigWrite (header, argv[2]); + + // attempt to write image with NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create a test image + // psImage *tmpimage = psImageAlloc (10, 10, PS_DATA_F32); + + // create a test table + table = psArrayAllocEmpty (10); + + for (int i = 0; i < 10; i++) { + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "ROW", PS_DATA_S32, "", i); + psMetadataAdd (row, PS_LIST_TAIL, "FROW", PS_TYPE_F32, "", 0.1*i); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "test line"); + + table->data[i] = row; + } + table->n = 10; + + psMetadata *theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); + + psFits *fits = psFitsOpen (argv[3], "w"); + // psFitsWriteImage (fits, header, tmpimage, 0); + psFitsWriteHeader (header, fits); + psFitsWriteTable (fits, theader, table); + +# endif Index: /tags/RHL_v1_2branch_0/psphot/src/psphotTestArguments.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotTestArguments.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotTestArguments.c (revision 10106) @@ -0,0 +1,19 @@ +# include "psphot.h" +static int usage (); + +void psphotTestArguments (int *argc, char **argv) { + + // basic pslib options + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + if (*argc != 7) usage (); + + return; +} + +static int usage () { + + fprintf (stderr, "USAGE: psphotTest (input.fits) (output.jpg) (zero) (scale) (colormap) (rebin)\n"); + exit (2); +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotTestPSF.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotTestPSF.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotTestPSF.c (revision 10106) @@ -0,0 +1,131 @@ +# include "psphot.h" + +bool psphotTestPSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + char *modelName; + pmPSF *psf = NULL; + psArray *stars = NULL; + + psTimerStart ("psphot"); + + psphotSaveImage (NULL, readout->image, "image.fits"); + + // check if a PSF model is supplied by the user + psf = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.PSF"); + if (psf != NULL) return psf; + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // array to store candidate PSF stars + int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS"); + if (!status) { + NSTARS = PS_MIN (sources->n, 200); + psWarning("PSF_MAX_NSTARS is not set in the recipe --- defaulting to %d\n", NSTARS); + } + + // use poissonian errors or local-sky errors + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); + if (!status) { + POISSON_ERRORS = true; + psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n"); + } + pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS); + + // how to model the PSF variations across the field + // XXX make a default value? or not? + psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK"); + psPolynomial2D *psfTrendMask; + if (!status || !md) { + psWarning("PSF.TREND.MASK is not set in the recipe --- defaulting to use zeroth order.\n"); + psfTrendMask = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0); + } else { + psfTrendMask = psPolynomial2DfromMetadata (md); + if (!psfTrendMask) { + psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); + return false; + } + } + + stars = psArrayAllocEmpty (sources->n); + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { + pmSource *source = sources->data[i]; + // if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source); + psArrayAdd (stars, 200, source); + } + psLogMsg ("psphot.pspsf", 4, "selected candidate %ld PSF objects\n", stars->n); + + if (stars->n == 0) { + psError(PSPHOT_ERR_PSF, true, "Failed to find any PSF candidates"); + return NULL; + } + + // get the fixed PSF fit radius + // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS + float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + if (!status) { + psWarning("PSF_FIT_RADIUS is not set in the recipe --- defaulting to 20.0\n"); + RADIUS = 20.0; + } + + // for this test, require a single model + psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + if (mdi->type != PS_DATA_STRING) psAbort ("psphotChoosePSF", "choose a single PSF_MODEL"); + modelName = mdi->data.V; + + pmPSFtestModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask); + + psphotSaveImage (NULL, readout->image, "resid.fits"); + psphotSaveImage (NULL, readout->mask, "mask.fits"); + psphotSaveImage (NULL, readout->weight, "weight.fits"); + + return true; +} + +bool pmPSFtestModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask) +{ + bool status; + float x; + float y; + + pmModelType type = pmModelSetType (modelName); + pmPSF *psf = pmPSFAlloc (type, poissonErrors, psfTrendMask); + if (psf == NULL) psAbort ("psphotTestPSF", "unknown model"); + + FILE *f = fopen ("params.dat", "w"); + + // stage 1: fit an independent model (freeModel) to all sources + psTimerStart ("fit"); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmModel *model = pmSourceModelGuess (source, psf->type); + x = source->peak->x; + y = source->peak->y; + + // set temporary object mask and fit object + // fit model as EXT, not PSF + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK); + status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // write fitted parameters to file + fprintf (f, "%f ", model->params->data.F32[PM_PAR_XPOS]); + fprintf (f, "%f ", model->params->data.F32[PM_PAR_YPOS]); + + fprintf (f, "%f ", model->params->data.F32[PM_PAR_SXX]); + fprintf (f, "%f ", model->params->data.F32[PM_PAR_SYY]); + fprintf (f, "%f ", model->params->data.F32[PM_PAR_SXY]); + + fprintf (f, "%f %d\n", model->chisq, model->nIter); + + // subtract model flux + pmModelSub (source->pixels, source->mask, model, false, false); + } + fclose (f); + psLogMsg ("psphot.psftest", 4, "fit ext: %f sec for %ld sources\n", psTimerMark ("fit"), sources->n); + return true; +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotVersion.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotVersion.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotVersion.c (revision 10106) @@ -0,0 +1,6 @@ +# include "psphot.h" + +const char *psphotCVSName(void) { + static char *name = "$Name: not supported by cvs2svn $"; + return(name); +} Index: /tags/RHL_v1_2branch_0/psphot/src/psphotWeightBias.c =================================================================== --- /tags/RHL_v1_2branch_0/psphot/src/psphotWeightBias.c (revision 10106) +++ /tags/RHL_v1_2branch_0/psphot/src/psphotWeightBias.c (revision 10106) @@ -0,0 +1,89 @@ +# include "psphot.h" + +// select objects fitted with PSF model +// re-fit all of them with the non-poisson errors +// only allow the normalization to vary +// XXX eventually, we should be able to do this with linear fitting... +bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { + + int Nfit = 0; + bool status; + psF32 *PARp, *PARc; + psF32 x, y; + + psTimerStart ("psphot"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // set fitting method to use non-poisson errors (local sky error) + float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV"); + if (!status) { + SKY_STDEV = 1.0; + psWarning("SKY_STDEV is not set --- defaulting to %f\n", SKY_STDEV); + } + // use poissonian errors or local-sky errors + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); + if (!status) { + POISSON_ERRORS = true; + psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n"); + } + pmSourceFitModelInit (15, 0.1, PS_SQR(SKY_STDEV), POISSON_ERRORS); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined"); + + FILE *f = fopen ("bias.dat", "w"); + if (f == NULL) psAbort ("psphot", "can't open output file"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip lower-quality objects + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (source->mode & PM_SOURCE_MODE_POOR) continue; + if (source->mode & PM_SOURCE_MODE_FAIL) continue; + if (source->mode & PM_SOURCE_MODE_PAIR) continue; + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + // if model is NULL, we don't have a starting guess + if (source->modelPSF == NULL) continue; + + if (source->moments->x < AnalysisRegion.x0) continue; + if (source->moments->y < AnalysisRegion.y0) continue; + if (source->moments->x > AnalysisRegion.x1) continue; + if (source->moments->y > AnalysisRegion.y1) continue; + + // replace object in image + pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + + // make a temporary model (we don't keep the result of this analysis) + pmModel *PSF = pmModelCopy (source->modelPSF); + + // extend source radius as needed + psphotCheckRadiusPSF (readout, source, PSF); + + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); + pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM); + psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // re-subtract PSF for object, leave local sky + pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + + PARp = source->modelPSF->params->data.F32; + PARc = PSF->params->data.F32; + fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky); + Nfit ++; + } + fclose (f); + psLogMsg ("psphot", 3, "measure PSF weighting bias for %d objects: %f sec\n", Nfit, psTimerMark ("psphot")); + + return (true); +}