Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c (revision 28439) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c (revision 28440) @@ -174,5 +174,15 @@ PS_ASSERT_PTR_NON_NULL(chip, false); - if (!pmPSFmodelWrite(chip->analysis, view, file, config)) { + // We need the readout as well, because that has the PSF analysis data (e.g., clumps) + // There is only one, because photometry is done on chip-mosaicked data. + pmFPAview *roView = pmFPAviewAlloc(0); // View to readout + *roView = *view; + roView->cell = 0; + roView->readout = 0; + pmReadout *ro = pmFPAviewThisReadout(roView, chip->parent); // Readout with analysis data + psFree(roView); + PM_ASSERT_READOUT_NON_NULL(ro, false); + + if (!pmPSFmodelWrite(chip->analysis, ro->analysis, view, file, config)) { psError(psErrorCodeLast(), false, "Failed to write PSF for chip"); return false; @@ -189,5 +199,5 @@ // else // - psf table (+header) : FITS Table -bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, +bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config) { @@ -198,6 +208,10 @@ char *headName, *tableName, *residName; - if (!analysis) { + if (!chipAnalysis) { psError(PM_ERR_PROG, true, "No analysis metadata for chip."); + return false; + } + if (!roAnalysis) { + psError(PM_ERR_PROG, true, "No analysis metadata for readout."); return false; } @@ -314,5 +328,5 @@ // select the psf of interest - pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF"); + pmPSF *psf = psMetadataLookupPtr (&status, chipAnalysis, "PSPHOT.PSF"); if (!psf) { psError(PM_ERR_PROG, true, "missing PSF for this analysis metadata"); @@ -346,10 +360,10 @@ // we now save clump parameters for each region : need to save all of those - int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + int nRegions = psMetadataLookupS32 (&status, roAnalysis, "PSF.CLUMP.NREGIONS"); psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions); for (int i = 0; i < nRegions; i++) { char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); assert (status); @@ -492,23 +506,23 @@ // write the residuals as planes of the image - psArray *images = psArrayAllocEmpty (1); - psArrayAdd (images, 1, psf->residuals->Ro); // z = 0 is Ro + psArray *images = psArrayAllocEmpty (1); + psArrayAdd (images, 1, psf->residuals->Ro); // z = 0 is Ro if (psf->residuals->Rx) { psArrayAdd (images, 1, psf->residuals->Rx); psArrayAdd (images, 1, psf->residuals->Ry); - } - - // note that all N plane are implicitly of the same type, so we convert the mask - if (psf->residuals->mask) { - psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type); - psArrayAdd (images, 1, mask); - psFree (mask); - } - - // psFitsWriteImageCube (file->fits, header, images, residName); - // psFree (images); - - if (!psFitsWriteImageCube (file->fits, header, images, residName)) { + } + + // note that all N plane are implicitly of the same type, so we convert the mask + if (psf->residuals->mask) { + psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type); + psArrayAdd (images, 1, mask); + psFree (mask); + } + + // psFitsWriteImageCube (file->fits, header, images, residName); + // psFree (images); + + if (!psFitsWriteImageCube (file->fits, header, images, residName)) { psError(psErrorCodeLast(), false, "Unable to write PSF residuals."); psFree(images); @@ -728,5 +742,29 @@ PS_ASSERT_PTR_NON_NULL(file->fpa, false); - if (!pmPSFmodelRead (chip->analysis, view, file, config)) { + // We need the readout as well, because that has the PSF analysis data (e.g., clumps) + // There may be only one, because photometry is done on chip-mosaicked data. + if (chip->cells->n != 1) { + psError(PM_ERR_PROG, true, "Chip to receive PSF has %ld cells (should be only one)", + chip->cells->n); + return false; + } + pmCell *cell = chip->cells->data[0]; // Cell to receive PSF + pmReadout *ro = NULL; // Readout to receive PSF + if (cell->readouts->n == 0) { + ro = pmReadoutAlloc(cell); + psFree(ro); // Drop reference + } else if (cell->readouts->n != 1) { + psError(PM_ERR_PROG, true, "Cell to receive PSF has %ld readouts (should be only one)", + cell->readouts->n); + return false; + } else { + ro = cell->readouts->data[0]; + } + PM_ASSERT_READOUT_NON_NULL(ro, false); + if (!ro->analysis) { + ro->analysis = psMetadataAlloc(); + } + + if (!pmPSFmodelRead(chip->analysis, ro->analysis, view, file, config)) { psError(psErrorCodeLast(), false, "Failed to write PSF for chip"); return false; @@ -737,6 +775,8 @@ // for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, // and header + image for the PSF residual images -bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { + PS_ASSERT_METADATA_NON_NULL(chipAnalysis, false); + PS_ASSERT_METADATA_NON_NULL(roAnalysis, false); PS_ASSERT_PTR_NON_NULL(view, false); PS_ASSERT_PTR_NON_NULL(file, false); @@ -809,10 +849,10 @@ char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.000"); - psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); - - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); + + psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); if (!regionMD) { regionMD = psMetadataAlloc(); - psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); psFree (regionMD); } @@ -831,5 +871,5 @@ psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); } else { - psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); + psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); for (int i = 0; i < nRegions; i++) { @@ -838,8 +878,8 @@ snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); if (!regionMD) { regionMD = psMetadataAlloc(); - psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); psFree (regionMD); } @@ -1019,20 +1059,20 @@ } - // note that all N plane are implicitly of the same type, so we convert the mask - psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type); - psImageInit (psf->residuals->mask, 0); - psImageInit (psf->residuals->Rx, 0.0); - psImageInit (psf->residuals->Ry, 0.0); - switch (Nz) { - case 1: // Ro only - break; - case 2: // Ro and mask + // note that all N plane are implicitly of the same type, so we convert the mask + psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type); + psImageInit (psf->residuals->mask, 0); + psImageInit (psf->residuals->Rx, 0.0); + psImageInit (psf->residuals->Ry, 0.0); + switch (Nz) { + case 1: // Ro only + break; + case 2: // Ro and mask if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) { psError(psErrorCodeLast(), false, "Unable to read PSF residual image."); return false; } - psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); - break; - case 3: // Ro, Rx and Ry, no mask + psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); + break; + case 3: // Ro, Rx and Ry, no mask if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) { psError(psErrorCodeLast(), false, "Unable to read PSF residual image."); @@ -1043,6 +1083,6 @@ return false; } - break; - case 4: // Ro, Rx, Ry, and mask: + break; + case 4: // Ro, Rx, Ry, and mask: if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) { psError(psErrorCodeLast(), false, "Unable to read PSF residual image."); @@ -1057,11 +1097,11 @@ return false; } - psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); - break; - } - psFree (mask); - } - - psMetadataAdd (analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); + psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); + break; + } + psFree (mask); + } + + psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); psFree (psf); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.h (revision 28439) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.h (revision 28440) @@ -20,5 +20,5 @@ bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, pmConfig *config); bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, pmConfig *config); -bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config); +bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config); bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); @@ -27,5 +27,5 @@ bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); -bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c (revision 28439) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c (revision 28440) @@ -107,12 +107,12 @@ // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y); - psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); - psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); - source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]; - source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0]; - source->psfMag = -2.5*log10(source->psfFlux); + psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); + psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); + source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]; + source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0]; + source->psfMag = -2.5*log10(source->psfFlux); } else { status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF); - source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]); + source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]); } @@ -120,18 +120,18 @@ if (source->modelFits) { bool foundEXT = false; - for (int i = 0; i < source->modelFits->n; i++) { - pmModel *model = source->modelFits->data[i]; - status = pmSourcePhotometryModel (&model->mag, NULL, model); - if (model == source->modelEXT) foundEXT = true; - } - if (foundEXT) { - source->extMag = source->modelEXT->mag; - } else { - status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); - } + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + status = pmSourcePhotometryModel (&model->mag, NULL, model); + if (model == source->modelEXT) foundEXT = true; + } + if (foundEXT) { + source->extMag = source->modelEXT->mag; + } else { + status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); + } } else { - if (source->modelEXT) { - status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); - } + if (source->modelEXT) { + status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); + } } @@ -204,5 +204,5 @@ } if (mode & PM_SOURCE_PHOT_APCORR) { - // XXX this should be removed -- we no longer fit for the 'sky bias' + // XXX this should be removed -- we no longer fit for the 'sky bias' rflux = pow (10.0, 0.4*source->psfMag); source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux; @@ -236,11 +236,11 @@ flux = model->modelFlux (model->params); if (flux > 0) { - mag = -2.5*log10(flux); + mag = -2.5*log10(flux); } if (fitMag) { - *fitMag = mag; + *fitMag = mag; } if (fitFlux) { - *fitFlux = flux; + *fitFlux = flux; } @@ -380,5 +380,5 @@ if (source->diffStats == NULL) { - source->diffStats = pmSourceDiffStatsAlloc(); + source->diffStats = pmSourceDiffStatsAlloc(); } @@ -388,5 +388,5 @@ int nMask = 0; int nBad = 0; - + psImage *flux = source->pixels; psImage *variance = source->variance; @@ -394,28 +394,28 @@ for (int iy = 0; iy < flux->numRows; iy++) { - for (int ix = 0; ix < flux->numCols; ix++) { + for (int ix = 0; ix < flux->numCols; ix++) { if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) { - nMask ++; - continue; - } - - float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - - if (SN > +FLUX_LIMIT) { - nGood ++; - fGood += fabs(flux->data.F32[iy][ix]); - } - - if (SN < -FLUX_LIMIT) { - nBad ++; - fBad += fabs(flux->data.F32[iy][ix]); - } - } + nMask ++; + continue; + } + + float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + + if (SN > +FLUX_LIMIT) { + nGood ++; + fGood += fabs(flux->data.F32[iy][ix]); + } + + if (SN < -FLUX_LIMIT) { + nBad ++; + fBad += fabs(flux->data.F32[iy][ix]); + } + } } source->diffStats->nGood = nGood; - source->diffStats->fRatio = (fGood + fBad == 0.0) ? NAN : fGood / (fGood + fBad); - source->diffStats->nRatioBad = (nGood + nBad == 0) ? NAN : nGood / (float) (nGood + nBad); - source->diffStats->nRatioMask = (nGood + nMask == 0) ? NAN : nGood / (float) (nGood + nMask); + source->diffStats->fRatio = (fGood + fBad == 0.0) ? NAN : fGood / (fGood + fBad); + source->diffStats->nRatioBad = (nGood + nBad == 0) ? NAN : nGood / (float) (nGood + nBad); + source->diffStats->nRatioMask = (nGood + nMask == 0) ? NAN : nGood / (float) (nGood + nMask); source->diffStats->nRatioAll = (nGood + nMask + nBad == 0) ? NAN : nGood / (float) (nGood + nMask + nBad); @@ -628,5 +628,5 @@ } } - model->nPix = Npix; + model->nPix = Npix; model->nDOF = Npix - 1; model->chisq = dC; @@ -636,5 +636,5 @@ -double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor) +double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(Mi, NAN); @@ -652,5 +652,5 @@ for (int yi = 0; yi < Pi->numRows; yi++) { for (int xi = 0; xi < Pi->numCols; xi++) { - if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi]) + if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal) continue; if (!unweighted_sum) { @@ -684,5 +684,5 @@ } -double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor) +double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(Mi, NAN); @@ -727,7 +727,7 @@ for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi]) - continue; - if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj]) + if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal) + continue; + if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal) continue; @@ -746,5 +746,5 @@ } -double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor) +double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(Mi, NAN); @@ -789,7 +789,7 @@ for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi]) - continue; - if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj]) + if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal) + continue; + if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal) continue; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h (revision 28439) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h (revision 28440) @@ -48,5 +48,5 @@ psImage *image, ///< image pixels to be used psImage *mask, ///< mask of pixels to ignore - psImageMaskType maskVal ///< Value to mask + psImageMaskType maskVal ///< Value to mask ); @@ -58,7 +58,7 @@ bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal); -double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor); -double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor); -double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor); +double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); +double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); +double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); // retire these: Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c (revision 28439) +++ /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 28439) +++ /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 28439) +++ /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 28439) +++ /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 28439) +++ /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 28439) +++ /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];