Index: branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c =================================================================== --- branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c (revision 28435) +++ 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 28435) +++ 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 28435) +++ 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 28435) +++ 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: