IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31154


Ignore:
Timestamp:
Apr 4, 2011, 1:12:39 PM (15 years ago)
Author:
eugene
Message:

use the smoothed image for footprint culling; use a more stringent culling for saturated stars; more tweaks to get sat stars to be fitted; various updates to psphotStack; unify psphotImageLoop varients; be a bit careful about number of stars in a PSF clump

Location:
trunk/psphot
Files:
68 edited
3 copied

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/doc/stack.txt

    r30624 r31154  
     1
     220110303
     3
     4  psphotStack sample timing:
     5
     6  * convolve images & stamps: 78 sec
     7  * standard photometry : 250?
     8  * petrosians : 6 sec
     9  * ext fits : 350 sec
     10  * radial aps : 140 sec
     11
     1220110302
     13
     14  * total number of images generated by this process:
     15
     16 load:
     17 raw      : Nfilter x (image + variance + mask) (20M)
     18 cnv      : Nfilter x (image + variance + mask) (20M)
     19 (8 * 4B) + (4 * 2B)
     20
     21 match: (added)
     22 match    : Nfilter x (image + variance + mask) (20M)
     23(4 * 4B) + (2 * 2B)
     24
     25 phot: (added)
     26 chisq    :       1 x (image + variance + mask) (10M) (freed at end of psphotStackReadout?)
     27 backgnd  :       1 x (image) (10M) or Nfilter (20M) [unsure]
     28
     29 load:
     30     1089174         4000000    psFitsImage.c:467
     31     1092246         4000000    psFitsImage.c:467
     32     1095902         4000000    psFitsImage.c:467
     33     1098974         4000000    psFitsImage.c:467
     34     1102619         4000000    psFitsImage.c:467
     35     1105669         4000000    psFitsImage.c:467
     36     1109303         4000000    psFitsImage.c:467
     37     1112353         4000000    psFitsImage.c:467
     38     1090766         2000000    psFitsImage.c:467
     39     1097494         2000000    psFitsImage.c:467
     40     1104200         2000000    psFitsImage.c:467
     41     1110884         2000000    psFitsImage.c:467
     42
     43  match:
     44     1128976         4000000    pmSubtractionMatch.c:178
     45     1128979         4000000    pmSubtractionMatch.c:183
     46     1652553         4000000    pmSubtractionMatch.c:178
     47     1652556         4000000    pmSubtractionMatch.c:183
     48     1128982         2000000    pmSubtractionMatch.c:189
     49     1652559         2000000    pmSubtractionMatch.c:189
     50
     51  chisq:
     52     2133004         4000000    pmFPACopy.c:71
     53     2133010         4000000    pmFPACopy.c:71
     54     2133007         2000000    pmFPACopy.c:71
     55
     56  backmdl: 
     57     2132413         4000000    pmFPAfileDefine.c:1275
     58
     59  ???
     60     8614548         2000000    psBinaryOp.c:502
     61     8617313         2000000    psBinaryOp.c:502
     62
     63  pmSource (modelFlx):
     64     2775 * 6k = 16.9M
     65
     66  pmSource (maskObj):
     67     2775 * 3k =  8.4M
     68
     69  (span + footprints account for ~5 - 10 M, rest in little things)
    170
    27120101221

  • trunk/psphot/src/Makefile.am

    r30624 r31154  
    6969        psphotForced.c             \
    7070        psphotForcedArguments.c    \
    71         psphotForcedImageLoop.c    \
    72         psphotForcedReadout.c      \
    7371        psphotParseCamera.c        \
     72        psphotImageLoop.c          \
    7473        psphotMosaicChip.c         \
    7574        psphotCleanup.c
    … …  
    7978        psphotMakePSF.c            \
    8079        psphotMakePSFArguments.c   \
    81         psphotMakePSFImageLoop.c   \
    82         psphotMakePSFReadout.c     \
    8380        psphotParseCamera.c        \
     81        psphotImageLoop.c   \
    8482        psphotMosaicChip.c         \
    8583        psphotCleanup.c
    … …  
    128126        psphotReadoutForcedKnownSources.c    \
    129127        psphotReadoutMinimal.c         \
     128        psphotForcedReadout.c          \
     129        psphotMakePSFReadout.c         \
    130130        psphotModelBackground.c        \
     131        psphotMaskBackground.c         \
    131132        psphotSubtractBackground.c     \
    132133        psphotFindDetections.c         \

  • trunk/psphot/src/models/pmModel_STRAIL.c

    r19881 r31154  
    496496
    497497    params[PM_PAR_SKY]  = Smoments->Sky;
    498     params[PM_PAR_I0]   = peak->flux;
     498    params[PM_PAR_I0]   = peak->rawFlux;
    499499    params[PM_PAR_XPOS] = peak->xf;
    500500    params[PM_PAR_YPOS] = peak->yf;

  • trunk/psphot/src/models/pmModel_TEST1.c

    r19881 r31154  
    139139
    140140    PAR[PM_PAR_SKY] = moments->Sky;
    141     PAR[PM_PAR_I0]   = peak->flux;
     141    PAR[PM_PAR_I0]   = peak->rawFlux;
    142142    PAR[PM_PAR_XPOS] = peak->xf;
    143143    PAR[PM_PAR_YPOS] = peak->yf;

  • trunk/psphot/src/psphot.c

    r24144 r31154  
    11# include "psphotStandAlone.h"
     2# define FORCED_PHOTOMETRY 0
    23
    34int main (int argc, char **argv) {
    … …  
    2021
    2122    // call psphot for each readout
    22     if (!psphotImageLoop (config)) {
     23    if (!psphotImageLoop (config, PSPHOT_SINGLE)) {
    2324        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
    2425        exit (psphotGetExitStatus());

  • trunk/psphot/src/psphot.h

    r30624 r31154  
    1212
    1313#define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs
     14
     15# define READOUT_OR_INTERNAL(VIEW,FILE)((FILE)->mode == PM_FPA_MODE_INTERNAL) ? (FILE)->readout : pmFPAviewThisReadout((VIEW), (FILE)->fpa)
     16
     17typedef enum {
     18    PSPHOT_SINGLE,
     19    PSPHOT_FORCED,
     20    PSPHOT_MAKE_PSF,
     21} psphotImageLoopMode;
    1422
    1523// top-level psphot functions
    … …  
    2230void            psphotVersionPrint(void);
    2331
     32bool            psphotImageLoop (pmConfig *config, psphotImageLoopMode mode);
     33
    2434bool            psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
    2535bool            psphotInit (void);
    … …  
    5363bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    5464
     65bool            psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
     66bool            psphotMaskBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     67
    5568bool            psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
    5669bool            psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    … …  
    168181
    169182// used by psphotFindDetections
    170 psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
    171 psArray        *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
    172 bool            psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
    173 psErrorCode     psphotCullPeaks(const pmReadout *readout, const psMetadata *recipe, psArray *footprints);
     183pmReadout      *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
     184psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
     185bool            psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
     186psErrorCode     psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints);
    174187
    175188// in psphotApResid.c:
    … …  
    248261bool            psphotVisualShowLogSignificance (psImage *image, float min, float max);
    249262bool            psphotVisualShowPeaks (pmDetections *detections);
     263bool            psphotVisualShowSources (psArray *sources);
    250264bool            psphotVisualShowFootprints (pmDetections *detections);
    251265bool            psphotVisualShowMoments (psArray *sources);
    … …  
    305319
    306320pmConfig *psphotForcedArguments(int argc, char **argv);
    307 bool psphotForcedImageLoop (pmConfig *config);
    308321bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
    309322
    310323pmConfig *psphotMakePSFArguments(int argc, char **argv);
    311 bool psphotMakePSFImageLoop (pmConfig *config);
    312324bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
    313325
    … …  
    430442bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
    431443bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    432 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry);
     444bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
    433445
    434446bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);

  • trunk/psphot/src/psphotApResid.c

    r30624 r31154  
    11# include "psphotInternal.h"
     2// # define DEBUG
    23
    34# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
    … …  
    116117
    117118    // set limits on the aperture magnitudes
    118     pmSourceMagnitudesInit (recipe);
     119    pmSourceMagnitudesInit (config, recipe);
    119120
    120121    // threaded measurement of the source magnitudes
    … …  
    465466        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
    466467
    467         bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     468        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    468469
    469470        // clear the mask bit
    … …  
    481482        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    482483            Nfail ++;
    483             psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
     484            psTrace ("psphot", 3, "fail : %f, %f : nan mags : %f %f", source->peak->xf, source->peak->yf, source->apMag, source->psfMag);
    484485            continue;
    485486        }

  • trunk/psphot/src/psphotArguments.c

    r29004 r31154  
    109109    }
    110110
    111     PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv );
     111    // generic arguments (version, dumpconfig)
     112    PS_ARGUMENTS_GENERIC( psphot, config, argc, argv );
     113
     114    // thread arguments
     115    PS_ARGUMENTS_THREADS( psphot, config, argc, argv )
    112116
    113117    // save the following additional recipe values based on command-line options
    114118    // these options override the PSPHOT recipe values loaded from recipe files
    115119    psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE);
    116 
    117     // Number of threads is handled
    118     PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
    119120
    120121    // run the test model (requires X,Y coordinate)

  • trunk/psphot/src/psphotBasicDeblend.c

    r29936 r31154  
    6363    for (int i = 0; i < SN->n; i++) {
    6464        source = sources->data[i];
    65         SN->data.F32[i] = source->peak->SN;
     65        SN->data.F32[i] = source->peak->rawFlux;
    6666    }
    6767    psVector *index = psVectorSortIndex (NULL, SN);
    … …  
    115115        // threshold is fraction of the source peak flux
    116116        // image is background subtracted; source->moments->Sky should always be 0.0
    117         threshold = FRACTION * source->peak->SN;
     117        threshold = FRACTION * sqrt(source->peak->detValue);
    118118        // threshold is no less than NSIGMA
    119119        threshold = PS_MAX (threshold, NSIGMA);
    … …  
    138138            testSource = overlap->data[k];
    139139            if (testSource->mode & PM_SOURCE_MODE_BLEND) continue;
    140             if (testSource->peak->value > source->peak->value) continue;
     140            if (testSource->peak->rawFlux > source->peak->rawFlux) continue;
    141141            for (int j = 0; j < xv->n; j+=2) {
    142142                if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;

  • trunk/psphot/src/psphotBlendFit.c

    r30624 r31154  
    9494    fitOptions->weight        = PS_SQR(skySig);
    9595    fitOptions->mode          = PM_SOURCE_FIT_PSF;
     96    fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
    9697
    9798    psphotInitLimitsPSF (recipe, readout);
    … …  
    103104
    104105    // source analysis is done in S/N order (brightest first)
    105     sources = psArraySort (sources, pmSourceSortBySN);
     106    sources = psArraySort (sources, pmSourceSortByFlux);
    106107    if (!sources->n) {
    107108        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
    … …  
    176177            }
    177178            psFree(job);
    178             }
     179        }
    179180    }
    180181    psFree (cellGroups);
    … …  
    242243
    243244        // limit selection to some SN limit
    244         if (source->peak->SN < FIT_SN_LIM) continue;
     245        if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
    245246
    246247        // exclude sources outside optional analysis region

  • trunk/psphot/src/psphotChoosePSF.c

    r30624 r31154  
    7070
    7171    // examine PSF sources in S/N order (brightest first)
    72     sources = psArraySort (sources, pmSourceSortBySN);
     72    sources = psArraySort (sources, pmSourceSortByFlux);
    7373
    7474    // structure to store user options defining the psf
    … …  
    160160    options->fitOptions->weight        = PS_SQR(SKY_SIG);
    161161    options->fitOptions->mode          = PM_SOURCE_FIT_PSF;
     162    options->fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
     163
    162164   
    163165    psArray *stars = psArrayAllocEmpty (sources->n);
    … …  
    178180    psphotCheckStarDistribution (stars, sources, options);
    179181
    180     psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
     182    psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected %ld candidate PSF objects\n", stars->n);
    181183
    182184    if (stars->n < 50) {
    … …  
    336338    }
    337339
     340    // XXX does this work here?
     341    psphotVisualShowPSFStars (recipe, try->psf, try->sources);
     342
    338343    // build the flux-to-magnitude conversion information
    339344    if (!psphotMakeFluxScale (readout->image, recipe, try->psf)) {
    … …  
    443448            }
    444449            psFree (modelPSF);
     450
     451            // float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
     452            // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
    445453        }
    446454    }

  • trunk/psphot/src/psphotCullPeaks.c

    r27673 r31154  
    66psErrorCode
    77psphotCullPeaks(const pmReadout *readout,   // the image wherein lives the footprint
     8                const pmReadout *signifR, // smoothed image and significance
    89                const psMetadata *recipe,
    910                psArray *footprints) {  // array of pmFootprints
    … …  
    2627    psAssert (status, "cannot find SKY_STDEV in readout->analysis");
    2728
    28     const float min_threshold = nsigma_min*skyStdev;
     29    // the sky has been smoothed, so we need to scale down the raw sky stdev by this amount:
     30    float scaleFactor = psMetadataLookupF32(&status, readout->analysis, "SIGNIFICANCE_SCALE_FACTOR");
     31    if (!status) scaleFactor = 1.0;
     32
     33    const float MIN_THRESHOLD = nsigma_min*skyStdev / sqrt(scaleFactor);
    2934   
     35    // for saturated stars, we should be somewhat more agressive about culling: instead of
     36    // letting the height of the intervening col (saddle point) be set by the S/N of the image
     37    // pixel, it should be set to a fraction of the saturation value.  example for a
     38    // near-saturation pixel:
     39    // flux = 40k, sigma = 200
     40    // nsigma_delta = 4.0, nsigma_min = 1.0, fPadding = 0.01,
     41    // fStdev = 0.005, stdev_pad = 0.011*40k = 400
     42    // threshold = 40k - 4*400 = 38400
     43    // this gives too tight of a tolerance on the bright stars
     44    // in practice, we should make the threshold much lower. 
     45    // below I am using 0.05 * saturation (eg, 2000 DN above sky in a GPC1 image)
     46
     47    float SATURATION = NAN;
     48
     49    // do not completely trust the values in the header...
     50    float CELL_SATURATION = psMetadataLookupF32 (&status, readout->parent->concepts, "CELL.SATURATION");
     51    float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION");
     52    if (!status || !isfinite(MIN_SATURATION)) {
     53        MIN_SATURATION = 40000.0;
     54    }
     55    if (!isfinite(CELL_SATURATION)) {
     56        SATURATION = MIN_SATURATION;
     57    } else {
     58        SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION);
     59    }
     60    float SAT_TEST_LEVEL = 0.50*SATURATION;
     61    float SAT_THRESHOLD  = 0.05*SATURATION;
     62
     63# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE)
     64    psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the smoothed image");
     65# else
     66    psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the raw (unsmoothed) image");
     67# endif
     68
    3069    for (int i = 0; i < footprints->n; i++) {
    31         // if (i % 50 == 0) fprintf (stderr, "cull %d\n", i);
    3270        pmFootprint *fp = footprints->data[i];
     71        if (fp->peaks == NULL) continue;
     72        if (fp->peaks->n == 0) continue;
     73
    3374        if (fp->npix > 30000) {
    3475            fprintf (stderr, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix);
    3576        }
    3677        psTimerStart ("psphot.cull.footprints");
    37         if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) {
     78
     79        pmPeak *brightPeak = fp->peaks->data[0];
     80        float max_threshold = SAT_TEST_LEVEL;
     81        if (brightPeak->rawFlux > SAT_TEST_LEVEL) {
     82            max_threshold = SAT_THRESHOLD;
     83            brightPeak->type = PM_PEAK_SUSPECT_SATURATION;
     84        }
     85
     86# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE)
     87        // New (post r30869) style of culling using the smoothed image and variance (S/N)
     88        // if we cull using the significance image, then the definition of variance is different (thus the bool in arg 8)
     89        if (pmFootprintCullPeaks(signifR->image, signifR->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) {
    3890            return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
    3991        }
     92# else
     93        // Old (pre r30869) style of culling using the raw image and variance
     94        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, false) != PS_ERR_NONE) {
     95             return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
     96        }
     97# endif
     98
    4099        float dtime = psTimerMark ("psphot.cull.footprints");
    41100        if (dtime > 1.0) {

  • trunk/psphot/src/psphotDeblendSatstars.c

    r29936 r31154  
    7474    for (int i = 0; i < SN->n; i++) {
    7575        source = sources->data[i];
    76         SN->data.F32[i] = source->peak->SN;
     76        SN->data.F32[i] = source->peak->rawFlux;
    7777    }
    7878    psVector *index = psVectorSortIndex (NULL, SN);
    … …  
    8686        // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    8787        if (source->mode & PM_SOURCE_MODE_BLEND) continue;
    88         if (source->peak->flux < SATURATION) continue;
     88        if (source->peak->rawFlux < SAT_TEST_LEVEL) continue;
    8989
    9090        // save these for reference below

  • trunk/psphot/src/psphotEllipticalProfile.c

    r27819 r31154  
    8282    // }
    8383
    84     psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, source->peak->flux, 0.0);
     84    psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, source->peak->rawFlux, 0.0);
    8585    // psphotPetrosianVisualProfileByAngle (radius, flux);
    8686

  • trunk/psphot/src/psphotExtendedSourceAnalysis.c

    r30624 r31154  
    11# include "psphotInternal.h"
    22
    3 // ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis
    4 // XXX option to choose a consistent position
    5 // XXX option to choose a consistent elliptical contour
    6 // XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands
    7 // XXX consistent choice of extendedness...
     3// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
     4// XXX this function needs to be threaded
    85
    96// for now, let's store the detections on the readout->analysis for each readout
    … …  
    4037    int Next = 0;
    4138    int Npetro = 0;
    42     int Nisophot = 0;
    4339    int Nannuli = 0;
    44     int Nkron = 0;
    4540
    4641    psTimerStart ("psphot.extended");
    … …  
    6863    assert (maskVal);
    6964
    70     // XXX require petrosian analysis for non-linear fits?
    71 
    72     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    73     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
     65    // get the sky noise from the background analysis; if this is missing, get the user-supplied value
     66    float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV");
     67    if (!status) {
     68        skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
     69        psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise);
     70    }
    7471
    7572    // S/N limit to perform full non-linear fits
    … …  
    7875    // which extended source analyses should we perform?
    7976    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    80     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    8177    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    82     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    83 
    84 # if (0)
    85     // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    86     // XXX use this to set skynoise
    87     pmReadout *backModel = psphotSelectBackground (config, view);
    88     pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
    89 # endif
     78    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
    9079
    9180    // source analysis is done in S/N order (brightest first)
    92     sources = psArraySort (sources, pmSourceSortBySN);
     81    sources = psArraySort (sources, pmSourceSortByFlux);
    9382
    9483    // option to limit analysis to a specific region
    … …  
    10392
    10493        // skip PSF-like and non-astronomical objects
    105         if (source->type == PM_SOURCE_TYPE_STAR) continue;
    10694        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    10795        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    10896        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
    10997        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    110         if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     98
     99        // optionally allow non-extended objects to get petrosians as well
     100        if (!doPetroStars) {
     101            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     102            if (source->type == PM_SOURCE_TYPE_STAR) continue;
     103        }
    111104
    112105        // limit selection to some SN limit
    113106        assert (source->peak); // how can a source not have a peak?
    114         if (source->peak->SN < SN_LIM) continue;
     107        if (sqrt(source->peak->detValue) < SN_LIM) continue;
    115108
    116109        // limit selection by analysis region
    … …  
    119112        if (source->peak->x > AnalysisRegion.x1) continue;
    120113        if (source->peak->y > AnalysisRegion.y1) continue;
    121 
    122         // fprintf (stderr, "xsrc: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN);
    123114
    124115        // replace object in image
    … …  
    136127
    137128        // if we request any of these measurements, we require the radial profile
    138         if (doPetrosian || doIsophotal || doAnnuli || doKron) {
     129        if (doPetrosian || doAnnuli) {
    139130            if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
    140131                // all measurements below require the radial profile; skip them all
    … …  
    144135                continue;
    145136            }
     137            Nannuli ++;
    146138            source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    147139        }
    … …  
    169161    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    170162    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    171     psLogMsg ("psphot", PS_LOG_INFO, "  %d isophotal\n", Nisophot);
    172163    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    173     psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
    174164
    175165    psphotVisualShowResidualImage (readout, false);

  • trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c

    r30624 r31154  
    6161
    6262    // source analysis is done in S/N order (brightest first)
    63     objects = psArraySort (objects, pmPhotObjSortBySN);
     63    objects = psArraySort (objects, pmPhotObjSortByFlux);
    6464
    6565    // process the objects in order. 
    … …  
    8888            // limit selection to some SN limit
    8989            assert (source->peak); // how can a source not have a peak?
    90             if (source->peak->SN < SN_LIM) continue;
     90            if (sqrt(source->peak->detValue) < SN_LIM) continue;
    9191            measureSource = true;
    9292        }
    … …  
    148148                psFree(source->extpars->radFlux);
    149149                psFree(source->extpars->ellipticalFlux);
     150                psFree(source->extpars->petProfile);
    150151            }
    151152        }

  • trunk/psphot/src/psphotExtendedSourceFits.c

    r30624 r31154  
    4343    int NplainPass = 0;
    4444    int Nfaint = 0;
     45    int Nfail = 0;
    4546
    4647    psTimerStart ("psphot.extended");
    … …  
    145146
    146147    // source analysis is done in S/N order (brightest first)
    147     sources = psArraySort (sources, pmSourceSortBySN);
     148    sources = psArraySort (sources, pmSourceSortByFlux);
    148149
    149150    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    … …  
    176177            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
    177178            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
    178             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
    179 
    180             if (false && !psThreadJobAddPending(job)) {
     179            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
     180            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     181
     182// set this to 0 to run without threading
     183# if (1)           
     184            if (!psThreadJobAddPending(job)) {
    181185                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    182186                psFree(AnalysisRegion);
    183187                return false;
    184             } else {
    185                 if (!psphotExtendedSourceFits_Threaded(job)) {
    186                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    187                     psFree(AnalysisRegion);
    188                     return false;
    189                 }
    190                 psScalar *scalar = NULL;
    191                 scalar = job->args->data[7];
    192                 Next += scalar->data.S32;
    193                 scalar = job->args->data[8];
    194                 Nconvolve += scalar->data.S32;
    195                 scalar = job->args->data[9];
    196                 NconvolvePass += scalar->data.S32;
    197                 scalar = job->args->data[10];
    198                 Nplain += scalar->data.S32;
    199                 scalar = job->args->data[11];
    200                 NplainPass += scalar->data.S32;
    201                 scalar = job->args->data[12];
    202                 Nfaint += scalar->data.S32;
    203                 psFree(job);
     188            }
     189# else
     190            if (!psphotExtendedSourceFits_Threaded(job)) {
     191                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     192                psFree(AnalysisRegion);
     193                return false;
    204194            }
    205         }
     195            psScalar *scalar = NULL;
     196            scalar = job->args->data[7];
     197            Next += scalar->data.S32;
     198            scalar = job->args->data[8];
     199            Nconvolve += scalar->data.S32;
     200            scalar = job->args->data[9];
     201            NconvolvePass += scalar->data.S32;
     202            scalar = job->args->data[10];
     203            Nplain += scalar->data.S32;
     204            scalar = job->args->data[11];
     205            NplainPass += scalar->data.S32;
     206            scalar = job->args->data[12];
     207            Nfaint += scalar->data.S32;
     208            scalar = job->args->data[13];
     209            Nfail += scalar->data.S32;
     210            psFree(job);
     211# endif
     212        }
    206213
    207214        // wait for the threads to finish and manage results
    … …  
    231238                scalar = job->args->data[12];
    232239                Nfaint += scalar->data.S32;
     240                scalar = job->args->data[13];
     241                Nfail += scalar->data.S32;
    233242            }
    234243            psFree(job);
    … …  
    238247    psFree(AnalysisRegion);
    239248
    240     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     249    psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    241250    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
    242251    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
    243     psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit\n", Nfaint);
     252    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit, %d failed\n", Nfaint, Nfail);
    244253    return true;
    245254}
    … …  
    250259    bool status;
    251260    int Next = 0;
     261    int Nfaint = 0;
     262    int Nfail = 0;
    252263    int Nconvolve = 0;
    253264    int NconvolvePass = 0;
    254265    int Nplain = 0;
    255     int Nfaint = 0;
    256266    int NplainPass = 0;
    257267    bool savePics = false;
    … …  
    268278    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    269279
     280    // pthread_t tid = pthread_self();     // Thread identifier
     281
    270282    // Define source fitting parameters for extended source fits
    271283    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
    272     fitOptions->mode          = PM_SOURCE_FIT_EXT;
     284    fitOptions->mode           = PM_SOURCE_FIT_EXT;
     285    fitOptions->saveCovariance = true;  // XXX make this a user option?
     286    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
     287
    273288    // XXX for now, use the defaults for the rest:
    274289    // fitOptions->nIter         = fitIter;
    … …  
    296311        // XXX use the parameters guessed from moments
    297312        // if (source->modelEXT == NULL) continue;
     313
     314        // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
    298315
    299316        // replace object in image
    … …  
    366383          // limit selection to some SN limit
    367384          assert (source->peak); // how can a source not have a peak?
    368           if (source->peak->SN < SNlim) {
     385          if (sqrt(source->peak->detValue) < SNlim) {
    369386              Nfaint ++;
    370387              continue;
    … …  
    388405              if (!modelFit) {
    389406                  psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
     407                  Nfail ++;
    390408                  continue;
    391409              }
    … …  
    402420              if (!modelFit) {
    403421                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
     422                  Nfail ++;
    404423                  continue;
    405424              }
    … …  
    519538    scalar->data.S32 = Nfaint;
    520539
     540    scalar = job->args->data[13];
     541    scalar->data.S32 = Nfail;
     542
    521543    return true;
    522544}

  • trunk/psphot/src/psphotFindDetections.c

    r30624 r31154  
    8484
    8585    // generate the smoothed significance image
    86     psImage *significance = psphotSignificanceImage (readout, recipe, maskVal);
     86    pmReadout *significance = psphotSignificanceImage (readout, recipe, maskVal);
    8787
    8888    // display the log significance image
    89     psphotVisualShowLogSignificance (significance, 0.0, 4.5);
     89    psphotVisualShowLogSignificance (significance->variance, 0.0, 4.5);
    9090
    9191    // display the significance image
    92     psphotVisualShowSignificance (significance, 0.98*threshold, 1.02*threshold);
     92    psphotVisualShowSignificance (significance->variance, 0.98*threshold, 1.02*threshold);
    9393
    9494    // detect the peaks in the significance image

  • trunk/psphot/src/psphotFindFootprints.c

    r30624 r31154  
    11# include "psphotInternal.h"
    22
    3 bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) {
     3bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) {
    44
    55    bool status;
    … …  
    1818    PS_ASSERT (status, false);
    1919
    20     // find the raw footprints & assign the peaks to those footprints
    21     psArray *footprints = pmFootprintsFind (significance, threshold, npixMin);
     20    // find the raw footprints in the smoothed significance image & assign the peaks to those footprints
     21    psArray *footprints = pmFootprintsFind (significance->variance, threshold, npixMin);
    2222
    2323    if (pmFootprintsAssignPeaks(footprints, detections->peaks) != PS_ERR_NONE) {
    … …  
    5656    }
    5757
    58     psphotCullPeaks(readout, recipe, detections->footprints);
     58    psphotCullPeaks(readout, significance, recipe, detections->footprints);
    5959    detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
    6060    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));

  • trunk/psphot/src/psphotFindPeaks.c

    r26894 r31154  
    44// image must be constructed to represent (S/N)^2.  If nMax is non-zero, only return a maximum
    55// of nMax peaks
    6 psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) {
     6psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) {
    77
    88    bool status = false;
    … …  
    1111
    1212    // find the peaks in the smoothed image
    13     psArray *peaks = pmPeaksInImage (significance, threshold);
     13    // NOTE : significance->variance actually carries the detection S/N image
     14    psArray *peaks = pmPeaksInImage (significance->variance, threshold);
    1415    if (peaks == NULL) {
    1516        // we only get a NULL peaks array due to a programming or config error.
    … …  
    3435    for (int i = 0; i < peaks->n; i++) {
    3536        pmPeak *peak = peaks->data[i];
    36         peak->SN = sqrt(peak->value);
    37         peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
    38         // if (peak->flux / peak->value > 5.0/12.0) {
    39         //     psWarning ("odd peak levels (1)");
    40         // }
    41         // if (peak->value / peak->flux > 5*12.0) {
    42         //     psWarning ("odd peak levels (2)");
    43         // }
     37        peak->rawFlux = readout->image->data.F32[peak->y-row0][peak->x-col0];
     38        peak->rawFluxStdev = sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
     39        peak->smoothFlux = significance->image->data.F32[peak->y-row0][peak->x-col0];
     40        peak->smoothFluxStdev = peak->smoothFlux / sqrt(significance->variance->data.F32[peak->y-row0][peak->x-col0]);
     41        // NOTE smoothFluxStdev is actually (sqrt(variance) / covar_factor)
     42
     43        // do we need this or not?
     44        // peak->SN = sqrt(peak->detValue);
     45
    4446        if (readout->variance && isfinite (peak->dx)) {
    4547            peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
    … …  
    5153
    5254    // limit the total number of returned peaks as specified
    53     psArraySort (peaks, pmPeakSortBySN);
     55    psArraySort (peaks, pmPeaksSortByRawFluxDescend);
    5456    if (nMax && (peaks->n > nMax)) {
    5557        psArray *tmpPeaks = psArrayAllocEmpty (nMax);

  • trunk/psphot/src/psphotFitSet.c

    r29004 r31154  
    2727    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
    2828    fitOptions->mode          = PM_SOURCE_FIT_EXT;
     29    fitOptions->covarFactor   = 1.0;
    2930    // XXX for now, use the defaults for the rest:
    3031    // fitOptions->nIter         = fitIter;

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r30624 r31154  
    121121    // covarFactor = 1.0;
    122122
     123    int Nsat = 0;
     124
    123125    // select the sources which will be used for the fitting analysis
    124126    for (int i = 0; i < sources->n; i++) {
    … …  
    134136        // do not include CRs in the full ensemble fit
    135137        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
     138
     139        // XXX count saturated stars
     140        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     141            Nsat ++;
     142        }
    136143
    137144        if (final) {
    … …  
    167174        if (modelSum < 0.8) {
    168175            fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n",
    169                      source->peak->xf, source->peak->yf, modelSum, source->peak->flux);
     176                     source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux);
    170177        }
    171178
    … …  
    179186    }
    180187    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
     188
     189    // fprintf (stderr, "****** Nsat : %d ********\n", Nsat);
    181190
    182191    if (fitSources->n == 0) {
    … …  
    303312    for (int i = 0; i < fitSources->n; i++) {
    304313        pmSource *source = fitSources->data[i];
    305         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
    306314        pmModel *model = pmSourceGetModel (NULL, source);
    307         pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor);
     315        if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
     316            model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
     317        }
     318        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
    308319    }
    309320    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    … …  
    324335    // We have to place this visualization here because the models are not realized until
    325336    // psphotGuessModels or fitted until psphotFitSourcesLinear.
    326     psphotVisualShowPSFStars (recipe, psf, sources);
     337    // psphotVisualShowPSFStars (recipe, psf, sources);
    327338
    328339    return true;

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r30624 r31154  
    163163    for (int i = 0; i < fitSources->n; i++) {
    164164        pmSource *source = fitSources->data[i];
    165         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
    166165        pmModel *model = pmSourceGetModel (NULL, source);
    167         pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR);
     166        if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
     167            model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
     168        }
     169        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
    168170    }
    169171    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);

  • trunk/psphot/src/psphotForced.c

    r25981 r31154  
    11# include "psphotStandAlone.h"
     2# define FORCED_PHOTOMETRY 1
    23
    34int main (int argc, char **argv) {
    … …  
    2021
    2122    // call psphot for each readout
    22     if (!psphotForcedImageLoop (config)) {
     23    if (!psphotImageLoop (config, PSPHOT_FORCED)) {
    2324        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
    2425        exit (psphotGetExitStatus());

  • trunk/psphot/src/psphotForcedArguments.c

    r25981 r31154  
    103103    }
    104104
    105     PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv );
     105    PS_ARGUMENTS_GENERIC( psphot, config, argc, argv );
    106106
    107107    // save the following additional recipe values based on command-line options
    … …  
    110110
    111111    // Number of threads is handled
    112     PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
     112    PS_ARGUMENTS_THREADS( psphot, config, argc, argv )
    113113
     114    if (psArgumentGet(argc, argv, "-help") ||
     115        psArgumentGet(argc, argv, "-h"))
     116      writeHelpInfo(argv[0], config, stdout);
     117     
    114118    // visual : interactive display mode
    115119    if ((N = psArgumentGet (argc, argv, "-visual"))) {
    … …  
    176180    }
    177181
    178     if (psArgumentGet(argc, argv, "-help") ||
    179         psArgumentGet(argc, argv, "-h"))
    180       writeHelpInfo(argv[0], config, stdout);
    181      
    182182    // the input file is a required argument; if not found, we will exit
    183183    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");

  • trunk/psphot/src/psphotForcedImageLoop.c

    r29936 r31154  
    11# include "psphotStandAlone.h"
    22
    3 # define ESCAPE(MESSAGE) { \
    4   psError(PSPHOT_ERR_DATA, false, MESSAGE); \
    5   psFree (view); \
    6   return false; \
    7 }
     3# define ESCAPE(MESSAGE) {                              \
     4        psError(PSPHOT_ERR_DATA, false, MESSAGE);       \
     5        psFree (view);                                  \
     6        return false;                                   \
     7    }
    88
    99bool psphotForcedImageLoop (pmConfig *config) {
    … …  
    3131    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
    3232
     33    // select the appropriate recipe information
     34    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     35    psAssert (recipe, "missing recipe?");
     36
     37    psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     38
    3339    // for psphot, we force data to be read at the chip level
    3440    while ((chip = pmFPAviewNextChip (view, load->fpa, 1)) != NULL) {
    … …  
    4551        // mosaic the cells of a chip into a single contiguous (trimmed) chip
    4652        if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip.");
     53
     54        // Read WCS if easy.
     55        // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell.
     56        {
     57            pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip
     58            pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL);
     59            if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) {
     60                psWarning("Unable to read WCS astrometry from header.");
     61                psErrorClear();
     62                pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL);
     63                if (hduHigh && hduHigh != hduLow &&
     64                    !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) {
     65                    psWarning("Unable to read WCS astrometry from primary header.");
     66                    psErrorClear();
     67                }
     68            }
     69        }
    4770
    4871        // try to load other supporting data (PSF, SRC, etc).
    … …  
    7699                if (readout->mask) {
    77100                    psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
    78                     if (!pmReadoutMaskNonfinite(readout, maskSat)) {
     101                    if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) {
    79102                        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
    80103                        psFree(view);
    … …  
    91114            }
    92115
     116            // drop all versions of the internal files
    93117            status = true;
    94118            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");

  • trunk/psphot/src/psphotGuessModels.c

    r30038 r31154  
    160160        pmSource *source = sources->data[i];
    161161
     162        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
     163            fprintf (stderr, "moment failure\n");
     164        }
     165
    162166        // this is used to mark sources for which the model is measured. We check later that
    163167        // all are used.
    … …  
    172176
    173177        // the guess central intensity comes from the peak:
    174         float Io = source->peak->flux;
     178        float Io = source->peak->rawFlux;
     179        if (!isfinite(Io) && source->moments) {
     180          Io = source->moments->Peak;
     181        }
    175182
    176183        // We have two options to get a guess for the object position: the position from the
    … …  
    178185        // moments
    179186
    180 # if (1)
    181         bool useMoments = true;
    182 # else
    183         bool useMoments = false;
    184         useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR);  // we only want to try if SATSTAR is set, but..
    185 # endif
    186 
    187         useMoments = (useMoments && source->moments);          // can't if there are no moments
    188         useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
    189         useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
     187        bool useMoments = pmSourcePositionUseMoments(source);
    190188
    191189        float Xo, Yo;
    … …  
    198196        }
    199197
     198        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     199            fprintf (stderr, "satstar: %f,%f vs %f,%f : %c\n",
     200                     source->moments->Mx, source->moments->My,
     201                     source->peak->xf, source->peak->yf,
     202                     (useMoments ? 'T' : 'F'));
     203        }
     204
    200205        // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
    201206        pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);

  • trunk/psphot/src/psphotImageLoop.c

    r29936 r31154  
    11# include "psphotStandAlone.h"
    22
    3 # define ESCAPE(MESSAGE) { \
    4   psError(PSPHOT_ERR_DATA, false, MESSAGE); \
    5   psFree (view); \
    6   return false; \
    7 }
     3# define ESCAPE(MESSAGE) {                              \
     4        psError(PSPHOT_ERR_DATA, false, MESSAGE);       \
     5        psFree (view);                                  \
     6        return false;                                   \
     7    }
    88
    9 bool psphotImageLoop (pmConfig *config) {
     9bool psphotImageLoop (pmConfig *config, psphotImageLoopMode mode) {
    1010
    1111    bool status;
    … …  
    9090
    9191                // Update the header
    92                 {
    93                     pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
    94                     if (hdu && hdu != lastHDU) {
    95                         psphotVersionHeaderFull(hdu->header);
    96                         lastHDU = hdu;
    97                     }
     92                pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
     93                if (hdu && hdu != lastHDU) {
     94                    psphotVersionHeaderFull(hdu->header);
     95                    lastHDU = hdu;
    9896                }
    9997
    … …  
    109107
    110108                // run the actual photometry analysis on this chip/cell/readout
    111                 if (!psphotReadout (config, view, "PSPHOT.INPUT")) {
    112                     psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
    113                     psFree (view);
    114                     return false;
    115                 }
     109                switch (mode) {
     110                  case PSPHOT_SINGLE:
     111                    if (!psphotReadout (config, view, "PSPHOT.INPUT")) {
     112                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
     113                        psFree (view);
     114                        return false;
     115                    }
     116                    break;
     117                  case PSPHOT_FORCED:
     118                    if (!psphotForcedReadout (config, view, "PSPHOT.INPUT")) {
     119                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
     120                        psFree (view);
     121                        return false;
     122                    }
     123                    break;
     124                  case PSPHOT_MAKE_PSF:
     125                    if (!psphotMakePSFReadout (config, view, "PSPHOT.INPUT")) {
     126                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
     127                        psFree (view);
     128                        return false;
     129                    }
     130                    break;
     131                }
    116132            }
    117133

  • trunk/psphot/src/psphotLoadSRCTEXT.c

    r29004 r31154  
    7878
    7979            source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    80             source->peak->flux = peakFlux;
    8180            source->peak->dx   = dPAR[PM_PAR_XPOS];
    8281            source->peak->dy   = dPAR[PM_PAR_YPOS];

  • trunk/psphot/src/psphotMagnitudes.c

    r29936 r31154  
    7676    maskVal |= markVal;
    7777
    78     pmSourceMagnitudesInit (recipe);
     78    pmSourceMagnitudesInit (config, recipe);
    7979
    8080    // the binning details are saved on the analysis metadata
    … …  
    176176        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
    177177
    178         status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     178        status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    179179        if (status && isfinite(source->apMag)) Nap ++;
    180180
    … …  
    295295    psArray *sources                = job->args->data[0];
    296296    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[1],PS_TYPE_IMAGE_MASK_DATA);
    297     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
    298297
    299298    for (int i = 0; i < sources->n; i++) {
    … …  
    309308        }
    310309
    311         status = pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
     310        status = pmSourcePixelWeight (source, model, source->maskObj, maskVal, source->apRadius);
    312311        if (!status) {
    313312          psTrace ("psphot", 3, "fail to measure pixel weight");

  • trunk/psphot/src/psphotMakeGrowthCurve.c

    r21183 r31154  
    88
    99    // set limits on the aperture magnitudes
    10     pmSourceMagnitudesInit (recipe);
     10    pmSourceMagnitudesInit (NULL, recipe);
    1111
    1212    // bit-masks to test for good/bad pixels

  • trunk/psphot/src/psphotMakePSF.c

    r25982 r31154  
    2020
    2121    // call psphot for each readout
    22     if (!psphotMakePSFImageLoop (config)) {
     22    if (!psphotImageLoop (config, PSPHOT_MAKE_PSF)) {
    2323        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
    2424        exit (psphotGetExitStatus());

  • trunk/psphot/src/psphotMakePSFArguments.c

    r25982 r31154  
    103103    }
    104104
    105     PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv );
     105    PS_ARGUMENTS_GENERIC( psphot, config, argc, argv );
    106106
    107107    // save the following additional recipe values based on command-line options
    … …  
    110110
    111111    // Number of threads is handled
    112     PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
     112    PS_ARGUMENTS_THREADS( psphot, config, argc, argv )
    113113
    114114    // visual : interactive display mode

  • trunk/psphot/src/psphotMergeSources.c

    r30624 r31154  
    11# include "psphotInternal.h"
    22
     3// Mask to apply for PSF sources : only exclude bad sources -- we will re-test for extendedness
    34#define PSF_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
    45                         PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    5                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply for PSF sources
     6                         PM_SOURCE_MODE_CR_LIMIT)
    67
    78// for now, let's store the detections on the readout->analysis for each readout
    … …  
    110111
    111112                // the supplied peak flux needs to be re-normalized
    112                 source->peak->flux = 1.0;
    113                 source->peak->value = 1.0;
     113                source->peak->rawFlux = 1.0;
     114                source->peak->smoothFlux = 1.0;
     115                source->peak->detValue = 1.0;
    114116
    115117                // drop the loaded source modelPSF
    … …  
    281283    for (int i = 0; i < sources->n; i++) {
    282284      pmSource *source = sources->data[i];
    283       source->peak->flux = source->moments->Peak;
     285      source->peak->rawFlux = source->moments->Peak;
     286      source->peak->smoothFlux = source->moments->Peak;
    284287    }
    285288
    … …  
    331334        float ypos = model->params->data.F32[PM_PAR_YPOS];
    332335
    333         pmPeak *peak = pmPeakAlloc(xpos, ypos, 1.0, PM_PEAK_LONE);
     336        pmPeak *peak = pmPeakAlloc(xpos, ypos, flux, PM_PEAK_LONE);
    334337        peak->xf = xpos;
    335338        peak->yf = ypos;
    336         peak->flux = flux; // this are being set wrong, but does it matter?
    337 
    338         if (isfinite (source->errMag) && (source->errMag > 0.0)) {
    339           peak->SN = 1.0 / source->errMag;
    340         } else {
    341           peak->SN = 0.0;
    342         }
    343339
    344340        psArrayAdd (detections->peaks, 100, peak);
    … …  
    347343
    348344    psLogMsg ("psphot", 3, "%ld PSF sources loaded", detections->peaks->n);
     345    psphotVisualShowSources (sources);
     346    psphotVisualShowPeaks (detections);
    349347
    350348    // save detections on the readout->analysis
    … …  
    354352    }
    355353    psFree (detections);
     354
    356355    return true;
    357356}
    358357
    359358// generate the detection structure for the supplied array of sources
     359// XXX this function is currently unused
    360360bool psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf) {
    361361
    … …  
    378378        float ypos = model->params->data.F32[PM_PAR_YPOS];
    379379
    380         pmPeak *peak = pmPeakAlloc(xpos, ypos, 1.0, PM_PEAK_LONE);
     380        pmPeak *peak = pmPeakAlloc(xpos, ypos, flux, PM_PEAK_LONE);
    381381        peak->xf = xpos;
    382382        peak->yf = ypos;
    383         peak->flux = flux; // this are being set wrong, but does it matter?
    384 
    385         if (isfinite (source->errMag) && (source->errMag > 0.0)) {
    386           peak->SN = 1.0 / source->errMag;
    387         } else {
    388           peak->SN = 0.0;
    389         }
     383        peak->rawFlux = flux; // this are being set wrong, but does it matter?
     384        peak->smoothFlux = flux; // this are being set wrong, but does it matter?
    390385
    391386        source->peak = peak;
    … …  
    613608
    614609      // allocate image, weight, mask for the new image for each peak
    615       pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
     610      if (sourceOut->modelPSF) {
     611        pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
     612      }
    616613
    617614      // child sources have not been subtracted in this image, but this flag may be raised if
    … …  
    679676        objectsOut->data[k] = objectOut;
    680677
    681         objectOut->SN = objectSrc->SN;
    682         objectOut->x  = objectSrc->x;
    683         objectOut->y  = objectSrc->y;
     678        objectOut->flux = objectSrc->flux;
     679        objectOut->x    = objectSrc->x;
     680        objectOut->y    = objectSrc->y;
    684681       
    685682        objectOut->sources = psArrayAlloc(objectSrc->sources->n);
    … …  
    697694            sourceOut->parent = sourceSrc;
    698695
    699             // keep the original source flags
     696            // keep the original source flags and sequence ID (if set)
     697            sourceOut->seq      = sourceSrc->seq;
    700698            sourceOut->type     = sourceSrc->type;
    701699            sourceOut->mode     = sourceSrc->mode;
    … …  
    723721
    724722            // allocate image, weight, mask for the new image for each peak
    725             pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
     723            if (sourceOut->modelPSF) {
     724              pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
     725            }
    726726
    727727            // child sources have not been subtracted in this image, but this flag may be raised if

  • trunk/psphot/src/psphotModelBackground.c

    r29936 r31154  
    5959    assert (maskVal);
    6060
     61    // mark pixel may be used for optional iteration -- is not required by all processes
     62    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
     63
     64    // apply both MASK and MARK to make background model
     65    maskVal |= markVal;
     66
    6167    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     68
     69    float dXsample = psMetadataLookupF32(&status, recipe, "BACKGROUND.XSAMPLE");
     70    if (!status) dXsample = 1.0;
     71    float dYsample = psMetadataLookupF32(&status, recipe, "BACKGROUND.YSAMPLE");
     72    if (!status) dYsample = 1.0;
    6273
    6374    // subtract this amount extra from the sky
    … …  
    7889                                 PS_STAT_ROBUST_QUARTILE |
    7990                                 PS_STAT_CLIPPED_MEAN |
    80                                  PS_STAT_FITTED_MEAN |
    81                                  PS_STAT_FITTED_MEAN_V2 |
    82                                  PS_STAT_FITTED_MEAN_V3 |
    83                                  PS_STAT_FITTED_MEAN_V4))) {
     91                                 PS_STAT_FITTED_MEAN))) {
    8492        statsOptionLocation = PS_STAT_FITTED_MEAN;
    8593    }
    … …  
    8997        statsOptionWidth = PS_STAT_SAMPLE_STDEV;
    9098    } else if (statsOptionLocation & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE)) {
    91 #if 1
    92         statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN
    93 #else
    94         statsOptionWidth = PS_STAT_FITTED_STDEV;
    95 #endif
     99        statsOptionWidth = PS_STAT_ROBUST_STDEV;
    96100    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN) {
    97101        statsOptionWidth = PS_STAT_FITTED_STDEV;
    98102    } else if (statsOptionLocation & PS_STAT_CLIPPED_MEAN) {
    99103        statsOptionWidth = PS_STAT_CLIPPED_STDEV;
    100     } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V2) {
    101         statsOptionWidth = PS_STAT_FITTED_STDEV_V2;
    102     } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V3) {
    103         statsOptionWidth = PS_STAT_FITTED_STDEV_V3;
    104     } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V4) {
    105         statsOptionWidth = PS_STAT_FITTED_STDEV_V4;
    106104    } else {
    107105        psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation);
    … …  
    132130    stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
    133131    if (!status) {
    134         if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) {
     132        if (stats->options & PS_STAT_FITTED_MEAN) {
    135133            stats->clipSigma = 1.0;
    136134        } else {
    … …  
    154152    int nFailures = 0;
    155153    psImageBackgroundInit();
     154
     155    // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than
     156    // 1 superpixel.  If we have a window of size dXsample * dYsample, then the regions run from:
     157    // (ix + 0.5) - dX to (ix + 0.5) + dX
    156158
    157159    // measure clipped median for subimages
    … …  
    162164
    163165            // convert the ruff grid cell to the equivalent fine grid cell
    164             // XXX we need to watch out for row0,col0
    165             ruffRegion = psRegionSet (ix, ix + 1, iy, iy + 1);
     166            ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample,
     167                                      iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
    166168            fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
    167169            fineRegion = psRegionForImage (image, fineRegion);
    … …  
    199201
    200202            if (gotX && gotY) {
    201                 (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
     203                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
    202204                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
    203205            } else {
    204                 (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
     206                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
    205207                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
    206208            }
    … …  
    264266    for (int iy = 0; iy < model->numRows; iy++) {
    265267        for (int ix = 0; ix < model->numCols; ix++) {
    266             if (!isnan(modelData[iy][ix])) {
     268            if (isfinite(modelData[iy][ix])) {
    267269                Value += modelData[iy][ix];
    268270                ValueStdev += modelStdevData[iy][ix];
    … …  
    280282                    if (jx   <   0) continue;
    281283                    if (jx   >= model->numCols) continue;
     284                    if (!isfinite(modelData[jy][jx])) continue;
    282285                    value += modelData[jy][jx];
    283286                    count += 1.0;
    284287                }
    285288            }
    286             if (count > 0) modelData[iy][ix] = value / count;
     289            if (count > 0) {
     290                psLogMsg ("psphot", PS_LOG_DETAIL, "patching background %d, %d: %f (%d pts)\n", ix, iy, (value / count), (int) count);
     291                modelData[iy][ix] = value / count;
     292            }
    287293        }
    288294    }
    … …  
    306312            modelStdevData[iy][ix] = ValueStdev;
    307313        }
     314    }
     315
     316    if (psTraceGetLevel("psphot") > 5) {
     317        char name[256];
     318        sprintf (name, "backfill.%02d.fits", npass);
     319        psphotSaveImage (NULL, model, name);
     320        sprintf (name, "backfist.%02d.fits", npass);
     321        psphotSaveImage (NULL, modelStdev, name);
    308322    }
    309323

  • trunk/psphot/src/psphotMosaicSubimage.c

    r21253 r31154  
    2929    inRegion = psRegionForImage (inImage, inRegion);
    3030
    31     float peak = source->peak->flux;
     31    float peak = source->peak->rawFlux;
    3232
    3333    psImage *subImage = psImageSubset (inImage, inRegion);

  • trunk/psphot/src/psphotOutput.c

    r30624 r31154  
    2929}
    3030
    31 pmReadout *psphotSelectBackground (pmConfig *config,
    32                                    const pmFPAview *view) {
     31pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view) {
    3332
    3433    bool status;
    35     pmReadout *background;
     34   
    3635
    3736    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
    3837    if (!file) return NULL;
    39     if (file->mode == PM_FPA_MODE_INTERNAL) {
    40         background = file->readout;
    41     } else {
    42         background = pmFPAviewThisReadout (view, file->fpa);
    43     }
     38
     39    pmReadout *background = READOUT_OR_INTERNAL(view, file);
    4440    return background;
    4541}
    4642
    47 pmReadout *psphotSelectBackgroundStdev (pmConfig *config,
    48                                         const pmFPAview *view) {
     43pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view) {
    4944
    5045    bool status;
    51     pmReadout *background;
    5246
    5347    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL.STDEV");
    5448    if (!file) return NULL;
    55     if (file->mode == PM_FPA_MODE_INTERNAL) {
    56         background = file->readout;
    57     } else {
    58         background = pmFPAviewThisReadout (view, file->fpa);
    59     }
     49
     50    pmReadout *background = READOUT_OR_INTERNAL(view, file);
    6051    return background;
    6152}
    … …  
    7970        // float mpeak = model ? model->params->data.F32[PM_PAR_I0] : NAN;
    8071        // bool subtracted = source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED;
    81         // fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->flux, source->pixels->data.F32[yc][xc], mcore, mpeak);
     72        // fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->rawFlux, source->pixels->data.F32[yc][xc], mcore, mpeak);
    8273
    8374        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
    … …  
    445436        if (!source->moments) continue;
    446437
    447         float Io = source->peak->flux;
     438        float Io = source->peak->rawFlux;
    448439
    449440        fprintf (f, "%f %f : %f %f :: %f\n",

  • trunk/psphot/src/psphotPetrosianAnalysis.c

    r30624 r31154  
    2929
    3030    // source analysis is done in S/N order (brightest first)
    31     sources = psArraySort (sources, pmSourceSortBySN);
     31    sources = psArraySort (sources, pmSourceSortByFlux);
    3232
    3333    // choose the sources of interest

  • trunk/psphot/src/psphotPetrosianRadialBins.c

    r28013 r31154  
    5151    psVector *binRad     = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin
    5252    psVector *binArea    = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping)
     53    psVector *binFill    = psVectorAllocEmpty(nMax, PS_TYPE_F32); // fraction of radial bin with valid pixels
    5354
    5455    psVectorInit (binSB, 0.0);
    … …  
    144145            binSB->data.F32[nOut] = value;
    145146            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
     147            binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut];
    146148
    147149            // error in the SB is the stdev per bin / sqrt (number of pixels)
    … …  
    176178        psFree(binRad);
    177179        psFree(binArea);
     180        psFree(binFill);
    178181        psFree(radMin);
    179182        psFree(radMax);
    … …  
    202205    psFree(profile->radialBins);
    203206    psFree(profile->area);
     207    psFree(profile->binFill);
    204208
    205209    // save the vectors
    206210    profile->radialBins = binRad;
    207211    profile->area       = binArea;
     212    profile->binFill    = binFill;
    208213    profile->binSB      = binSB;
    209214    profile->binSBstdev = binSBstdev;

  • trunk/psphot/src/psphotPetrosianStats.c

    r28013 r31154  
    66// generate the Petrosian radius and flux from the mean surface brightness (r_i)
    77
    8 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
     8float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
     9float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
     10float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1);
    911
    1012bool psphotPetrosianStats (pmSource *source) {
    … …  
    2527    psVector *binRad     = profile->radialBins;
    2628    psVector *area       = profile->area;
     29    psVector *binFill    = profile->binFill;
    2730
    2831    psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    … …  
    3336    psVector *meanSB      = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    3437    psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     38    psVector *apixSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    3539
    3640    float petRadius = NAN;
     41    float petRadiusErr = NAN;
    3742    float petFlux = NAN;
     43    float petFluxErr = NAN;
     44    float petArea = NAN;
     45    float petApix = NAN;
    3846
    3947    bool anyPetro = false;
    … …  
    4149    bool above = true;
    4250    float Asum = 0.0;
     51    float Psum = 0.0;
    4352    float Fsum = 0.0;
    4453    float dFsum2 = 0.0;
    … …  
    5665
    5766        float Area = area->data.F32[i];
    58         Asum += Area;
     67        Asum += Area;                   // Asum is the cumulative area interior to this bin
    5968        Fsum += binSB->data.F32[i] * Area;
    6069        dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area);
     70        Psum += Area*binFill->data.F32[i]; // Psum is the cumulative number of pixels interior to this bin
    6171
    6272        float areaInner = 0.5 * Area;
    … …  
    8797
    8898        psVectorAppend(areaSum, Asum);
     99        psVectorAppend(apixSum, Psum);
    89100        psVectorAppend(fluxSum, Fsum);
    90101        psVectorAppend(fluxSumErr2, dFsum2);
    91102        psVectorAppend(refRadius, binRad->data.F32[i]);
    92103
    93         psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
     104        psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f  %5.1f\n",
    94105                 i, refRadius->data.F32[nOut],
    95106                 binSB->data.F32[i], binSBstdev->data.F32[i],
    96107                 meanSB->data.F32[nOut], meanSBerr,
    97108                 petRatio->data.F32[nOut], petRatioErr->data.F32[nOut],
    98                  fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
     109                 fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], apixSum->data.F32[nOut], areaInner);
    99110   
    100111        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
    … …  
    104115            if (i == 0) {
    105116                // assume Fmax @ R = 0.0
    106                 petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
    107             } else {
    108                 petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     117                petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     118                petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
     119            } else {
     120                petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     121                petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
    109122            }
    110123            above = false;
    … …  
    128141    }
    129142
     143    // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!)
    130144    if (!anyPetro) {
    131145        // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
    132146        if (lowestSignificantRadius == 0) {
    133147            // assume Fmax @ R = 0.0
    134             petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     148            petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     149            petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]);
     150
    135151        } else {
    136             petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     152            int n0 = lowestSignificantRadius-1;
     153            int n1 = lowestSignificantRadius;
     154            petRadius    = InterpolateValues     (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO);
     155            petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]);
    137156        }
    138157    }
    … …  
    147166                continue;
    148167            } else {
    149                 petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     168                petFlux    = InterpolateValues     (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     169                petFluxErr = InterpolateValuesErrY (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i]));
     170                petArea    = InterpolateValues     (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius);
     171                petApix    = InterpolateValues     (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius);
    150172                break;
    151173            }
    … …  
    158180    float R50 = NAN;
    159181    float R90 = NAN;
     182    float R50err = NAN;
     183    float R90err = NAN;
    160184    bool found50 = false;
    161185    bool found90 = false;
    … …  
    167191                continue;
    168192            } else {
    169                 R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     193                R50    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     194                R50err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i]));
    170195                found50 = true;
    171196            }
    … …  
    176201                continue;
    177202            } else {
    178                 R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     203                R90    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     204                R90err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i]));
    179205                found90 = true;
    180206            }
    … …  
    188214    source->extpars->petrosianR50    = R50;
    189215    source->extpars->petrosianR90    = R90;
     216    source->extpars->petrosianFill      = petApix / petArea;
    190217   
    191218    // XXX add the errors
    192     source->extpars->petrosianRadiusErr = NAN;
    193     source->extpars->petrosianFluxErr   = NAN;
    194     source->extpars->petrosianR50Err    = NAN;
    195     source->extpars->petrosianR90Err    = NAN;
     219    source->extpars->petrosianRadiusErr = petRadiusErr;
     220    source->extpars->petrosianFluxErr   = petFluxErr;
     221    source->extpars->petrosianR50Err    = R50err;
     222    source->extpars->petrosianR90Err    = R90err;
    196223
    197224    // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
    … …  
    205232    psFree(meanSB);
    206233    psFree(areaSum);
     234    psFree(apixSum);
    207235
    208236    return true;
    … …  
    210238
    211239float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
    212     float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0);
     240    float dydx = (Y1 - Y0) / (X1 - X0);
     241    float Y = Y0 + dydx * (X - X0);
    213242    return Y;
    214243}
    215244
     245float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1) {
     246
     247    float dydx = (Y1 - Y0) / (X1 - X0);
     248    float dxdx = (X  - X0) / (X1 - X0);
     249   
     250    float dY = sqrt(PS_SQR(dX1*dydx*dxdx) + PS_SQR(dX0*dydx*(dxdx - 1.0)));
     251    return dY;
     252}
     253
     254float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1) {
     255
     256    float dxdx = (X  - X0) / (X1 - X0);
     257   
     258    float dY = sqrt(PS_SQR(dY1*dxdx) + PS_SQR(dY0*(1.0 - dxdx)));
     259    return dY;
     260}
     261

  • trunk/psphot/src/psphotRadialApertures.c

    r30624 r31154  
    3737    int Nradial = 0;
    3838
     39    // perform full non-linear fits / extended source analysis?
     40    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
     41        psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n");
     42        return true;
     43    }
     44
    3945    psTimerStart ("psphot.radial");
    4046
    … …  
    6268    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    6369    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
     70    float outerRadius = radMax->data.F32[radMax->n - 1];
    6471
    6572    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    6774    assert (maskVal);
    6875
    69     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    70     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
    71 
    7276    // S/N limit to perform full non-linear fits
    7377    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    … …  
    7579    // source analysis is done in S/N order (brightest first)
    7680    // XXX are we getting the objects out of order? does it matter?
    77     sources = psArraySort (sources, pmSourceSortBySN);
     81    sources = psArraySort (sources, pmSourceSortByFlux);
    7882
    7983    // option to limit analysis to a specific region
    … …  
    9599        // limit selection to some SN limit
    96100        assert (source->peak); // how can a source not have a peak?
    97         if (source->peak->SN < SN_LIM) continue;
     101        if (sqrt(source->peak->detValue) < SN_LIM) continue;
    98102
    99103        // limit selection by analysis region
    … …  
    112116            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    113117        }
     118
     119        // we need to change the view for the radial aperture analysis, but we want to recover exactly
     120        // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
     121        psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
     122        psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
     123        psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
     124        psRegion oldRegion    = source->region;
     125
    114126        Nradial ++;
    115127
    116128        // force source image to be a bit larger...
    117         float radius = source->peak->xf - source->pixels->col0;
    118         radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
    119         radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    120         radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    121         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
    122 
    123         if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
     129        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
     130
     131        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
    124132            psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    125133        } else {
    … …  
    127135        }
    128136
     137        pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
     138        psFree(source->maskObj);   source->maskObj   = oldMaskObj;
     139        psFree(source->modelFlux); source->modelFlux = oldModelFlux;
     140        psFree(source->psfImage);  source->psfImage  = oldPSFimage;
     141       
    129142        // re-subtract the object, leave local sky
    130143        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    135148}
    136149
    137 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
     150bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
    138151
    139152    // if we are a child source, save the results to the parent source radial aperture array
    … …  
    163176    }
    164177
    165     // center of the apertures
    166     float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
    167     float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
     178    float xCM = NAN, yCM = NAN;
     179    if (pmSourcePositionUseMoments(source)) {
     180        xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
     181        yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
     182    } else {
     183        xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
     184        yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
     185    }
    168186
    169187    // one pass through the pixels to select the valid pixels and calculate R^2
    … …  
    203221    psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    204222    psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     223    psVector *fluxStd = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    205224    psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    206225
    207226    psVectorInit (flux,    0.0);
     227    psVectorInit (fluxStd, 0.0);
    208228    psVectorInit (fluxErr, 0.0);
    209229    psVectorInit (fill,    0.0);
    … …  
    212232    for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
    213233
     234        int j = 0;
    214235        float *aRad2 = aperRadii2->data.F32;
    215         for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
     236        for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
     237        for (; j < aperRadii2->n; j++, aRad2++) {
    216238            flux->data.F32[j]    += pixFlux->data.F32[i];
     239            fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
    217240            fluxErr->data.F32[j] += pixVar->data.F32[i];
    218241            fill->data.F32[j]    += 1.0;
    219242        }
    220243    }
     244
     245    /* for each radial bin, R(i), we measure:
     246       1) the flux within that aperture: F(i) = \sum_{r_j<R_i}(F_j)
     247       2) the fractional fill factor (count of valid pixels / effective area of the aperture
     248       3) the error on the flux within that aperture
     249     */
    221250
    222251    for (int i = 0; i < flux->n; i++) {
    223252        // calculate the total flux for bin 'nOut'
    224253        float Area = M_PI*aperRadii2->data.F32[i];
    225         fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
    226         fill->data.F32[i] /= Area;
    227         psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n",
    228                  i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area);
     254
     255        int nPix = fill->data.F32[i];
     256        float SBmean = flux->data.F32[i] / nPix;
     257        float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean));
     258
     259        // flux->data.F32[i]    = SBmean * Area;
     260        fluxStd->data.F32[i] = SBstdv * Area;
     261        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
     262
     263        // fill->data.F32[i] /= Area;
     264        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n",
     265                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area);
    229266    }
    230267   
    231268    radialAper->flux = flux;
     269    radialAper->fluxStdev = fluxStd;
    232270    radialAper->fluxErr = fluxErr;
    233271    radialAper->fill = fill;
    … …  
    245283static int nPix = 0;
    246284
    247 bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
     285bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry) {
    248286
    249287    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");

  • trunk/psphot/src/psphotRadialAperturesByObject.c

    r30624 r31154  
    4343    assert (maskVal);
    4444
    45     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    46     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
    47 
    4845    // S/N limit to perform full non-linear fits
    4946    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    … …  
    6360   
    6461    // source analysis is done in S/N order (brightest first)
    65     objects = psArraySort (objects, pmPhotObjSortBySN);
     62    objects = psArraySort (objects, pmPhotObjSortByFlux);
    6663
    6764    // generate look-up arrays for readouts
    … …  
    109106            // limit selection to some SN limit
    110107            assert (source->peak); // how can a source not have a peak?
    111             if (source->peak->SN < SN_LIM) continue;
     108            if (sqrt(source->peak->detValue) < SN_LIM) continue;
    112109
    113110            int index = source->imageID;
    … …  
    147144
    148145            // force source image to be a bit larger...
    149             // float radius = source->peak->xf - source->pixels->col0;
    150             // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
    151             // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    152             // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    153146            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
    154147
    155             if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) {
     148            if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, nMatchedPSF)) {
    156149                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    157150            } else {
    158151                source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
     152                if (source->parent) {
     153                    source->parent->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
     154                }
    159155            }
    160156

  • trunk/psphot/src/psphotRadialProfile.c

    r27819 r31154  
    1414    float Rmax = 200;
    1515    float fluxMin = 0.0;
    16     float fluxMax = source->peak->flux;
     16    float fluxMax = source->peak->rawFlux;
    1717
    1818    bool RAW_RADIUS = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_RAW_RADIUS");
    … …  
    2525        return false;
    2626    }
     27    // allocate: extpars->radFlux->radii,fluxes,theta
    2728
    2829    // use the radial profiles to determine the radius of a given isophote.  this isophote
    … …  
    3334        return false;
    3435    }
     36    // allocate : extpars->radFlux->isophotalRadii (use profile->radii,fluxes)
     37
    3538
    3639    // convert the isophotal radius vs angle measurements to an elliptical contour
    … …  
    3942        return false;
    4043    }
    41  
     44    // use extpars->radFlux->isophotalRadii,theta (result in extpars->axes)
     45
    4246    // generate a single, normalized radial profile following the elliptical contours.
    4347    // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
    … …  
    4650        return false;
    4751    }
     52    // allocate extpars->ellipticalFlux->radiusElliptical,fluxElliptical (use axes to scale raw pixels)
    4853 
    4954    // generated profile in averaged bins
    … …  
    5257        return false;
    5358    }
     59    // allocate extpars->radProfile->binSB, binSBstdv, binSum, binFill, radialBins, area (small lengths)
     60    // use radiusElliptical, fluxElliptical,
    5461 
    5562    return true;

  • trunk/psphot/src/psphotReadout.c

    r30624 r31154  
    99}
    1010
     11// for now, let's store the detections on the readout->analysis for each readout
     12bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule)
     13{
     14    static int npass = 0;
     15    char filename[64];
     16
     17    return true;
     18
     19    bool status = true;
     20
     21    int num = psphotFileruleCount(config, filerule);
     22
     23    snprintf (filename, 64, "chisq.%02d.dat", npass);
     24    FILE *f = fopen (filename, "w");
     25
     26    // loop over the available readouts
     27    for (int i = 0; i < num; i++) {
     28
     29        // find the currently selected readout
     30        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     31        psAssert (file, "missing file?");
     32
     33        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     34        psAssert (readout, "missing readout?");
     35
     36        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     37        psAssert (detections, "missing detections?");
     38
     39        psArray *sources = detections->allSources;
     40        psAssert (sources, "missing sources?");
     41
     42        for (int i = 0; i < sources->n; i++) {
     43            pmSource *source = sources->data[i];
     44            if (!source) continue;
     45
     46            pmModel *model = pmSourceGetModel (NULL, source);
     47            if (!model) continue;
     48       
     49            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
     50                fprintf (f, "%f %f %f %d %d %f  1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
     51            } else {
     52                fprintf (f, "%f %f %f %d %d %f  0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
     53            }
     54        }
     55    }
     56    fclose (f);
     57    npass ++;
     58
     59    return true;
     60}
     61
    1162bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
    1263
    … …  
    52103    }
    53104
     105# if (0)
     106    // XXX test to mask outliers, not very helpful
     107    // mask the high values in the image (with MARK)
     108    if (!psphotMaskBackground (config, view, filerule)) {
     109        return psphotReadoutCleanup (config, view, filerule);
     110    }
     111
     112    // generate a background model (median, smoothed image)
     113    if (!psphotModelBackground (config, view, filerule)) {
     114        return psphotReadoutCleanup (config, view, filerule);
     115    }
     116# endif
     117
    54118    if (!psphotSubtractBackground (config, view, filerule)) {
    55119        return psphotReadoutCleanup (config, view, filerule);
    … …  
    83147
    84148    // find blended neighbors of very saturated stars (detections->newSources)
    85     if (!psphotDeblendSatstars (config, view, filerule)) {
    86         psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    87         return psphotReadoutCleanup (config, view, filerule);
    88     }
     149    // if (!psphotDeblendSatstars (config, view, filerule)) {
     150    // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
     151    // return psphotReadoutCleanup (config, view, filerule);
     152    // }
    89153
    90154    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
    … …  
    134198    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    135199    psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
     200    psphotDumpChisqs (config, view, filerule);
    136201
    137202    // identify CRs and extended sources (only unmeasured sources are measured)
    … …  
    144209    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    145210    psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources)
     211    psphotDumpChisqs (config, view, filerule);
    146212
    147213    // replace all sources
    … …  
    151217    // NOTE : apply to ALL sources (extended + psf)
    152218    psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources)
     219    psphotDumpChisqs (config, view, filerule);
    153220
    154221    // if we only do one pass, skip to extended source analysis
    … …  
    157224    // NOTE: possibly re-measure background model here with objects subtracted / or masked
    158225
    159     // add noise for subtracted objects
    160     psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources)
    161 
    162     // find fainter sources
    163     // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
    164     psphotFindDetections (config, view, filerule, false); // pass 2 (detections->peaks, detections->footprints)
    165 
    166     // remove noise for subtracted objects (ie, return to normal noise level)
    167     // NOTE: this needs to operate only on the OLD sources
    168     psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
    169 
    170     // define new sources based on only the new peaks
    171     // NOTE: new sources are saved on detections->newSources
    172     psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
    173 
    174     // set source type
    175     // NOTE: apply only to detections->newSources
    176     if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources)
    177         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    178         return psphotReadoutCleanup (config, view, filerule);
    179     }
    180 
    181     // create full input models, set the radius to fitRadius, set circular fit mask
    182     // NOTE: apply only to detections->newSources
    183     psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources)
    184 
    185     // replace all sources so fit below applies to all at once
    186     // NOTE: apply only to OLD sources (which have been subtracted)
    187     psphotReplaceAllSources (config, view, filerule); // pass 2
    188 
    189     // merge the newly selected sources into the existing list
    190     // NOTE: merge OLD and NEW
    191     // XXX check on free of sources...
    192     psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
    193 
    194     // NOTE: apply to ALL sources
    195     psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
     226    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     227    {
     228        // add noise for subtracted objects
     229        psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources)
     230
     231        // find fainter sources
     232        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
     233        psphotFindDetections (config, view, filerule, false); // pass 2 (detections->peaks, detections->footprints)
     234
     235        // remove noise for subtracted objects (ie, return to normal noise level)
     236        // NOTE: this needs to operate only on the OLD sources
     237        psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
     238
     239        // define new sources based on only the new peaks
     240        // NOTE: new sources are saved on detections->newSources
     241        psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
     242
     243        // set source type
     244        // NOTE: apply only to detections->newSources
     245        if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources)
     246            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
     247            return psphotReadoutCleanup (config, view, filerule);
     248        }
     249
     250        // create full input models, set the radius to fitRadius, set circular fit mask
     251        // NOTE: apply only to detections->newSources
     252        psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources)
     253
     254        // replace all sources so fit below applies to all at once
     255        // NOTE: apply only to OLD sources (which have been subtracted)
     256        psphotReplaceAllSources (config, view, filerule); // pass 2
     257
     258        // merge the newly selected sources into the existing list
     259        // NOTE: merge OLD and NEW
     260        // XXX check on free of sources...
     261        psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
     262
     263        // NOTE: apply to ALL sources
     264        psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
     265        psphotDumpChisqs (config, view, filerule);
     266    }
     267
     268    // NOTE: this block performs the 2nd pass low-significance EXT detection stage (smooth or rebin by NxN times PSF size)
     269    if (0) {
     270        // add noise for subtracted objects
     271        psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources)
     272
     273        // find fainter sources
     274        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
     275        psphotFindDetections (config, view, filerule, false); // pass 2 (detections->peaks, detections->footprints)
     276
     277        // remove noise for subtracted objects (ie, return to normal noise level)
     278        // NOTE: this needs to operate only on the OLD sources
     279        psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
     280
     281        // define new sources based on only the new peaks
     282        // NOTE: new sources are saved on detections->newSources
     283        psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
     284
     285        // set source type
     286        // NOTE: apply only to detections->newSources
     287        if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources)
     288            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
     289            return psphotReadoutCleanup (config, view, filerule);
     290        }
     291
     292        // create full input models, set the radius to fitRadius, set circular fit mask
     293        // NOTE: apply only to detections->newSources
     294        psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources)
     295
     296        // replace all sources so fit below applies to all at once
     297        // NOTE: apply only to OLD sources (which have been subtracted)
     298        psphotReplaceAllSources (config, view, filerule); // pass 2
     299
     300        // merge the newly selected sources into the existing list
     301        // NOTE: merge OLD and NEW
     302        // XXX check on free of sources...
     303        psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
     304
     305        // NOTE: apply to ALL sources
     306        psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
     307    }
    196308
    197309pass1finish:
    … …  
    203315    psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
    204316    psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
     317    psphotRadialApertures(config, view, filerule);
    205318
    206319finish:

  • trunk/psphot/src/psphotReplaceUnfit.c

    r30624 r31154  
    155155      // sources have not yet been subtracted in this image (but this flag may be raised)
    156156      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     157      if (!source->modelPSF) continue;
    157158
    158159      float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    … …  
    231232      bool isPSF = false;
    232233      pmModel *model = pmSourceGetModel(&isPSF, source);
     234      if (!model) continue;
     235
    233236      float radius = model->fitRadius; // save for future use below
    234237
    … …  
    236239      if (isPSF || model->isPCM) {
    237240          // the guess central intensity comes from the peak:
    238           float Io = source->peak->flux;
     241          float Io = source->peak->rawFlux;
    239242          float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    240243          float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];

  • trunk/psphot/src/psphotRoughClass.c

    r29936 r31154  
    6363    }
    6464
     65    int NstarsInClump = psMetadataLookupS32 (&status, readout->analysis, "PSF_CLUMP_NSTARS");
     66    // if NstarsInClump is not defined, use the user-selected option:
     67    if (!status) {
     68        NstarsInClump = 1000;
     69    }
     70
     71    int ScaleForClump = 1;
     72    if (NstarsInClump >= 12) ScaleForClump = 2; // 4 cells
     73    if (NstarsInClump >= 27) ScaleForClump = 3; // 9 cells
     74    if (NstarsInClump >= 48) ScaleForClump = 4; // 16 cells
     75    if (NstarsInClump >  75) ScaleForClump = 5; // 25 cells
     76
    6577    // we make this measurement on a NxM grid of regions across the readout
    6678    int NX  = psMetadataLookupS32 (&status, recipe, "PSF_CLUMP_NX");  CHECK_STATUS (status, "PSF_CLUMP_NX");
    6779    int NY  = psMetadataLookupS32 (&status, recipe, "PSF_CLUMP_NY");  CHECK_STATUS (status, "PSF_CLUMP_NY");
    68     int dX  = readout->image->numCols / NX;
    69     int dY  = readout->image->numRows / NY;
     80
     81    int NXuse, NYuse;
     82
     83    ScaleForClump = PS_MIN(ScaleForClump, PS_MAX(NX, NY));
     84    if (NX > NY) {
     85        NXuse = ScaleForClump;
     86        NYuse = (int) (ScaleForClump * (NY / NX) + 0.5);
     87    } else {
     88        NYuse = ScaleForClump;
     89        NXuse = (int) (ScaleForClump * (NY / NX) + 0.5);
     90    }
     91
     92    psLogMsg ("psphot", 4, "With %d stars, using %d x %d grid for PSF clump\n", NstarsInClump, NXuse, NYuse);
     93
     94    int dX  = readout->image->numCols / NXuse;
     95    int dY  = readout->image->numRows / NYuse;
    7096
    7197    int nRegion = 0;
    72     for (int ix = 0; ix < NX; ix ++) {
    73         for (int iy = 0; iy < NY; iy ++) {
     98    for (int ix = 0; ix < NXuse; ix ++) {
     99        for (int iy = 0; iy < NYuse; iy ++) {
    74100
    75101            psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY);
    … …  
    181207        return false;
    182208    }
    183     psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
    184     psLogMsg ("psphot", 3, "psf clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
     209    psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal);
    185210
    186211    // get basic parameters, or set defaults

  • trunk/psphot/src/psphotSavePSFStars.c

    r19869 r31154  
    1717
    1818    // examine PSF sources in S/N order (brightest first)
    19     sources = psArraySort (sources, pmSourceSortBySN);
     19    sources = psArraySort (sources, pmSourceSortByFlux);
    2020
    2121    // counters to track the size of the image and area used in a row

  • trunk/psphot/src/psphotSetThreads.c

    r30624 r31154  
    3535    psFree(task);
    3636
    37     task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 13);
     37    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 14);
    3838    task->function = &psphotExtendedSourceFits_Threaded;
    3939    psThreadTaskAdd(task);

  • trunk/psphot/src/psphotSignificanceImage.c

    r29548 r31154  
    44// (S/N)^2.  If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the
    55// smoothing kernel.
    6 psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) {
     6pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) {
    77
    88    float SIGMA_SMTH, NSIGMA_SMTH;
    … …  
    2121        minGauss = 0.5;
    2222    }
     23
     24    // NOTE: for a faint extended-source detection pass, we over-smooth by SOMETHING
    2325
    2426    // if we have already determined the PSF model, then we have a better idea how to smooth this image
    … …  
    9294    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
    9395
    94     // XXX multithread this if needed
     96    // we are going to return both the image and the weight here: the image contains the signal
     97    // while the 'weight' will contain the significance (NOTE the deviation from the usual
     98    // definition)
     99
     100    // save the smoothed significance image in the weight array
    95101    for (int j = 0; j < smooth_im->numRows; j++) {
    96102        for (int i = 0; i < smooth_im->numCols; i++) {
    97103            float value = smooth_im->data.F32[j][i];
    98104            if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)) {
    99                 smooth_im->data.F32[j][i] = 0.0;
     105                smooth_wt->data.F32[j][i] = 0.0;
    100106            } else {
    101                 float v2 = value + PS_SQR(value/1000.0);
    102                 smooth_im->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
     107                // XXX the value of 100 here (or 1000 before) must depend on the FWHM of the smoothing kernel, right??
     108                float v2 = value + PS_SQR(value/100.0);
     109                smooth_wt->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
    103110            }
    104111        }
    … …  
    112119        static int pass = 0;
    113120        sprintf (name, "snsmooth.v%d.fits", pass);
    114         psphotSaveImage (NULL, smooth_im, name);
     121        psphotSaveImage (NULL, smooth_wt, name);
    115122        pass ++;
    116123    }
    117 
    118     psFree(smooth_wt);
    119 
    120124    psImageConvolveSetThreads(oldThreads);
    121125
    122     return smooth_im;
     126    pmReadout *significanceRO = pmReadoutAlloc(NULL);
     127    significanceRO->variance = smooth_wt;   
     128    significanceRO->image = smooth_im;   
     129
     130    return significanceRO;
    123131}
    124132

  • trunk/psphot/src/psphotSourceFits.c

    r30624 r31154  
    5151    psArrayAdd (sourceSet, 16, source);
    5252
    53     psTrace ("psphot", 4, "fitting blended source at %f %f : %f\n", source->peak->xf, source->peak->yf, source->peak->flux);
     53    psTrace ("psphot", 4, "fitting blended source at %f %f : %f\n", source->peak->xf, source->peak->yf, source->peak->rawFlux);
    5454
    5555    // we need to include all blends in the fit (unless primary is saturated?)
    … …  
    6969
    7070        // XXX assume local sky is 0.0?
    71         model->params->data.F32[PM_PAR_I0] = blend->peak->flux;
     71        model->params->data.F32[PM_PAR_I0] = blend->peak->rawFlux;
    7272        model->params->data.F32[PM_PAR_XPOS] = blend->peak->xf;
    7373        model->params->data.F32[PM_PAR_YPOS] = blend->peak->yf;
    … …  
    8383        psArrayAdd (sourceSet, 16, blend);
    8484
    85         psTrace ("psphot", 5, "adding source at %f %f : %f\n", blend->peak->xf, blend->peak->yf, blend->peak->flux);
     85        psTrace ("psphot", 5, "adding source at %f %f : %f\n", blend->peak->xf, blend->peak->yf, blend->peak->rawFlux);
    8686
    8787        // free to avoid double counting model
    … …  
    100100
    101101    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    102 
    103     // correct model chisq for flux trend
    104     double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]);
    105     PSF->chisqNorm = PSF->chisq / chiTrend;
    106102
    107103    // evaluate the blend objects, subtract if good, free otherwise
    … …  
    111107
    112108        if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    113 
    114         // correct model chisq for flux trend
    115         chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[PM_PAR_I0]);
    116         model->chisqNorm = model->chisq / chiTrend;
    117109
    118110        // if this one failed, skip it
    … …  
    159151bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    160152
    161     double chiTrend;
    162153    pmSourceFitOptions options = *fitOptions;
    163154
    … …  
    182173    // clear the circular mask
    183174    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    184 
    185     // correct model chisq for flux trend
    186     chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]);
    187     PSF->chisqNorm = PSF->chisq / chiTrend;
    188175
    189176    // does the PSF model succeed?
    … …  
    224211    float radius;
    225212    bool okEXT, okDBL;
    226     float chiEXT, chiDBL;
    227     double chiTrend;
    228213    pmModel *ONE = NULL;
    229214    pmSource *tmpSrc = NULL;
    … …  
    258243    maskVal |= markVal;
    259244
     245    float chiEXT = NAN;
     246    float chiDBL = NAN;
     247
    260248    // this temporary source is used as a place-holder by the psphotEval functions below
    261249    tmpSrc = pmSourceAlloc ();
    … …  
    271259
    272260        // correct first model chisqs for flux trend
    273         chiDBL = NAN;
    274261        ONE = DBL->data[0];
    275262        if (ONE) {
    276263            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    277             chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    278             ONE->chisqNorm = ONE->chisq / chiTrend;
    279             chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
     264            chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
    280265            ONE->fitRadius = radius;
    281266        }
    … …  
    285270        if (ONE) {
    286271            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    287             chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    288             ONE->chisqNorm = ONE->chisq / chiTrend;
    289272            ONE->fitRadius = radius;
    290273        }
    … …  
    298281
    299282        okEXT = psphotEvalEXT (tmpSrc, EXT);
    300         chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
     283        chiEXT = EXT ? EXT->chisqNorm : NAN;
    301284    }
    302285

  • trunk/psphot/src/psphotSourceMatch.c

    r30624 r31154  
    8686        if (!src) NEXT1;
    8787        if (!src->peak) NEXT1;
    88         if (!finite(src->peak->xf)) NEXT1;
    89         if (!finite(src->peak->yf)) NEXT1;
     88        if (!isfinite(src->peak->xf)) NEXT1;
     89        if (!isfinite(src->peak->yf)) NEXT1;
    9090 
    9191        if (!obj) NEXT2;
    92         if (!finite(obj->x)) NEXT2;
    93         if (!finite(obj->y)) NEXT2;
     92        if (!isfinite(obj->x)) NEXT2;
     93        if (!isfinite(obj->y)) NEXT2;
    9494 
    9595        dx = src->peak->xf - obj->x;
    … …  
    229229            int col0 = readout->image->col0;
    230230
    231             // XXX the peak type is not really used in psphot
    232             // PM_PEAK_LONE is certainly not true, but irrelevant
     231            // The peak type is not used in psphot. PM_PEAK_LONE may be wrong, but irrelevant
    233232            float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)];
    234233            pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE);
    235             peak->flux = peakFlux;
    236             peak->SN = 1.0;
    237234            peak->xf = obj->x;
    238235            peak->yf = obj->y;

  • trunk/psphot/src/psphotSourcePlots.c

    r25755 r31154  
    1818
    1919    // examine PSF sources in S/N order (brightest first)
    20     sources = psArraySort (sources, pmSourceSortBySN);
     20    sources = psArraySort (sources, pmSourceSortByFlux);
    2121
    2222    // counters to track the size of the image and area used in a row

  • trunk/psphot/src/psphotSourceSize.c

    r30818 r31154  
    207207        num++;
    208208
    209         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     209        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    210210
    211211        float kMag = -2.5*log10(source->moments->KronFlux);
    … …  
    327327
    328328        // XXX can we test if psfMag is set and calculate only if needed?
    329         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     329        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    330330
    331331        // convert to Mmaj, Mmin:
    … …  
    501501        // psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    502502        // psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    503         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     503        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    504504
    505505        // clear the mask bit
    … …  
    12681268
    12691269        // Soften variances (add systematic error)
    1270         float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
     1270        float softening = options->soft * PS_SQR(source->peak->rawFlux); // Softening for variances
    12711271
    12721272        // Across the middle: y = 0

  • trunk/psphot/src/psphotSourceStats.c

    r29936 r31154  
    122122        // allocate space for moments
    123123        source->moments = pmMomentsAlloc();
     124
     125        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
     126            fprintf (stderr, "moment failure\n");
     127        }
    124128
    125129        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
    … …  
    361365
    362366    bool status = false;
    363     float BIG_RADIUS;
    364367    psScalar *scalar = NULL;
    365368
    … …  
    373376    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    374377    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     378
     379    // if no valid pixels, massive swing or very large Mrf, likely saturated source,
     380    // try a much larger box
     381    float BIG_RADIUS = 3.0*RADIUS;
     382    float BIG_SIGMA  = 3.0*SIGMA;
    375383
    376384    // maskVal is used to test for rejected pixels, and must include markVal
    … …  
    390398
    391399        // skip faint sources for moments measurement
    392         if (source->peak->SN < MIN_SN) {
     400        if (sqrt(source->peak->detValue) < MIN_SN) {
    393401            source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    394402            Nfaint++;
    … …  
    416424        }
    417425
    418         // measure basic source moments (no S/N clipping on input pixels)
    419         status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
    420         // XXX moments / aperture test:
    421         if (0) {
    422             // clear the mask bit and set the circular mask pixels
    423             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    424             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, 4.0*source->moments->Mrf, "OR", markVal);
    425            
    426             float apMag = NAN;
    427             pmSourcePhotometryAper (&apMag, NULL, NULL, NULL, source->pixels, source->variance, source->maskObj, maskVal);
    428             fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux));
    429 
    430             // clear the mask bit
    431             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    432         }           
    433         if (status) {
    434             Nmoments ++;
    435             continue;
    436         }
    437 
    438         // if no valid pixels, or massive swing, likely saturated source,
    439         // try a much larger box
    440         BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    441         psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    442         status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0, maskVal);
    443         if (status) {
    444             source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    445             Nmoments ++;
    446             continue;
    447         }
    448 
    449         source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
    450         Nfail ++;
    451         psErrorClear();
    452         continue;
     426        // XXX how can this be set: this is raised in pmSourceRoughClass, called after this function?
     427        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     428            fprintf (stderr, "satstar: %f,%f\n", source->peak->xf, source->peak->yf);
     429        }
     430
     431        if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) {
     432            // measure basic source moments (no S/N clipping on input pixels)
     433            status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
     434        } else {
     435            // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but
     436            // only slightly: pmSourceRedefinePixels uses the readout to pass the pointers to
     437            // the parent image data.  I guess the API could be simplified: we could recover
     438            // this from the source in the function
     439
     440            pmReadout tmpReadout;
     441            tmpReadout.image    = (psImage *)source->pixels->parent;
     442            tmpReadout.mask     = (psImage *)source->maskView->parent;
     443            tmpReadout.variance = (psImage *)source->variance->parent;
     444
     445            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
     446            pmSourceRedefinePixels (source, &tmpReadout, source->peak->x, source->peak->y, BIG_RADIUS + 2);
     447
     448            psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
     449            status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal);
     450            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
     451        }
     452        if (!status) {
     453            source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
     454            Nfail ++;
     455            psErrorClear();
     456            continue;
     457        }
     458        Nmoments ++;
     459
     460        // re-try big sources or not??
     461        // if (source->moments->Mrf < 2.0*SIGMA)
    453462    }
    454463
    … …  
    473482
    474483    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
    475     float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
    476484    psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
     485
     486    // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit
     487    // since we are using all sources (true?)
     488    float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
    477489
    478490    // XXX move this to a config file?
    … …  
    480492    float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
    481493    float Sout[NSIGMA];
    482 
    483     // this sorts by peak->SN
    484     sources = psArraySort (sources, pmSourceSortBySN);
     494    int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
     495
     496    // this sorts by peak->rawFlux
     497    sources = psArraySort (sources, pmSourceSortByFlux);
    485498
    486499    // loop over radii:
    … …  
    495508
    496509            // skip faint sources for moments measurement
    497             if (source->peak->SN < MIN_SN) {
     510            if (sqrt(source->peak->detValue) < MIN_SN) {
    498511                continue;
    499512            }
    … …  
    510523        // determine the PSF parameters from the source moment values
    511524        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
    512         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
     525        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    513526
    514527#if 0
    … …  
    531544
    532545        Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     546        Nout[i] = psfClump.nStars;
    533547    }
    534548
    535549    // we are looking for sigma for which Sout = 0.65 (or some other value)
     550    int Nstars = 0;
    536551    float Sigma = NAN;
    537552    float minS = Sout[0];
    … …  
    549564        if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
    550565        Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
     566        Nstars = 0.5*(Nout[i] + Nout[i+1]);
    551567    }
    552568    psAssert (isfinite(Sigma), "did we miss a case?");
    … …  
    558574    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
    559575    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
     576    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars);
    560577
    561578    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);

  • trunk/psphot/src/psphotStack.c

    r28960 r31154  
    22
    33int main (int argc, char **argv) {
     4
     5    // uncomment to turn on memory dumps (move this to an option)
     6    // psMemDumpSetState(true);
    47
    58    psTimerStart ("complete");
    … …  
    1215
    1316    psphotVersionPrint();
     17
     18    psMemDump("start");
    1419
    1520    // load input data (config and images (signal, noise, mask)

  • trunk/psphot/src/psphotStackArguments.c

    r30624 r31154  
    2222    }
    2323
    24     // -version and -dumpconfig arguments
    25     PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv );
     24    // generic arguments (version, dumpconfig)
     25    PS_ARGUMENTS_GENERIC( psphot, config, argc, argv );
     26
     27    // thread arguments
     28    PS_ARGUMENTS_THREADS( psphot, config, argc, argv )
    2629
    2730    // save the following additional recipe values based on command-line options
    … …  
    2932    psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE);
    3033
    31     // Number of threads is handled
    32     PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
    33 
    3434    // visual : interactive display mode
    3535    if ((N = psArgumentGet (argc, argv, "-visual"))) {
    3636        psArgumentRemove (N, &argc, argv);
    3737        pmVisualSetVisual(true);
     38    }
     39
     40    // memdump : enable memory spot checks
     41    if ((N = psArgumentGet (argc, argv, "-memdump"))) {
     42        psArgumentRemove (N, &argc, argv);
     43        psMemDumpSetState(true);
    3844    }
    3945

  • trunk/psphot/src/psphotStackChisqImage.c

    r30624 r31154  
    101101        for (int ix = 0; ix < inImage->numCols; ix++) {
    102102            if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
    103             chiImage->data.F32[iy][ix] += PS_SQR(inImage->data.F32[iy][ix]) / inVariance->data.F32[iy][ix];
     103            // XXX TEST chiImage->data.F32[iy][ix] += PS_SQR(inImage->data.F32[iy][ix]) / inVariance->data.F32[iy][ix];
     104            chiImage->data.F32[iy][ix] = 0.0;
    104105            chiVariance->data.F32[iy][ix] = 1.0; // ?? what is the right value?  just init to this?
    105             chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x00; // we have valid data so unmask this pixel
     106            // XXX TEST chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x00; // we have valid data so unmask this pixel
     107            chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x01; // we have valid data so unmask this pixel
    106108        }
    107109    }

  • trunk/psphot/src/psphotStackImageLoop.c

    r30624 r31154  
    11# include "psphotStandAlone.h"
    2 #define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
    32
    43# define ESCAPE(MESSAGE) {                              \
    … …  
    87    }
    98
    10 bool GetAstrometryFPA (pmConfig *config, pmFPAview *view);
    11 bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
    12 bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
     9// XXX this implementation is not smart about multi-level astrometry headers
     10bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view);
     11bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view);
    1312
    1413bool psphotStackImageLoop (pmConfig *config) {
    … …  
    1817    pmCell *cell;
    1918    pmReadout *readout;
     19
     20    psMemDump("startloop");
    2021
    2122    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
    … …  
    3233    // XXX for now, just load the full set of images up front
    3334    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
    34 
    35     bool bilevelAstrometry = GetAstrometryFPA(config, view);
    3635
    3736    // for psphot, we force data to be read at the chip level
    … …  
    5150                if (! readout->data_exists) { continue; }
    5251
     52                psMemDump("load");
     53
    5354                // PSF matching
    5455                if (!psphotStackMatchPSFs (config, view)) {
    … …  
    5758                    return false;
    5859                }
     60                psMemDump("stackmatch");
    5961
    6062                // XXX for now, we assume there is only a single chip in the PHU:
    … …  
    6466                    return false;
    6567                }
    66 
     68                psMemDump("psphot");
    6769            }
    6870            // drop all versions of the internal files
    … …  
    7880        }
    7981
    80         GetAstrometryChip(config, view, bilevelAstrometry);
     82        UpdateHeadersForChip(config, view);
    8183
    8284        // save output which is saved at the chip level
    8385        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
    8486    }
     87    psMemDump("doneloop");
    8588
    86     SetAstrometryFPA(config, view, bilevelAstrometry);
    87    
     89    UpdateHeadersForFPA(config, view);
     90
    8891    // save output which is saved at the fpa level
    8992    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
    … …  
    9396
    9497    psFree (view);
     98
     99    psMemDump("doneoutput");
    95100    return true;
    96101}
    97102
    98 /*
    99    the easiest way to implement this is to assume we can pre-load the full set of images up front.
    100    with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad.
    101 */
    102 
    103 # define UPDATE_HEADER 1
    104 
    105 bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) {
    106 
    107     bool foundAstrometry = false;
    108     bool bilevelAstrometry = false;
     103bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view) {
    109104
    110105    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
    … …  
    122117        psAssert (input, "missing input file");
    123118
    124         // find the FPA phu
    125         pmHDU *phu = pmFPAviewThisPHU(view, input->fpa);
    126         if (!phu) {
    127             psWarning("Unable to read bilevel mosaic astrometry for input FPA entry %d", i);
    128             continue;
     119        // just copy the input headers to the output headers, then update version info
     120        pmChip *inChip = pmFPAviewThisChip(view, input->fpa); ///< Chip in the input
     121        pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa);
     122
     123        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
     124        pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
     125        if (!outHDU) {
     126            pmFPAAddSourceFromView(output->fpa, view, output->format);
     127            outHDU = pmFPAviewThisHDU (view, output->fpa);
     128            psAssert (outHDU, "failed to make HDU");
    129129        }
    130 
    131         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    132         if (!ctype) {
    133             psWarning("Error in WCS keywords for input FPA entry %d", i);
    134             continue;
     130        if (!outHDU->header) {
     131            outHDU->header = psMetadataCopy(NULL, inHDU->header);
    135132        }
    136 
    137         if (!foundAstrometry) {
    138             bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    139         } else {
    140             if (bilevelAstrometry != !strcmp (&ctype[4], "-DIS")) {
    141                 psAbort("astrometry type mis-match");
    142             }
    143         }
    144 
    145         if (bilevelAstrometry) {
    146             // update the output structures
    147             if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
    148                 psWarning("Unable to read bilevel mosaic astrometry for input FPA.");
    149             }
    150         }
     133        outChip->toFPA = psMemIncrRefCounter(inChip->toFPA);
     134        outChip->fromFPA = psMemIncrRefCounter(inChip->fromFPA);
    151135    }
    152     return bilevelAstrometry;
     136    return true;
    153137}
    154138
    155 bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
     139bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view) {
    156140
    157141    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
    … …  
    169153        psAssert (input, "missing input file");
    170154
    171         // Need to update the output for astrometry.  Read WCS data from the corresponding
    172         // header and save in the output fpa
    173         pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa);
    174         pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
     155        output->fpa->toTPA = psMemIncrRefCounter(input->fpa->toTPA);
     156        output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA);
     157        output->fpa->toSky = psMemIncrRefCounter(input->fpa->toSky);
    175158
    176         pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
    177         if (!outHDU) {
    178             pmFPAAddSourceFromView(output->fpa, view, output->format);
    179             outHDU = pmFPAviewThisHDU (view, output->fpa);
    180             psAssert (outHDU, "failed to make HDU");
    181         }
    182         if (!outHDU->header) {
    183           outHDU->header = psMetadataAlloc();
    184         }
     159        pmConceptsCopyFPA(output->fpa, input->fpa, true, true);
    185160
    186         if (bilevelAstrometry) {
    187             if (!pmAstromReadBilevelChip (outChip, inHDU->header)) {
    188                 psWarning("Unable to read bilevel chip astrometry for input FPA.");
    189                 continue;
    190             }
    191             if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) {
    192                 psWarning("Unable to generate WCS header.");
    193                 continue;
    194             }
    195         } else {
    196             // we use a default FPA pixel scale of 1.0
    197             if (!pmAstromReadWCS (output->fpa, outChip, inHDU->header, 1.0)) {
    198                 psWarning("Unable to read WCS astrometry for input FPA.");
    199                 continue;
    200             }
    201             if (!pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
    202                 psWarning("Unable to generate WCS header.");
    203                 continue;
    204             }
    205         }
     161        // XXX TEST
     162        // pmFPASetFileStatus(output->fpa, true);
     163        // pmFPASetDataStatus(output->fpa, true);
     164        // pmChip *chip = output->fpa->chips->data[0];
     165        // pmCell *cell = chip->cells->data[0];
     166        // pmReadout *readout = cell->readouts->data[0];
    206167    }
    207 
    208168    return true;
    209169}
    210170
    211 bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
     171/*
     172   the easiest way to implement this is to assume we can pre-load the full set of images up front.
     173   with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad.
     174*/
    212175
    213     if (!bilevelAstrometry) return true;
    214 
    215     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
    216 
    217     // loop over the available readouts
    218     for (int i = 0; i < num; i++) {
    219 
    220         // find the currently selected readout
    221         pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
    222         psAssert (output, "missing file?");
    223 
    224         pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa);
    225         if (!PHU) {
    226             pmFPAAddSourceFromView(output->fpa, view, output->format);
    227             PHU = pmFPAviewThisPHU (view, output->fpa);
    228             psAssert (PHU, "failed to make PHU");
    229         }
    230         if (!PHU->header) {
    231           PHU->header = psMetadataAlloc();
    232         }
    233 
    234         if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) {
    235             psWarning("Unable to generate WCS header.");
    236         }
    237     }
    238 
    239     return true;
    240 }

  • trunk/psphot/src/psphotStackMatchPSFs.c

    r30624 r31154  
    9292    }
    9393
     94    // copy the header data from Src to Out
     95    // pmHDU *hduSrc = pmHDUFromReadout(readoutSrc);
     96    // psAssert (hduSrc, "input missing hdu?");
     97
     98
    9499    // set NAN pixels to 'SAT'
    95100    psImageMaskType maskSat = pmConfigMaskGet("SAT", config);

  • trunk/psphot/src/psphotStackMatchPSFsUtils.c

    r30624 r31154  
    6464        }
    6565        if (!source->peak) continue;
    66         if (source->peak->SN < SN_MIN) continue;
     66        if (sqrt(source->peak->detValue) < SN_MIN) continue;
    6767        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
    6868        numGood++;
    … …  
    8181        }
    8282        if (!source->peak) continue;
    83         if (source->peak->SN < SN_MIN) continue;
     83        if (sqrt(source->peak->detValue) < SN_MIN) continue;
    8484        float xSource, ySource;         // Coordinates of source
    8585        coordsFromSource(&xSource, &ySource, source);
    … …  
    287287
    288288    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
    289     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     289    int threads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads
    290290
    291291    int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order
    … …  
    363363
    364364        // we need to register the FWHM values for use downstream
    365         pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing);
     365        pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]);
    366366
    367367        pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);

  • trunk/psphot/src/psphotStackPSF.c

    r30624 r31154  
    5656        }
    5757
    58         float Sxx = sqrt(2.0)*targetFWHM / 2.35;
     58        // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
     59        // GAUSS  : 1.000
     60        // PGAUSS : 1.006
     61        // QGAUSS : 1.151
     62        // RGAUSS : 0.883
     63        // PS1_V1 : 1.134
     64       
     65        float scaleFactor = NAN;
     66        if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
     67            scaleFactor = 1.000;
     68        }
     69        if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
     70            scaleFactor = 1.0006;
     71        }
     72        if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
     73            scaleFactor = 1.151;
     74        }
     75        if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
     76            scaleFactor = 0.883;
     77        }
     78        if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
     79            scaleFactor = 1.134;
     80        }
     81        psAssert (isfinite(scaleFactor), "invalid model for PSF");
     82
     83        float Sxx = sqrt(2.0)*targetFWHM / 2.35 / scaleFactor;
    5984
    6085        // XXX probably should make the model type (and par 7) optional from recipe
    61         psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
     86        // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
     87        psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
    6288        if (!psf) {
    6389            psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");

  • trunk/psphot/src/psphotStackParseCamera.c

    r30624 r31154  
    172172    }
    173173
     174    // select the appropriate recipe information
     175    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     176    bool saveCnv = psMetadataLookupBool(&status, recipe, "SAVE.CNV");
     177    bool saveChisq = psMetadataLookupBool(&status, recipe, "SAVE.CHISQ");
     178
     179    // loop over the available readouts
     180    for (int i = 0; i < nInputs; i++) {
     181        pmFPAfile *file = NULL;
     182
     183        file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i);
     184        file->save = saveCnv;
     185
     186        file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.MASK", i);
     187        file->save = saveCnv;
     188
     189        file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.VARIANCE", i);
     190        file->save = saveCnv;
     191    }
     192
    174193    // generate an pmFPAimage for the chisqImage
    175194    // XXX output of these files should be optional
    … …  
    180199            return false;
    181200        }
    182         chisqImage->save = true;
     201        chisqImage->save = saveChisq;
    183202
    184203        pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK");
    … …  
    191210            return NULL;
    192211        }
    193         chisqMask->save = true;
     212        chisqMask->save = saveChisq;
    194213
    195214        pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE");
    … …  
    202221            return NULL;
    203222        }
    204         chisqVariance->save = true;
     223        chisqVariance->save = saveChisq;
    205224    }
    206225

  • trunk/psphot/src/psphotStackReadout.c

    r30624 r31154  
    127127    }
    128128
     129    psMemDump("sourcestats");
     130
    129131    // generate the objects (object unify the sources from the different images)
    130132    // XXX this could just match the detections for the chisq image, and not bother measuring the
    … …  
    132134    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
    133135
     136    psMemDump("matchsources");
     137
    134138    if (!strcasecmp (breakPt, "TEST2")) {
    135139        psFree(objects);
    … …  
    143147        return psphotReadoutCleanup (config, view, STACK_SRC);
    144148    }
     149
     150    psMemDump("sourcestats");
    145151
    146152    // find blended neighbors of very saturated stars (detections->newSources)
    … …  
    211217    }
    212218
     219    psMemDump("psfstats");
     220
    213221    psphotStackObjectsUnifyPosition (objects);
    214222
    … …  
    230238    }
    231239
     240    psMemDump("extmeas");
     241
    232242    bool smoothAgain = true;
    233243    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
    … …  
    250260        // smooth to the next FWHM, or set 'smoothAgain' to false if no more
    251261        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
     262        psMemDump("matched");
    252263    }
    253264
    … …  
    264275
    265276    // drop the references to the image pixels held by each source
    266     // psphotSourceFreePixels (config, view, STACK_OUT);
     277    psphotSourceFreePixels (config, view, STACK_OUT);
    267278    psphotSourceFreePixels (config, view, STACK_SRC);
    268279

  • trunk/psphot/src/psphotStandAlone.h

    r23487 r31154  
    1414pmConfig       *psphotArguments (int argc, char **argv);
    1515bool            psphotParseCamera (pmConfig *config);
    16 bool            psphotImageLoop (pmConfig *config);
    1716bool            psphotMosaicChip(pmConfig *config, const pmFPAview *view, char *outFile, char *inFile);
    1817void            psphotCleanup (pmConfig *config);

  • trunk/psphot/src/psphotSubtractBackground.c

    r29936 r31154  
    2424    assert (modelFile);
    2525
    26     pmReadout *model = NULL;
    27     if (modelFile->mode == PM_FPA_MODE_INTERNAL) {
    28         model = modelFile->readout;
    29     } else {
    30         model = pmFPAviewThisReadout (view, modelFile->fpa);
    31     }
     26    pmReadout *model = READOUT_OR_INTERNAL(view, modelFile);
    3227    assert (model);
    3328
    … …  
    4439    if (file) {
    4540        // we are using PSPHOT.BACKGND as an I/O file: select readout or create
    46         if (file->mode == PM_FPA_MODE_INTERNAL) {
    47             background = file->readout;
    48         } else {
    49             background = pmFPAviewThisReadout (view, file->fpa);
    50         }
     41        background = READOUT_OR_INTERNAL(view, file);
    5142        if (background == NULL) {
    5243            // readout does not yet exist: create from input

  • trunk/psphot/src/psphotTestPSF.c

    r21183 r31154  
    1717
    1818    // examine PSF sources in S/N order (brightest first)
    19     sources = psArraySort (sources, pmSourceSortBySN);
     19    sources = psArraySort (sources, pmSourceSortByFlux);
    2020
    2121    // array to store candidate PSF stars

  • trunk/psphot/src/psphotVisual.c

    r30624 r31154  
    3939    switch (channel) {
    4040      case 1:
    41         if (kapa1 == -1) {
    42             kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images");
    43             if (kapa1 == -1) {
    44                 fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    45                 pmVisualSetVisual(false);
    46             }
    47         }
    48         return kapa1;
     41        pmVisualInitWindow (&kapa1, "psphot:images");
     42        return kapa1;
    4943      case 2:
    50         if (kapa2 == -1) {
    51             kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    52             if (kapa2 == -1) {
    53                 fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    54                 pmVisualSetVisual(false);
    55             }
    56         }
     44        pmVisualInitWindow (&kapa2, "psphot:plots");
    5745        return kapa2;
    5846      case 3:
    59         if (kapa3 == -1) {
    60             kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps");
    61             if (kapa3 == -1) {
    62                 fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    63                 pmVisualSetVisual(false);
    64             }
    65         }
     47        pmVisualInitWindow (&kapa3, "psphot:stamps");
    6648        return kapa3;
    6749      default:
    … …  
    7456
    7557    int myKapa = psphotKapaChannel (channel);
     58    if (myKapa == -1) return false;
     59
    7660    if (!(strcasecmp (overlay, "all"))) {
    7761        KiiEraseOverlay (myKapa, "red");
    … …  
    218202}
    219203
    220 bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, int channel) {
     204bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, float factor, int channel) {
    221205
    222206    KiiImage image;
    … …  
    239223    strcpy (data.name, name);
    240224    strcpy (data.file, name);
    241     data.zero = stats->robustMedian - stats->robustStdev;
    242     data.range = 5*stats->robustStdev;
     225    data.zero = stats->robustMedian - factor*stats->robustStdev;
     226
     227    // XXX I we have a smoothed image, this make a much-too-tight display range
     228    data.range = 5*factor*stats->robustStdev;
    243229    data.logflux = 0;
    244230
    … …  
    283269    if (kapa == -1) return false;
    284270
     271    float factor = 1.0;
     272    if (readout->covariance) {
     273        factor = psImageCovarianceFactorForAperture(readout->covariance, 10.0);
     274    }
     275
    285276    psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    286     psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1);
    287     psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", 0);
     277    psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1.0, 1);
     278    psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", sqrt(factor), 0);
    288279
    289280    pmVisualAskUser(NULL);
    … …  
    292283
    293284bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) {
    294 
    295     pmReadout *backgnd;
    296285
    297286    if (!pmVisualTestLevel("psphot.image.backgnd", 2)) return true;
    … …  
    303292    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
    304293
    305     if (file->mode == PM_FPA_MODE_INTERNAL) {
    306         backgnd = file->readout;
    307     } else {
    308         backgnd = pmFPAviewThisReadout (view, file->fpa);
    309     }
    310 
    311     psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 2);
    312     psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", 0);
     294    pmReadout *backgnd = READOUT_OR_INTERNAL(view, file);
     295
     296    float factor = 1.0;
     297    if (readout->covariance) {
     298        factor = psImageCovarianceFactorForAperture(readout->covariance, 10.0);
     299    }
     300
     301    psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 1.0, 2);
     302    psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", sqrt(factor), 0);
    313303
    314304    pmVisualAskUser(NULL);
    … …  
    339329    psphotVisualRangeImage (kapa, lsig, "log-signif", 2, min, max);
    340330    psFree (lsig);
     331
     332    pmVisualAskUser(NULL);
     333    return true;
     334}
     335
     336// requires psphotVisualShowImage
     337bool psphotVisualShowSources (psArray *sources) {
     338
     339    int Noverlay;
     340    KiiOverlay *overlay;
     341
     342    if (!pmVisualTestLevel("psphot.objects.sources", 1)) return true;
     343
     344    int kapa = psphotKapaChannel (1);
     345    if (kapa == -1) return false;
     346
     347    // note: this uses the Ohana allocation tools:
     348    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
     349    ALLOCATE (overlay, KiiOverlay, sources->n + 2);
     350
     351    Noverlay = 0;
     352    for (int i = 0; i < sources->n; i++) {
     353
     354        pmSource *source = sources->data[i];
     355        if (!source) continue;
     356
     357        pmPeak *peak = source->peak;
     358        if (!peak) continue;
     359
     360        overlay[Noverlay].type = KII_OVERLAY_BOX;
     361        overlay[Noverlay].x = peak->xf;
     362        overlay[Noverlay].y = peak->yf;
     363        overlay[Noverlay].dx = 4.0;
     364        overlay[Noverlay].dy = 4.0;
     365        overlay[Noverlay].angle = 0.0;
     366        overlay[Noverlay].text = NULL;
     367        Noverlay ++;
     368    }
     369
     370    KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
     371    FREE (overlay);
    341372
    342373    pmVisualAskUser(NULL);
    … …  
    839870
    840871        if (source->type != type) continue;
    841         if (mode && !(source->mode & mode)) continue;
     872
     873        if (mode == PM_SOURCE_MODE_PSFSTAR) {
     874            bool keep = false;
     875            keep |= (source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR);
     876            keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
     877            if (!keep) continue;
     878        } else {
     879            if (mode && !(source->mode & mode)) continue;
     880        }
    842881
    843882        pmMoments *moments = source->moments;
    … …  
    9741013
    9751014    // examine PSF sources in S/N order (brightest first)
    976     sources = psArraySort (sources, pmSourceSortBySN);
     1015    sources = psArraySort (sources, pmSourceSortByFlux);
    9771016
    9781017    // counters to track the size of the image and area used in a row
    … …  
    11261165
    11271166    // examine PSF sources in S/N order (brightest first)
    1128     sources = psArraySort (sources, pmSourceSortBySN);
     1167    sources = psArraySort (sources, pmSourceSortByFlux);
    11291168
    11301169    // counters to track the size of the image and area used in a row
    … …  
    12271266    }
    12281267
    1229     psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 2);
     1268    psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 1.0, 2);
    12301269
    12311270    pmVisualAskUser(NULL);
    … …  
    13551394    graphdata.xmax = +30.05;
    13561395    graphdata.ymin = -0.05;
    1357     graphdata.ymax = +5.05;
     1396    graphdata.ymax = +8.05;
    13581397    KapaSetLimits (myKapa, &graphdata);
    13591398
    … …  
    14151454    graphdata.xmax = +1.51;
    14161455    graphdata.ymin = -0.05;
    1417     graphdata.ymax = +5.05;
     1456    graphdata.ymax = +8.05;
    14181457    graphdata.color = KapaColorByName ("black");
    14191458    KapaSetLimits (myKapa, &graphdata);
    … …  
    24962535    if (myKapa == -1) return false;
    24972536
     2537    float factor = 1.0;
     2538    if (readout->covariance) {
     2539        factor = psImageCovarianceFactorForAperture(readout->covariance, 10.0);
     2540    }
     2541
    24982542    if (reshow) {
    24992543        psphotVisualShowMask (myKapa, readout->mask, "mask", 2);
    2500         psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1);
    2501     }
    2502     psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 0);
     2544        psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1.0, 1);
     2545    }
     2546    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 0);
    25032547
    25042548    pmVisualAskUser(NULL);
    … …  
    26532697    graphdata.ymax = -32.0;
    26542698
    2655     FILE *f = fopen ("chisq.dat", "w");
    2656 
    26572699    // construct the plot vectors
    26582700    int n = 0;
    … …  
    26722714        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
    26732715        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
    2674 
    2675         fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]);
    2676 
    26772716        n++;
    26782717    }
    26792718    x->n = y->n = n;
    2680     fclose (f);
    26812719
    26822720    float range;
Note: See TracChangeset for help on using the changeset viewer.