IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28440


Ignore:
Timestamp:
Jun 23, 2010, 2:36:55 PM (16 years ago)
Author:
eugene
Message:

merge changes from trunk (pauls fixes for loading pre-calculated psfs)

Location:
branches/eam_branches/ipp-20100621
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c

    r27531 r28440  
    174174    PS_ASSERT_PTR_NON_NULL(chip, false);
    175175
    176     if (!pmPSFmodelWrite(chip->analysis, view, file, config)) {
     176    // We need the readout as well, because that has the PSF analysis data (e.g., clumps)
     177    // There is only one, because photometry is done on chip-mosaicked data.
     178    pmFPAview *roView = pmFPAviewAlloc(0); // View to readout
     179    *roView = *view;
     180    roView->cell = 0;
     181    roView->readout = 0;
     182    pmReadout *ro = pmFPAviewThisReadout(roView, chip->parent); // Readout with analysis data
     183    psFree(roView);
     184    PM_ASSERT_READOUT_NON_NULL(ro, false);
     185
     186    if (!pmPSFmodelWrite(chip->analysis, ro->analysis, view, file, config)) {
    177187        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
    178188        return false;
    … …  
    189199// else
    190200//   - psf table (+header) : FITS Table
    191 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view,
     201bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view,
    192202                      pmFPAfile *file, pmConfig *config)
    193203{
    … …  
    198208    char *headName, *tableName, *residName;
    199209
    200     if (!analysis) {
     210    if (!chipAnalysis) {
    201211        psError(PM_ERR_PROG, true, "No analysis metadata for chip.");
     212        return false;
     213    }
     214    if (!roAnalysis) {
     215        psError(PM_ERR_PROG, true, "No analysis metadata for readout.");
    202216        return false;
    203217    }
    … …  
    314328
    315329    // select the psf of interest
    316     pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
     330    pmPSF *psf = psMetadataLookupPtr (&status, chipAnalysis, "PSPHOT.PSF");
    317331    if (!psf) {
    318332        psError(PM_ERR_PROG, true, "missing PSF for this analysis metadata");
    … …  
    346360
    347361        // we now save clump parameters for each region : need to save all of those
    348         int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
     362        int nRegions = psMetadataLookupS32 (&status, roAnalysis, "PSF.CLUMP.NREGIONS");
    349363        psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions);
    350364        for (int i = 0; i < nRegions; i++) {
    351365            char regionName[64];
    352366            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    353             psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     367            psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
    354368
    355369            psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   assert (status);
    … …  
    492506
    493507        // write the residuals as planes of the image
    494         psArray *images = psArrayAllocEmpty (1);
    495         psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
     508        psArray *images = psArrayAllocEmpty (1);
     509        psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
    496510
    497511        if (psf->residuals->Rx) {
    498512            psArrayAdd (images, 1, psf->residuals->Rx);
    499513            psArrayAdd (images, 1, psf->residuals->Ry);
    500         }
    501 
    502         // note that all N plane are implicitly of the same type, so we convert the mask
    503         if (psf->residuals->mask) {
    504             psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
    505             psArrayAdd (images, 1, mask);
    506             psFree (mask);
    507         }
    508 
    509         // psFitsWriteImageCube (file->fits, header, images, residName);
    510         // psFree (images);
    511 
    512         if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
     514        }
     515
     516        // note that all N plane are implicitly of the same type, so we convert the mask
     517        if (psf->residuals->mask) {
     518            psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
     519            psArrayAdd (images, 1, mask);
     520            psFree (mask);
     521        }
     522
     523        // psFitsWriteImageCube (file->fits, header, images, residName);
     524        // psFree (images);
     525
     526        if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
    513527            psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
    514528            psFree(images);
    … …  
    728742    PS_ASSERT_PTR_NON_NULL(file->fpa, false);
    729743
    730     if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
     744    // We need the readout as well, because that has the PSF analysis data (e.g., clumps)
     745    // There may be only one, because photometry is done on chip-mosaicked data.
     746    if (chip->cells->n != 1) {
     747        psError(PM_ERR_PROG, true, "Chip to receive PSF has %ld cells (should be only one)",
     748                chip->cells->n);
     749        return false;
     750    }
     751    pmCell *cell = chip->cells->data[0]; // Cell to receive PSF
     752    pmReadout *ro = NULL;                // Readout to receive PSF
     753    if (cell->readouts->n == 0) {
     754        ro = pmReadoutAlloc(cell);
     755        psFree(ro);                     // Drop reference
     756    } else if (cell->readouts->n != 1) {
     757        psError(PM_ERR_PROG, true, "Cell to receive PSF has %ld readouts (should be only one)",
     758                cell->readouts->n);
     759        return false;
     760    } else {
     761        ro = cell->readouts->data[0];
     762    }
     763    PM_ASSERT_READOUT_NON_NULL(ro, false);
     764    if (!ro->analysis) {
     765        ro->analysis = psMetadataAlloc();
     766    }
     767
     768    if (!pmPSFmodelRead(chip->analysis, ro->analysis, view, file, config)) {
    731769        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
    732770        return false;
    … …  
    737775// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters,
    738776// and header + image for the PSF residual images
    739 bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
     777bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    740778{
     779    PS_ASSERT_METADATA_NON_NULL(chipAnalysis, false);
     780    PS_ASSERT_METADATA_NON_NULL(roAnalysis, false);
    741781    PS_ASSERT_PTR_NON_NULL(view, false);
    742782    PS_ASSERT_PTR_NON_NULL(file, false);
    … …  
    809849        char regionName[64];
    810850        snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
    811         psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
    812 
    813         psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     851        psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
     852
     853        psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
    814854        if (!regionMD) {
    815855            regionMD = psMetadataAlloc();
    816             psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
     856            psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
    817857            psFree (regionMD);
    818858        }
    … …  
    831871        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
    832872    } else {
    833         psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
     873        psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
    834874
    835875        for (int i = 0; i < nRegions; i++) {
    … …  
    838878            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    839879
    840             psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     880            psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
    841881            if (!regionMD) {
    842882                regionMD = psMetadataAlloc();
    843                 psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
     883                psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
    844884                psFree (regionMD);
    845885            }
    … …  
    10191059        }
    10201060
    1021         // note that all N plane are implicitly of the same type, so we convert the mask
    1022         psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
    1023         psImageInit (psf->residuals->mask, 0);
    1024         psImageInit (psf->residuals->Rx, 0.0);
    1025         psImageInit (psf->residuals->Ry, 0.0);
    1026         switch (Nz) {
    1027           case 1: // Ro only
    1028             break;
    1029           case 2: // Ro and mask
     1061        // note that all N plane are implicitly of the same type, so we convert the mask
     1062        psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
     1063        psImageInit (psf->residuals->mask, 0);
     1064        psImageInit (psf->residuals->Rx, 0.0);
     1065        psImageInit (psf->residuals->Ry, 0.0);
     1066        switch (Nz) {
     1067          case 1: // Ro only
     1068            break;
     1069          case 2: // Ro and mask
    10301070            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) {
    10311071                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    10321072                return false;
    10331073            }
    1034             psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
    1035             break;
    1036           case 3: // Ro, Rx and Ry, no mask
     1074            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
     1075            break;
     1076          case 3: // Ro, Rx and Ry, no mask
    10371077            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
    10381078                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    … …  
    10431083                return false;
    10441084            }
    1045             break;
    1046           case 4: // Ro, Rx, Ry, and mask:
     1085            break;
     1086          case 4: // Ro, Rx, Ry, and mask:
    10471087            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
    10481088                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    … …  
    10571097                return false;
    10581098            }
    1059             psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
    1060             break;
    1061         }
    1062         psFree (mask);
    1063     }
    1064 
    1065     psMetadataAdd (analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     1099            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
     1100            break;
     1101        }
     1102        psFree (mask);
     1103    }
     1104
     1105    psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
    10661106    psFree (psf);
    10671107

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.h

    r18601 r28440  
    2020bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    2121bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    22 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
     22bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    2323
    2424bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    … …  
    2727bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
    2828bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
    29 bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
     29bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
    3030
    3131bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file);

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c

    r27531 r28440  
    107107        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    108108        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
    109         psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
    110         psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
    111         source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
    112         source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
    113         source->psfMag = -2.5*log10(source->psfFlux);
     109        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
     110        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
     111        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
     112        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
     113        source->psfMag = -2.5*log10(source->psfFlux);
    114114    } else {
    115115        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
    116         source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     116        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
    117117    }
    118118
    … …  
    120120    if (source->modelFits) {
    121121        bool foundEXT = false;
    122         for (int i = 0; i < source->modelFits->n; i++) {
    123             pmModel *model = source->modelFits->data[i];
    124             status = pmSourcePhotometryModel (&model->mag, NULL, model);
    125             if (model == source->modelEXT) foundEXT = true;
    126         }
    127         if (foundEXT) {
    128             source->extMag = source->modelEXT->mag;
    129         } else {
    130             status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
    131         }
     122        for (int i = 0; i < source->modelFits->n; i++) {
     123            pmModel *model = source->modelFits->data[i];
     124            status = pmSourcePhotometryModel (&model->mag, NULL, model);
     125            if (model == source->modelEXT) foundEXT = true;
     126        }
     127        if (foundEXT) {
     128            source->extMag = source->modelEXT->mag;
     129        } else {
     130            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     131        }
    132132    } else {
    133         if (source->modelEXT) {
    134             status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
    135         }
     133        if (source->modelEXT) {
     134            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     135        }
    136136    }
    137137
    … …  
    204204        }
    205205        if (mode & PM_SOURCE_PHOT_APCORR) {
    206             // XXX this should be removed -- we no longer fit for the 'sky bias'
     206            // XXX this should be removed -- we no longer fit for the 'sky bias'
    207207            rflux   = pow (10.0, 0.4*source->psfMag);
    208208            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    … …  
    236236    flux = model->modelFlux (model->params);
    237237    if (flux > 0) {
    238         mag = -2.5*log10(flux);
     238        mag = -2.5*log10(flux);
    239239    }
    240240    if (fitMag) {
    241         *fitMag = mag;
     241        *fitMag = mag;
    242242    }
    243243    if (fitFlux) {
    244         *fitFlux = flux;
     244        *fitFlux = flux;
    245245    }
    246246
    … …  
    380380
    381381    if (source->diffStats == NULL) {
    382         source->diffStats = pmSourceDiffStatsAlloc();
     382        source->diffStats = pmSourceDiffStatsAlloc();
    383383    }
    384384
    … …  
    388388    int   nMask = 0;
    389389    int   nBad  = 0;
    390    
     390
    391391    psImage *flux     = source->pixels;
    392392    psImage *variance = source->variance;
    … …  
    394394
    395395    for (int iy = 0; iy < flux->numRows; iy++) {
    396         for (int ix = 0; ix < flux->numCols; ix++) {
     396        for (int ix = 0; ix < flux->numCols; ix++) {
    397397            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
    398                 nMask ++;
    399                 continue;
    400             }
    401 
    402             float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    403 
    404             if (SN > +FLUX_LIMIT) {
    405                 nGood ++;
    406                 fGood += fabs(flux->data.F32[iy][ix]);
    407             }
    408 
    409             if (SN < -FLUX_LIMIT) {
    410                 nBad ++;
    411                 fBad += fabs(flux->data.F32[iy][ix]);
    412             }
    413         }
     398                nMask ++;
     399                continue;
     400            }
     401
     402            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     403
     404            if (SN > +FLUX_LIMIT) {
     405                nGood ++;
     406                fGood += fabs(flux->data.F32[iy][ix]);
     407            }
     408
     409            if (SN < -FLUX_LIMIT) {
     410                nBad ++;
     411                fBad += fabs(flux->data.F32[iy][ix]);
     412            }
     413        }
    414414    }
    415415
    416416    source->diffStats->nGood      = nGood;
    417     source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);         
    418     source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);         
    419     source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);         
     417    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);
     418    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);
     419    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);
    420420    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
    421421
    … …  
    628628        }
    629629    }
    630     model->nPix = Npix; 
     630    model->nPix = Npix;
    631631    model->nDOF = Npix - 1;
    632632    model->chisq = dC;
    … …  
    636636
    637637
    638 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
     638double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    639639{
    640640    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    652652    for (int yi = 0; yi < Pi->numRows; yi++) {
    653653        for (int xi = 0; xi < Pi->numCols; xi++) {
    654             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
     654            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
    655655                continue;
    656656            if (!unweighted_sum) {
    … …  
    684684}
    685685
    686 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
     686double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    687687{
    688688    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    727727    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    728728        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    729             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
    730                 continue;
    731             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
     729            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     730                continue;
     731            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    732732                continue;
    733733
    … …  
    746746}
    747747
    748 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
     748double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    749749{
    750750    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    789789    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    790790        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    791             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
    792                 continue;
    793             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
     791            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     792                continue;
     793            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    794794                continue;
    795795

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h

    r27531 r28440  
    4848    psImage *image,                     ///< image pixels to be used
    4949    psImage *mask,                      ///< mask of pixels to ignore
    50     psImageMaskType maskVal             ///< Value to mask
     50    psImageMaskType maskVal             ///< Value to mask
    5151);
    5252
    … …  
    5858bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
    5959
    60 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
    61 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
    62 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor);
     60double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     61double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     62double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    6363
    6464// retire these:

  • branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c

    r28399 r28440  
    171171
    172172        // diagonal elements of the sparse matrix (auto-cross-product)
    173         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     173        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    174174        psSparseMatrixElement (sparse, i, i, f);
    175175
    176176        // the formal error depends on the weighting scheme
    177177        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    178             float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
     178            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
    179179            errors->data.F32[i] = 1.0 / sqrt(var);
    180180        } else {
    … …  
    184184
    185185        // find the image x model value
    186         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     186        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    187187        psSparseVectorElement (sparse, i, f);
    188188
    … …  
    190190        switch (SKY_FIT_ORDER) {
    191191          case 1:
    192             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     192            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    193193            psSparseBorderElementB (border, i, 1, f);
    194             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     194            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    195195            psSparseBorderElementB (border, i, 2, f);
    196196
    197197          case 0:
    198             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     198            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    199199            psSparseBorderElementB (border, i, 0, f);
    200200            break;
    … …  
    216216
    217217            // got an overlap; calculate cross-product and add to output array
    218             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     218            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    219219            psSparseMatrixElement (sparse, j, i, f);
    220220        }

  • branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinearStack.c

    r28013 r28440  
    4343    for (int i = 0; i < objects->n; i++) {
    4444        pmPhotObj *object = objects->data[i];
    45         if (!object) continue;
    46         if (!object->sources) continue;
     45        if (!object) continue;
     46        if (!object->sources) continue;
    4747
    48         // XXX check an element of the group to see if we should use it
    49         // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
     48        // XXX check an element of the group to see if we should use it
     49        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
    5050
    51         for (int j = 0; j < object->sources->n; j++) {
    52           pmSource *source = object->sources->data[j];
    53           if (!source) continue;
     51        for (int j = 0; j < object->sources->n; j++) {
     52          pmSource *source = object->sources->data[j];
     53          if (!source) continue;
    5454
    55           // turn this bit off and turn it on again if we keep this source
    56           source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
     55          // turn this bit off and turn it on again if we keep this source
     56          source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
    5757
    58           // generate model for sources without, or skip if we can't
    59           if (!source->modelFlux) {
     58          // generate model for sources without, or skip if we can't
     59          if (!source->modelFlux) {
    6060            if (!pmSourceCacheModel (source, maskVal)) continue;
    61           }
     61          }
    6262
    63           source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
    64           psArrayAdd (fitSources, 100, source);
    65         }
     63          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
     64          psArrayAdd (fitSources, 100, source);
     65        }
    6666    }
    6767    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n);
    … …  
    8585
    8686        // diagonal elements of the sparse matrix (auto-cross-product)
    87         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
     87        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
    8888        psSparseMatrixElement (sparse, i, i, f);
    8989
    9090        // the formal error depends on the weighting scheme
    9191        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    92             float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR);
     92            float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR, maskVal);
    9393            errors->data.F32[i] = 1.0 / sqrt(var);
    9494        } else {
    … …  
    9797
    9898        // find the image x model value
    99         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
     99        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
    100100        psSparseVectorElement (sparse, i, f);
    101101
    … …  
    104104            pmSource *SRCj = fitSources->data[j];
    105105
    106             // we only need to generate dot terms for source on the same image
    107             if (SRCj->imageID != SRCi->imageID) { continue; }
     106            // we only need to generate dot terms for source on the same image
     107            if (SRCj->imageID != SRCi->imageID) { continue; }
    108108
    109109            // skip over disjoint source images, break after last possible overlap
    … …  
    114114
    115115            // got an overlap; calculate cross-product and add to output array
    116             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
     116            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
    117117            psSparseMatrixElement (sparse, j, i, f);
    118118        }

  • branches/eam_branches/ipp-20100621/psphot/src/psphotImageLoop.c

    r27657 r28440  
    4646        if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip.");
    4747
     48        // Read WCS if easy.
     49        // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell.
     50        {
     51            pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip
     52            pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL);
     53            if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) {
     54                psWarning("Unable to read WCS astrometry from header.");
     55                psErrorClear();
     56                pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL);
     57                if (hduHigh && hduHigh != hduLow &&
     58                    !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) {
     59                    psWarning("Unable to read WCS astrometry from primary header.");
     60                    psErrorClear();
     61                }
     62            }
     63        }
     64
    4865        // try to load other supporting data (PSF, SRC, etc).
    4966        // do not re-load the following three files
    … …  
    6784
    6885                // Update the header
    69                 pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
    70                 if (hdu && hdu != lastHDU) {
    71                     psphotVersionHeaderFull(hdu->header);
    72                     lastHDU = hdu;
     86                {
     87                    pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
     88                    if (hdu && hdu != lastHDU) {
     89                        psphotVersionHeaderFull(hdu->header);
     90                        lastHDU = hdu;
     91                    }
    7392                }
    7493
    75                 // if an external mask is supplied, ensure that NAN pixels are also masked
    76                 if (readout->mask) {
    77                     psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
    78                     if (!pmReadoutMaskNonfinite(readout, maskSat)) {
    79                         psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
    80                         psFree(view);
    81                         return false;
    82                     }
    83                 }
     94                // if an external mask is supplied, ensure that NAN pixels are also masked
     95                if (readout->mask) {
     96                    psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
     97                    if (!pmReadoutMaskNonfinite(readout, maskSat)) {
     98                        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
     99                        psFree(view);
     100                        return false;
     101                    }
     102                }
    84103
    85104                // run the actual photometry analysis on this chip/cell/readout
    … …  
    91110            }
    92111
    93             // drop all versions of the internal files
    94             status = true;
    95             status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
    96             status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
    97             status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
    98             if (!status) {
    99                 psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
    100                 psFree (view);
    101                 return false;
    102             }
     112            // drop all versions of the internal files
     113            status = true;
     114            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
     115            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
     116            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
     117            if (!status) {
     118                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
     119                psFree (view);
     120                return false;
     121            }
    103122        }
    104123        // save output which is saved at the chip level

  • branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c

    r26894 r28440  
    44static float PSF_FIT_NSIGMA;
    55static float PSF_FIT_PADDING;
    6 static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
    7 static float PSF_FIT_RADIUS = 0;        // radius to use in fitting (ignored if <= 0,
    8                                         // and a per-object radius is calculated)
     6static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
     7static float PSF_FIT_RADIUS = 0;        // radius to use in fitting (ignored if <= 0,
     8                                        // and a per-object radius is calculated)
    99
    1010bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {
    … …  
    1717    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
    1818    if (!status) {
    19         PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
     19        PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
    2020    }
    2121
    2222    PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
    2323    if (!status) {
    24         PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
     24        PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
     25    }
     26
     27    // The PSF_FIT_RADIUS and PSF_APERTURE may not be set if the PSF was loaded and not chosen
     28
     29    if (PSF_FIT_RADIUS == 0.0) {
     30        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
     31        if (!status) {
     32            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
     33        }
     34        float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE");
     35        PSF_FIT_RADIUS = (int)(fitScale*gaussSigma);
     36    }
     37
     38    if (PSF_APERTURE == 0.0) {
     39        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
     40        if (!status) {
     41            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
     42        }
     43        float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
     44        PSF_APERTURE = (int)(apScale*gaussSigma);
    2545    }
    2646
    … …  
    3858    // set the fit radius based on the object flux limit and the model
    3959    float radiusFit = PSF_FIT_RADIUS;
    40     if (radiusFit <= 0) {               // use fixed radius
    41         if (moments == NULL) {
    42             radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
    43         } else {
    44             radiusFit = model->modelRadius(model->params, 1.0);
    45         }
    46         model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
     60    if (radiusFit <= 0) {               // use fixed radius
     61        if (moments == NULL) {
     62            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
     63        } else {
     64            radiusFit = model->modelRadius(model->params, 1.0);
     65        }
     66        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
    4767    } else {
    48         model->fitRadius = radiusFit;
     68        model->fitRadius = radiusFit;
    4969    }
    5070    if (isnan(model->fitRadius)) psAbort("error in radius");
    51        
     71
    5272    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    53         model->fitRadius *= 2;
     73        model->fitRadius *= 2;
    5474    }
    5575
    … …  
    7393    // set the fit radius based on the object flux limit and the model
    7494    float radiusFit = PSF_FIT_RADIUS;
    75     if (radiusFit <= 0) {               // use fixed radius
    76         if (moments == NULL) {
    77             radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
    78         } else {
    79             radiusFit = model->modelRadius(model->params, 1.0);
    80         }
    81         model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
     95    if (radiusFit <= 0) {               // use fixed radius
     96        if (moments == NULL) {
     97            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
     98        } else {
     99            radiusFit = model->modelRadius(model->params, 1.0);
     100        }
     101        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
    82102    } else {
    83         model->fitRadius = radiusFit;
     103        model->fitRadius = radiusFit;
    84104    }
    85105    if (isnan(model->fitRadius)) psAbort("error in radius");
    … …  
    89109
    90110    if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
    91         model->fitRadius *= 2;
     111        model->fitRadius *= 2;
    92112    }
    93113
    … …  
    134154    float radius = 0.0;
    135155    for (int j = 0; j < footprint->spans->n; j++) {
    136         pmSpan *span = footprint->spans->data[j];
    137 
    138         float dY  = span->y  - peak->yf;
    139         float dX0 = span->x0 - peak->xf;
    140         float dX1 = span->x1 - peak->xf;
    141 
    142         radius = PS_MAX (radius, hypot(dY, dX0));
    143         radius = PS_MAX (radius, hypot(dY, dX1));
     156        pmSpan *span = footprint->spans->data[j];
     157
     158        float dY  = span->y  - peak->yf;
     159        float dX0 = span->x0 - peak->xf;
     160        float dX1 = span->x1 - peak->xf;
     161
     162        radius = PS_MAX (radius, hypot(dY, dX0));
     163        radius = PS_MAX (radius, hypot(dY, dX1));
    144164    }
    145165

  • branches/eam_branches/ipp-20100621/psphot/src/psphotRoughClass.c

    r28013 r28440  
    2525    // loop over the available readouts
    2626    for (int i = 0; i < num; i++) {
    27         if (i == chisqNum) continue; // skip chisq image
    28         if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) {
     27        if (i == chisqNum) continue; // skip chisq image
     28        if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) {
    2929            psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i);
    30             return false;
    31         }
     30            return false;
     31        }
    3232    }
    3333    return true;
    … …  
    5050    bool havePSF = false;
    5151    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
    52         havePSF = true;
     52        havePSF = true;
    5353    }
    5454
    … …  
    6060
    6161    if (!sources->n) {
    62         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");
    63         return true;
     62        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");
     63        return true;
    6464    }
    6565
    … …  
    7878                psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n",
    7979                         region->x0, region->y0, region->x1, region->y1);
     80
     81                // If in doubt, it's a PSF
     82                for (int i = 0; i < sources->n; i++) {
     83                    pmSource *source = sources->data[i]; // Source of interest
     84                    if (!source || !source->peak) {
     85                        continue;
     86                    }
     87                    if (source->peak->x <  region->x0) continue;
     88                    if (source->peak->x >= region->x1) continue;
     89                    if (source->peak->y <  region->y0) continue;
     90                    if (source->peak->y >= region->y1) continue;
     91                    source->type = PM_SOURCE_TYPE_STAR;
     92                }
    8093                psFree (region);
    8194                continue;
    … …  
    124137        // XXX why not save the psfClump as a PTR?
    125138
    126         float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
    127         float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
     139        float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     140        float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
    128141
    129         float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE");
    130         if (!status) {
    131             PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
    132             psAssert (status, "missing PSF_CLUMP_GRID_SCALE");
    133         }
    134         float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX");
    135         if (!status) {
    136             MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
    137             psAssert (status, "missing MOMENTS_SX_MAX");
    138         }
    139         float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX");
    140         if (!status) {
    141             MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
    142             psAssert (status, "missing MOMENTS_SY_MAX");
    143         }
     142        float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE");
     143        if (!status) {
     144            PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
     145            psAssert (status, "missing PSF_CLUMP_GRID_SCALE");
     146        }
     147        float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX");
     148        if (!status) {
     149            MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
     150            psAssert (status, "missing MOMENTS_SX_MAX");
     151        }
     152        float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX");
     153        if (!status) {
     154            MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
     155            psAssert (status, "missing MOMENTS_SY_MAX");
     156        }
    144157
    145158        psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);

  • branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c

    r28435 r28440  
    489489        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    490490
     491#if 0
     492        // Modifying clump parameters without restoring!
    491493        psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
    492494        psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000");
    … …  
    500502        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
    501503        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
    502 
    503504        if (pmVisualTestLevel("psphot.moments.full", 2)) {
    504505            psphotVisualPlotMoments (recipe, analysis, sources);
    505506        }
     507#endif
    506508
    507509        Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
Note: See TracChangeset for help on using the changeset viewer.