IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 19, 2012, 5:24:19 PM (14 years ago)
Author:
mhuber
Message:

merging latest r34040 trunk changes to branch

Location:
branches/meh_branches/ppstack_test
Files:
35 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/meh_branches/ppstack_test

  • branches/meh_branches/ppstack_test/psphot

  • branches/meh_branches/ppstack_test/psphot/src

  • branches/meh_branches/ppstack_test/psphot/src/Makefile.am

    r33415 r34041  
    206206        psphotPetrosianStats.c         \
    207207        psphotPetrosianVisual.c        \
    208         psphotEfficiency.c
     208        psphotEfficiency.c             \
     209        psphotSetNFrames.c
    209210
    210211# re-instate these

  • branches/meh_branches/ppstack_test/psphot/src/psphot.h

    r33415 r34041  
    9999
    100100bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
     101
     102# if (HAVE_MODEL_VAR)
     103bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
     104# else
    101105bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
     106# endif
    102107
    103108bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
    … …  
    189194// used by psphotFindDetections
    190195pmReadout      *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
    191 psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
     196psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks);
    192197bool            psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
    193198psErrorCode     psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints);
    … …  
    472477
    473478bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
    474 bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
     479bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index);
    475480bool psphotKronIterate_Threaded (psThreadJob *job);
    476481
    … …  
    481486bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
    482487
     488bool psphotSetNFrames (pmConfig *config, const pmFPAview *view, const char *filerule);
     489bool psphotSetNFramesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     490
    483491#endif

  • branches/meh_branches/ppstack_test/psphot/src/psphotAddNoise.c

    r33415 r34041  
    4747
    4848    psArray *sources = detections->allSources;
    49     psAssert (sources, "missing sources?");
     49    //psAssert (sources, "missing sources?");
     50    // if no work, should just return true
     51    if (!sources) return true;
    5052
    5153    psTimerStart ("psphot.noise");

  • branches/meh_branches/ppstack_test/psphot/src/psphotArguments.c

    r31154 r34041  
    217217    pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf",      "-psflist");
    218218    pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC",        "-src",      "-srclist");
     219    pmConfigFileSetsMD (config->arguments, &argc, argv, "EXPNUM",     "-expnum",   "-expnumlist");
    219220
    220221    if (argc == 1) {

  • branches/meh_branches/ppstack_test/psphot/src/psphotDefineFiles.c

    r25983 r34041  
    8585        output->save = true;
    8686    }
     87//    // optionally save the smoothed variance model (small FITS image)
     88//    if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) {
     89//        int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
     90//        int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
     91//        pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL");
     92//        if (!output) {
     93//            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL");
     94//            return false;
     95//        }
     96//        output->save = true;
     97//    }
     98//    // optionally save the smoothed variance model's standard deviation (small FITS image)
     99//    if (psMetadataLookupBool(NULL, recipe, "SAVE.VARMDL.STDEV")) {
     100//        int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
     101//        int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
     102//        pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL.STDEV");
     103//        if (!output) {
     104//            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL.STDEV");
     105//            return false;
     106//        }
     107//        output->save = true;
     108//    }
    87109    // optionally save the PSF Model
    88110    if (psMetadataLookupBool(NULL, recipe, "SAVE.PSF")) {

  • branches/meh_branches/ppstack_test/psphot/src/psphotEfficiency.c

    r33415 r34041  
    420420
    421421    // psphotFitSourcesLinearReadout subtracts the model fits
     422# if (HAVE_MODEL_VAR)
     423    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
     424# else
    422425    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
     426# endif
    423427        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    424428        psFree(fakeSources);

  • branches/meh_branches/ppstack_test/psphot/src/psphotFindDetections.c

    r33415 r34041  
    4848    // Use the new pmFootprints approach?
    4949    const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS");
     50    const bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
    5051
    5152    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    … …  
    6465    }
    6566
     67    bool replaceSourcesForFootprints = false;
    6668    if (firstPass) {
    6769        pass = 1;
    … …  
    7072    } else {
    7173        pass = 2;
     74        replaceSourcesForFootprints = footprintUseUnsubtracted;   
    7275        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); PS_ASSERT (status, NULL);
    7376        NMAX = 0; // unlimited number of peaks in final pass: allow a limit (PEAKS_NMAX_2) ?
    … …  
    9699
    97100    // detect the peaks in the significance image
    98     detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX);
     101    int totalPeaks = 0;
     102    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks);
    99103    psMetadataAddF32  (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
    100104    if (!detections->peaks) {
    … …  
    105109        return false;
    106110    }
     111    // hard limit on number of peaks we will accept. (To avoid memory overload in psphotStack)
     112    int maxPeaks = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX_TOTAL"); PS_ASSERT (status, NULL);
     113    if (maxPeaks && (totalPeaks > maxPeaks)) {
     114        psFree (detections);
     115        psError (PSPHOT_ERR_DATA, true, "Too many peaks %d found PEAKS_NMAX_TOTAL: %d", totalPeaks, maxPeaks);
     116        return false;
     117    }
    107118
    108119    // optionally merge peaks into footprints
    109120    if (useFootprints) {
     121        if (replaceSourcesForFootprints) {
     122            psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
     123            psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
     124            psFree (significance);
     125            significance = psphotSignificanceImage (readout, recipe, maskVal);
     126        }
     127
    110128        psphotFindFootprints (detections, significance, readout, recipe, threshold, pass, maskVal);
     129
     130        if (replaceSourcesForFootprints) {
     131            psphotRemoveAllSourcesReadout (config, view, filerule, index, recipe, false);
     132        }
    111133    }
    112134

  • branches/meh_branches/ppstack_test/psphot/src/psphotFindPeaks.c

    r33415 r34041  
    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 (pmReadout *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, int *totalPeaks) {
    77
    88    bool status = false;
    … …  
    1818        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
    1919        return NULL;
     20    }
     21    // return the total number of peaks found before the nMax limit is applied
     22    if (totalPeaks) {
     23        *totalPeaks = peaks->n;
    2024    }
    2125

  • branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinear.c

    r33415 r34041  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
     14# if (HAVE_MODEL_VAR)
     15bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
     16pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
     17bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
     18# endif
     19
    1420// for now, let's store the detections on the readout->analysis for each readout
    1521bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
    … …  
    2430    assert (recipe);
    2531
     32# if (HAVE_MODEL_VAR)
     33    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
     34    if (!fitVarMode) {
     35        psError (PSPHOT_ERR_CONFIG, true, "failed to get LINEAR_FIT_VARIANCE_MODE");
     36        return false;
     37    }
     38    // MODEL_VAR requires 2 passes -- in the first, we get the rough fluxes; in the second, we
     39    // use the flux to define the model variance before fitting the objects.  Other modes only
     40    // do a single pass.
     41    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
     42# endif
     43
    2644    int num = psphotFileruleCount(config, filerule);
    2745
    … …  
    5068        psAssert (psf, "missing psf?");
    5169
    52         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
     70# if (HAVE_MODEL_VAR)
     71        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
     72# else
     73        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
     74# endif
     75        {
    5376            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    5477            return false;
    5578        }
     79
     80# if (HAVE_MODEL_VAR)
     81        // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
     82        // after we have determined the initial set of fits, then we can generate the variance image and
     83        // re-run the fit against that variance.
     84        if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) {
     85            // generate the model variance image & source pointers
     86            if (!psphotGenerateModelVariance (config, view, file, i, recipe, readout, sources)) {
     87                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
     88                return false;
     89            }
     90
     91            // replace all sources (use TMPF_SUBTRACTED as test)
     92            psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, false);
     93
     94            // rerun fit with correct fitVarMode
     95            if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) {
     96                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
     97                return false;
     98            }
     99
     100            // free the model variance image & source pointers
     101            if (!psphotFreeModelVariance (readout, sources)) {
     102                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
     103                return false;
     104            }
     105        }
     106# endif
    56107
    57108        psphotVisualShowResidualImage (readout, (num > 0));
    … …  
    62113}
    63114
    64 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) {
    65 
     115# if (HAVE_MODEL_VAR)
     116// look up the fit variance mode from the recipe; older recipes do not have the value
     117// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
     118// set the value to either CONST or IMAGE_VAR
     119pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe) {
     120
     121    bool status = false;
     122
     123    char *fitVarModeString = psMetadataLookupStr(&status, recipe, "LINEAR_FIT_VARIANCE_MODE");
     124    if (!status) {
     125        bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     126        if (!status) {
     127            psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS");
     128        }
     129        pmSourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
     130        return fitVarMode;
     131    }
     132    if (!strcasecmp(fitVarModeString, "CONSTANT") || !strcasecmp(fitVarModeString, "CONST")) {
     133        return PM_SOURCE_PHOTFIT_CONST;
     134    }
     135    if (!strcasecmp(fitVarModeString, "IMAGE") || !strcasecmp(fitVarModeString, "IMAGE_VAR")) {
     136        return PM_SOURCE_PHOTFIT_IMAGE_VAR;
     137    }
     138    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
     139        return PM_SOURCE_PHOTFIT_MODEL_VAR;
     140    }
     141    psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString);
     142    return PM_SOURCE_PHOTFIT_NONE;
     143}
     144# endif
     145
     146# if (HAVE_MODEL_VAR)
     147bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
     148# else
     149bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
     150# endif
     151{
    66152    bool status;
    67153    float x;
    … …  
    99185    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    100186
     187# if (!HAVE_MODEL_VAR)
    101188    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    102189        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    … …  
    104191        psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    105192    }
     193# endif
    106194    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    107195    if (!status) {
    … …  
    231319    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
    232320
     321# if (HAVE_MODEL_VAR)
     322    // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
     323    // XXX we have two possibilities here:
     324
     325    // 1) do 2 passes, where in the first case we use the CONST weighting, and in the second
     326    // use the fitted model values to define the model
     327
     328    // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
     329# endif
     330
    233331    // fill out the sparse matrix elements and border elements (B)
    234332    // SRCi is the current source of interest
    … …  
    238336
    239337        // diagonal elements of the sparse matrix (auto-cross-product)
     338# if (HAVE_MODEL_VAR)
     339        f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
     340# else
    240341        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
     342# endif
    241343        psSparseMatrixElement (sparse, i, i, f);
    242344
     345# if (HAVE_MODEL_VAR)
     346        // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
     347        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
     348            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
     349            errors->data.F32[i] = 1.0 / sqrt(var);
     350        } else {
     351            errors->data.F32[i] = 1.0 / sqrt(f);
     352        }
     353# else
    243354        // the formal error depends on the weighting scheme
    244355        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    … …  
    248359            errors->data.F32[i] = 1.0 / sqrt(f);
    249360        }
    250 
     361# endif
    251362
    252363        // find the image x model value
     364# if (HAVE_MODEL_VAR)
     365        f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
     366# else
    253367        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
     368# endif
    254369        psSparseVectorElement (sparse, i, f);
    255370
    … …  
    257372        switch (SKY_FIT_ORDER) {
    258373          case 1:
     374# if (HAVE_MODEL_VAR)
     375            f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
     376            psSparseBorderElementB (border, i, 1, f);
     377            f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
     378            psSparseBorderElementB (border, i, 2, f);
     379# else
    259380            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    260381            psSparseBorderElementB (border, i, 1, f);
    261382            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    262383            psSparseBorderElementB (border, i, 2, f);
     384# endif
    263385
    264386          case 0:
     387# if (HAVE_MODEL_VAR)
     388            f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
     389# else
    265390            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
     391# endif
    266392            psSparseBorderElementB (border, i, 0, f);
    267393            break;
    … …  
    283409
    284410            // got an overlap; calculate cross-product and add to output array
     411# if (HAVE_MODEL_VAR)
     412            f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
     413# else
    285414            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
     415# endif
    286416            psSparseMatrixElement (sparse, j, i, f);
    287417        }
    … …  
    321451
    322452    // set the sky, sky_x, sky_y components of border matrix
     453# if (HAVE_MODEL_VAR)
     454    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
     455# else
    323456    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
     457# endif
    324458
    325459    psSparseConstraint constraint;
    … …  
    479613    return true;
    480614}
     615
     616# if (HAVE_MODEL_VAR)
     617bool psphotModelBackgroundReadout(psImage *model,  // Model image
     618                                  psImage *modelStdev, // Model stdev image
     619                                  psMetadata *analysis, // Analysis metadata for outputs
     620                                  pmReadout *readout, // Readout for which to generate a background model
     621                                  psImageBinning *binning, // Binning parameters
     622                                  const pmConfig *config,// Configuration
     623                                  bool useVarianceImage
     624    );
     625
     626bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
     627
     628    bool status = false;
     629    psRegion fullRegion = psRegionSet (0, 0, 0, 0);
     630
     631    // bit-masks to test for good/bad pixels
     632    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     633    assert (maskVal);
     634
     635    // bit-mask to mark pixels not used in analysis
     636    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     637    assert (markVal);
     638
     639    // maskVal is used to test for rejected pixels, and must include markVal
     640    maskVal |= markVal;
     641
     642    // create a model variance image
     643    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
     644
     645    // find the binning information
     646    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
     647    assert (backBinning);
     648   
     649    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     650    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     651
     652    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
     653        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
     654        psFree (varModel);
     655        psFree (varModelStdev);
     656        return false;
     657    }
     658
     659    // linear interpolation to full-scale
     660    if (!psImageUnbin (modelVar, varModel, backBinning)) {
     661        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
     662        psFree (varModel);
     663        psFree (varModelStdev);
     664        return false;
     665    }
     666
     667    psFree (varModel);
     668    psFree (varModelStdev);
     669
     670    // XXX for a test:
     671    psphotSaveImage (NULL, modelVar, "model.bck.fits");
     672
     673    // insert all of the source models
     674    for (int i = 0; i < sources->n; i++) {
     675
     676        // source of interest
     677        pmSource *source = sources->data[i];
     678
     679        // skip sources which were not fitted already
     680        if (!(source->mode & PM_SOURCE_MODE_LINEAR_FIT)) continue;
     681
     682        // pixel region appropriate for the source
     683        psRegion region = psRegionForImage (source->pixels, fullRegion);
     684
     685        // define the source->modelVar pixels (view on modelVar image)
     686        psAssert (!source->modelVar, "programming error : modelVar should be NULL here");
     687        psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here");
     688        psAssert (source->modelFlux->data.F32, "programming error : modelFlux should not be NULL here");
     689        source->modelVar = psImageSubset(modelVar, region);
     690
     691        // add the source model to the model variance image
     692        pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
     693    }
     694
     695    // XXX for a test:
     696    psphotSaveImage (NULL, modelVar, "model.var.fits");
     697    psphotSaveImage (NULL, readout->variance, "image.var.fits");
     698
     699    // we save the model variance for future reference
     700    psMetadataAddImage(readout->analysis, PS_LIST_TAIL, "PSPHOT.MODEL.VAR", PS_META_REPLACE, "model variance", modelVar);
     701    psFree (modelVar);
     702
     703    return true;
     704}
     705
     706bool psphotFreeModelVariance (pmReadout *readout, psArray *sources) {
     707
     708    bool status = false;
     709
     710    // find the binning information
     711    psImage *modelVar = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.MODEL.VAR");
     712    assert (modelVar);
     713
     714    psMetadataRemoveKey (readout->analysis, "PSPHOT.MODEL.VAR");
     715
     716    // clear modelVar pointers for all of the source models
     717    for (int i = 0; i < sources->n; i++) {
     718
     719        // source of interest
     720        pmSource *source = sources->data[i];
     721        psFree (source->modelVar);
     722    }
     723
     724    return true;
     725}
     726# endif

  • branches/meh_branches/ppstack_test/psphot/src/psphotImageLoop.c

    r33415 r34041  
    7575        pmFPAfileActivate (config->files, false, "PSPHOT.MASK");
    7676        pmFPAfileActivate (config->files, false, "PSPHOT.VARIANCE");
     77        // do not load the exposure number file yet file until we are done with the psphot analysis
     78        pmFPAfileActivate (config->files, false, "PSPHOT.EXPNUM");
    7779        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphot.");
    7880
    7981        // re-activate files so they will be closed and freed below
    80         pmFPAfileActivate (config->files, true, NULL);
     82//        XXX: Defer this
     83//        pmFPAfileActivate (config->files, true, NULL);
    8184
    8285        // there is now only a single chip (multiple readouts?). loop over it and process
    … …  
    150153            }
    151154        }
    152         // save output which is saved at the chip level
    153         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
     155        // Defer output and closing of files until we've (possibly) done the NFrames analysis below
     156        pmFPAfileActivate (config->files, false, NULL);
     157
     158        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psphot.");
    154159    }
    155     // save output which is saved at the fpa level
     160    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot.");
     161
     162    // activate the EXPNUM image, cause it to be read in (if it exists), use it to set N Frames for the sources
     163    // and finally iterate up to trigger writing of the output files
     164    pmFPAfileActivate (config->files, true, "PSPHOT.EXPNUM");
     165    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     166        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed attempting to load EXPNUM input for Chip in psphot.");
     167        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
     168            while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
     169                if (! readout->data_exists) { continue; }
     170                if (!psphotSetNFrames (config, view, "PSPHOT.INPUT")) ESCAPE ("failed to setNFrames.");
     171            }
     172        }
     173        // now activate all files to trigger the output
     174        pmFPAfileActivate (config->files, true, NULL);
     175        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed to iterate up Chip in psphot.");
     176    }
    156177    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
    157178

  • branches/meh_branches/ppstack_test/psphot/src/psphotKronIterate.c

    r33415 r34041  
    11# include "psphotInternal.h"
    2 # ifndef ROUND
    3 # define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
    4 # endif
    5 
    6 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert);
    7 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal);
     2
     3bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
     4bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);
     5
    86
    97bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
    … …  
    119    bool status = true;
    1210
    13     // return true;
    14 
    1511    fprintf (stdout, "\n");
    1612    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Iterate ---");
    … …  
    4642        // psAssert (psf, "missing psf?");
    4743
    48         if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) {
     44        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
    4945            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    5046            return false;
    … …  
    5854bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
    5955
    60 bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
     56bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
    6157
    6258    bool status = false;
    … …  
    6965    psTimerStart ("psphot.kron");
    7066
     67
    7168    // determine the number of allowed threads
    7269    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    … …  
    8481        MIN_KRON_RADIUS = 0.25*RADIUS;
    8582    }
     83
     84    int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS");
     85    if (!status) {
     86        KRON_ITERATIONS = 1;
     87    }
     88
     89    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
     90    if (!status) {
     91        KRON_APPLY_WEIGHT = true;
     92    }
     93
     94    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
     95    if (!status) {
     96        KRON_APPLY_WINDOW = false;
     97    }
     98    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
     99    if (!status) {
     100        KRON_SMOOTH = false;
     101    }
     102    float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA");
     103    if (!status) {
     104        KRON_SMOOTH_SIGMA = 1.7;
     105    }
     106    float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA");
     107    if (!status) {
     108        KRON_SMOOTH_NSIGMA = 2;
     109    }
     110    /*
     111     *  Parameter for calculating maximum integration radius based on source's surface
     112     *  brightness
     113     *  Given minimum surface brightness SBmin and a flux the maximum radius is found from
     114     *
     115     *  SBmin = source->flux / (pi * Rmax**2)
     116     *  Rmax = sqrt (source->flux ) / sqrt (SBmin * pi)
     117     *
     118     *  Now what do we use for SBmin?
     119     *
     120     *  SBmin = ( some flux ) / (some area)
     121     *  some flux ~ flux of Ns sigma PSF source
     122     *  some area ~ K times area of a PSF
     123     * 
     124     * flux of Ns sigma source ~ Ns * SKY_STDEV * PSF_EFFECTIVE_AREA
     125     * PSF_EFFECTIVE_AREA = 4 pi sigma_PSF^2
     126     * (the 4 accounts for the fact that 1 sigma is not the total area, it is
     127     * actually a larger region).
     128     *
     129     *   SBmin = Ns * SKY_STDEV * 4 * pi * sigma_PSF^2 / (K * pi * sigma_PSF^2)
     130     *         = Ns * SKY_STDEV * 4 / K
     131     *
     132     * We combine the two parameters Ns and K, and the constant 4 into a single recipe value
     133     * KRON_SB_MIN_FACTOR with Ns = 5 and K = 2 the corresponding value
     134     * for KRON_SB_MIN_FACTOR is 5 * 4 / 2 = 10
     135     */
     136    float KRON_SB_MIN_FACTOR = psMetadataLookupF32 (&status, recipe, "KRON_SB_MIN_FACTOR");
     137    if (!status) {
     138        KRON_SB_MIN_FACTOR = 10;
     139    }
     140    float SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_DEV");
     141    float KRON_SB_MIN_DIVISOR = sqrt ( M_PI * KRON_SB_MIN_FACTOR * SKY_STDEV );
    86142
    87143    // bit-masks to test for good/bad pixels
    … …  
    114170        // set a window function for each source based on the moments
    115171        // (this skips really bad sources (no peak, no moments, DEFECT)
    116         psphotKronWindowSetSource (source, kronWindow, false, true);
     172        psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
     173    }
     174
     175    // We measure the Kron Radius on a smoothed copy of the readout image
     176    psImage *smoothedImage = NULL;
     177    if (KRON_SMOOTH) {
     178        // Build the smoothed source image
     179        // Replace the subtracted sources
     180        psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false);
     181        // Copy the image and smooth
     182        psTimerStart ("psphot.kron.smooth");
     183        smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32);
     184        psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA);
     185        psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth"));
     186
     187        // remove the sources
     188        psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false );
     189        // Now subtract smooth versions of the sources from the smoothed image
     190        psTimerStart ("psphot.kron.smooth.sources");
     191        for (int i=0; i< sources->n; i++) {
     192            pmSource *source = sources->data[i];
     193            // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage
     194            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     195                // cache copy of smoothedPixels in the source
     196                // tmpPtr is for use by a single "module" and must be null otherwise
     197                psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!");
     198
     199                psImage *smoothedPixels = psImageSubset(smoothedImage, source->region);
     200                source->tmpPtr = (psPtr) smoothedPixels;
     201                pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
     202            }
     203        }
     204        psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources"));
     205
    117206    }
    118207
    … …  
    136225            psArrayAdd(job->args, 1, kronWindow);
    137226            psArrayAdd(job->args, 1, cells->data[j]); // sources
     227            psArrayAdd(job->args, 1, smoothedImage);
    138228            PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
    139229            PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
    140230            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,             PS_TYPE_F32);
    141231            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
     232            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
     233            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
     234            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
     235            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
     236            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
    142237
    143238// set this to 0 to run without threading
    … …  
    150245            if (!psphotKronIterate_Threaded(job)) {
    151246                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    152                 psFree(AnalysisRegion);
     247                // psFree(AnalysisRegion);
    153248                return false;
    154249            }
    … …  
    174269    psFree (cellGroups);
    175270    psFree (kronWindow);
     271    if (KRON_SMOOTH) {
     272        for (int i = 0; i < sources->n; i++) {
     273            pmSource *source = sources->data[i];
     274            psFree(source->tmpPtr);
     275            source->tmpPtr = NULL;
     276        }
     277    }
     278    psFree (smoothedImage);
    176279
    177280    psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    … …  
    184287    psImage *kronWindow             = job->args->data[1];
    185288    psArray *sources                = job->args->data[2];
    186     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
    187     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    188     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
    189     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
    190 
    191     // XXX TEST : set iteration to 1
    192     for (int j = 0; j < 1; j++) {
     289    psImage *smoothedImage          = job->args->data[3];
     290    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     291    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     292    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
     293    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
     294    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
     295    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
     296    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
     297    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
     298    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
     299
     300    for (int j = 0; j < KRON_ITERATIONS; j++) {
    193301        for (int i = 0; i < sources->n; i++) {
    194302
    … …  
    196304            if (!source->peak) continue; // XXX how can we have a peak-less source?
    197305
    198             // allocate space for moments
     306            // check status of this source's moments
    199307            if (!source->moments) continue;
     308            if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
     309            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
    200310
    201311            // replace object in image
    202312            bool reSubtract = false;
     313            psImage *smoothedPixels = NULL;
    203314            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    204                 pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     315                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     316                smoothedPixels = (psImage *) source->tmpPtr;
     317                if (smoothedPixels) {
     318                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0);
     319                }
    205320                reSubtract = true;
    206321            }
    207322
    208             // use S/N to control max window size
    209             // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
    210 
    211             // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux)
    212             // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2)
    213             // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0)
    214             // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
    215 
    216             // iterate to the window radius
    217             // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
    218 
    219             // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
    220             // if we lack the skyRadius (eg MATCHED sources), go to the default value
    221             float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     323            // On first iteration set window radius to sky radius (if valid) on second iteration
     324            // use a factor times the previous radial moment value up to a maximum value that
     325            // depends on the surface brightness of the source
     326            float maxWindow;
     327            if (j == 0) {
     328                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     329            } else {
     330                if (KRON_SB_MIN_DIVISOR) {
     331                    if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) {
     332                        // Limit window radius based on surface brightness
     333                        float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;
     334
     335                        if (source->moments->Mrf > 0) {
     336                            maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
     337                        } else {
     338                            maxWindow = Rmax;
     339                        }
     340                    } else {
     341                        maxWindow = RADIUS;
     342                    }
     343                } else {
     344                    // old code
     345                    maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
     346                }
     347            }
    222348            float windowRadius = PS_MAX(RADIUS, maxWindow);
    223349
    224350            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    225             pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     351            bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    226352            psAssert (source->pixels, "WTF?");
     353            if (extend && smoothedPixels) {
     354                psFree(source->tmpPtr);
     355                smoothedPixels = psImageSubset(smoothedImage, source->region);
     356                psAssert (smoothedPixels, "WTF?");
     357                source->tmpPtr = (psPtr) smoothedPixels ;
     358            }
     359
    227360
    228361            // clear the window function for this source based on the moments
    229             psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
    230 
    231             // this function populates moments->Mrf,KronFlux,KronFluxErr
    232             psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    233             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    234 
    235             // set a window function for each source based on the moments
    236             psphotKronWindowSetSource (source, kronWindow, true, true);
     362            // Note: this function also applies cuts on the source and returns false if it
     363            // does not meet the requirements for measuring the Kron Radius or Magnitude.
     364            // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
     365            if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
     366
     367                // this function populates moments->Mrf,KronFlux,KronFluxErr
     368                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
     369                psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     370
     371                // set a window function for each source based on the moments
     372                psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
     373            }
    237374
    238375            // if we subtracted it above, re-subtract the object, leave local sky
    239376            if (reSubtract) {
    240                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     377                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     378                if (smoothedPixels) {
     379                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
     380                }
    241381            }
    242382        }
    243383    }
     384
    244385    return true;
    245386}
    246387
    247 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) {
     388bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
     389    bool applyWeight, psImage *smoothedPixels) {
    248390
    249391    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    253395
    254396    psF32 R2 = PS_SQR(radius);
    255     float rsigma2 = 0.5 / PS_SQR(radius/2.0);
     397    float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
    256398
    257399    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    266408    // Xn  = SUM (x - xc)^n * (z - sky)
    267409
     410
    268411    psF32 RF = 0.0;
    269412    psF32 RS = 0.0;
    … …  
    284427    int Ywo = source->pixels->row0;
    285428
    286     psF32 **vPix = source->pixels->data.F32;
     429    psF32 **vPix;
     430    if (smoothedPixels) {
     431        vPix = smoothedPixels->data.F32;
     432    } else {
     433        vPix = source->pixels->data.F32;
     434    }
    287435    psF32 **vWin = kronWindow->data.F32;
    288436    psF32 **vWgt = source->variance->data.F32;
    … …  
    357505    float Var = 0.0;
    358506    float Win = 0.0;
     507
     508    // set vPix to the source pixels (it may have been set to the
     509    // smoothed image above)
     510    vPix = source->pixels->data.F32;
     511
    359512
    360513    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    395548}
    396549
    397 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) {
     550bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
    398551
    399552    if (!source) return false;
    … …  
    402555    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    403556    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
     557    if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
     558    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    404559    psAssert(kronWindow, "need a window");
     560
     561    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
     562
     563    if (!applyWindow) {
     564        return true;
     565    }
    405566
    406567    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has
    … …  
    413574    float Yo = source->moments->My;
    414575
    415     float Mxx = source->moments->Mxx;
    416     float Mxy = source->moments->Mxy;
    417     float Myy = source->moments->Myy;
    418 
    419     float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    420     float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    421 
    422     // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width
    423     float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
    424 
    425     float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
    426     float Sxy = Mxy / (scale * Mmajor * Mminor);
    427     float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2
    428 
    429     float Smajor = sqrt(Mmajor);
     576    psEllipseMoments moments;
     577    moments.x2 = source->moments->Mxx;
     578    moments.y2 = source->moments->Myy;
     579    moments.xy = source->moments->Mxy;
     580
     581    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
     582    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
     583        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
     584        // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
     585        return false;
     586    }
     587
     588    // Why this factor of 0.5 ?
     589    float scale = 0.5 * source->moments->Mrf / axes.major;
     590    axes.major *= scale;
     591    axes.minor *= scale;
     592
     593    psEllipseShape shape = psEllipseAxesToShape(axes);
     594    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
     595        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
     596        // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
     597        return false;
     598    }
     599
     600    float Smajor = axes.major;
    430601
    431602    int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    … …  
    434605    int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
    435606
    436     float rMxx = 0.5 / Sxx;
    437     float rMyy = 0.5 / Syy;
    438    
     607    float rMxx = 0.5 / PS_SQR(shape.sx);
     608    float rMyy = 0.5 / PS_SQR(shape.sy);
     609    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
     610                                    // implementation. XXX: Is this correct?
     611
    439612    for (int iy = minY; iy < maxY; iy++) {
    440613        for (int ix = minX; ix < maxX; ix++) {
    … …  
    447620            float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
    448621
    449             kronWindow->data.F32[iy][ix] *= f;
     622            kronWindow->data.F32[iy][ix] *= f;
    450623        }
    451624    }
     625
    452626    return true;
    453627}

  • branches/meh_branches/ppstack_test/psphot/src/psphotMergeSources.c

    r33415 r34041  
    677677      if (sourceOut->modelPSF) {
    678678        pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
    679       }
     679      } else {
     680        // if we have no pixels we can't use it to determine the psf so make sure this bit is off
     681        sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
     682       }
    680683
    681684      // child sources have not been subtracted in this image, but this flag may be raised if
    … …  
    741744    psAssert (readoutOut, "missing readout?");
    742745
     746    pmDetections *detectionsOutOld = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS");
     747    psArray *oldFootprints = detectionsOutOld ? detectionsOutOld->footprints : NULL;
     748
    743749    // replace any existing DETECTION container on readoutOut->analysis with the new one
    744750    pmDetections *detectionsOut = pmDetectionsAlloc();
     751    if (oldFootprints) {
     752        // ... but hang on to any existing footprints so that they can be merged with new footprints in pass 2
     753        detectionsOut->footprints = psMemIncrRefCounter(oldFootprints);
     754    }
    745755    detectionsOut->allSources = psArrayAllocEmpty (100);
    746756    if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) {
    … …  
    900910            // allocate image, weight, mask for the new image for each peak
    901911            if (sourceOut->modelPSF) {
    902               pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
    903             }
     912                pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y,
     913                                                                        sourceOut->modelPSF->fitRadius);
     914            } else {
     915                // if we have no pixels we can't use it to determine the psf so make sure this bit is off
     916                sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
     917            }
    904918
    905919            // child sources have not been subtracted in this image, but this flag may be raised if

  • branches/meh_branches/ppstack_test/psphot/src/psphotModelBackground.c

    r33415 r34041  
    3737// corresponding to the model.  Other information about the background model is saved on the
    3838// readout->analysis
    39 static bool psphotModelBackgroundReadout(psImage *model,  // Model image
    40                                   psImage *modelStdev, // Model stdev image
    41                                   psMetadata *analysis, // Analysis metadata for outputs
    42                                   pmReadout *readout, // Readout for which to generate a background model
    43                                   psImageBinning *binning, // Binning parameters
    44                                   const pmConfig *config // Configuration
     39bool psphotModelBackgroundReadout(psImage *model,  // Model image
     40                                  psImage *modelStdev, // Model stdev image
     41                                  psMetadata *analysis, // Analysis metadata for outputs
     42                                  pmReadout *readout, // Readout for which to generate a background model
     43                                  psImageBinning *binning, // Binning parameters
     44                                  const pmConfig *config,// Configuration
     45                                  bool useVarianceImage
    4546    )
    4647{
    … …  
    4950    bool status = true;
    5051
    51     psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout
     52    psImage *image = useVarianceImage ? readout->variance : readout->image;
     53    psImage *mask = readout->mask; // Image and mask for readout
    5254
    5355    // select the appropriate recipe information
    … …  
    143145
    144146    // we save the binning structure for use in psphotMagnitudes
    145     psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     147    if (!useVarianceImage) {
     148        psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     149    }
    146150
    147151    psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    159163
    160164    // measure clipped median for subimages
    161     psRegion ruffRegion = {0,0,0,0};
    162     psRegion fineRegion = {0,0,0,0};
     165    psRegion ruffRegion = psRegionSet (0,0,0,0);
     166    psRegion fineRegion = psRegionSet (0,0,0,0);
    163167    for (int iy = 0; iy < model->numRows; iy++) {
    164168        for (int ix = 0; ix < model->numCols; ix++) {
    … …  
    330334    psVectorStats (statsDQ, dQ, NULL, NULL, 0);
    331335
    332     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN",  PS_META_REPLACE, "sky mean", Value);
    333     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
    334     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ",    PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian);
    335     psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
    336 
    337     // measure image and background stats and save for later output
    338     psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
    339                                       PS_STAT_SAMPLE_STDEV |
    340                                       PS_STAT_MIN |
    341                                       PS_STAT_MAX);
    342     psImageStats (statsBck, model, NULL, 0);
    343     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN",  PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    344     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
    345     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev",               ValueStdev);
    346     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ",  PS_META_REPLACE, "sky quartile slope",      statsDQ->sampleMedian);
    347     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
    348     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
    349     psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)",      model->numCols);
    350     psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)",      model->numRows);
    351     psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
    352               statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
     336    if (!useVarianceImage) {
     337        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN",  PS_META_REPLACE, "sky mean", Value);
     338        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
     339        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ",    PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian);
     340        psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
     341
     342        // measure image and background stats and save for later output
     343        psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
     344                                          PS_STAT_SAMPLE_STDEV |
     345                                          PS_STAT_MIN |
     346                                          PS_STAT_MAX);
     347        psImageStats (statsBck, model, NULL, 0);
     348        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN",  PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     349        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
     350        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev",               ValueStdev);
     351        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ",  PS_META_REPLACE, "sky quartile slope",      statsDQ->sampleMedian);
     352        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
     353        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
     354        psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)",      model->numCols);
     355        psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)",      model->numRows);
     356        psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
     357                  statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
     358        psFree(statsBck);
     359    } else {
     360        psLogMsg ("psphot", PS_LOG_INFO, "variance data : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
     361    }
    353362
    354363    psFree(statsDQ);
    … …  
    356365
    357366    psFree(stats);
    358     psFree(statsBck);
    359367    psFree(statsDefaults);
    360368    psFree(binning);
    … …  
    375383    psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation
    376384
    377     if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) {
     385    if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config, false)) {
    378386        psFree(model);
    379387        psFree(modelStdev);
    … …  
    400408    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index);
    401409
    402     if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) {
     410    if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) {
    403411        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    404412        return false;

  • branches/meh_branches/ppstack_test/psphot/src/psphotParseCamera.c

    r29936 r34041  
    22
    33// define the needed / desired I/O files
     4
     5
    46bool psphotParseCamera (pmConfig *config) {
    57
    … …  
    4042    }
    4143
     44    pmFPAfileBindFromArgs (&status, input, config, "PSPHOT.EXPNUM", "EXPNUM");
     45    if (!status) {
     46        psError (PS_ERR_UNKNOWN, false, "failed to load find definition");
     47        return NULL;
     48    }
     49
     50
    4251    // define the additional input/output files associated with psphot
    4352    if (!psphotDefineFiles (config, input)) {

  • branches/meh_branches/ppstack_test/psphot/src/psphotRadialProfileWings.c

    r33415 r34041  
    5555static float MIN_RADIUS = NAN;
    5656static float SKY_STDEV  = NAN;
     57static float SKY_SLOPE_MIN = NAN;
    5758// static FILE *file = NULL;
    5859
    … …  
    8182
    8283    MIN_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
    83     if (!status) {
     84    if (!status || (MIN_RADIUS > MAX_RADIUS)) {
    8485        MIN_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    8586    }
    … …  
    8990    if (!status) {
    9091        SKY_STDEV = 1.0; // a crude default value (why would this not exist?)
     92    }
     93
     94    // SKY_SLOPE_MIN is the sigma of the sky model (ie, smoothed on large scales)
     95    SKY_SLOPE_MIN = psMetadataLookupF32 (&status, recipe, "SKY_SLOPE_MIN");
     96    if (!status) {
     97        SKY_SLOPE_MIN = 3.0;
    9198    }
    9299
    … …  
    302309            limitSlope = slope;
    303310        }
    304         if (!limit && isfinite(slope) && (fabs(slope) < 3.0)) {
     311        if (!limit && isfinite(slope) && (fabs(slope) < SKY_SLOPE_MIN)) {
    305312            // SB no longer changing.       
    306313            limit = true;
    307314            // linearly interpolate to the radius at which we hit the sky, using the last flux and the limiting slope
    308315            if (isfinite(lastFlux)) {
    309                 limitRadius = lastRadius + lastFlux / 3.0;
     316                float interpolatedRadius = lastRadius + lastFlux / SKY_SLOPE_MIN;
     317                if (interpolatedRadius < MAX_RADIUS) {
     318                    limitRadius = interpolatedRadius;
     319                } else {
     320                    // XXX should we keep going in this case?
     321                    limitRadius = meanRadius;
     322                }
    310323            } else {
    311324                limitRadius = meanRadius;

  • branches/meh_branches/ppstack_test/psphot/src/psphotReadout.c

    r33415 r34041  
    185185        // remove noise for subtracted objects (ie, return to normal noise level)
    186186        // NOTE: this needs to operate only on the OLD sources
    187         psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
     187        bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
     188        // Note: if footprintUseUnsubtracted is true the noise was already subtracted in psphotFindDetections()
     189        if (!footprintUseUnsubtracted) {
     190            psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
     191        }
    188192
    189193        // define new sources based on only the new peaks & measure moments

  • branches/meh_branches/ppstack_test/psphot/src/psphotReadoutKnownSources.c

    r33415 r34041  
    7171    psphotSourceFreePixels (config, view, filerule);
    7272
     73    // set source->nFrames
     74    psphotSetNFrames (config, view, filerule);
     75
    7376    // create the exported-metadata and free local data
    7477    return psphotReadoutCleanup (config, view, filerule);

  • branches/meh_branches/ppstack_test/psphot/src/psphotReadoutMinimal.c

    r33415 r34041  
    6464    // Construct an initial model for each object, set the radius to fitRadius, set circular
    6565    // fit mask.  NOTE: only applied to sources without guess models
    66     psphotGuessModels (config, view, filerule);
     66    if (!psphotGuessModels (config, view, filerule)) {
     67        psLogMsg ("psphot", 3, "failure to Guess Model");
     68        return psphotReadoutCleanup (config, view, filerule);
     69    }
    6770
    6871    // linear PSF fit to source peaks

  • branches/meh_branches/ppstack_test/psphot/src/psphotReplaceUnfit.c

    r33415 r34041  
    6060
    6161    psArray *sources = detections->allSources;
    62     psAssert (sources, "missing sources?");
     62    //psAssert (sources, "missing sources?");
     63    if (!sources) return true;
    6364
    6465    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    139140
    140141    psArray *sources = detections->allSources;
    141     psAssert (sources, "missing sources?");
    142 
     142    //psAssert (sources, "missing sources?");
     143    if (!sources) return true;
     144   
    143145    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    144146    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels

  • branches/meh_branches/ppstack_test/psphot/src/psphotSersicModelClass.c

    r33415 r34041  
    546546    // set the normalization by linear fit between model and data
    547547    psphotSersicModelNorm (pcm, source);
     548    if (!isfinite(PAR[PM_PAR_I0])) {
     549        fprintf(stderr, "psphotSersicModelClassGuessPCM: psphotSerisicModelNorm set PM_PAR_I0 to NAN\n");
     550        return false;
     551    }
    548552
    549553    if (TIMING) { t5 = psTimerMark ("SersicGuess"); }

  • branches/meh_branches/ppstack_test/psphot/src/psphotSetThreads.c

    r33415 r34041  
    3030    psFree(task);
    3131
    32     task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 7);
     32    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 13);
    3333    task->function = &psphotKronIterate_Threaded;
    3434    psThreadTaskAdd(task);

  • branches/meh_branches/ppstack_test/psphot/src/psphotSourceMatch.c

    r33415 r34041  
    33bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
    44bool psphotMatchSourcesSetIDs (psArray *objects);
     5
     6static psArray     *psphotMatchFootprintCacheAlloc (int nImages);
     7static pmFootprint *psphotMatchLookupFootprint (psArray *cache, int footprintID, int imageID);
     8static pmFootprint *psphotMatchCopyFootprint (psArray *cache, pmFootprint *footprint, int imageID, psImage *image, psArray *footprints);
    59 
    610psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule)
    … …  
    193197    // vector to track if source for an image is found
    194198    psVector *found = psVectorAlloc(nImages, PS_TYPE_U8);
     199    psArray *footprintCache = psphotMatchFootprintCacheAlloc(nImages);
    195200
    196201    for (int i = 0; i < objects->n; i++) {
    … …  
    219224
    220225        // we make a copy of the largest footprint; this will be used for all new sources associated with this object
    221         pmFootprint *footprint = NULL;
     226        pmFootprint *largestFootprint = NULL;
    222227        if (iSpansMax != -1) { // copy the footprint info
    223228            pmSource *src = obj->sources->data[iSpansMax];
    … …  
    226231            psAssert(src->peak->footprint->nspans == nSpansMax, "wrong footprint?");
    227232           
    228             // we only care about the spans, do not worry about the image of this footprint
    229             footprint = pmFootprintCopyData(src->peak->footprint, NULL);
     233            largestFootprint = src->peak->footprint;
    230234        }
    231235
    … …  
    248252           
    249253            // assign to a footprint on this readout->image
    250             if (footprint) {
    251                 peak->footprint = pmFootprintCopyData(footprint, readout->image);
    252 
    253                 // the peak does not claim ownership of the footprint (it does not free it). save a copy of this
    254                 // footprint on detections->footprints so we can free it later
    255                 psArrayAdd(detections->footprints, 100, peak->footprint);
    256                 psFree (peak->footprint);
     254            if (largestFootprint) {
     255                // we save the copies that we make of the of the footprints in a hash so that we can reuse them
     256                // for all sources that share the fooprint. Without this we had serious memory explosion in
     257                // dense fields with lots of footprints (spans are small but when you have enough of them ...)
     258                peak->footprint = psphotMatchLookupFootprint(footprintCache, largestFootprint->id, index);
     259                if (!peak->footprint) {
     260                    // the peak does not claim ownership of the footprint (it does not free it).
     261                    // psphotMatchCopyFootprint saves a copy of this
     262                    // footprint on detections->footprints so we can free it later
     263                    peak->footprint = psphotMatchCopyFootprint(footprintCache, largestFootprint,
     264                                            index, readout->image, detections->footprints);
     265                }
    257266            }
    258267           
    … …  
    272281            pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    273282
     283#if (0)
     284            fprintf(stderr, "Add mising source for obj: %5d %5d image: %d flux: %f size: %4d %4d\n",
     285                                                      i, obj->id, index, peakFlux, source->pixels->numRows, source->pixels->numCols);
     286#endif
     287
    274288            peak->assigned = true;
    275289            pmPhotObjAddSource(obj, source);
    … …  
    277291            psFree (source);
    278292        }
    279         psFree (footprint);
    280     }
     293    }
     294    psFree(footprintCache);
    281295
    282296    // how many sources do we have now?
    … …  
    309323    return true;
    310324}
     325
     326// Cache of footprints created for unmatched sources.
     327static psArray * psphotMatchFootprintCacheAlloc (int nImages) {
     328    psArray *cache = psArrayAlloc(nImages);
     329    for (int i = 0; i < nImages; i++) {
     330        psHash *hash = psHashAlloc(5000);
     331        cache->data[i] = hash;
     332    }
     333    return cache;
     334}
     335
     336// Find copy of footprint with given ID made for a given image
     337static pmFootprint *psphotMatchLookupFootprint (psArray *cache, int footprintID, int imageID) {
     338    psHash *hash = (psHash *) cache->data[imageID];
     339
     340    psAssert(hash != NULL, "missing hash for image %d", imageID);
     341
     342    // footprintID is the id of the original footprint.
     343    char key[32];
     344    sprintf(key, "%d", footprintID);
     345
     346    // The footprints in our hashes are the copies of that footprint for the respective images
     347    pmFootprint *copy = psHashLookup(hash, key);
     348
     349    return copy;
     350}
     351
     352// Create a copy of a given footprint for a given image
     353static pmFootprint *psphotMatchCopyFootprint (psArray *cache, pmFootprint *footprint, int imageID, psImage *image, psArray *footprints) {
     354
     355    psAssert((imageID >= 0 && imageID < cache->n), "invalid imageID %d", imageID);
     356
     357    psHash *hash = (psHash *) cache->data[imageID];
     358
     359    psAssert(hash != NULL, "missing hash for image %d", imageID);
     360
     361    char key[32];
     362    sprintf(key, "%d", footprint->id);
     363
     364    pmFootprint *copy = pmFootprintCopyData (footprint, image);
     365
     366    // save a copy in this image's hash
     367    psHashAdd(hash, key, copy);
     368
     369    // save a copy of this footprint on the passed in footprints Array so we can free it later
     370    psArrayAdd(footprints, 100, copy);
     371
     372    // drop our ref
     373    psFree(copy);
     374
     375    return copy;
     376}

  • branches/meh_branches/ppstack_test/psphot/src/psphotSourceStats.c

    r33415 r34041  
    605605    }
    606606    psAssert (isfinite(Sigma), "did we miss a case?");
     607    if (!isfinite(minKronRadius)) {
     608        minKronRadius = Sigma;
     609    }
    607610
    608611    // choose a grid scale that is a fixed fraction of the psf sigma^2

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackArguments.c

    r31154 r34041  
    5353        psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]);
    5454        psArgumentRemove (N, &argc, argv);
     55    }
     56    if ((N = psArgumentGet (argc, argv, "-ds9regions"))) {
     57        psArgumentRemove (N, &argc, argv);
     58        pmSubtractionRegions(true);
    5559    }
    5660

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackImageLoop.c

    r33415 r34041  
    3030    }
    3131
     32    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     33    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
     34    if (useRaw && inputRaw == NULL) {
     35        psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input.");
     36        useRaw = false;
     37    }
     38
    3239    pmFPAview *view = pmFPAviewAlloc (0);
    3340
    34     // XXX for now, just load the full set of images up front
     41
     42    // XXX for now, just load the full set of images up front except for EXPNUM which we defer
     43    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
     44    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
    3545    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
    3646
    … …  
    8696        UpdateHeadersForChip(config, view);
    8797
    88         // save output which is saved at the chip level
    89         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
     98        // Defer output until we have performed psphotSetNFrames()
     99        pmFPAfileActivate (config->files, false, NULL);
     100
     101        // Iterate up
     102        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psphot.");
    90103    }
    91104    psMemDump("doneloop");
    … …  
    93106    UpdateHeadersForFPA(config, view);
    94107
    95     // save output which is saved at the fpa level
     108    // Iterate up output which is saved at the fpa level
     109    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput pmFPAfileIOChecks FPA in psphot.");
     110
     111    // Load the appropriate EXPNUM image
     112    pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV");
     113
     114    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot.");
     115
     116    // for psphot, we force data to be read at the chip level
     117    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     118        psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     119        if (! chip->process || ! chip->file_exists) { continue; }
     120        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip EXPNUM in psphotStack.");
     121
     122        // there is now only a single chip (multiple readouts?). loop over it and process
     123        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
     124            psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     125            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack.");
     126
     127            // process each of the readouts
     128            while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
     129                psLogMsg ("psphot", 6, "Readout %d: %x %x\n", view->readout, cell->file_exists, cell->process);
     130                if (! readout->data_exists) { continue; }
     131
     132                if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames.");
     133            }
     134        }
     135        // now activate all files to trigger final output
     136        pmFPAfileActivate (config->files, true, NULL);
     137
     138        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
     139    }
    96140    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
    97141

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsPrepare.c

    r33415 r34041  
    3232        pmDetections *detections = psMetadataLookupPtr(NULL, readout->analysis, "PSPHOT.DETECTIONS"); // Sources
    3333        if (!detections || !detections->allSources) {
    34             psWarning("No detections found for image %d --- rejecting.", index);
    35             options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0x01;
    36             return true;  // XXX not an error: we continue processing other images
     34            // An input must have sources. If we got here without any the input file was likely corrupt
     35            psError(PSPHOT_ERR_IO, true, "Input %d has no sources.\n\tCheck all instances of: %s", index, inputSrc->origname);
     36            return false;
    3737        }
    3838        psAssert (detections->allSources, "missing sources?");

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsUtils.c

    r33415 r34041  
    105105    psAssert(psphotRecipe, "Need PSPHOT recipe");
    106106
     107#if (0)
     108    // DON'T OVERWRITE PSPHOT recipe value
    107109    psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad
    108110    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    109111
    110112    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    111     psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
     113    psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
     114#endif
    112115
    113116    psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackObjects.c

    r33415 r34041  
    119119            }
    120120            if (!skipSourcePetro) keepObjectPetro = true;
    121 
    122             keepObjectPetro = true;
    123121        }
    124122

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackParseCamera.c

    r31154 r34041  
    5959                }
    6060            }
     61            psString expnum = psMetadataLookupStr(&status, input, "RAW:EXPNUM"); // Name of expnum image
     62            if (expnum && strlen(expnum) > 0) {
     63                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.EXPNUM.RAW", expnum, PM_FPA_FILE_EXPNUM)) {
     64                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from expnum %d (%s)", i, expnum);
     65                    return false;
     66                }
     67            }
    6168            nRaw ++;
    6269        }
    … …  
    8491                }
    8592            }
     93            psString expnum = psMetadataLookupStr(&status, input, "CNV:EXPNUM"); // Name of EXPNUM image
     94            if (expnum && strlen(expnum) > 0) {
     95                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.EXPNUM.CNV", expnum, PM_FPA_FILE_EXPNUM)) {
     96                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from expnum %d (%s)", i, expnum);
     97                    return false;
     98                }
     99            }
    86100            psString psf = psMetadataLookupStr(&status, input, "CNV:PSF"); // Name of mask
    87101            if (psf && strlen(psf) > 0) {
    … …  
    115129        }
    116130
     131        psS64 stack_id = psMetadataLookupS64(&status, input, "STACK_ID");
     132        if (!status) {
     133            psError(PS_ERR_UNKNOWN, false, "Unable to find STACK_ID from sources %d", i);
     134            return false;
     135        }
    117136        // generate an pmFPAimage for the output convolved image
    118137        // XXX output of these files should be optional
    … …  
    125144            }
    126145            outputImage->save = true;
    127             outputImage->fileID = i;            // this is used to generate output names
     146            outputImage->fileID = stack_id;             // this is used to generate output names
    128147
    129148            pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK");
    … …  
    137156            }
    138157            outputMask->save = true;
    139             outputMask->fileID = i;             // this is used to generate output names
     158            outputMask->fileID = stack_id;              // this is used to generate output names
    140159
    141160            pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE");
    … …  
    149168            }
    150169            outputVariance->save = true;
    151             outputVariance->fileID = i;         // this is used to generate output names
     170            outputVariance->fileID = stack_id;          // this is used to generate output names
    152171
    153172            // the output sources are carried on the outputImage->fpa structures
    … …  
    158177            }
    159178            outsources->save = true;
    160             outsources->fileID = i;             // this is used to generate output names
     179            outsources->fileID = stack_id;              // this is used to generate output names
    161180        }
    162181    }

  • branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c

    r33415 r34041  
    8989        return psphotReadoutCleanup (config, view, STACK_SRC);
    9090    }
    91     if (!psphotSubtractBackground (config, view, STACK_SRC)) {
    92         return psphotReadoutCleanup (config, view, STACK_SRC);
     91    if (strcmp(STACK_SRC, STACK_DET)) {
     92#define MODEL_BACKGROUND_SRC 1
     93#ifdef MODEL_BACKGROUND_SRC
     94        // work around the fact that the background levels on the convolved
     95        // and unconvolved stacks can be different
     96        if (!psphotModelBackground (config, view, STACK_SRC)) {
     97            return psphotReadoutCleanup (config, view, STACK_SRC);
     98        }
     99#endif
     100        if (!psphotSubtractBackground (config, view, STACK_SRC)) {
     101            return psphotReadoutCleanup (config, view, STACK_SRC);
     102        }
    93103    }
    94104    if (!strcasecmp (breakPt, "BACKMDL")) {
    … …  
    190200    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
    191201
     202
    192203    // if we only do one pass, skip to extended source analysis
    193204    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
    … …  
    211222            //  subtract all sources from DET (this will subtract using the psf model for SRC, which
    212223            //  will somewhat oversubtract the sources -- this is OK
    213             psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
     224            psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
    214225        }
    215226
    … …  
    223234        // remove noise for subtracted objects (ie, return to normal noise level)
    224235        // NOTE: this needs to operate only on the OLD sources
    225         psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     236        // NOTE: if fooprintsUseUnsubtracted, the noise has already been removed by psphotFindDetections
     237        bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
     238        if (!footprintsUseUnsubtracted) {
     239            psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     240        }
    226241
    227242        // if DET and SRC are different images, copy the detections from DET to SRC
    … …  
    302317    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
    303318
    304     // create source children for the OUT filerule (for radial aperture photometry)
    305     psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
    306     if (!objectsRadial) {
     319    // create source children for the OUT filerule (for radial aperture photometry and output)
     320    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
     321    if (!objectsOut) {
    307322        psFree(objects);
    308323        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    … …  
    310325    }
    311326
    312     // measure circular, radial apertures (objects sorted by S/N)
    313     // this forces photometry on the undetected sources from other images
    314     psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
    315     psMemDump("extmeas");
    316 
    317     int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
    318 
    319     for (int entry = 1; entry < nRadialEntries; entry++) {
    320         // NOTE: entry 0 is the unmatched image set
    321 
    322         // re-measure the PSF for the smoothed image (using entries in 'allSources')
    323         psphotChoosePSF (config, view, STACK_OUT, false);
    324 
    325         // this is necessary to update the models based on the new PSF
    326         psphotResetModels (config, view, STACK_OUT);
    327 
    328         // this is necessary to get the right normalization for the new models
    329         psphotFitSourcesLinear (config, view, STACK_OUT, false);
    330 
    331         // measure circular, radial apertures (objects sorted by S/N)
    332         // entry 0 == unmatched? pass entry + 1?
    333         psphotRadialApertures (config, view, STACK_OUT, entry);
    334 
    335         // replace the flux in the image so it is returned to its original state
    336         psphotReplaceAllSources (config, view, STACK_OUT, false);
    337 
    338         // smooth to the next FWHM, or set 'smoothAgain' to false if no more
    339         psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
    340         psMemDump("matched");
     327
     328    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
     329    if (radial_apertures) {
     330        // measure circular, radial apertures (objects sorted by S/N)
     331        // this forces photometry on the undetected sources from other images
     332
     333        // NOTE: we always do the radial apertures analysis on the convolved image since
     334        // those are the ones that are psf matched and the source of STACK_OUT's pixels
     335        // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true.
     336        // but in that case we don't get past the psf matching step because there is no
     337        // target psf for the RAW inputs
     338
     339        // If useRaw copy the sources to the convolved readout
     340        if (strcmp(STACK_SRC, STACK_CNV)) {
     341            if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) {
     342                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     343                return psphotReadoutCleanup (config, view, STACK_SRC);
     344            }
     345        }
     346        psphotRadialApertures (config, view, STACK_CNV, 0); // XXX entry 0 == unmatched?
     347        psMemDump("extmeas");
     348
     349        int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
     350
     351        for (int entry = 1; entry < nRadialEntries; entry++) {
     352            // NOTE: entry 0 is the unmatched image set
     353
     354            // re-measure the PSF for the smoothed image (using entries in 'allSources')
     355            psphotChoosePSF (config, view, STACK_OUT, false);
     356
     357            // this is necessary to update the models based on the new PSF
     358            psphotResetModels (config, view, STACK_OUT);
     359
     360            // this is necessary to get the right normalization for the new models
     361            psphotFitSourcesLinear (config, view, STACK_OUT, false);
     362
     363            // measure circular, radial apertures (objects sorted by S/N)
     364            // entry 0 == unmatched? pass entry + 1?
     365            psphotRadialApertures (config, view, STACK_OUT, entry);
     366
     367            // replace the flux in the image so it is returned to its original state
     368            psphotReplaceAllSources (config, view, STACK_OUT, false);
     369
     370            // smooth to the next FWHM, or set 'smoothAgain' to false if no more
     371            psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
     372            psMemDump("matched");
     373        }
    341374    }
    342375
    … …  
    344377    if (!psphotApResid (config, view, STACK_SRC)) {
    345378        psFree (objects);
     379        psFree (objectsOut);
    346380        psLogMsg ("psphot", 3, "failed on psphotApResid");
    347381        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    375409
    376410    psFree (objects);
    377     psFree (objectsRadial);
     411    psFree (objectsOut);
    378412
    379413    // create the exported-metadata and free local data

  • branches/meh_branches/ppstack_test/psphot/src/psphotTest.c

    r33415 r34041  
    103103# if (0)
    104104
    105     psRegion region = {0,0,0,0};        // a region representing the entire array
     105psRegion region = psRegionSet (0,0,0,0);        // a region representing the entire array
    106106    psphotTestArguments (&argc, argv);
    107107

  • branches/meh_branches/ppstack_test/psphot/src/psphotTestSourceOutput.c

    r21183 r34041  
    136136
    137137            // generate working image for this source
    138             psRegion region = {ix - dx, ix + dx, iy - dy, iy + dy};
     138            psRegion region = psRegionSet(ix - dx, ix + dx, iy - dy, iy + dy);
    139139
    140140            psImage *vM = psImageSubset (imMo, region);
Note: See TracChangeset for help on using the changeset viewer.