Index: branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c (revision 28435) +++ branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c (revision 28440) @@ -171,10 +171,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor); + float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -184,5 +184,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseVectorElement (sparse, i, f); @@ -190,11 +190,11 @@ switch (SKY_FIT_ORDER) { case 1: - f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 1, f); - f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 2, f); case 0: - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 0, f); break; @@ -216,5 +216,5 @@ // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseMatrixElement (sparse, j, i, f); } Index: branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinearStack.c (revision 28435) +++ branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinearStack.c (revision 28440) @@ -43,25 +43,25 @@ for (int i = 0; i < objects->n; i++) { pmPhotObj *object = objects->data[i]; - if (!object) continue; - if (!object->sources) continue; + if (!object) continue; + if (!object->sources) continue; - // XXX check an element of the group to see if we should use it - // if (!object->flags & PM_PHOT_OBJ_BAD) continue; + // XXX check an element of the group to see if we should use it + // if (!object->flags & PM_PHOT_OBJ_BAD) continue; - for (int j = 0; j < object->sources->n; j++) { - pmSource *source = object->sources->data[j]; - if (!source) continue; + for (int j = 0; j < object->sources->n; j++) { + pmSource *source = object->sources->data[j]; + if (!source) continue; - // turn this bit off and turn it on again if we keep this source - source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; + // turn this bit off and turn it on again if we keep this source + source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; - // generate model for sources without, or skip if we can't - if (!source->modelFlux) { + // generate model for sources without, or skip if we can't + if (!source->modelFlux) { if (!pmSourceCacheModel (source, maskVal)) continue; - } + } - source->mode |= PM_SOURCE_MODE_LINEAR_FIT; - psArrayAdd (fitSources, 100, source); - } + source->mode |= PM_SOURCE_MODE_LINEAR_FIT; + psArrayAdd (fitSources, 100, source); + } } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n); @@ -85,10 +85,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR); + float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR, maskVal); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -97,5 +97,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); psSparseVectorElement (sparse, i, f); @@ -104,6 +104,6 @@ pmSource *SRCj = fitSources->data[j]; - // we only need to generate dot terms for source on the same image - if (SRCj->imageID != SRCi->imageID) { continue; } + // we only need to generate dot terms for source on the same image + if (SRCj->imageID != SRCi->imageID) { continue; } // skip over disjoint source images, break after last possible overlap @@ -114,5 +114,5 @@ // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); psSparseMatrixElement (sparse, j, i, f); } Index: branches/eam_branches/ipp-20100621/psphot/src/psphotImageLoop.c =================================================================== --- branches/eam_branches/ipp-20100621/psphot/src/psphotImageLoop.c (revision 28435) +++ branches/eam_branches/ipp-20100621/psphot/src/psphotImageLoop.c (revision 28440) @@ -46,4 +46,21 @@ if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip."); + // Read WCS if easy. + // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell. + { + pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip + pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } + // try to load other supporting data (PSF, SRC, etc). // do not re-load the following three files @@ -67,19 +84,21 @@ // Update the header - pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); - if (hdu && hdu != lastHDU) { - psphotVersionHeaderFull(hdu->header); - lastHDU = hdu; + { + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; + } } - // if an external mask is supplied, ensure that NAN pixels are also masked - if (readout->mask) { - psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels - if (!pmReadoutMaskNonfinite(readout, maskSat)) { - psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); - psFree(view); - return false; - } - } + // if an external mask is supplied, ensure that NAN pixels are also masked + if (readout->mask) { + psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels + if (!pmReadoutMaskNonfinite(readout, maskSat)) { + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); + psFree(view); + return false; + } + } // run the actual photometry analysis on this chip/cell/readout @@ -91,14 +110,14 @@ } - // drop all versions of the internal files - status = true; - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); - if (!status) { - psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); - psFree (view); - return false; - } + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } } // save output which is saved at the chip level Index: branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c =================================================================== --- branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c (revision 28435) +++ branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c (revision 28440) @@ -4,7 +4,7 @@ static float PSF_FIT_NSIGMA; static float PSF_FIT_PADDING; -static float PSF_APERTURE = 0; // radius to use in PSF aperture mags -static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0, - // and a per-object radius is calculated) +static float PSF_APERTURE = 0; // radius to use in PSF aperture mags +static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0, + // and a per-object radius is calculated) bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) { @@ -17,10 +17,30 @@ PSF_FIT_RADIUS = psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS"); if (!status) { - PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); + PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); } PSF_APERTURE = psMetadataLookupF32(&status, analysis, "PSF_APERTURE"); if (!status) { - PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); + PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); + } + + // The PSF_FIT_RADIUS and PSF_APERTURE may not be set if the PSF was loaded and not chosen + + if (PSF_FIT_RADIUS == 0.0) { + float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE"); + PSF_FIT_RADIUS = (int)(fitScale*gaussSigma); + } + + if (PSF_APERTURE == 0.0) { + float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE"); + PSF_APERTURE = (int)(apScale*gaussSigma); } @@ -38,18 +58,18 @@ // set the fit radius based on the object flux limit and the model float radiusFit = PSF_FIT_RADIUS; - if (radiusFit <= 0) { // use fixed radius - if (moments == NULL) { - radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); - } else { - radiusFit = model->modelRadius(model->params, 1.0); - } - model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); + if (radiusFit <= 0) { // use fixed radius + if (moments == NULL) { + radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + } else { + radiusFit = model->modelRadius(model->params, 1.0); + } + model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); } else { - model->fitRadius = radiusFit; + model->fitRadius = radiusFit; } if (isnan(model->fitRadius)) psAbort("error in radius"); - + if (source->mode & PM_SOURCE_MODE_SATSTAR) { - model->fitRadius *= 2; + model->fitRadius *= 2; } @@ -73,13 +93,13 @@ // set the fit radius based on the object flux limit and the model float radiusFit = PSF_FIT_RADIUS; - if (radiusFit <= 0) { // use fixed radius - if (moments == NULL) { - radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); - } else { - radiusFit = model->modelRadius(model->params, 1.0); - } - model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); + if (radiusFit <= 0) { // use fixed radius + if (moments == NULL) { + radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + } else { + radiusFit = model->modelRadius(model->params, 1.0); + } + model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); } else { - model->fitRadius = radiusFit; + model->fitRadius = radiusFit; } if (isnan(model->fitRadius)) psAbort("error in radius"); @@ -89,5 +109,5 @@ if (source->mode & PM_SOURCE_MODE_SATSTAR) { - model->fitRadius *= 2; + model->fitRadius *= 2; } @@ -134,12 +154,12 @@ float radius = 0.0; for (int j = 0; j < footprint->spans->n; j++) { - pmSpan *span = footprint->spans->data[j]; - - float dY = span->y - peak->yf; - float dX0 = span->x0 - peak->xf; - float dX1 = span->x1 - peak->xf; - - radius = PS_MAX (radius, hypot(dY, dX0)); - radius = PS_MAX (radius, hypot(dY, dX1)); + pmSpan *span = footprint->spans->data[j]; + + float dY = span->y - peak->yf; + float dX0 = span->x0 - peak->xf; + float dX1 = span->x1 - peak->xf; + + radius = PS_MAX (radius, hypot(dY, dX0)); + radius = PS_MAX (radius, hypot(dY, dX1)); } Index: branches/eam_branches/ipp-20100621/psphot/src/psphotRoughClass.c =================================================================== --- branches/eam_branches/ipp-20100621/psphot/src/psphotRoughClass.c (revision 28435) +++ branches/eam_branches/ipp-20100621/psphot/src/psphotRoughClass.c (revision 28440) @@ -25,9 +25,9 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (i == chisqNum) continue; // skip chisq image - if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i); - return false; - } + return false; + } } return true; @@ -50,5 +50,5 @@ bool havePSF = false; if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { - havePSF = true; + havePSF = true; } @@ -60,6 +60,6 @@ if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); + return true; } @@ -78,4 +78,17 @@ psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + + // If in doubt, it's a PSF + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source || !source->peak) { + continue; + } + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + source->type = PM_SOURCE_TYPE_STAR; + } psFree (region); continue; @@ -124,22 +137,22 @@ // XXX why not save the psfClump as a PTR? - float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); - float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); + float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); + float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); - float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE"); - if (!status) { - PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); - psAssert (status, "missing PSF_CLUMP_GRID_SCALE"); - } - float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX"); - if (!status) { - MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); - psAssert (status, "missing MOMENTS_SX_MAX"); - } - float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX"); - if (!status) { - MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); - psAssert (status, "missing MOMENTS_SY_MAX"); - } + float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE"); + if (!status) { + PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); + psAssert (status, "missing PSF_CLUMP_GRID_SCALE"); + } + float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX"); + if (!status) { + MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); + psAssert (status, "missing MOMENTS_SX_MAX"); + } + float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX"); + if (!status) { + MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); + psAssert (status, "missing MOMENTS_SY_MAX"); + } psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); Index: branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c =================================================================== --- branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c (revision 28435) +++ branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c (revision 28440) @@ -489,4 +489,6 @@ 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]); +#if 0 + // Modifying clump parameters without restoring! psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000"); @@ -500,8 +502,8 @@ psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - if (pmVisualTestLevel("psphot.moments.full", 2)) { psphotVisualPlotMoments (recipe, analysis, sources); } +#endif Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];