IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34086 for trunk/psphot/src


Ignore:
Timestamp:
Jun 26, 2012, 11:33:10 AM (14 years ago)
Author:
eugene
Message:

re-enable MODEL_VAR option for linear photometry fit

Location:
trunk/psphot
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src

  • trunk/psphot/src/Makefile.am

    r33690 r34086  
    183183        psphotSourcePlots.c            \
    184184        psphotRadialPlot.c             \
    185         psphotKronMasked.c             \
    186185        psphotKronIterate.c            \
    187186        psphotRadialProfileWings.c     \
    … …  
    209208        psphotSetNFrames.c
    210209
    211 # re-instate these
     210# not currently used
    212211#       psphotIsophotal.c              \
    213212#       psphotAnnuli.c                 \
    214213#       psphotKron.c                   \
     214#       psphotKronMasked.c             \
    215215#
    216216

  • trunk/psphot/src/psphot.h

    r33994 r34086  
    6161bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    6262
     63bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
     64bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     65
    6366bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
    6467bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    … …  
    99102
    100103bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
    101 
    102 # if (HAVE_MODEL_VAR)
    103104bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
    104 # else
    105 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
    106 # endif
    107105
    108106bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
    … …  
    489487bool psphotSetNFramesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    490488
     489bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
     490pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
     491bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
     492
     493bool psphotModelBackgroundReadout(psImage *model,  // Model image
     494                                  psImage *modelStdev, // Model stdev image
     495                                  psMetadata *analysis, // Analysis metadata for outputs
     496                                  pmReadout *readout, // Readout for which to generate a background model
     497                                  psImageBinning *binning, // Binning parameters
     498                                  const pmConfig *config,// Configuration
     499                                  bool useVarianceImage
     500    );
     501
    491502#endif

  • trunk/psphot/src/psphotCleanup.c

    r29548 r34086  
    2929psExit psphotGetExitStatus (void) {
    3030
    31     psErrorCode err = psErrorCodeLast ();
     31    // gcc -Wswitch complains here if err is declared as type psErrorCode
     32    // the collection of ps*ErrorCode values are enums defined separately for
     33    // each module (psphot, pswarp, etc).  the lowest type, psErrorCode is only the base set and does
     34    // not include the possible psphot values
     35
     36    // for now, to get around this, we just use an int for the switch
     37
     38    // psErrorCode err = psErrorCodeLast ();
     39    int err = psErrorCodeLast ();
    3240    switch (err) {
    3341      case PS_ERR_NONE:

  • trunk/psphot/src/psphotEfficiency.c

    r33993 r34086  
    420420
    421421    // psphotFitSourcesLinearReadout subtracts the model fits
    422 # if (HAVE_MODEL_VAR)
    423422    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
    424 # else
    425     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    426 # endif
    427423        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    428424        psFree(fakeSources);

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r33994 r34086  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 # if (HAVE_MODEL_VAR)
    15 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
    16 pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
    17 bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
    18 # endif
    19 
    2014// for now, let's store the detections on the readout->analysis for each readout
    2115bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
    … …  
    3024    assert (recipe);
    3125
    32 # if (HAVE_MODEL_VAR)
    3326    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
    3427    if (!fitVarMode) {
    … …  
    4033    // do a single pass.
    4134    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
    42 # endif
    4335
    4436    int num = psphotFileruleCount(config, filerule);
    … …  
    6860        psAssert (psf, "missing psf?");
    6961
    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         {
     62        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
    7663            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    7764            return false;
    7865        }
    7966
    80 # if (HAVE_MODEL_VAR)
    8167        // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
    8268        // after we have determined the initial set of fits, then we can generate the variance image and
    … …  
    10490            }
    10591        }
    106 # endif
    10792
    10893        psphotVisualShowResidualImage (readout, (num > 0));
    … …  
    11398}
    11499
    115 # if (HAVE_MODEL_VAR)
    116100// look up the fit variance mode from the recipe; older recipes do not have the value
    117101// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
    … …  
    136120        return PM_SOURCE_PHOTFIT_IMAGE_VAR;
    137121    }
     122    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
     123        return PM_SOURCE_PHOTFIT_MODEL_SKY;
     124    }
    138125    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
    139126        return PM_SOURCE_PHOTFIT_MODEL_VAR;
    … …  
    142129    return PM_SOURCE_PHOTFIT_NONE;
    143130}
    144 # endif
    145 
    146 # if (HAVE_MODEL_VAR)
    147 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
    148 # else
    149 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
    150 # endif
    151 {
     131
     132bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
    152133    bool status;
    153134    float x;
    … …  
    185166    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    186167
    187 # if (!HAVE_MODEL_VAR)
    188     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    189         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    190     if (!status) {
    191         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    192     }
    193 # endif
    194168    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    195169    if (!status) {
    … …  
    299273            source->mode |= PM_SOURCE_MODE_PSFMODEL;
    300274        }           
    301        
     275
    302276        psArrayAdd (fitSources, 100, source);
    303277    }
    … …  
    319293    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
    320294
    321 # if (HAVE_MODEL_VAR)
    322295    // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
    323296    // XXX we have two possibilities here:
    … …  
    327300
    328301    // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
    329 # endif
    330302
    331303    // fill out the sparse matrix elements and border elements (B)
    … …  
    336308
    337309        // diagonal elements of the sparse matrix (auto-cross-product)
    338 # if (HAVE_MODEL_VAR)
    339         f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    340 # else
    341         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    342 # endif
    343         psSparseMatrixElement (sparse, i, i, f);
    344 
    345 # if (HAVE_MODEL_VAR)
    346         // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
     310        float MM = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
     311        psSparseMatrixElement (sparse, i, i, MM);
     312
     313        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
    347314        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
    348315            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
    349316            errors->data.F32[i] = 1.0 / sqrt(var);
    350317        } else {
    351             errors->data.F32[i] = 1.0 / sqrt(f);
    352         }
    353 # else
    354         // the formal error depends on the weighting scheme
    355         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    356             float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
    357             errors->data.F32[i] = 1.0 / sqrt(var);
    358         } else {
    359             errors->data.F32[i] = 1.0 / sqrt(f);
    360         }
    361 # endif
     318            errors->data.F32[i] = 1.0 / sqrt(MM);
     319        }
    362320
    363321        // find the image x model value
    364 # if (HAVE_MODEL_VAR)
    365         f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    366 # else
    367         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    368 # endif
    369         psSparseVectorElement (sparse, i, f);
     322        float FM = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
     323        psSparseVectorElement (sparse, i, FM);
    370324
    371325        // add the per-source variances (border region)
    372326        switch (SKY_FIT_ORDER) {
    373327          case 1:
    374 # if (HAVE_MODEL_VAR)
    375328            f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
    376329            psSparseBorderElementB (border, i, 1, f);
    377330            f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
    378331            psSparseBorderElementB (border, i, 2, f);
    379 # else
    380             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    381             psSparseBorderElementB (border, i, 1, f);
    382             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    383             psSparseBorderElementB (border, i, 2, f);
    384 # endif
    385332
    386333          case 0:
    387 # if (HAVE_MODEL_VAR)
    388334            f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
    389 # else
    390             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    391 # endif
    392335            psSparseBorderElementB (border, i, 0, f);
    393336            break;
    … …  
    409352
    410353            // got an overlap; calculate cross-product and add to output array
    411 # if (HAVE_MODEL_VAR)
    412354            f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
    413 # else
    414             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    415 # endif
    416355            psSparseMatrixElement (sparse, j, i, f);
    417356        }
    … …  
    451390
    452391    // set the sky, sky_x, sky_y components of border matrix
    453 # if (HAVE_MODEL_VAR)
    454392    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
    455 # else
    456     SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
    457 # endif
    458393
    459394    psSparseConstraint constraint;
    … …  
    500435    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    501436
     437    // mean stats on fit windows
     438    float sumRadius = 0.0;
     439    float sumPixels = 0.0;
     440    float sumSource = 0.0;
     441
    502442    // measure chisq for each source
    503443    // for (int i = 0; final && (i < fitSources->n); i++) {
    … …  
    505445        pmSource *source = fitSources->data[i];
    506446        pmModel *model = pmSourceGetModel (NULL, source);
     447
     448        // accumulate fit windows statistics
     449        sumRadius += model->fitRadius;
     450        sumPixels += M_PI*PS_SQR(model->fitRadius);
     451        sumSource += 1.0;
     452
    507453        if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
    508454            model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
    … …  
    511457    }
    512458    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
     459
     460    float meanRadius = sumRadius / sumSource;
     461    float meanPixels = sumPixels / sumSource;
    513462
    514463    // psFree (index);
    … …  
    520469    psFree (border);
    521470
    522     psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
     471    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs (mean radius = %f pixels, mean area = %f pixels: %f sec\n", meanRadius, meanPixels, psTimerMark ("psphot.linear"));
    523472
    524473    psphotVisualPlotChisq (sources);
    … …  
    531480    return true;
    532481}
    533 
    534 // XXX do we need this?
    535 // XXX disallow the simultaneous sky fit and remove this code...
    536482
    537483// Calculate the weight terms for the sky fit component of the matrix.  This function operates
    … …  
    614560}
    615561
    616 # if (HAVE_MODEL_VAR)
    617562bool psphotModelBackgroundReadout(psImage *model,  // Model image
    618563                                  psImage *modelStdev, // Model stdev image
    … …  
    624569    );
    625570
    626 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
     571 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
    627572
    628573    bool status = false;
    … …  
    650595    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
    651596
     597    // generate an image of the mean variance image in DN
    652598    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
    653599        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    … …  
    664610        return false;
    665611    }
    666 
    667612    psFree (varModel);
    668613    psFree (varModelStdev);
    669614
    670     // XXX for a test:
    671     psphotSaveImage (NULL, modelVar, "model.bck.fits");
     615    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
     616    pmCell *cell = readout->parent; // The parent cell
     617    if (cell) {
     618      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
     619      if (!status) {
     620        gain = 1.0;           // set note above
     621      }
     622    }
     623    if (gain > 2.0) { /* warn? */ }
     624    // XXX we are not actually using the gain, but need to test it to avoid gcc pedantic warnings
    672625
    673626    // insert all of the source models
    … …  
    690643
    691644        // add the source model to the model variance image
     645        // XXX note that this should be added with gain applied
     646        // var_DN = flux_DN / gain [e/DN]
     647        // to do this requires an API upgrade...
    692648        pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
    693649    }
    694 
    695     // XXX for a test:
    696     psphotSaveImage (NULL, modelVar, "model.var.fits");
    697     psphotSaveImage (NULL, readout->variance, "image.var.fits");
    698650
    699651    // we save the model variance for future reference
    … …  
    724676    return true;
    725677}
    726 # endif

  • trunk/psphot/src/psphotImageLoop.c

    r33691 r34086  
    164164    pmFPAfileActivate (config->files, true, "PSPHOT.EXPNUM");
    165165    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     166        if (! chip->process || ! chip->file_exists) { continue; }
    166167        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed attempting to load EXPNUM input for Chip in psphot.");
    167168        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {

  • trunk/psphot/src/psphotMakeResiduals.c

    r30624 r34086  
    177177
    178178                mflux = 0;
    179                 bool offImage = false;
    180179                if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
    181180                    // This pixel is off the image
    182                     offImage = true;
    183181                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
    184182                } else {

  • trunk/psphot/src/psphotMaskReadout.c

    r29936 r34086  
    9494
    9595    // test output of files at this stage
    96     if (psTraceGetLevel("psphot") >= 5) {
     96    if (psTraceGetLevel("psphot.imsave") >= 5) {
    9797        psphotSaveImage (NULL, readout->image,  "image.fits");
    9898        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    … …  
    105105    return true;
    106106}
     107
     108// XXX this function and support below was created to test the theory that the faint-end
     109// bias results from the Poisson variation of the background pixels.  This is NOT the
     110// case.  Using the code below maintains the faint-end bias.
     111bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
     112
     113    bool status = false;
     114
     115    // select the appropriate recipe information
     116    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     117    psAssert (recipe, "missing recipe?");
     118
     119    int num = psphotFileruleCount(config, filerule);
     120
     121    // loop over the available readouts
     122    for (int i = 0; i < num; i++) {
     123
     124        // Generate the mask and weight images, including the user-defined analysis region of interest
     125        if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
     126            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
     127            return false;
     128        }
     129    }
     130    return true;
     131}
     132
     133// determine the mean variance image (equivalent to the background model, but for the variance image)
     134// set the variance image to the MAX(input, mean)
     135bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     136
     137    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     138    psAssert (file, "missing file?");
     139
     140    // find the currently selected readout
     141    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     142    psAssert (readout, "missing readout?");
     143
     144    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
     145    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
     146      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
     147      return false;
     148    }
     149
     150    // make this an option via the recipe
     151    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
     152
     153    // create a model variance image (full-scale image to take result of psImageUnbin below)
     154    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
     155
     156    // find the binning information
     157    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
     158    assert (backBinning);
     159   
     160    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     161    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     162
     163    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
     164        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
     165        psFree (varModel);
     166        psFree (varModelStdev);
     167        return false;
     168    }
     169
     170    // linear interpolation to full-scale
     171    if (!psImageUnbin (modelVar, varModel, backBinning)) {
     172        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
     173        psFree (varModel);
     174        psFree (varModelStdev);
     175        return false;
     176    }
     177
     178    // XXX save these?
     179    psFree (varModel);
     180    psFree (varModelStdev);
     181
     182    psImage *im = readout->image;
     183    psImage *wt = readout->variance;
     184    for (int j = 0; j < im->numRows; j++) {
     185      for (int i = 0; i < im->numCols; i++) {
     186        if (!isfinite(im->data.F32[j][i])) continue;
     187        if (!isfinite(wt->data.F32[j][i])) continue;
     188        // XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
     189        wt->data.F32[j][i] = modelVar->data.F32[j][i];
     190      }
     191    }
     192
     193    // test output of files at this stage
     194    if (psTraceGetLevel("psphot.imsave") >= 5) {
     195        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
     196        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
     197        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
     198    }
     199
     200    psFree (modelVar);
     201
     202    return true;
     203}

  • trunk/psphot/src/psphotPetrosianStats.c

    r32348 r34086  
    4646
    4747    bool anyPetro = false;
    48     bool manyPetro = false;
     48    // bool manyPetro = false;  XXX not used
    4949    bool above = true;
    5050    float Asum = 0.0;
    … …  
    122122            }
    123123            above = false;
    124             if (anyPetro) manyPetro = true;
     124            // if (anyPetro) manyPetro = true;
    125125            anyPetro = true;
    126126        }
    … …  
    212212    source->extpars->petrosianR50    = R50;
    213213    source->extpars->petrosianR90    = R90;
    214     source->extpars->petrosianFill      = petApix / petArea;
     214    source->extpars->petrosianFill   = petApix / petArea;
    215215   
    216216    // XXX add the errors

  • trunk/psphot/src/psphotReadout.c

    r33915 r34086  
    99}
    1010
     11# if (0)
     12// TEST CODE, can be removed
     13bool psphotDumpFlux (pmConfig *config, const pmFPAview *view, const char *filerule) {
     14
     15    bool status = false;
     16
     17    // find the currently selected readout
     18    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
     19    psAssert (file, "missing file?");
     20
     21    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     22    psAssert (readout, "missing readout?");
     23
     24    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     25    psAssert (detections, "missing detections?");
     26
     27    psArray *sources = detections->allSources;
     28    psAssert (sources, "missing sources?");
     29
     30    static int npass;
     31    char filename[64];
     32    snprintf (filename, 64, "mags.%d.dat", npass);
     33    FILE *ftest = fopen (filename, "w");
     34    for (int j = 0; j < sources->n; j++) {
     35        pmSource *source = sources->data[j];
     36
     37        float psfMag;
     38        status = pmSourcePhotometryModel (&psfMag, NULL, source->modelPSF);
     39
     40        float psfMagNorm;
     41        float Io = source->modelPSF->params->data.F32[PM_PAR_I0];
     42        source->modelPSF->params->data.F32[PM_PAR_I0] = 1.0;
     43        status = pmSourcePhotometryModel (&psfMagNorm, NULL, source->modelPSF);
     44        source->modelPSF->params->data.F32[PM_PAR_I0] = Io;
     45
     46        // double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
     47        fprintf (ftest, "%d %d %d  %f %f %f %f  %f %f\n", j, source->peak->x, source->peak->y, source->modelPSF->params->data.F32[PM_PAR_I0], source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_SXY], psfMag, psfMagNorm);
     48    }
     49    fclose (ftest);
     50    npass++;
     51
     52    return true;
     53}
     54# endif
     55
    1156bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
    1257

  • trunk/psphot/src/psphotReadoutCleanup.c

    r29936 r34086  
    77
    88    // remove internal pmFPAfiles, if created
    9     if (psErrorCodeLast() == PSPHOT_ERR_DATA) {
     9    if (psErrorCodeLast() == (psErrorCode) PSPHOT_ERR_DATA) {
    1010        psErrorStackPrint(stderr, "Error in the psphot readout analysis");
    1111        psErrorClear();

  • trunk/psphot/src/psphotSignificanceImage.c

    r31154 r34086  
    88    float SIGMA_SMTH, NSIGMA_SMTH;
    99    bool status = false;
    10     bool guess = false;
    1110
    1211    // smooth the image and variance map
    … …  
    3534        SIGMA_SMTH  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0)));
    3635        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    37         guess = false;
    3836    } else {
    3937        // if we do not know the FWHM, use the guess smoothing kernel supplied.
    … …  
    4341        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    4442        PS_ASSERT (status, NULL);
    45         guess = true;
    4643    }
    4744    // record the actual smoothing sigma

  • trunk/psphot/src/psphotSourceFits.c

    r32744 r34086  
    224224    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    225225    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    226 
    227 # define TEST_X -420.0
    228 # define TEST_Y 300.0
    229    
    230     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    231         fprintf (stderr, "test galaxy\n");
    232     }
    233 
    234 # undef TEST_X
    235 # undef TEST_Y
    236226
    237227    // set the radius based on the footprint (also sets the mask pixels)
    … …  
    513503    }
    514504
    515 # define TEST_X -540.0
    516 # define TEST_Y 540.0
    517    
    518     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    519         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    520     }
    521 
    522505    float t1, t2, t4, t5;
    523506    if (TIMING) { psTimerStart ("psphotFitPCM"); }
    … …  
    569552    }
    570553
    571     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    572         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    573     }
    574 
    575554    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    576555    psFree (pcm);
    … …  
    578557    return model;
    579558}
    580 
    581 # undef TEST_X
    582 # undef TEST_Y
    583559
    584560// note that these should be 1/2n of the standard sersic index
    … …  
    603579    float xMin = NAN;
    604580    float chiSquare[N_INDEX_GUESS];
    605 
    606 # define TEST_X -540.0
    607 # define TEST_Y 540.0
    608    
    609     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    610         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    611     }
    612581
    613582    for (int i = 0; i < N_INDEX_GUESS; i++) {
    … …  
    635604    assert (iMin >= 0);
    636605
    637     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    638         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    639     }
    640 
    641606    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    642607    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
    … …  
    666631    float xMin = NAN;
    667632    float chiSquare[N_INDEX_GUESS];
    668 
    669     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    670         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    671     }
    672633
    673634    for (int i = 0; i < N_INDEX_GUESS; i++) {
    … …  
    701662    }
    702663    assert (iMin >= 0);
    703    
    704     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    705         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    706     }
    707664
    708665    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it

  • trunk/psphot/src/psphotStackImageLoop.c

Note: See TracChangeset for help on using the changeset viewer.