IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 20, 2010, 1:14:11 PM (16 years ago)
Author:
eugene
Message:

Many changes:

  • psphot-related issues

1) added a new feature to psLib/src/fft/psImageFFT to allow for a pre-FFTed kernel generated for a specific image dimension which can then be applied to an arbitrary number of images for convolution -- this avoids the kernel FFT for every image convolution.
2) updated the psMinimizeLMM to include 2 levels of tolerance: minTol and maxTol: if the minimization reaches minTol, it stops. if it only reaches maxTol within maxIter, it stops and is considered successful, if it fails to reach maxTol, then it fails. This allows us to accept fits that are actually acceptable (within error) even if they are not as close as ideally possible.
3) add new stat: psfWeightNotPoor (vs psfWeightNotBad) : the first gives the fraction of psf-weighted unmasked pixels considering any mask bits (except the internal 'mark' bit), while the former considers only 'bad' mask bits -- these are written to QF_PSF and QF_PSF_PERFECT in the CMF files.
4) define user-set parameters for max and min valid flux in the linear fit analysis: Note this was the cause of the non-negative fluxes in forced photometry -- the min limit was hard-wired to 0.0.
5) significance image is now constructed as (image + (image/1000)2) / sqrt(variance) so that bright sources have well-defined peaks (other wise, they become somewhat flat-topped as image = sqrt(variance), leading to ill-defined peaks).
6) modification of the visualization functions to accept facility and level values akin to psTrace
7) modification of the source fitting APIs to allow thread- and source- independent fit options (iterations, tolerances, etc)
8) use Kron magnitude as test for source size (CR vs EXT) instead of PSF-based aperture
9) set a min systematic error in the aperture mags when used for size classification significance
10) threaded the psf-convolved model (PCM) fitting process (extended source fits)
11) for extended sources, adjust the radius based on the footprint and re-calculate moments for guess; save the psf-based moments for output in psf table
12) for Sersic fitting, choose the best index with a grid search using only a few iterations; iterate fully on the selected value.
13) same for Sersic fitting using the PCM fitting process
14) fixed a bug in which the PSF candidate stars which failed in the psf fitting were not correctly unmarked as being PSF stars
15) define galaxy fit radius based on sky stdev (global, not local)
16) move the PCM code to psModules
17) save and report the raw aperture mag in addition to the curve-of-growth corrected one (PS1_V3)
18) save and report the Kron parameters and higher-order moments
19) some psModule header reorganization for code clarity
20) fixed a bug in which the diff stats were counting 'marked' pixels (outside area of interest) as 'masked'.
21) save the radial profile aperture sizes in the headers
22) better guess for Sersic parameters based on moments (requires setting the functional form so that the scale length is right)

  • ppSub-related:

1) ensure masked pixels are NANed in output diff image
2) add code to flag detections if they have bright positive neighbors (+ output of these in PS1_DV2)
3) define separate penalties for each image (based on their fwhm values) (this requires the penalties to be calculated later in the code).
4) define separate apertures for each image for flux normalization
5) choose aperture based on curve-of-growth (was based on fixed fraction of full aperture flux, and thus noisy)
6) some fine tuning of the penalty factor (this still seems arbitrary, and results are somewhat sensitive to the right value)

Location:
trunk/psModules
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

    r28013 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4655
    4756// NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes
    48 // XXX currently in the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes)
     57// XXX currently, the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes)
    4958
    5059bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    … …  
    6271    psF32 errMag, chisq, apRadius;
    6372    psS32 nPix, nDOF;
     73    char keyword1[80], keyword2[80];
    6474
    6575    pmChip *chip = readout->parent->parent;
    … …  
    95105    // we use this just to define the output vectors (which must be present for all objects)
    96106    bool status = false;
     107    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    97108    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    98109    psAssert (radMax, "this must have been defined and tested earlier!");
    99110    psAssert (radMax->n, "this must have been defined and tested earlier!");
     111    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
     112
     113    // write the radial profile apertures to header
     114    for (int i = 0; i < radMax->n; i++) {
     115      sprintf (keyword1, "RMIN_%02d", i);
     116      sprintf (keyword2, "RMAX_%02d", i);
     117      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     118      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     119    }
     120
    100121
    101122    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    193214        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    194215        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     216        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
     217        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    196218        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    197219        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    198220
    199221        // distinguish moments measure from window vs S/N > XX ??
    200         float mxx = source->moments ? source->moments->Mxx : NAN;
    201         float mxy = source->moments ? source->moments->Mxy : NAN;
    202         float myy = source->moments ? source->moments->Myy : NAN;
    203         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    204         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    205         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
    206 
    207         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
     222        float Mxx = source->moments ? source->moments->Mxx : NAN;
     223        float Mxy = source->moments ? source->moments->Mxy : NAN;
     224        float Myy = source->moments ? source->moments->Myy : NAN;
     225
     226        float Mrf  = source->moments ? source->moments->Mrf : NAN;
     227        float Mrh  = source->moments ? source->moments->Mrh : NAN;
     228        float Krf  = source->moments ? source->moments->KronFlux : NAN;
     229        float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
     230
     231        float Kinner = source->moments ? source->moments->KronFinner : NAN;
     232        float Kouter = source->moments ? source->moments->KronFouter : NAN;
     233
     234        float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
     235        float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
     236        float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
     237        float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
     238
     239        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
     240        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
     241        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
     242
     243        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
     244        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
     245        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
     246        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
     247
     248        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
     249        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
     250        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
     251        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
     252
     253        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
     254        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
     255
     256        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
     257        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
    208258
    209259        psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    … …  
    352402        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    353403        source->peak->flux = peakFlux;
     404        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     405        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    354406        source->peak->dx   = dPAR[PM_PAR_XPOS];
    355407        source->peak->dy   = dPAR[PM_PAR_YPOS];
    356         source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
    357 
    358         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     408        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     409          source->peak->SN = 1.0 / source->errMag;
     410        } else {
     411          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
     412        }
     413
     414        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
     415        source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
    359416        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    360417        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    … …  
    371428        source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
    372429
     430        source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
     431        source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
     432        source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
     433        source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
     434
     435        source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
     436        source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
     437
     438        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
     439        // we are storing enough information so the output will be consistent with the input
     440        source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
     441        source->moments->Mxxy = 0.0;
     442        source->moments->Mxyy = 0.0;
     443        source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
     444
     445        source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
     446        source->moments->Mxxxy = 0.0;
     447        source->moments->Mxxyy = 0.0;
     448        source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
     449        source->moments->Myyyy = 0.0;
     450
    373451        source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
     452        source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
    374453        assert (status);
    375454
    … …  
    391470    psF32 xErr, yErr;
    392471    int nRow = -1;
     472    char keyword1[80], keyword2[80];
    393473
    394474    // create a header to hold the output data
    … …  
    422502    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    423503    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    424     // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    425     // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    426504
    427505    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    429507    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
    430508
    431     // int nRadialBins = 0;
    432     // if (doAnnuli) {
    433     //  // get the max count of radial bins
    434     //  for (int i = 0; i < sources->n; i++) {
    435     //      pmSource *source = sources->data[i];
    436     //      if (!source->extpars) continue;
    437     //      if (!source->extpars->radProfile ) continue;
    438     //         if (!source->extpars->radProfile->binSB) continue;
    439     //      nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n);
    440     //  }
    441     // }
     509    // write the radial profile apertures to header
     510    for (int i = 0; i < radMax->n; i++) {
     511      sprintf (keyword1, "RMIN_%02d", i);
     512      sprintf (keyword2, "RMAX_%02d", i);
     513      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     514      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     515    }
    442516
    443517    // we write out all sources, regardless of quality.  the source flags tell us the state
    444518    for (int i = 0; i < sources->n; i++) {
    445         // skip source if it is not a ext sourc
    446         // XXX we have two places that extended source parameters are measured:
    447         // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
    448         // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
    449         // should we require both?
    450 
    451519        pmSource *source = sources->data[i];
    452520
    … …  
    514582        }
    515583
    516 # if (0)
    517         // Kron measurements
    518         if (doKron) {
    519             pmSourceKronValues *kron = source->extpars->kron;
    520             if (kron) {
    521                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
    522                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
    523                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
    524                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
    525             } else {
    526                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
    527                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
    528                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
    529                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
    530             }
    531         }
    532 
    533         // Isophot measurements
    534         // XXX insert header data: isophotal level
    535         if (doIsophotal) {
    536             pmSourceIsophotalValues *isophot = source->extpars->isophot;
    537             if (isophot) {
    538                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
    539                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
    540                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
    541                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
    542             } else {
    543                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
    544                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
    545                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
    546                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
    547             }
    548         }
    549 # endif
    550 
    551584        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
    552585        if (doAnnuli) {
    … …  
    616649
    617650// XXX this layout is still the same as PS1_DEV_1
    618 bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
     651bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    619652{
    620653
    … …  
    665698            assert (model);
    666699
     700            // skip models which were not actually fitted
     701            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     702
    667703            PAR = model->params->data.F32;
    668704            dPAR = model->dparams->data.F32;
Note: See TracChangeset for help on using the changeset viewer.