IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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/psphot/src
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • 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.