Changeset 32633

Timestamp:

Nov 8, 2011, 2:56:56 PM (15 years ago)

Author:

eugene

Message:

add unofficial PS1_V4 output format (not yet standard); add clean radial profile to find sky limit; use sky limit to control kron analysis window; fix thread / memory problem for radial photometry; fix object matching for forced psf photometry; fix object matching to get the right matches

Location:

trunk

Files:

• psModules/src/objects/Makefile.am (modified) (3 diffs)
• psModules/src/objects/mksource.pl (modified) (1 diff)
• psModules/src/objects/pmSource.c (modified) (1 diff)
• psModules/src/objects/pmSource.h (modified) (2 diffs)
• psModules/src/objects/pmSourceIO.c (modified) (2 diffs)
• psModules/src/objects/pmSourceIO.h (modified) (2 diffs)
• psModules/src/objects/pmSourceIO_CMF.c.in (modified) (7 diffs)
• psphot (modified) (1 prop)
• psphot/src/Makefile.am (modified) (2 diffs)
• psphot/src/psphot.h (modified) (2 diffs)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (1 diff)
• psphot/src/psphotFitSourcesLinear.c (modified) (1 diff)
• psphot/src/psphotGuessModels.c (modified) (1 diff)
• psphot/src/psphotKronIterate.c (modified) (3 diffs)
• psphot/src/psphotRadialApertures.c (modified) (20 diffs)
• psphot/src/psphotRadialProfileWings.c (copied) (copied from branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c )
• psphot/src/psphotReadout.c (modified) (2 diffs)
• psphot/src/psphotSetThreads.c (modified) (1 diff)
• psphot/src/psphotSourceFits.c (modified) (1 diff)
• psphot/src/psphotSourceMatch.c (modified) (8 diffs)
• psphot/src/psphotSourceSize.c (modified) (1 diff)
• psphot/src/psphotSourceStats.c (modified) (1 diff)
• psphot/src/psphotStackObjects.c (modified) (5 diffs)
• psphot/src/psphotStackReadout.c (modified) (7 diffs)

trunk/psModules/src/objects/Makefile.am

@@ -45 +45 @@
         pmSourceIO_CMF_PS1_V2.c \
         pmSourceIO_CMF_PS1_V3.c \
+        pmSourceIO_CMF_PS1_V4.c \
         pmSourceIO_CMF_PS1_SV1.c \
         pmSourceIO_CMF_PS1_DV1.c \
@@ -130 +131 @@
 
 # pmSourceID_CMF_* functions use a common framework
-BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c
+BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c pmSourceIO_CMF_PS1_V4.c
 
 pmSourceIO_CMF_PS1_V1.c : pmSourceIO_CMF.c.in mksource.pl
@@ -141 +142 @@
         mksource.pl pmSourceIO_CMF.c.in PS1_V3 pmSourceIO_CMF_PS1_V3.c
 
+pmSourceIO_CMF_PS1_V4.c : pmSourceIO_CMF.c.in mksource.pl
+        mksource.pl pmSourceIO_CMF.c.in PS1_V4 pmSourceIO_CMF_PS1_V4.c
+
 # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in

trunk/psModules/src/objects/mksource.pl

@@ -16 +16 @@
 %cmfmodes = ("PS1_V1", 1,
              "PS1_V2", 2,
-             "PS1_V3", 3);
+             "PS1_V3", 3,
+             "PS1_V4", 4);
 
 print "1: $cmfmodes{1}\n";

trunk/psModules/src/objects/pmSource.c

@@ -137 +137 @@
     source->apFlux           = NAN;
     source->apFluxErr        = NAN; 
+
+    source->skyRadius        = NAN;
+    source->skyFlux          = NAN;
+    source->skySlope         = NAN;
 
     source->pixWeightNotBad  = NAN;

trunk/psModules/src/objects/pmSource.h

@@ -36 +36 @@
     PM_SOURCE_TMPF_MOMENTS_MEASURED  = 0x0010,
     PM_SOURCE_TMPF_CANDIDATE_PSFSTAR = 0x0020,
-    PM_SOURCE_TMPF_STACK_KEEP        = 0x0040,
-    PM_SOURCE_TMPF_STACK_SKIP        = 0x0080,
+    PM_SOURCE_TMPF_RADIAL_KEEP       = 0x0040,
+    PM_SOURCE_TMPF_RADIAL_SKIP       = 0x0080,
+    PM_SOURCE_TMPF_PETRO_KEEP        = 0x0100,
+    PM_SOURCE_TMPF_PETRO_SKIP        = 0x0200,
 } pmSourceTmpF;
 
@@ -95 +97 @@
     float apFluxErr;                    ///< apFluxErr corresponding to psfMag or extMag (depending on type)
 
+    float skyRadius;                    ///< radius at which profile hits local sky (or goes flat)
+    float skyFlux;                      ///< mean flux per pixel in aperture at which profile hits local sky (or goes flat)
+    float skySlope;                     ///< mean flux slope at which profile hits local sky (or goes flat)
+
     float pixWeightNotBad;              ///< PSF-weighted coverage of unmasked (not BAD) pixels
     float pixWeightNotPoor;             ///< PSF-weighted coverage of unmasked (not POOR) pixels

trunk/psModules/src/objects/pmSourceIO.c

@@ -569 +569 @@
             PM_SOURCES_WRITE("PS1_V2",    CMF_PS1_V2);
             PM_SOURCES_WRITE("PS1_V3",    CMF_PS1_V3);
+            PM_SOURCES_WRITE("PS1_V4",    CMF_PS1_V4);
             PM_SOURCES_WRITE("PS1_SV1",   CMF_PS1_SV1);
             PM_SOURCES_WRITE("PS1_DV1",   CMF_PS1_DV1);
@@ -1025 +1026 @@
                 sources = pmSourcesRead_CMF_PS1_V3 (file->fits, hdu->header);
             }
+            if (!strcmp (exttype, "PS1_V4")) {
+                sources = pmSourcesRead_CMF_PS1_V4 (file->fits, hdu->header);
+            }
             if (!strcmp (exttype, "PS1_SV1")) {
                 sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header);

trunk/psModules/src/objects/pmSourceIO.h

@@ -62 +62 @@
 bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 
+bool pmSourcesWrite_CMF_PS1_V4(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_V4_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_V4_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_V4_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
+
 bool pmSourcesWrite_CMF_PS1_SV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_SV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
@@ -86 +91 @@
 psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header);
+psArray *pmSourcesRead_CMF_PS1_V4 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);

trunk/psModules/src/objects/pmSourceIO_CMF.c.in

@@ -127 +127 @@
 
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
         @<PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
@@ -138 +138 @@
         @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
 
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
@@ -150 +150 @@
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
@@ -158 +158 @@
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
 
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
-
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
-
-        // XXX do not keep this long term, just a TEST:
-        // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF",    PS_DATA_F32, "Kron Flux",                                  moments.KronPSF);
-        // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF_SIG",PS_DATA_F32, "Kron Flux",                                  moments.KronPSFErr);
-        // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d
-        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX",   PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.KronCore);
-        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
+
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
+
+        @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
+        @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
+        @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
 
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
 
         // XXX not sure how to get this : need to load Nimages with weight?
@@ -286 +283 @@
         @ALL@     source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
-        @=PS1_V3@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
+        @>PS1_V2@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
 
         // XXX use these to determine PAR[PM_PAR_I0] if they exist?
-        @=PS1_V3@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
-        @=PS1_V3@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
+        @>PS1_V2@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
+        @>PS1_V2@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
@@ -311 +308 @@
 
         @ALL@     source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
-        @=PS1_V3@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
+        @>PS1_V2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
         @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
         @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
@@ -329 +326 @@
         @ALL@     source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
 
-        @=PS1_V3@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
-        @=PS1_V3@ source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
-        @=PS1_V3@ source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
-        @=PS1_V3@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
-
-        @=PS1_V3@ source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
-        @=PS1_V3@ source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
+        @>PS1_V2@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
+        @>PS1_V2@ source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
+        @>PS1_V2@ source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
+        @>PS1_V2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
+
+        @>PS1_V2@ source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
+        @>PS1_V2@ source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
+
+        @>PS1_V3@ source->skyRadius            = psMetadataLookupF32 (&status, row, "SKY_LIMIT_RAD");
+        @>PS1_V3@ source->skyFlux              = psMetadataLookupF32 (&status, row, "SKY_LIMIT_FLUX");
+        @>PS1_V3@ source->skySlope             = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE");
 
         // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
         // we are storing enough information so the output will be consistent with the input
-        @=PS1_V3@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
-        @=PS1_V3@ source->moments->Mxxy = 0.0;
-        @=PS1_V3@ source->moments->Mxyy = 0.0;
-        @=PS1_V3@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
-
-        @=PS1_V3@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
-        @=PS1_V3@ source->moments->Mxxxy = 0.0;
-        @=PS1_V3@ source->moments->Mxxyy = 0.0;
-        @=PS1_V3@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
-        @=PS1_V3@ source->moments->Myyyy = 0.0;
+        @>PS1_V2@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
+        @>PS1_V2@ source->moments->Mxxy = 0.0;
+        @>PS1_V2@ source->moments->Mxyy = 0.0;
+        @>PS1_V2@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
+
+        @>PS1_V2@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
+        @>PS1_V2@ source->moments->Mxxxy = 0.0;
+        @>PS1_V2@ source->moments->Mxxyy = 0.0;
+        @>PS1_V2@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
+        @>PS1_V2@ source->moments->Myyyy = 0.0;
 
         @ALL@     source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
-        @=PS1_V3@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
+        @>PS1_V2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
         assert (status);
 

trunk/psphot/src/Makefile.am

@@ -185 +185 @@
         psphotKronMasked.c             \
         psphotKronIterate.c            \
+        psphotRadialProfileWings.c     \
         psphotDeblendSatstars.c        \
         psphotMosaicSubimage.c         \
@@ -201 +202 @@
         psphotRadialBins.c             \
         psphotRadialApertures.c        \
-        psphotRadialAperturesByObject.c \
         psphotPetrosian.c              \
         psphotPetrosianRadialBins.c    \

trunk/psphot/src/psphot.h

@@ -428 +428 @@
 bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry);
 bool psphotRadialApertures_Threaded (psThreadJob *job);
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVer);
+// bool psphotRadialApertureSource (pmSource *source, int entry);
 
 bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
@@ -466 +467 @@
 bool psphotKronIterate_Threaded (psThreadJob *job);
 
+bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources);
+bool psphotRadialProfileWings_Threaded (psThreadJob *job);
+
 bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
 

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -209 +209 @@
         // if we have checked the source validity on the basis of the object set, then 
         // we either skip these tests below or we skip the source completely
-        if (source->tmpFlags & PM_SOURCE_TMPF_STACK_SKIP) continue;
-        if (source->tmpFlags & PM_SOURCE_TMPF_STACK_KEEP) goto keepSource;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -134 +134 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+
+        // psAssert (source->peak, "source without peak??");
+        // psAssert (source->peak->footprint, "peak without footprint??");
 
         // turn this bit off and turn it on again if we pass this test

trunk/psphot/src/psphotGuessModels.c

@@ -170 +170 @@
         if (!source->peak) continue;
 
+        // psAssert (source->peak->footprint, "peak without footprint??");
+
         nSrc ++;
 

trunk/psphot/src/psphotKronIterate.c

@@ -151 +151 @@
             }
 # else
-            if (!psphotExtendedSourceFits_Threaded(job)) {
+            if (!psphotKronIterate_Threaded(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
                 psFree(AnalysisRegion);
@@ -192 +192 @@
     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
 
-    for (int j = 0; j < 5; j++) {
+    // XXX TEST : set iteration to 1
+    for (int j = 0; j < 1; j++) {
         for (int i = 0; i < sources->n; i++) {
 
@@ -209 +210 @@
 
             // use S/N to control max window size
-            float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
+            // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
 
             // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux)
             // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2)
             // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0)
-            float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
+            // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
 
             // iterate to the window radius
-            float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
+            // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
+
+            // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
+            float maxWindow = source->skyRadius;
+            float windowRadius = PS_MAX(RADIUS, maxWindow);
 
             // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS

trunk/psphot/src/psphotRadialApertures.c

@@ -26 +26 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
         if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
@@ -35 +41 @@
     return true;
 }
+
+// these values are used by all threads repeatedly (and are not modified)
+static psVector *aperRadii = NULL;
+static psVector *aperRadii2 = NULL;
+static float outerRadius = NAN;
+static float SN_LIM = NAN;
+static psImageMaskType maskVal = 0;
 
 // aperture-like measurements for extended sources
@@ -75 +88 @@
         nThreads = 0;
     }
+
+    // aperRadii stores the upper bounds of the annuli
+    // XXX keep the same name here as for the petrosian / elliptical apertures?
+    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
+    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+
+    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
+
+    // save the R^2 values as well for quicker comparison
+    aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
+    for (int i = 0; i < aperRadii->n; i++) {
+        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // S/N limit to perform full non-linear fits
+    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     // source analysis is done in S/N order (brightest first)
@@ -117 +151 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, AnalysisRegion);
-            psArrayAdd(job->args, 1, recipe);
-
             PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
-
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
 
-// set this to 0 to run without threading
+            // set this to 0 to run without threading
 # if (1)            
             if (!psThreadJobAddPending(job)) {
@@ -138 +169 @@
             }
             psScalar *scalar = NULL;
-            scalar = job->args->data[6];
+            scalar = job->args->data[5];
             Nradial += scalar->data.S32;
             psFree(job);
@@ -158 +189 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[6];
+                scalar = job->args->data[5];
                 Nradial += scalar->data.S32;
             }
@@ -166 +197 @@
     psFree (cellGroups);
     psFree(AnalysisRegion);
+    psFree (aperRadii2);
 
     psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
@@ -173 +205 @@
 bool psphotRadialApertures_Threaded (psThreadJob *job) {
 
-    bool status;
     int Nradial = 0;
 
@@ -180 +211 @@
     psArray *sources        = job->args->data[1];
     psRegion *region        = job->args->data[2];
-    psMetadata *recipe      = job->args->data[3];
-
-    int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
-    int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
-
-    // radMax stores the upper bounds of the annuli
-    // XXX keep the same name here as for the petrosian / elliptical apertures?
-    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
-    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
-
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    assert (maskVal);
-
-    // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
-
-    float outerRadius = radMax->data.F32[radMax->n - 1];
+    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
+
+    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
 
     // choose the sources of interest
@@ -204 +223 @@
 
         pmSource *source = sources->data[i];
+
+        // if we have checked the source validity on the basis of the object set, then 
+        // we either skip these tests below or we skip the source completely
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects
@@ -222 +246 @@
         if (source->peak->x > region->x1) continue;
         if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
         // allocate pmSourceExtendedParameters, if not already defined
@@ -240 +266 @@
         }
 
-        // we need to change the view for the radial aperture analysis, but we want to recover exactly
-        // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
-        psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
-        psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
-        psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
-        psRegion oldRegion    = source->region;
-
         Nradial ++;
 
-        // force source image to be a bit larger...
-        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
-
-        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
+        if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
             psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
         } else {
@@ -258 +274 @@
         }
 
-        pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
-        psFree(source->maskObj);   source->maskObj   = oldMaskObj;
-        psFree(source->modelFlux); source->modelFlux = oldModelFlux;
-        psFree(source->psfImage);  source->psfImage  = oldPSFimage;
-        
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
-    psScalar *scalar = job->args->data[6];
+    psScalar *scalar = job->args->data[5];
     scalar->data.S32 = Nradial;
+
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
 
     return true;
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
-
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
+                                            
     // if we are a child source, save the results to the parent source radial aperture array
     psArray *radialAperSet = source->radialAper;
@@ -285 +300 @@
     radialAperSet->data[entry] = radialAper;
 
-    // storage for the derived pixel values
-    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
+    // find the largest aperture of interest (use only apertures with inner radii <=
+    // source->skyRadius)
+    int lastAp = aperRadii->n;
+    for (int i = 0; i < aperRadii->n; i++) {
+        if (aperRadii->data.F32[i] < source->skyRadius) continue;
+        lastAp = i + 1;
+        break;
+    }
 
     // outer-most radius for initial truncation
-    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax  = aperRadii->data.F32[lastAp - 1];
     float Rmax2 = PS_SQR(Rmax);
 
-    // store the R^2 values for the apertures
-    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
-    for (int i = 0; i < aperRadii->n; i++) {
-        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
-    }
+    // in this function, the operatins are relative to the full image (readout->image, etc)
 
     float xCM = NAN, yCM = NAN;
     if (pmSourcePositionUseMoments(source)) {
-        xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-        yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+        xCM = source->moments->Mx; // index coord of peak in readout
+        yCM = source->moments->My; // index coord of peak in readout
     } else {
-        xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
-        yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
-    }
+        xCM = source->peak->xf; // index coord of peak in readout
+        yCM = source->peak->yf; // index coord of peak in readout
+    }
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    pixRadius2->n = 0;
+    pixFlux->n = 0;
+    pixVar->n = 0;
 
     // one pass through the pixels to select the valid pixels and calculate R^2
-    for (int iy = 0; iy < source->pixels->numRows; iy++) {
-
-        float yDiff = iy - yCM;
-        if (fabs(yDiff) > Rmax) continue;
-
-        float *vPix = source->pixels->data.F32[iy];
-        float *vWgt = source->variance->data.F32[iy];
-        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
-
-        for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
-
-            if (vMsk) {
-                if (*vMsk & maskVal) {
-                    vMsk++;
-                    continue;
-                }
-                vMsk++;
-            }
-            if (isnan(*vPix)) continue;
-
-            float xDiff = ix - xCM;
-            if (fabs(xDiff) > Rmax) continue;
+    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
+
+        float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
+        int yPix = (int) yDiff;
+
+        if (yPix < 0) continue;
+        if (yPix > Ny - 1) continue;
+        if (fabs(iy) > Rmax) continue;
+
+        float *vPix = readout->image->data.F32[yPix];
+        float *vWgt = readout->variance->data.F32[yPix];
+        psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
+
+        for (int ix = -Rmax; ix < Rmax + 1; ix++) {
+
+            float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
+            int xPix = (int) xDiff;
+            
+            if (xPix < 0) continue;
+            if (xPix > Nx - 1) continue;
+            if (fabs(ix) > Rmax) continue;
+            
+            if (vMsk[xPix] & maskVal) continue;
+            if (isnan(vPix[xPix])) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(ix) + PS_SQR(iy);
             if (r2 > Rmax2) continue;
 
             psVectorAppend(pixRadius2, r2);
-            psVectorAppend(pixFlux, *vPix);
-            psVectorAppend(pixVar, *vWgt);
+            psVectorAppend(pixFlux, vPix[xPix]);
+            psVectorAppend(pixVar, vWgt[xPix]);
         }
     }
@@ -348 +372 @@
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
+    // init the apertures of interest to 0.0, the rest go to NAN
     psVectorInit (flux,    0.0);
     psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
+    for (int i = lastAp; i < flux->n; i++) {
+        flux->data.F32[i] = NAN;
+        fluxStd->data.F32[i] = NAN;
+        fluxErr->data.F32[i] = NAN;
+        fill->data.F32[i] = NAN;
+    }
 
     float *rPix2 = pixRadius2->data.F32;
@@ -358 +389 @@
         int j = 0;
         float *aRad2 = aperRadii2->data.F32;
-        for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
-        for (; j < aperRadii2->n; j++, aRad2++) {
+        for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
+
+        // XXX I can speed this up by only saving this single aperture
+        for (; j < lastAp; j++, aRad2++) {
             flux->data.F32[j]    += pixFlux->data.F32[i];
             fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
@@ -371 +404 @@
        2) the fractional fill factor (count of valid pixels / effective area of the aperture
        3) the error on the flux within that aperture
-     */
-
-    for (int i = 0; i < flux->n; i++) {
+    */
+
+    for (int i = 0; i < lastAp; i++) {
         // calculate the total flux for bin 'nOut'
         float Area = M_PI*aperRadii2->data.F32[i];
@@ -383 +416 @@
         // XXX report the total flux or the mask-corrected flux?
         // flux->data.F32[i]    = SBmean * Area;
-        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
 
         fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
@@ -393 +426 @@
     }
     
+# if (1)
     radialAper->flux = flux;
     radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
-
-    psFree (aperRadii2);
-    psFree (pixRadius2);
-    psFree (pixFlux);
-    psFree (pixVar);
+# else
+    // XXX TEST
+    psFree(flux);
+    psFree(fluxStd);
+    psFree(fluxErr);
+    psFree(fill);
+# endif
 
     return true;
 }
+
+/*** below is a test to use a sort to speed this up, not very successfully ***/
 
 static int nCalls = 0;
@@ -495 +533 @@
 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
 # define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_VECT(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B];  \
-    radius->data.F32[B] = tmp; \
-    tmp = pixFlux->data.F32[A]; \
-    pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
-    pixFlux->data.F32[B] = tmp; \
-    tmp = pixVar->data.F32[A]; \
-    pixVar->data.F32[A] = pixVar->data.F32[B]; \
-    pixVar->data.F32[B] = tmp; \
-  } \
-}
+# define SWAP_VECT(TYPE,A,B) {                                  \
+        float tmp;                                              \
+        if (A != B) {                                           \
+            tmp = radius->data.F32[A];                          \
+            radius->data.F32[A] = radius->data.F32[B];          \
+            radius->data.F32[B] = tmp;                          \
+            tmp = pixFlux->data.F32[A];                         \
+            pixFlux->data.F32[A] = pixFlux->data.F32[B];        \
+            pixFlux->data.F32[B] = tmp;                         \
+            tmp = pixVar->data.F32[A];                          \
+            pixVar->data.F32[A] = pixVar->data.F32[B];          \
+            pixVar->data.F32[B] = tmp;                          \
+        }                                                       \
+    }
 
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {

trunk/psphot/src/psphotReadout.c

@@ -190 +190 @@
     psphotDumpChisqs (config, view, filerule);
 
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, filerule);
+
     // re-measure the kron mags with models subtracted.  this pass uses a circular window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     
@@ -305 +308 @@
 pass1finish:
 
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, filerule);
+
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, filerule);

trunk/psphot/src/psphotSetThreads.c

@@ -50 +50 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 7);
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 6);
     task->function = &psphotRadialApertures_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_PROFILE_WINGS", 3);
+    task->function = &psphotRadialProfileWings_Threaded;
     psThreadTaskAdd(task);
     psFree(task);

trunk/psphot/src/psphotSourceFits.c

@@ -357 +357 @@
     pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
+    newSrc->peak->footprint = source->peak->footprint; // just a reference; the peak does not own the footprint
+    psArrayAdd(newSrc->peak->footprint->peaks, 1, newSrc->peak); // the footprint owns the peak
 
     // build cached models and subtract

trunk/psphot/src/psphotSourceMatch.c

@@ -67 +67 @@
     objects = psArraySort (objects, pmPhotObjSortByX); 
  
-    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
-    psVectorInit (found, 0);
+    psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundSrc, 0);
+
+    psVector *foundObj = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundObj, 0);
 
     // match sources to existing objects
@@ -93 +96 @@
         if (dx > +1.02*RADIUS) NEXT2; 
  
+        /* this block will match a given detection to the closest object within range of that detection.
+           XXX note that this matches ALL detections within range of the single object to that same object 
+           this is bad, but I cannot just go in linear order (ie, mark off each object as they are
+           used).  I should make a list of all Nobj * Ndet pairs in range and choose the matches
+           based on their separations.  UGH
+        */
+    
         // we are within match range, look for matches: 
         int Jmin = -1;
@@ -98 +108 @@
         for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
  
+            // skip objects that are already assigned:
+            if (foundObj->data.U8[J]) continue;
             obj = objects->data[J]; 
             
@@ -117 +129 @@
         }
         obj = objects->data[Jmin]; 
+        foundObj->data.U8[Jmin] = 1;
 
         // add to object
         pmPhotObjAddSource (obj, src);
-        found->data.U8[i] = 1;
+        foundSrc->data.U8[i] = 1;
         i++; 
     } 
@@ -128 +141 @@
     for (i = 0; i < sources->n; i++) {
 
-        if (found->data.U8[i]) continue;
+        if (foundSrc->data.U8[i]) continue;
 
         pmSource *src = sources->data[i]; 
@@ -139 +152 @@
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
 
-    psFree (found);
+    psFree (foundSrc);
+    psFree (foundObj);
     return true;
 } 
@@ -234 +248 @@
             
             // assign to a footprint on this readout->image
-            peak->footprint = pmFootprintCopyData(footprint, readout->image);
-
-            // the peak does not claim ownership of the footprint (it does not free it). save a copy of this 
-            // footprint on detections->footprints so we can free it later
-            psArrayAdd(detections->footprints, 100, peak->footprint); 
-            psFree (peak->footprint);
+            if (footprint) {
+                peak->footprint = pmFootprintCopyData(footprint, readout->image);
+
+                // the peak does not claim ownership of the footprint (it does not free it). save a copy of this 
+                // footprint on detections->footprints so we can free it later
+                psArrayAdd(detections->footprints, 100, peak->footprint); 
+                psFree (peak->footprint);
+            }
             
             // create a new source
@@ -268 +284 @@
         pmPhotObj *obj = objects->data[i]; 
         nSources += obj->sources->n;
+        psAssert (obj->sources->n == nImages, "failed to match sources?");
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);

trunk/psphot/src/psphotSourceSize.c

@@ -156 +156 @@
 
     // XXX fix this (was source->n  - first)
-    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
+    psLogMsg ("psphot.size", PS_LOG_WARN, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
 
     psphotVisualPlotSourceSize (recipe, readout->analysis, sources);

trunk/psphot/src/psphotSourceStats.c

@@ -124 +124 @@
         // add the peak
         source->peak = psMemIncrRefCounter(peak);
+
+        // psAssert (source->peak->footprint, "peak without footprint??");
 
         // allocate space for moments

trunk/psphot/src/psphotStackObjects.c

@@ -70 +70 @@
 
     // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_PETRO  = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
@@ -81 +82 @@
         // we check each source for an object and keep the object if any source is valid
 
-        bool keepObject = false;
+        bool keepObjectRadial = false;
+        bool keepObjectPetro = false;
         for (int j = 0; j < object->sources->n; j++) {
 
@@ -94 +96 @@
             if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
             
-        // optionally allow non-extended objects to get petrosians as well
-        if (!doPetroStars) {
-            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-            if (source->type == PM_SOURCE_TYPE_STAR) continue;
-        }
-
-            // limit selection to some SN limit
-            // assert (source->peak); // how can a source not have a peak?
-            // limit selection to some SN limit
-            bool skipSource = false;
-            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-                skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
-            } else {
-                skipSource = (sqrt(source->peak->detValue) < SN_LIM);
-            }
-            if (skipSource) continue;
-
             // limit selection by analysis region (this automatically apply
             if (source->peak->x < AnalysisRegion.x0) continue;
@@ -117 +102 @@
             if (source->peak->y > AnalysisRegion.y1) continue;
             
-            keepObject = true;
+            // SN limit tests for RADIAL APERTURES:
+            bool skipSourceRadial = false;
+            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+                skipSourceRadial = (source->moments->KronFlux < SN_LIM_RADIAL * source->moments->KronFluxErr);
+            } else {
+                skipSourceRadial = (sqrt(source->peak->detValue) < SN_LIM_RADIAL);
+            }
+            if (!skipSourceRadial) keepObjectRadial = true;
+
+            // SN limit tests for PETRO
+            bool skipSourcePetro = false;
+            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+                skipSourcePetro = (source->moments->KronFlux < SN_LIM_PETRO * source->moments->KronFluxErr);
+            } else {
+                skipSourcePetro = doPetroStars ? (sqrt(source->peak->detValue) < SN_LIM_PETRO) : true;
+            }
+            if (!skipSourcePetro) keepObjectPetro = true;
+
+            keepObjectPetro = true;
         }
 
@@ -128 +131 @@
             // avoid the single-detection tests
 
-            if (keepObject) {
-                source->tmpFlags |=  PM_SOURCE_TMPF_STACK_KEEP;
-                source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_SKIP;
+            if (keepObjectPetro) {
+                source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_KEEP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_SKIP;
             } else {
-                source->tmpFlags |=  PM_SOURCE_TMPF_STACK_SKIP;
-                source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_KEEP;
+                source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_SKIP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_KEEP;
+            }       
+
+            if (keepObjectRadial) {
+                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_KEEP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP;
+            } else {
+                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_SKIP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP;
             }       
         }

trunk/psphot/src/psphotStackReadout.c

@@ -43 +43 @@
 
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
+
+    psArray *objects = NULL; // used below after 'pass1finish' label
 
     // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
@@ -169 +171 @@
     psphotStackVisualFilerule(config, view, STACK_SRC);
 
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
     // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
@@ -184 +189 @@
     psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
+    // if we only do one pass, skip to extended source analysis
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
-
-    // if we only do one pass, skip to extended source analysis
-    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -238 +243 @@
         // XXX check on free of sources...
         psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
-
-        // NOTE: apply to ALL sources
-        psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
     }
 
 pass1finish:
+
+    // generate the objects (objects unify the sources from the different images) NOTE: could
+    // this just match the detections for the chisq image, and not bother measuring the source
+    // stats in that case...?
+    objects = psphotMatchSources (config, view, STACK_SRC);
+    psMemDump("matchsources");
+
+    psphotStackObjectsUnifyPosition (objects);
+
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+
+    // NOTE: apply to ALL sources
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
+
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, STACK_SRC);
 
     // re-measure the kron mags with models subtracted
@@ -254 +272 @@
 
     psMemDump("psfstats");
-
-    // generate the objects (objects unify the sources from the different images)
-    // XXX this could just match the detections for the chisq image, and not bother measuring the
-    // source stats in that case...
-    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
-    psMemDump("matchsources");
-
-    psphotStackObjectsUnifyPosition (objects);
-
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
@@ -281 +289 @@
 
     // measure circular, radial apertures (objects sorted by S/N)
-    // XXX can we just use psphotRadialApertures
-    // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none')
-    // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF); 
+    // this forces photometry on the undetected sources from other images
     psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     psMemDump("extmeas");
@@ -323 +329 @@
     psphotMagnitudes(config, view, STACK_SRC);
 
-    if (0 && !psphotEfficiency(config, view, STACK_DET)) {
+    // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the 
+    // PSF of the SRC image.  can we fake it by generating the PSF for DET as well (up above)? 
+    if (false && !psphotEfficiency(config, view, STACK_DET)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();