Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 28010) +++ trunk/psModules/src/objects/Makefile.am (revision 28013) @@ -43,4 +43,5 @@ pmSourceIO_CMF_PS1_V1.c \ pmSourceIO_CMF_PS1_V2.c \ + pmSourceIO_CMF_PS1_SV1.c \ pmSourceIO_CMF_PS1_DV1.c \ pmSourceIO_MatchedRefs.c \ Index: trunk/psModules/src/objects/pmPhotObj.c =================================================================== --- trunk/psModules/src/objects/pmPhotObj.c (revision 28010) +++ trunk/psModules/src/objects/pmPhotObj.c (revision 28013) @@ -58,8 +58,26 @@ if (!object->sources) { object->sources = psArrayAllocEmpty(1); - object->x = source->peak->xf; - object->y = source->peak->yf; + object->x = source->peak->xf; + object->y = source->peak->yf; + object->SN = source->peak->SN; + } else { + object->SN = PS_MAX(object->SN, source->peak->SN); } psArrayAdd (object->sources, 1, source); return true; } + +// sort by SN (descending) +int pmPhotObjSortBySN (const void **a, const void **b) +{ + pmPhotObj *objA = *(pmPhotObj **)a; + pmPhotObj *objB = *(pmPhotObj **)b; + + psF32 fA = objA->SN; + psF32 fB = objB->SN; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} Index: trunk/psModules/src/objects/pmPhotObj.h =================================================================== --- trunk/psModules/src/objects/pmPhotObj.h (revision 28010) +++ trunk/psModules/src/objects/pmPhotObj.h (revision 28013) @@ -39,8 +39,12 @@ float x; float y; + float SN; // max of peak->SN for all matched sources } pmPhotObj; +pmPhotObj *pmPhotObjAlloc(void); + bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source); -pmPhotObj *pmPhotObjAlloc(void); + +int pmPhotObjSortBySN (const void **a, const void **b); /// @} Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 28010) +++ trunk/psModules/src/objects/pmSource.c (revision 28013) @@ -54,4 +54,5 @@ psFree(tmp->moments); psFree(tmp->diffStats); + psFree(tmp->radial); psFree(tmp->blends); psTrace("psModules.objects", 10, "---- end ----\n"); @@ -116,4 +117,5 @@ source->extpars = NULL; source->diffStats = NULL; + source->radial = NULL; source->region = psRegionSet(NAN, NAN, NAN, NAN); Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 28010) +++ trunk/psModules/src/objects/pmSource.h (revision 28013) @@ -91,4 +91,5 @@ pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections + pmSourceRadialApertures *radial; ///< radial flux in circular apertures int imageID; }; Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 28010) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 28013) @@ -65,4 +65,30 @@ } +// pmSourceRadialApertures carries the raw radial flux information, including angular bins +static void pmSourceRadialAperturesFree(pmSourceRadialApertures *radial) +{ + if (!radial) return; + psFree(radial->flux); + psFree(radial->fluxErr); + psFree(radial->fill); +} + +pmSourceRadialApertures *pmSourceRadialAperturesAlloc() +{ + pmSourceRadialApertures *radial = (pmSourceRadialApertures *)psAlloc(sizeof(pmSourceRadialApertures)); + psMemSetDeallocator(radial, (psFreeFunc) pmSourceRadialAperturesFree); + + radial->flux = NULL; + radial->fluxErr = NULL; + radial->fill = NULL; + return radial; +} + +bool psMemCheckSourceRadialApertures(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialAperturesFree); +} + // pmSourceEllipticalFlux carries the elliptical renormalized radial flux info static void pmSourceEllipticalFluxFree(pmSourceEllipticalFlux *flux) Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 28010) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 28013) @@ -20,4 +20,10 @@ psVector *isophotalRadii; // isophotal radius for the above angles } pmSourceRadialFlux; + +typedef struct { + psVector *flux; // fluxes measured at above radii + psVector *fluxErr; // formal error on the fluxes (sqrt\sum(variance)) + psVector *fill; // angles corresponding to above radial profiles +} pmSourceRadialApertures; typedef struct { @@ -62,4 +68,7 @@ bool psMemCheckSourceRadialFlux(psPtr ptr); +pmSourceRadialApertures *pmSourceRadialAperturesAlloc(); +bool psMemCheckSourceRadialApertures(psPtr ptr); + pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc(); bool psMemCheckSourceEllipticalFlux(psPtr ptr); Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 28010) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 28013) @@ -539,4 +539,7 @@ status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); } + if (!strcmp (exttype, "PS1_SV1")) { + status &= pmSourcesWrite_CMF_PS1_SV1 (file->fits, readout, sources, file->header, outhead, dataname, recipe); + } if (!strcmp (exttype, "PS1_DV1")) { status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname); @@ -556,4 +559,7 @@ status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } + if (!strcmp (exttype, "PS1_SV1")) { + status &= pmSourcesWrite_CMF_PS1_SV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); + } if (!strcmp (exttype, "PS1_DV1")) { status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); @@ -573,4 +579,7 @@ status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname); } + if (!strcmp (exttype, "PS1_SV1")) { + status &= pmSourcesWrite_CMF_PS1_SV1_XFIT (file->fits, readout, sources, xfitname); + } if (!strcmp (exttype, "PS1_DV1")) { status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname); @@ -981,5 +990,5 @@ // we need to find the corresponding table EXTNAME. // first check the header - char *extdata = psMetadataLookupStr (NULL, hdu->header, "EXTDATA"); + char *extdata = psMetadataLookupStr (&status, hdu->header, "EXTDATA"); if (extdata) { // if EXTDATA is defined in the header, use that value for 'dataname' @@ -1023,4 +1032,7 @@ if (!strcmp (exttype, "PS1_V2")) { sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); + } + if (!strcmp (exttype, "PS1_SV1")) { + sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header); } if (!strcmp (exttype, "PS1_DV1")) { Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 28010) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 28013) @@ -43,4 +43,8 @@ bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); + bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); @@ -57,4 +61,5 @@ psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28010) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28013) @@ -175,6 +175,6 @@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 28013) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 28013) @@ -0,0 +1,745 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 02:44:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceIO.h" + +// panstarrs-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data + +// NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes +// XXX currently in the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes) + +bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(extname, false); + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + psF32 errMag, chisq, apRadius; + psS32 nPix, nDOF; + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + bool status1 = false; + bool status2 = false; + float magOffset = NAN; + float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); + if (!isfinite(zeropt)) { + zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + } + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + + // if the sequence is defined, write these in seq order; otherwise + // write them in S/N order: + if (sources->n > 0) { + pmSource *source = (pmSource *) sources->data[0]; + if (source->seq == -1) { + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we use this just to define the output vectors (which must be present for all objects) + bool status = false; + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMax, "this must have been defined and tested earlier!"); + psAssert (radMax->n, "this must have been defined and tested earlier!"); + + // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state + // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // If source->seq is -1, source was generated in this analysis. If source->seq is + // not -1, source was read from elsewhere: in the latter case, preserve the source + // ID. source.seq is used instead of source.id since the latter is a const + // generated on Alloc, and would thus be wrong for read in sources. + if (source->seq == -1) { + source->seq = i; + } + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; + } + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + axes = pmPSF_ModelToAxes (PAR, 20.0); + } else { + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + } + chisq = model->chisq; + nDOF = model->nDOF; + nPix = model->nPix; + apRadius = source->apRadius; + errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + chisq = NAN; + nDOF = 0; + nPix = 0; + apRadius = NAN; + errMag = NAN; + } + + float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + psS16 nImageOverlap = 1; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); + + // distinguish moments measure from window vs S/N > XX ?? + float mxx = source->moments ? source->moments->Mxx : NAN; + float mxy = source->moments ? source->moments->Mxy : NAN; + float myy = source->moments ? source->moments->Myy : NAN; + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + + psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVectorInit (radFlux, NAN); + psVectorInit (radFluxErr, NAN); + psVectorInit (radFill, NAN); + if (!source->radial) goto empty_annuli; + if (!source->radial->flux) goto empty_annuli; + if (!source->radial->fill) goto empty_annuli; + psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths"); + + // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux + for (int j = 0; j < source->radial->flux->n; j++) { + radFlux->data.F32[j] = source->radial->flux->data.F32[j]; + radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j]; + radFill->data.F32[j] = source->radial->fill->data.F32[j]; + } + + empty_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFill); + + // XXX not sure how to get this : need to load Nimages with weight? + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + + psArrayAdd (table, 100, row); + psFree (row); + + // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not + // subtracted + + // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image + // edge; 3) any pixels in the 3x3 peak region are masked; + } + + // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? + psMetadata *header = psMetadataCopy(NULL, tableHeader); + pmSourceMasksHeader(header); + + if (table->n == 0) { + if (!psFitsWriteBlank(fits, header, extname)) { + psError(psErrorCodeLast(), false, "Unable to write blank sources file."); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable(fits, header, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // define PSF model type + // XXX need to carry the extra model parameters + int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + assert (modelType > -1); + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numSources = psFitsTableSize(fits); // Number of sources in table + psArray *sources = psArrayAlloc(numSources); // Array of sources, to return + + // convert the table to the pmSource entriesa + for (int i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); + psFree(sources); + return NULL; + } + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); + axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN; + + pmPSF_AxesToModel (PAR, axes); + + float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + source->peak->flux = peakFlux; + source->peak->dx = dPAR[PM_PAR_XPOS]; + source->peak->dy = dPAR[PM_PAR_YPOS]; + source->peak->SN = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N + + source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); + model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); + + source->moments = pmMomentsAlloc (); + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); + source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); + source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); + + source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + assert (status); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +{ + + bool status; + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psF32 xPos, yPos; + psF32 xErr, yErr; + int nRow = -1; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + // zero point corrections + bool status1 = false; + bool status2 = false; + float magOffset = 0.0; + float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); + if (!isfinite(zeropt)) { + zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + } + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // which extended source analyses should we perform? + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMin->n == radMax->n, "inconsistent annular bins"); + + // int nRadialBins = 0; + // if (doAnnuli) { + // // get the max count of radial bins + // for (int i = 0; i < sources->n; i++) { + // pmSource *source = sources->data[i]; + // if (!source->extpars) continue; + // if (!source->extpars->radProfile ) continue; + // if (!source->extpars->radProfile->binSB) continue; + // nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n); + // } + // } + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + // skip source if it is not a ext sourc + // XXX we have two places that extended source parameters are measured: + // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, + // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) + // should we require both? + + pmSource *source = sources->data[i]; + + // skip sources without measurements + if (source->extpars == NULL) continue; + + // we require a PSF model fit (ignore the real crud) + pmModel *model = source->modelPSF; + if (model == NULL) continue; + + // XXX I need to split the extended models from the extended aperture measurements + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); + + float AxialRatio = NAN; + float AxialTheta = NAN; + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + AxialRatio = extpars->axes.minor / extpars->axes.major; + AxialTheta = extpars->axes.theta; + } + psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); + psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); + + // Petrosian measurements + // XXX insert header data: petrosian ref radius, flux ratio + // XXX check flags to see if Pet was measured + if (doPetrosian) { + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + // XXX note that this mag is either calibrated or instrumental depending on existence of zero point + float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point + float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); + } + } + +# if (0) + // Kron measurements + if (doKron) { + pmSourceKronValues *kron = source->extpars->kron; + if (kron) { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); + } + } + + // Isophot measurements + // XXX insert header data: isophotal level + if (doIsophotal) { + pmSourceIsophotalValues *isophot = source->extpars->isophot; + if (isophot) { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); + } + } +# endif + + // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) + if (doAnnuli) { + psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVectorInit (radSB, NAN); + psVectorInit (radFlux, NAN); + psVectorInit (radFill, NAN); + if (!source->extpars) goto empty_annuli; + if (!source->extpars->radProfile) goto empty_annuli; + if (!source->extpars->radProfile->binSB) goto empty_annuli; + psAssert (source->extpars->radProfile->binSum, "programming error"); + psAssert (source->extpars->radProfile->binFill, "programming error"); + psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); + + // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux + for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { + radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; + radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; + radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; + } + + empty_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radSB); + psFree (radFlux); + psFree (radFill); + } + if (nRow < 0) { + nRow = row->list->n; + } else { + psAssert (nRow == row->list->n, "inconsistent row lengths"); + } + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + + return true; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) +{ + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + char name[64]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } + } + + // XXX other parameters which may be set. + // XXX flag / value to define the model + // XXX write out the model type, fit status flags + + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} Index: trunk/psphot/doc/stack.txt =================================================================== --- trunk/psphot/doc/stack.txt (revision 28010) +++ trunk/psphot/doc/stack.txt (revision 28013) @@ -1,2 +1,21 @@ + +20100506: + + we may load RAW and/or CNV images. + + * options->convolve: + * if only one of RAW or CNV exists, convolve it to match target PSF + * if both exist, convolve desired one + + * if matching PSF exists, use it (unless told to re-measure) + + * output goes to OUT (which is then used by psphot analysis routines) + + * detect + + * if (RAW) ? RAW : OUT + + * + 20100503: Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 28010) +++ trunk/psphot/src/Makefile.am (revision 28013) @@ -95,4 +95,10 @@ psphotFitSourcesLinearStack.c \ psphotSourceMatch.c \ + psphotStackMatchPSFs.c \ + psphotStackMatchPSFsUtils.c \ + psphotStackMatchPSFsPrepare.c \ + psphotStackOptions.c \ + psphotStackObjects.c \ + psphotStackPSF.c \ psphotCleanup.c @@ -159,4 +165,5 @@ psphotModelWithPSF.c \ psphotExtendedSourceAnalysis.c \ + psphotExtendedSourceAnalysisByObject.c \ psphotExtendedSourceFits.c \ psphotKernelFromPSF.c \ @@ -186,4 +193,6 @@ psphotEllipticalProfile.c \ psphotRadialBins.c \ + psphotRadialApertures.c \ + psphotRadialAperturesByObject.c \ psphotPetrosian.c \ psphotPetrosianRadialBins.c \ Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 28010) +++ trunk/psphot/src/psphot.h (revision 28013) @@ -37,6 +37,6 @@ bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view); -bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view); -bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotDefineFiles (pmConfig *config, pmFPAfile *input); @@ -51,86 +51,86 @@ // psphotReadout functions -bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view); +bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); -bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view); -bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotModelBackground (pmConfig *config, const pmFPAview *view); -bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filename, int index); - -bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view); -bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass); -bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass); - -bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow); -bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow); - -bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view); -bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); - -bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view); -bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); - -bool psphotRoughClass (pmConfig *config, const pmFPAview *view); -bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotFindDetections (pmConfig *config, const pmFPAview *view, const char *filerule, bool firstPass); +bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool firstPass); + +bool psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow); +bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow); + +bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotRoughClass (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *analysis, psMetadata *recipe, const bool havePSF); -bool psphotImageQuality (pmConfig *config, const pmFPAview *view); -bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view); -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotGuessModels (pmConfig *config, const pmFPAview *view); -bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); - -bool psphotMergeSources (pmConfig *config, const pmFPAview *view); -bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); - -bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final); +bool psphotImageQuality (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotMergeSources (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final); bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final); -bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize); -bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize); - -bool psphotBlendFit (pmConfig *config, const pmFPAview *view); -bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize); +bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool getPSFsize); + +bool psphotBlendFit (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotBlendFit_Threaded (psThreadJob *job); -bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view); -bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotAddNoise (pmConfig *config, const pmFPAview *view); -bool psphotSubNoise (pmConfig *config, const pmFPAview *view); -bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add); -bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add); - -bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view); -bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view); -bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotApResid (pmConfig *config, const pmFPAview *view); -bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); - -bool psphotMagnitudes (pmConfig *config, const pmFPAview *view); +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add); +bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add); + +bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotMagnitudes (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf); bool psphotMagnitudes_Threaded (psThreadJob *job); -bool psphotEfficiency (pmConfig *config, const pmFPAview *view); -bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotEfficiency (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources); bool psphotPSFWeights_Threaded (psThreadJob *job); -bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); -bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); - -bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view); -bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index); +bool psphotSkyReplace (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index); // in psphotSourceStats.c: @@ -147,11 +147,11 @@ // in psphotMergeSources.c: -bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view); +bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, const char *filerule); psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view); -bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view); -bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view); +bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view, const char *filerule); // generate the detection structure for the supplied array of sources -bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources); +bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *sources); // generate the detection structure for the supplied array of sources @@ -349,16 +349,16 @@ bool psphotStackImageLoop (pmConfig *config); bool psphotStackReadout (pmConfig *config, const pmFPAview *view); -bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view); +bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv); bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration const pmFPAview *view, + const char *filename, pmReadout **chiReadout, - char *filename, int index); -bool psphotStackRemoveChisqFromInputs (pmConfig *config); +bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule); bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num); -psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view); -bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index); +psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index); bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS); @@ -367,3 +367,98 @@ int pmPhotObjSortByX (const void **a, const void **b); +typedef enum { + PSPHOT_CNV_SRC_NONE, + PSPHOT_CNV_SRC_AUTO, + PSPHOT_CNV_SRC_CNV, + PSPHOT_CNV_SRC_RAW, +} psphotStackConvolveSource; + +/// Options for stacking process +typedef struct { + // Setup + + int numCols; // size of image (X) + int numRows; // size of image (Y) + + int num; // Number of inputs + bool convolve; // Convolve images? + psphotStackConvolveSource convolveSource; + // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images + + // bool matchZPs; // Adjust relative fluxes based on transparency analysis? + // bool photometry; // Perform photometry? + // psMetadata *stats; // Statistics for output + // FILE *statsFile; // File to which to write statistics + // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images + // psArray *origCovars; // Original covariances matrices + // int quality; // Bad data quality flag + + // Prepare + pmPSF *psf; // Target PSF + psVector *inputSeeing; // Input seeing FWHMs + psVector *inputMask; // Mask for inputs + + float targetSeeing; // Target seeing FWHM + psArray *sourceLists; // Individual lists of sources for matching + psVector *norm; // Normalisation for each image + psArray *psfs; + + // psVector *exposures; // Exposure times + // float sumExposure; // Sum of exposure times + // float zp; // Zero point for output + // psVector *inputMask; // Mask for inputs + // psArray *sources; // Matched sources + + // Convolve + psArray *kernels; // PSF-matching kernels --- required in the stacking + psArray *regions; // PSF-matching regions --- required in the stacking + psVector *matchChi2; // chi^2 for stamps from matching + psVector *weightings; // Combination weightings for images (1/noise^2) + // psArray *cells; // Cells for convolved images --- a handle for reading again + // int numCols, numRows; // Size of image + // psArray *convCovars; // Convolved covariance matrices + + // Combine initial + // pmReadout *outRO; // Output readout + // pmReadout *expRO; // Exposure readout + // psArray *inspect; // Array of arrays of pixels to inspect + + // Rejection + // psArray *rejected; // Rejected pixels +} psphotStackOptions; + +/*** psphotStackMatchPSF prototypes ***/ +bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view); +bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); +bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); + +// psphotStackMatchPSFsUtils +psVector *SetOptWidths (bool *optimum, psMetadata *recipe); +pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize); +bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index); +bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index); +bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index); +bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index); +bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname); +bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname); +bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index); + +pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask); + +psphotStackOptions *psphotStackOptionsAlloc (int num); +psphotStackConvolveSource psphotStackConvolveSourceFromString (const char *string); +pmFPAfile *psphotStackGetConvolveSource (pmConfig *config, psphotStackOptions *options, int index); + +bool psphotCopySources (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc); +bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index); + +bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax); + +bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); +bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); + +bool psphotStackObjectsUnifyPosition (psArray *objects); + #endif Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 28010) +++ trunk/psphot/src/psphotAddNoise.c (revision 28013) @@ -1,14 +1,14 @@ # include "psphotInternal.h" -bool psphotAddNoise (pmConfig *config, const pmFPAview *view) { - return psphotAddOrSubNoise (config, view, true); +bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule) { + return psphotAddOrSubNoise (config, view, filerule, true); } -bool psphotSubNoise (pmConfig *config, const pmFPAview *view) { - return psphotAddOrSubNoise (config, view, false); +bool psphotSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule) { + return psphotAddOrSubNoise (config, view, filerule, false); } // for now, let's store the detections on the readout->analysis for each readout -bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add) +bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add) { bool status = true; @@ -23,6 +23,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for PSPHOT.INPUT entry %d", i); + if (!psphotAddOrSubNoiseReadout (config, view, filerule, i, recipe, add)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for %s entry %d", filerule, i); return false; } @@ -31,5 +31,5 @@ } -bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) { +bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) { bool status = false; @@ -39,5 +39,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 28010) +++ trunk/psphot/src/psphotApResid.c (revision 28013) @@ -5,5 +5,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotApResid (pmConfig *config, const pmFPAview *view) +bool psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -23,6 +23,6 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i); + if (!psphotApResidReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for %s entry %d", filerule, i); return false; } @@ -31,5 +31,5 @@ } -bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) +bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { int Nfail = 0; @@ -43,5 +43,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- trunk/psphot/src/psphotBasicDeblend.c (revision 28010) +++ trunk/psphot/src/psphotBasicDeblend.c (revision 28013) @@ -2,5 +2,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view) +bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -11,6 +11,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on basic deblend analysis for PSPHOT.INPUT entry %d", i); + if (!psphotBasicDeblendReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on basic deblend analysis for %s entry %d", filerule, i); return false; } @@ -19,5 +19,5 @@ } -bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) { +bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { int N; @@ -31,5 +31,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 28010) +++ trunk/psphot/src/psphotBlendFit.c (revision 28013) @@ -2,5 +2,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotBlendFit (pmConfig *config, const pmFPAview *view) +bool psphotBlendFit (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -15,6 +15,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for PSPHOT.INPUT entry %d", i); + if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for %s entry %d", filerule, i); return false; } @@ -24,5 +24,5 @@ // XXX I don't like this name -bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { int Nfit = 0; @@ -35,5 +35,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 28010) +++ trunk/psphot/src/psphotChoosePSF.c (revision 28013) @@ -2,5 +2,5 @@ // generate a PSF model for inputs without PSF models already loaded -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view) +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -20,6 +20,6 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i); + if (!psphotChoosePSFReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for %s entry %d", filerule, i); return false; } @@ -29,5 +29,5 @@ // try PSF models and select best option -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status; @@ -36,5 +36,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotCleanup.c =================================================================== --- trunk/psphot/src/psphotCleanup.c (revision 28010) +++ trunk/psphot/src/psphotCleanup.c (revision 28013) @@ -19,4 +19,5 @@ pmConceptsDone (); pmConfigDone (); + psLibFinalize(); // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot"); Index: trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- trunk/psphot/src/psphotDeblendSatstars.c (revision 28010) +++ trunk/psphot/src/psphotDeblendSatstars.c (revision 28013) @@ -2,5 +2,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view) +bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -11,6 +11,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + if (!psphotDeblendSatstarsReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for %s entry %d", filerule, i); return false; } @@ -19,5 +19,5 @@ } -bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) { +bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { int N; @@ -31,5 +31,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotEfficiency.c =================================================================== --- trunk/psphot/src/psphotEfficiency.c (revision 28010) +++ trunk/psphot/src/psphotEfficiency.c (revision 28013) @@ -156,5 +156,5 @@ } -bool psphotEfficiency (pmConfig *config, const pmFPAview *view) +bool psphotEfficiency (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -173,7 +173,7 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (i == chisqNum) continue; // skip chisq image - if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i); + if (i == chisqNum) continue; // skip chisq image + if (!psphotEfficiencyReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for %s entry %d", filerule, i); return false; } @@ -183,5 +183,5 @@ // Determine detection efficiency -bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) +bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status = true; @@ -190,5 +190,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotErrorCodes.dat =================================================================== --- trunk/psphot/src/psphotErrorCodes.dat (revision 28010) +++ trunk/psphot/src/psphotErrorCodes.dat (revision 28013) @@ -11,4 +11,5 @@ APERTURE Problem with aperture photometry SKY Problem in sky determination +IO Problem in data I/O # these errors correspond to standard exit conditions ARGUMENTS Incorrect arguments Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 28010) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 28013) @@ -1,6 +1,12 @@ # include "psphotInternal.h" +// ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis +// XXX option to choose a consistent position +// XXX option to choose a consistent elliptical contour +// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands +// XXX consistent choice of extendedness... + // for now, let's store the detections on the readout->analysis for each readout -bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view) +bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -21,6 +27,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotExtendedSourceAnalysisReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for PSPHOT.INPUT entry %d", i); + if (!psphotExtendedSourceAnalysisReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i); return false; } @@ -30,5 +36,5 @@ // aperture-like measurements for extended sources -bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status; @@ -42,5 +48,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c (revision 28013) +++ trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c (revision 28013) @@ -0,0 +1,158 @@ +# include "psphotInternal.h" + +// XXX option to choose a consistent elliptical contour +// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands + +// aperture-like measurements for extended sources +bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) { + + bool status; + int Next = 0; + int Npetro = 0; + int Nannuli = 0; + + psTimerStart ("psphot.extended"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // XXX require petrosian analysis for non-linear fits? + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + + // which extended source analyses should we perform? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + // number of images used to define sources + int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // generate look-up arrays for readouts + psArray *readouts = psArrayAlloc(nImages); + for (int i = 0; i < nImages; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + readouts->data[i] = psMemIncrRefCounter(readout); + } + + // source analysis is done in S/N order (brightest first) + objects = psArraySort (objects, pmPhotObjSortBySN); + + // process the objects in order. + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + // we need to decide for an object if we are going to measure all sources or not + // simple rule : if *any* of the sources would be measured, measure the object + + // choose the sources of interest + bool measureSource = false; + for (int j = 0; !measureSource && (j < object->sources->n); j++) { + + pmSource *source = object->sources->data[j]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) continue; + measureSource = true; + } + if (!measureSource) continue; + + // choose the sources of interest + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + psAssert (source, "programming error"); // all entries in object->sources must exist, right? + psAssert (source->peak, "programming error"); // how can a source not have a peak? + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Next ++; + + int index = source->imageID; + pmReadout *readout = readouts->data[index]; + + // force source image to be a bit larger... + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + + // if we request any of these measurements, we require the radial profile + if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { + // all measurements below require the radial profile; skip them all + // re-subtract the object, leave local sky + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + continue; + } + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + } + + // Petrosian Mags + if (doPetrosian) { + if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { + psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + Npetro ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; + } + } + + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + if (source->extpars) { + psFree(source->extpars->radFlux); + psFree(source->extpars->ellipticalFlux); + } + } + } + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + + psFree(readouts); + return true; +} Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 28010) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 28013) @@ -2,5 +2,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view) +bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -21,6 +21,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotExtendedSourceFitsReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for PSPHOT.INPUT entry %d", i); + if (!psphotExtendedSourceFitsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for %s entry %d", filerule, i); return false; } @@ -30,5 +30,5 @@ // non-linear model fitting for extended sources -bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status; @@ -41,5 +41,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 28010) +++ trunk/psphot/src/psphotFindDetections.c (revision 28013) @@ -4,5 +4,5 @@ // peaks and new footprints. any old peaks are saved on oldPeaks. the resulting footprint set // contains all footprints (old and new) -bool psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass) +bool psphotFindDetections (pmConfig *config, const pmFPAview *view, const char *filerule, bool firstPass) { bool status = true; @@ -17,6 +17,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, recipe, firstPass)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); + if (!psphotFindDetectionsReadout (config, view, filerule, i, recipe, firstPass)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i); return false; } @@ -26,5 +26,5 @@ // smooth the image, search for peaks, optionally define footprints based on the peaks -bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass) { +bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool firstPass) { bool status; @@ -34,5 +34,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 28010) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 28013) @@ -13,5 +13,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final) +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final) { bool status = true; @@ -28,5 +28,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest psAssert (file, "missing file?"); @@ -44,5 +44,5 @@ if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i); + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; } Index: trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 28010) +++ trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 28013) @@ -28,5 +28,4 @@ // analysis is done in spatial order (to speed up overlap search) - // sort by first element in each source list objects = psArraySort (objects, pmPhotObjSortByX); Index: trunk/psphot/src/psphotForcedReadout.c =================================================================== --- trunk/psphot/src/psphotForcedReadout.c (revision 28010) +++ trunk/psphot/src/psphotForcedReadout.c (revision 28013) @@ -20,5 +20,5 @@ // set the photcode for this image - if (!psphotAddPhotcode (config, view)) { + if (!psphotAddPhotcode (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -30,18 +30,18 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view); + psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT"); if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // generate a background model (median, smoothed image) - if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, view); + if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } - if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, view); + if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } @@ -49,19 +49,19 @@ // this only happens if we had a programming error in psphotLoadPSF psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // include externally-supplied sources - psphotLoadExtSources (config, view); + psphotLoadExtSources (config, view, "PSPHOT.INPUT"); // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, view); + psphotGuessModels (config, view, "PSPHOT.INPUT"); // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - psphotMergeSources (config, view); + psphotMergeSources (config, view, "PSPHOT.INPUT"); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (config, view, false); + psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false); // identify CRs and extended sources @@ -71,5 +71,5 @@ // calculate source magnitudes - psphotMagnitudes(config, view); + psphotMagnitudes(config, view, "PSPHOT.INPUT"); // XXX do I want to do this? @@ -80,10 +80,10 @@ // replace background in residual image - psphotSkyReplace (config, view); + psphotSkyReplace (config, view, "PSPHOT.INPUT"); // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view); + psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 28010) +++ trunk/psphot/src/psphotGuessModels.c (revision 28013) @@ -8,5 +8,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotGuessModels (pmConfig *config, const pmFPAview *view) +bool psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -22,6 +22,6 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i); + if (!psphotGuessModelsReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for %s entry %d", filerule, i); return false; } @@ -31,5 +31,5 @@ // construct an initial PSF model for each object (new sources only) -bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { +bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) { bool status; @@ -38,5 +38,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotImageQuality.c =================================================================== --- trunk/psphot/src/psphotImageQuality.c (revision 28010) +++ trunk/psphot/src/psphotImageQuality.c (revision 28013) @@ -5,5 +5,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotImageQuality (pmConfig *config, const pmFPAview *view) +bool psphotImageQuality (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -23,6 +23,6 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i); + if (!psphotImageQualityReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for %s entry %d", filerule, i); return false; } @@ -32,10 +32,10 @@ // selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4 -bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) +bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status = true; // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 28010) +++ trunk/psphot/src/psphotMagnitudes.c (revision 28013) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotMagnitudes (pmConfig *config, const pmFPAview *view) +bool psphotMagnitudes (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -21,5 +21,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest psAssert (file, "missing file?"); @@ -37,5 +37,5 @@ if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i); + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; } Index: trunk/psphot/src/psphotMakePSFReadout.c =================================================================== --- trunk/psphot/src/psphotMakePSFReadout.c (revision 28010) +++ trunk/psphot/src/psphotMakePSFReadout.c (revision 28013) @@ -19,5 +19,5 @@ // set the photcode for this image - if (!psphotAddPhotcode (config, view)) { + if (!psphotAddPhotcode (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -29,21 +29,21 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view); + psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT"); if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // generate a background model (median, smoothed image) - if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, view); + if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } - if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, view); + if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } - psphotLoadExtSources (config, view); + psphotLoadExtSources (config, view, "PSPHOT.INPUT"); // If sources have been supplied, then these should be used to measure the PSF include @@ -53,24 +53,24 @@ // a text file have no valid SN values, but psphotChoosePSF needs to select the top // PSF_MAX_NSTARS to generate the PSF. - if (!psphotCheckExtSources (config, view)) { + if (!psphotCheckExtSources (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failure to select possible PSF sources (external or internal)"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // Use bright stellar objects to measure PSF. If we do not have enough stars to generate // the PSF, build one from the SEEING guess and model class - if (!psphotChoosePSF (config, view)) { + if (!psphotChoosePSF (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // measure aperture photometry corrections # if 0 - if (!psphotApResid (config, view)) { + if (!psphotApResid (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } # endif - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } Index: trunk/psphot/src/psphotMaskReadout.c =================================================================== --- trunk/psphot/src/psphotMaskReadout.c (revision 28010) +++ trunk/psphot/src/psphotMaskReadout.c (revision 28013) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view) { +bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = false; @@ -16,5 +16,5 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotSetMaskAndVarianceReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + if (!psphotSetMaskAndVarianceReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i); return false; @@ -25,9 +25,9 @@ // generate mask and variance if not defined, additional mask for restricted subregion -bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status; - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotMergeSources.c =================================================================== --- trunk/psphot/src/psphotMergeSources.c (revision 28010) +++ trunk/psphot/src/psphotMergeSources.c (revision 28013) @@ -6,5 +6,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotMergeSources (pmConfig *config, const pmFPAview *view) +bool psphotMergeSources (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -15,6 +15,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotMergeSourcesReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i); + if (!psphotMergeSourcesReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for %s entry %d", filerule, i); return false; } @@ -24,10 +24,10 @@ // add newly selected sources to the existing list of sources -bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { +bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) { bool status; // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); @@ -71,5 +71,5 @@ // only expect a single entry for PSPHOT.INPUT.CMF and PSPHOT.SOURCES.TEXT, so we can only // associate input sources with a single entry for PSPHOT.INPUT -bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view) { +bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status; @@ -79,5 +79,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); @@ -130,5 +130,5 @@ // load data from input TXT file: { - pmChip *chipTXT = pmFPAfileThisChip (config->files, view, "PSPHOT.INPUT"); + pmChip *chipTXT = pmFPAfileThisChip (config->files, view, filerule); if (!chipTXT) goto finish; @@ -167,5 +167,5 @@ // extract the input sources corresponding to this readout -// XXX this function needs to be updated to work with the new context of pshot inputs +// XXX this function needs to be updated to work with the new context of psphot inputs psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view) { @@ -197,8 +197,8 @@ // psphotDetectionsFromSources to psphotSourceStats and are now stored on // detections->newSources. -bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view) { - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest +bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view, const char *filerule) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest psAssert (file, "missing file?"); @@ -227,8 +227,8 @@ // generate the detection structure for the supplied array of sources // XXX this currently assumes there is a single input file -bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources) { - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest +bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *sources) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest psAssert (file, "missing file?"); @@ -335,5 +335,5 @@ } -bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view) { +bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status; @@ -343,5 +343,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest psAssert (file, "missing file?"); @@ -373,28 +373,84 @@ // find the detections (by peak and/or footprint) in the image. - if (!psphotFindDetections (config, view, true)) { + if (!psphotFindDetections (config, view, filerule, true)) { psError(PSPHOT_ERR_CONFIG, false, "unable to find detections in this image"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, filerule); } // construct sources and measure basic stats - psphotSourceStats (config, view, true); + psphotSourceStats (config, view, filerule, true); // find blended neighbors of very saturated stars - psphotDeblendSatstars (config, view); + psphotDeblendSatstars (config, view, filerule); // mark blended peaks PS_SOURCE_BLEND - if (!psphotBasicDeblend (config, view)) { + if (!psphotBasicDeblend (config, view, filerule)) { psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, filerule); } // classify sources based on moments, brightness - if (!psphotRoughClass (config, view)) { + if (!psphotRoughClass (config, view, filerule)) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, view); - } - } - - return true; -} + return psphotReadoutCleanup (config, view, filerule); + } + } + + return true; +} + +// copy the detections from one pmFPAfile to another +bool psphotCopySources (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // skip the chisq image because it is a duplicate of the detection version + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotCopySourcesReadout (config, view, ruleOut, ruleSrc, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to copy sources from %s to %s entry %d", ruleSrc, ruleOut, i); + return false; + } + } + return true; +} + +// add newly selected sources to the existing list of sources +bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { + + bool status; + + // find the currently selected readout + pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest + psAssert (fileSrc, "missing file?"); + + pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa); + psAssert (readoutSrc, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + // find the currently selected readout + pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest + psAssert (fileOut, "missing file?"); + + pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa); + psAssert (readoutOut, "missing readout?"); + + // save detections on the readout->analysis + // XXX this replaced any existing entry; allow this operation to merge? + if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + + return true; +} + Index: trunk/psphot/src/psphotModelBackground.c =================================================================== --- trunk/psphot/src/psphotModelBackground.c (revision 28010) +++ trunk/psphot/src/psphotModelBackground.c (revision 28013) @@ -384,5 +384,5 @@ // XXX supply filename or keep PSPHOT.INPUT fixed? -bool psphotModelBackground (pmConfig *config, const pmFPAview *view) +bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = false; @@ -393,5 +393,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotModelBackgroundReadoutFileIndex(config, view, "PSPHOT.INPUT", i)) { + if (!psphotModelBackgroundReadoutFileIndex(config, view, filerule, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); return false; Index: trunk/psphot/src/psphotModelTest.c =================================================================== --- trunk/psphot/src/psphotModelTest.c (revision 28010) +++ trunk/psphot/src/psphotModelTest.c (revision 28013) @@ -3,5 +3,5 @@ // XXX add more test information? -bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe) { +bool psphotModelTest (pmConfig *config, const pmFPAview *view, const char *filerule, psMetadata *recipe) { bool status; @@ -33,5 +33,5 @@ // use poissonian errors or local-sky errors - bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule); if (!status) POISSON_ERRORS = true; pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS); Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 28010) +++ trunk/psphot/src/psphotOutput.c (revision 28013) @@ -126,5 +126,5 @@ } -bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view) { +bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = false; @@ -135,5 +135,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotAddPhotcodeReadout (config, view, "PSPHOT.INPUT", i)) { + if (!psphotAddPhotcodeReadout (config, view, filerule, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to add photcode to PSPHOT.INPUT entry %d", i); return false; Index: trunk/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- trunk/psphot/src/psphotPetrosianRadialBins.c (revision 28010) +++ trunk/psphot/src/psphotPetrosianRadialBins.c (revision 28013) @@ -182,5 +182,5 @@ psFree(values); psFree(stats); - return false; + return true; } Index: trunk/psphot/src/psphotPetrosianStats.c =================================================================== --- trunk/psphot/src/psphotPetrosianStats.c (revision 28010) +++ trunk/psphot/src/psphotPetrosianStats.c (revision 28013) @@ -15,4 +15,9 @@ pmSourceRadialProfile *profile = source->extpars->petProfile; + + if (!profile->binSB) { + psLogMsg ("psphot", PS_LOG_DETAIL, "no petrosian profile, skipping source %f, %f", source->peak->xf, source->peak->yf); + return true; + } psVector *binSB = profile->binSB; @@ -190,5 +195,5 @@ source->extpars->petrosianR90Err = NAN; - fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf); + // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf); psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); Index: trunk/psphot/src/psphotRadialApertures.c =================================================================== --- trunk/psphot/src/psphotRadialApertures.c (revision 28013) +++ trunk/psphot/src/psphotRadialApertures.c (revision 28013) @@ -0,0 +1,236 @@ +# include "psphotInternal.h" + +bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar); + +// for now, let's store the detections on the readout->analysis for each readout +bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n"); + return true; + } + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +// aperture-like measurements for extended sources +// flux in simple, circular apertures +bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + int Nradial = 0; + + psTimerStart ("psphot.radial"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + // radMax stores the upper bounds of the annuli + // XXX keep the same name here as for the petrosian / elliptical apertures? + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe"); + psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); + + // source analysis is done in S/N order (brightest first) + // XXX are we getting the objects out of order? does it matter? + sources = psArraySort (sources, pmSourceSortBySN); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) continue; + + // limit selection by analysis region + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.y1) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Nradial ++; + + // force source image to be a bit larger... + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + + if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) { + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + } + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial); + return true; +} + +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax) { + + // allocate pmSourceExtendedParameters, if not already defined + if (!source->radial) { + source->radial = pmSourceRadialAperturesAlloc (); + } + + psVector *radius = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *pixVar = psVectorAllocEmpty(100, PS_TYPE_F32); + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + + // 0.5 PIX: get radius as a function of pixel coord + float x = ix + 0.5 - source->peak->xf + source->pixels->col0; + float y = iy + 0.5 - source->peak->yf + source->pixels->row0; + + float r = hypot(x, y); + + psVectorAppend(radius, r); + psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]); + psVectorAppend(pixVar, source->variance->data.F32[iy][ix]); + } + } + psphotRadialAperturesSortFlux(radius, pixFlux, pixVar); + + psVector *flux = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *fluxErr = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *fill = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin + + psVectorInit (flux, 0.0); + psVectorInit (fluxErr, 0.0); + psVectorInit (fill, 0.0); + + float fluxSum = 0.0; + float varSum = 0.0; + int nPixSum = 0; + + bool done = false; + int nOut = 0; + float Rmax = radMax->data.F32[nOut]; + + // XXX assume (or enforce) that the bins are contiguous and non-overlapping (Rmax[i] = Rmin[i+1]) + for (int i = 0; !done && (i < radius->n); i++) { + if (radius->data.F32[i] > Rmax) { + // calculate the total flux for bin 'nOut' + float Area = M_PI*PS_SQR(Rmax); + flux->data.F32[nOut] = fluxSum; + fluxErr->data.F32[nOut] = sqrt(varSum); + fill->data.F32[nOut] = nPixSum / Area; + + psTrace ("psphot", 5, "radial bins: %3d %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", + nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area); + + nOut ++; + if (nOut >= radMax->n) break; + Rmax = radMax->data.F32[nOut]; + } + fluxSum += pixFlux->data.F32[i]; + varSum += pixVar->data.F32[i]; + nPixSum ++; + } + flux->n = fluxErr->n = fill->n = nOut; + + psFree(source->radial->flux); + psFree(source->radial->fluxErr); + psFree(source->radial->fill); + + source->radial->flux = flux; + source->radial->fluxErr = fluxErr; + source->radial->fill = fill; + + psFree (radius); + psFree (pixFlux); + psFree (pixVar); + + return true; +} + +// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors +# define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B]) +# define SWAP_VECT(TYPE,A,B) { \ + float tmp; \ + if (A != B) { \ + tmp = radius->data.F32[A]; \ + radius->data.F32[A] = radius->data.F32[B]; \ + radius->data.F32[B] = tmp; \ + tmp = pixFlux->data.F32[A]; \ + pixFlux->data.F32[A] = pixFlux->data.F32[B]; \ + pixFlux->data.F32[B] = tmp; \ + tmp = pixVar->data.F32[A]; \ + pixVar->data.F32[A] = pixVar->data.F32[B]; \ + pixVar->data.F32[B] = tmp; \ + } \ +} + +bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) { + + // sort the vector set by the radius + PSSORT (radius->n, COMPARE_VECT, SWAP_VECT, NONE); + return true; +} + Index: trunk/psphot/src/psphotRadialAperturesByObject.c =================================================================== --- trunk/psphot/src/psphotRadialAperturesByObject.c (revision 28013) +++ trunk/psphot/src/psphotRadialAperturesByObject.c (revision 28013) @@ -0,0 +1,111 @@ +# include "psphotInternal.h" + +// aperture-like measurements for extended sources +// flux in simple, circular apertures +bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) { + + bool status; + int Nradial = 0; + + psTimerStart ("psphot.radial"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n"); + return true; + } + + // number of images used to define sources + int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // radMax stores the upper bounds of the annuli + // XXX keep the same name here as for the petrosian / elliptical apertures? + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe"); + psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); + + // source analysis is done in S/N order (brightest first) + objects = psArraySort (objects, pmPhotObjSortBySN); + + // generate look-up arrays for readouts + psArray *readouts = psArrayAlloc(nImages); + for (int i = 0; i < nImages; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + readouts->data[i] = psMemIncrRefCounter(readout); + } + + // process the objects in order. + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + // choose the sources of interest + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Nradial ++; + + int index = source->imageID; + pmReadout *readout = readouts->data[index]; + + // force source image to be a bit larger... + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + + if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) { + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + } + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + } + + psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial); + + psFree(readouts); + return true; +} Index: trunk/psphot/src/psphotRadialBins.c =================================================================== --- trunk/psphot/src/psphotRadialBins.c (revision 28010) +++ trunk/psphot/src/psphotRadialBins.c (revision 28013) @@ -44,6 +44,12 @@ psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - if (!radMin || !radMin->n) return false; - if (!radMax || !radMax->n) return false; + if (!radMin || !radMin->n) { + psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)"); + return false; + } + if (!radMax || !radMax->n) { + psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)"); + return false; + } psVector *binSB = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin @@ -139,5 +145,5 @@ } } - binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut; + binSB->n = binSBstdev->n = binSum->n = binRad->n = binArea->n = nOut; // interpolate any bins that were empty (extrapolate to center if needed) @@ -155,5 +161,5 @@ psFree(values); psFree(stats); - return false; + return true; } Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 28010) +++ trunk/psphot/src/psphotReadout.c (revision 28013) @@ -28,5 +28,5 @@ // set the photcode for this image - if (!psphotAddPhotcode (config, view)) { + if (!psphotAddPhotcode (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -34,99 +34,99 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotSetMaskAndVariance (config, view)) { - return psphotReadoutCleanup(config, view); + if (!psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup(config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // generate a background model (median, smoothed image) - if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, view); - } - if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, view); + if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); + } + if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are determined and saved on // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT - if (!psphotLoadPSF (config, view)) { + if (!psphotLoadPSF (config, view)) { // ??? need to supply 2 ? psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // find the detections (by peak and/or footprint) in the image. - if (!psphotFindDetections (config, view, true)) { // pass 1 + if (!psphotFindDetections (config, view, "PSPHOT.INPUT", true)) { // pass 1 // this only happens if we had an error in psphotFindDetections psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // construct sources and measure basic stats (saved on detections->newSources) - if (!psphotSourceStats (config, view, true)) { // pass 1 + if (!psphotSourceStats (config, view, "PSPHOT.INPUT", true)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // find blended neighbors of very saturated stars (detections->newSources) - if (!psphotDeblendSatstars (config, view)) { + if (!psphotDeblendSatstars (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // mark blended peaks PS_SOURCE_BLEND (detections->newSources) - if (!psphotBasicDeblend (config, view)) { + if (!psphotBasicDeblend (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // classify sources based on moments, brightness. if a PSF model has been loaded, the PSF // clump defined for it is used not measured (detections->newSources) - if (!psphotRoughClass (config, view)) { // pass 1 + if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { // pass 1 psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) - if (!psphotImageQuality (config, view)) { // pass 1 + if (!psphotImageQuality (config, view, "PSPHOT.INPUT")) { // pass 1 psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, // this step is skipped - if (!psphotChoosePSF (config, view)) { // pass 1 + if (!psphotChoosePSF (config, view, "PSPHOT.INPUT")) { // pass 1 psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // include externally-supplied sources // XXX fix this in the new multi-input context - // psphotLoadExtSources (config, view); // pass 1 + // psphotLoadExtSources (config, view, "PSPHOT.INPUT"); // pass 1 // construct an initial model for each object, set the radius to fitRadius, set circular // fit mask (detections->newSources) - psphotGuessModels (config, view); // pass 1 + psphotGuessModels (config, view, "PSPHOT.INPUT"); // pass 1 // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - psphotMergeSources (config, view); + psphotMergeSources (config, view, "PSPHOT.INPUT"); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (config, view, false); // pass 1 (detections->allSources) + psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false); // pass 1 (detections->allSources) // identify CRs and extended sources (only unmeasured sources are measured) - psphotSourceSize (config, view, true); // pass 1 (detections->allSources) + psphotSourceSize (config, view, "PSPHOT.INPUT", true); // pass 1 (detections->allSources) if (!strcasecmp (breakPt, "ENSEMBLE")) { goto finish; @@ -135,12 +135,12 @@ // non-linear PSF and EXT fit to brighter sources // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) - psphotBlendFit (config, view); // pass 1 (detections->allSources) + psphotBlendFit (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) // replace all sources - psphotReplaceAllSources (config, view); // pass 1 (detections->allSources) + psphotReplaceAllSources (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) // linear fit to include all sources (subtract again) // NOTE : apply to ALL sources (extended + psf) - psphotFitSourcesLinear (config, view, true); // pass 2 (detections->allSources) + psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", true); // pass 2 (detections->allSources) // if we only do one pass, skip to extended source analysis @@ -150,40 +150,40 @@ // add noise for subtracted objects - psphotAddNoise (config, view); // pass 1 (detections->allSources) + psphotAddNoise (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) // find fainter sources // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections - psphotFindDetections (config, view, false); // pass 2 (detections->peaks, detections->footprints) + psphotFindDetections (config, view, "PSPHOT.INPUT", false); // pass 2 (detections->peaks, detections->footprints) // remove noise for subtracted objects (ie, return to normal noise level) // NOTE: this needs to operate only on the OLD sources - psphotSubNoise (config, view); // pass 1 (detections->allSources) + psphotSubNoise (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) // define new sources based on only the new peaks // NOTE: new sources are saved on detections->newSources - psphotSourceStats (config, view, false); // pass 2 (detections->newSources) + psphotSourceStats (config, view, "PSPHOT.INPUT", false); // pass 2 (detections->newSources) // set source type // NOTE: apply only to detections->newSources - if (!psphotRoughClass (config, view)) { // pass 2 (detections->newSources) + if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { // pass 2 (detections->newSources) psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // create full input models, set the radius to fitRadius, set circular fit mask // NOTE: apply only to detections->newSources - psphotGuessModels (config, view); // pass 2 (detections->newSources) + psphotGuessModels (config, view, "PSPHOT.INPUT"); // pass 2 (detections->newSources) // replace all sources so fit below applies to all at once // NOTE: apply only to OLD sources (which have been subtracted) - psphotReplaceAllSources (config, view); // pass 2 + psphotReplaceAllSources (config, view, "PSPHOT.INPUT"); // pass 2 // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW // XXX check on free of sources... - psphotMergeSources (config, view); // (detections->newSources + detections->allSources -> detections->allSources) + psphotMergeSources (config, view, "PSPHOT.INPUT"); // (detections->newSources + detections->allSources -> detections->allSources) // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, true); // pass 3 (detections->allSources) + psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", true); // pass 3 (detections->allSources) pass1finish: @@ -191,8 +191,8 @@ // measure source size for the remaining sources // NOTE: applies only to NEW (unmeasured) sources - psphotSourceSize (config, view, false); // pass 2 (detections->allSources) - - psphotExtendedSourceAnalysis (config, view); // pass 1 (detections->allSources) - psphotExtendedSourceFits (config, view); // pass 1 (detections->allSources) + psphotSourceSize (config, view, "PSPHOT.INPUT", false); // pass 2 (detections->allSources) + + psphotExtendedSourceAnalysis (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) + psphotExtendedSourceFits (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) finish: @@ -202,13 +202,13 @@ // measure aperture photometry corrections - if (!psphotApResid (config, view)) { // pass 1 (detections->allSources) + if (!psphotApResid (config, view, "PSPHOT.INPUT")) { // pass 1 (detections->allSources) psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // calculate source magnitudes - psphotMagnitudes(config, view); // pass 1 (detections->allSources) - - if (!psphotEfficiency(config, view)) { // pass 1 + psphotMagnitudes(config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) + + if (!psphotEfficiency(config, view, "PSPHOT.INPUT")) { // pass 1 psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -219,10 +219,10 @@ // replace background in residual image - psphotSkyReplace (config, view); // pass 1 + psphotSkyReplace (config, view, "PSPHOT.INPUT"); // pass 1 // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view); // pass 1 + psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // pass 1 // create the exported-metadata and free local data - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup(config, view, "PSPHOT.INPUT"); } Index: trunk/psphot/src/psphotReadoutCleanup.c =================================================================== --- trunk/psphot/src/psphotReadoutCleanup.c (revision 28010) +++ trunk/psphot/src/psphotReadoutCleanup.c (revision 28013) @@ -2,5 +2,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view) +bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -24,6 +24,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotReadoutCleanupReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for PSPHOT.INPUT entry %d", i); + if (!psphotReadoutCleanupReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for %s entry %d", filerule, i); return false; } @@ -39,10 +39,10 @@ // not a DATA error, then there was a serious problem. Only in this case, or if the fail // on the stats measurement, do we return false -bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status = true; // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotReadoutFindPSF.c =================================================================== --- trunk/psphot/src/psphotReadoutFindPSF.c (revision 28010) +++ trunk/psphot/src/psphotReadoutFindPSF.c (revision 28013) @@ -8,5 +8,5 @@ // set the photcode for the PSPHOT.INPUT - if (!psphotAddPhotcode(config, view)) { + if (!psphotAddPhotcode(config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -14,5 +14,5 @@ // Generate the mask and variance images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view); + psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT"); // Note that in this implementation, we do NOT model the background and we do not @@ -21,11 +21,11 @@ // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF) // XXX we assume a single set of input sources is supplied - if (!psphotDetectionsFromSources (config, view, inSources)) { + if (!psphotDetectionsFromSources (config, view, "PSPHOT.INPUT", inSources)) { psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // construct detections->newSources and measure basic stats (moments, local sky) - if (!psphotSourceStats(config, view, true)) { + if (!psphotSourceStats(config, view, "PSPHOT.INPUT", true)) { psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); return false; @@ -33,5 +33,5 @@ // peak flux is wrong : use the peak measured in the moments analysis: - if (!psphotRepairLoadedSources(config, view)) { + if (!psphotRepairLoadedSources(config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources"); return false; @@ -39,17 +39,17 @@ // classify sources based on moments, brightness (psf is not known) - if (!psphotRoughClass (config, view)) { + if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } - if (!psphotImageQuality (config, view)) { + if (!psphotImageQuality (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } - if (!psphotChoosePSF(config, view)) { + if (!psphotChoosePSF(config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } @@ -59,23 +59,23 @@ // fits from that analysis, or run the linear PSF fit for all objects currently in hand // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, view); + psphotGuessModels (config, view, "PSPHOT.INPUT"); # endif // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - psphotMergeSources (config, view); + psphotMergeSources (config, view, "PSPHOT.INPUT"); # if 0 // measure aperture photometry corrections - if (!psphotApResid (config, view)) { + if (!psphotApResid (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } # endif // drop the references to the image pixels held by each source - psphotSourceFreePixels(config, view); + psphotSourceFreePixels(config, view, "PSPHOT.INPUT"); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } Index: trunk/psphot/src/psphotReadoutKnownSources.c =================================================================== --- trunk/psphot/src/psphotReadoutKnownSources.c (revision 28010) +++ trunk/psphot/src/psphotReadoutKnownSources.c (revision 28013) @@ -8,5 +8,5 @@ // set the photcode for this image - if (!psphotAddPhotcode(config, view)) { + if (!psphotAddPhotcode(config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -14,5 +14,5 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view); + psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT"); // Note that in this implementation, we do NOT model the background and we do not @@ -20,11 +20,11 @@ // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF) - if (!psphotDetectionsFromSources (config, view, inSources)) { + if (!psphotDetectionsFromSources (config, view, "PSPHOT.INPUT", inSources)) { psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // construct sources and measure basic stats - if (!psphotSourceStats (config, view, true)) { + if (!psphotSourceStats (config, view, "PSPHOT.INPUT", true)) { psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); return false; @@ -32,5 +32,5 @@ // peak flux is wrong : use the peak measured in the moments analysis: - if (!psphotRepairLoadedSources(config, view)) { + if (!psphotRepairLoadedSources(config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources"); return false; @@ -38,37 +38,37 @@ // classify sources based on moments, brightness (psf is not known) - if (!psphotRoughClass (config, view)) { + if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } - if (!psphotChoosePSF (config, view)) { + if (!psphotChoosePSF (config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // construct an initial model for each object - psphotGuessModels (config, view); + psphotGuessModels (config, view, "PSPHOT.INPUT"); // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - psphotMergeSources (config, view); + psphotMergeSources (config, view, "PSPHOT.INPUT"); // linear PSF fit to source peaks - psphotFitSourcesLinear (config, view, false); + psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false); // measure aperture photometry corrections - if (!psphotApResid (config, view)) { + if (!psphotApResid (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // calculate source magnitudes - psphotMagnitudes(config, view); + psphotMagnitudes(config, view, "PSPHOT.INPUT"); // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view); + psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } Index: trunk/psphot/src/psphotReadoutMinimal.c =================================================================== --- trunk/psphot/src/psphotReadoutMinimal.c (revision 28010) +++ trunk/psphot/src/psphotReadoutMinimal.c (revision 28013) @@ -18,5 +18,5 @@ // set the photcode for this image - if (!psphotAddPhotcode(config, view)) { + if (!psphotAddPhotcode(config, view, "PSPHOT.INPUT")) { psError(PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -24,63 +24,63 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view); + psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT"); // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved on readout->analysis if (!psphotLoadPSF (config, view)) { psError (PSPHOT_ERR_CONFIG, false, "missing psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // find the detections (by peak and/or footprint) in the image. (final pass) - if (!psphotFindDetections(config, view, false)) { + if (!psphotFindDetections(config, view, "PSPHOT.INPUT", false)) { psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // construct sources and measure basic stats (saved on detections->newSources) - if (!psphotSourceStats (config, view, false)) { // pass 1 + if (!psphotSourceStats (config, view, "PSPHOT.INPUT", false)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // find blended neighbors of very saturated stars - psphotDeblendSatstars (config, view); + psphotDeblendSatstars (config, view, "PSPHOT.INPUT"); // mark blended peaks PS_SOURCE_BLEND - if (!psphotBasicDeblend (config, view)) { + if (!psphotBasicDeblend (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // classify sources based on moments, brightness (use supplied psf shape parameters) - if (!psphotRoughClass (config, view)) { + if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } // construct an initial model for each object - psphotGuessModels (config, view); + psphotGuessModels (config, view, "PSPHOT.INPUT"); // merge the newly selected sources into the existing list - psphotMergeSources (config, view); + psphotMergeSources (config, view, "PSPHOT.INPUT"); // linear PSF fit to source peaks - psphotFitSourcesLinear (config, view, false); + psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false); // XXX eventually, add the extended source fits here # if (0) // measure source size for the remaining sources - psphotSourceSize (config, view); + psphotSourceSize (config, view, "PSPHOT.INPUT"); - psphotExtendedSourceAnalysis (config, view); + psphotExtendedSourceAnalysis (config, view, "PSPHOT.INPUT"); - psphotExtendedSourceFits (config, view); + psphotExtendedSourceFits (config, view, "PSPHOT.INPUT"); # endif // calculate source magnitudes - psphotMagnitudes(config, view); + psphotMagnitudes(config, view, "PSPHOT.INPUT"); // XXX ensure this is measured if the analysis succeeds (even if quality is low) - if (!psphotEfficiency(config, view)) { + if (!psphotEfficiency(config, view, "PSPHOT.INPUT")) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -88,7 +88,7 @@ // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view); + psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 28010) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 28013) @@ -23,5 +23,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view) +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -36,5 +36,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotReplaceAllSourcesReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for PSPHOT.INPUT entry %d", i); return false; Index: trunk/psphot/src/psphotRoughClass.c =================================================================== --- trunk/psphot/src/psphotRoughClass.c (revision 28010) +++ trunk/psphot/src/psphotRoughClass.c (revision 28013) @@ -8,5 +8,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotRoughClass (pmConfig *config, const pmFPAview *view) +bool psphotRoughClass (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -26,6 +26,6 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i); + if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i); return false; } @@ -34,5 +34,5 @@ } -bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status; @@ -41,5 +41,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotSetMaskBits.c =================================================================== --- trunk/psphot/src/psphotSetMaskBits.c (revision 28010) +++ trunk/psphot/src/psphotSetMaskBits.c (revision 28013) @@ -37,2 +37,4 @@ return true; } + +// XXX should these be in config->analysis or somewhere else besides 'recipe'? Index: trunk/psphot/src/psphotSkyReplace.c =================================================================== --- trunk/psphot/src/psphotSkyReplace.c (revision 28010) +++ trunk/psphot/src/psphotSkyReplace.c (revision 28013) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotSkyReplace (pmConfig *config, const pmFPAview *view) +bool psphotSkyReplace (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -15,6 +15,6 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i); + if (!psphotSkyReplaceReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for %s entry %d", filerule, i); return false; } @@ -25,10 +25,10 @@ // XXX make this an option? // in order to successfully replace the sky, we must define a corresponding file... -bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { +bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) { psTimerStart ("psphot.skyreplace"); // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotSourceFreePixels.c =================================================================== --- trunk/psphot/src/psphotSourceFreePixels.c (revision 28010) +++ trunk/psphot/src/psphotSourceFreePixels.c (revision 28013) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view) +bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -10,6 +10,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotSourceFreePixelsReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to free source pixels for PSPHOT.INPUT entry %d", i); + if (!psphotSourceFreePixelsReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to free source pixels for %s entry %d", filerule, i); return false; } @@ -18,10 +18,10 @@ } -bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index) { +bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index) { bool status; // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotSourceMatch.c =================================================================== --- trunk/psphot/src/psphotSourceMatch.c (revision 28010) +++ trunk/psphot/src/psphotSourceMatch.c (revision 28013) @@ -1,7 +1,8 @@ # include "psphotInternal.h" -bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects); - -psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view) +bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects); +bool psphotMatchSourcesSetIDs (psArray *objects); + +psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -14,5 +15,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) { + if (!psphotMatchSourcesReadout (objects, config, view, filerule, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i); psFree (objects); @@ -21,10 +22,14 @@ } - psphotMatchSourcesGenerate (config, view, objects); + // create sources for images where an object has been detected in the other images + psphotMatchSourcesAddMissing (config, view, filerule, objects); + + // choose a consistent position; set common sequence values + psphotMatchSourcesSetIDs (objects); return objects; } -bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { +bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) { bool status = false; @@ -38,5 +43,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); @@ -145,5 +150,5 @@ } -bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) { +bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) { bool status = false; @@ -167,5 +172,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest psAssert (file, "missing file?"); @@ -255,2 +260,16 @@ return true; } + +bool psphotMatchSourcesSetIDs (psArray *objects) { + + for (int i = 0; i < objects->n; i++) { + pmPhotObj *obj = objects->data[i]; + + // set the source->seq values + for (int j = 0; j < obj->sources->n; j++) { + pmSource *src = obj->sources->data[j]; + src->seq = i; + } + } + return true; +} Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 28010) +++ trunk/psphot/src/psphotSourceSize.c (revision 28013) @@ -33,5 +33,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize) +bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize) { bool status = true; @@ -51,5 +51,5 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) { + if (!psphotSourceSizeReadout (config, view, filerule, i, recipe, getPSFsize)) { psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i); return false; @@ -60,5 +60,5 @@ // this function use an internal flag to mark sources which have already been measured -bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize) +bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool getPSFsize) { bool status; @@ -68,5 +68,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 28010) +++ trunk/psphot/src/psphotSourceStats.c (revision 28013) @@ -5,5 +5,5 @@ // The new sources are added to any existing sources on detections->newSources. The sources // on detections->allSources are ignored. -bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow) +bool psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow) { bool status = true; @@ -16,8 +16,13 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + // int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + // if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotSourceStatsReadout (config, view, "PSPHOT.INPUT", i, recipe, setWindow)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); + // if (i == chisqNum) continue; // skip chisq image + if (!psphotSourceStatsReadout (config, view, filerule, i, recipe, setWindow)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i); return false; } @@ -26,5 +31,5 @@ } -bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow) { +bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow) { bool status = false; @@ -34,5 +39,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (file, "missing file?"); Index: trunk/psphot/src/psphotStackArguments.c =================================================================== --- trunk/psphot/src/psphotStackArguments.c (revision 28010) +++ trunk/psphot/src/psphotStackArguments.c (revision 28013) @@ -1,4 +1,5 @@ # include "psphotStandAlone.h" +static void dumpTemplate(void); static void usage(pmConfig *config, int exitCode); static void writeHelpInfo(FILE* ofile); @@ -11,4 +12,6 @@ if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout); if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout); + + if (psArgumentGet(argc, argv, "-template")) dumpTemplate(); // load config data from default locations @@ -84,7 +87,8 @@ "\n" "where INPUTS.mdc contains various METADATAs, each with:\n" - "\tIMAGE(STR): Image filename\n" - "\tMASK(STR): Mask filename\n" - "\tVARIANCE(STR): Variance map filename\n" + "\tIMAGE : Image filename\n" + "\tMASK : Mask filename\n" + "\tVARIANCE : Variance map filename\n" + "(use -template to generate a sample input.mdc file)\n" "OUTROOT is the 'root name' for output files\n" "\n" @@ -144,2 +148,30 @@ } +static void dumpTemplate(void) { + + fprintf (stdout, "# this line is required for multiple INPUT blocks to be accepted\n"); + fprintf (stdout, "INPUT MULTI\n\n"); + + fprintf (stdout, "# copy and repeat the following block as needed (one per input image set)\n"); + fprintf (stdout, "# either RAW (unconvolved) or CNV (convolved) input images are required\n"); + fprintf (stdout, "# if both are supplied, by default RAW is used for detection, CNV is convolved (further) to match target PSF\n"); + fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n"); + fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n"); + fprintf (stdout, "# PSF may be supplied for the convolution target\n"); + fprintf (stdout, "INPUT METADATA\n"); + fprintf (stdout, " RAW:IMAGE STR file.im.fits # signal image filename\n"); + fprintf (stdout, " RAW:MASK STR file.mk.fits # mask image filename\n"); + fprintf (stdout, " RAW:VARIANCE STR file.wt.fits # variance image filename\n"); + fprintf (stdout, " RAW:PSF STR file.psf.fits # psf from input unconvolved image\n"); + + fprintf (stdout, " CNV:IMAGE STR file.im.fits # signal image filename\n"); + fprintf (stdout, " CNV:MASK STR file.mk.fits # mask image filename\n"); + fprintf (stdout, " CNV:VARIANCE STR file.wt.fits # variance image filename\n"); + fprintf (stdout, " CNV:PSF STR file.psf.fits # psf from input convolved image\n"); + + fprintf (stdout, " SOURCES STR file.cmf # measured source positions\n"); + fprintf (stdout, "END\n"); + + psLibFinalize(); + exit(PS_EXIT_SUCCESS); +} Index: trunk/psphot/src/psphotStackChisqImage.c =================================================================== --- trunk/psphot/src/psphotStackChisqImage.c (revision 28010) +++ trunk/psphot/src/psphotStackChisqImage.c (revision 28013) @@ -6,5 +6,5 @@ // XXX supply filename or keep PSPHOT.INPUT fixed? -bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view) +bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv) { bool status = false; @@ -21,7 +21,8 @@ // loop over the available readouts + // generate the chisq image from the 'detection' images for (int i = 0; i < num; i++) { - if (!psphotStackChisqImageAddReadout(config, view, &chiReadout, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); + if (!psphotStackChisqImageAddReadout(config, view, ruleDet, &chiReadout, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to model background for %s entry %d", ruleDet, i); return false; } @@ -35,6 +36,12 @@ psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num); - // need to save the resulting image somewhere (PSPHOT.INPUT?) - if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) { + // save the resulting image in the 'detection' set + if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleDet, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) { + psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files"); + return false; + } + + // save the resulting image in the 'convolved' set + if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleCnv, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) { psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files"); return false; @@ -48,6 +55,6 @@ bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration const pmFPAview *view, + const char *filerule, pmReadout **chiReadout, - char *filename, int index) { @@ -55,5 +62,5 @@ // find the currently selected readout - pmFPAfile *input = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *input = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest psAssert (input, "missing file?"); @@ -111,5 +118,5 @@ } -bool psphotStackRemoveChisqFromInputs (pmConfig *config) { +bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule) { bool status = false; @@ -121,5 +128,5 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - pmFPAfileRemoveSingle (config->files, "PSPHOT.INPUT", chisqNum); + pmFPAfileRemoveSingle (config->files, filerule, chisqNum); inputNum --; Index: trunk/psphot/src/psphotStackImageLoop.c =================================================================== --- trunk/psphot/src/psphotStackImageLoop.c (revision 28010) +++ trunk/psphot/src/psphotStackImageLoop.c (revision 28013) @@ -1,9 +1,14 @@ # include "psphotStandAlone.h" - -# define ESCAPE(MESSAGE) { \ - psError(PSPHOT_ERR_DATA, false, MESSAGE); \ - psFree (view); \ - return false; \ -} +#define WCS_NONLIN_TOL 0.001 // Non-linear tolerance for header WCS + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + +bool GetAstrometryFPA (pmConfig *config, pmFPAview *view); +bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry); +bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry); bool psphotStackImageLoop (pmConfig *config) { @@ -14,6 +19,9 @@ pmReadout *readout; - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); - if (!status) { + pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW"); + pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV"); + pmFPAfile *input = inputRaw ? inputRaw : inputCnv; + + if (!input) { psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); return false; @@ -24,4 +32,6 @@ // XXX for now, just load the full set of images up front if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + bool bilevelAstrometry = GetAstrometryFPA(config, view); // for psphot, we force data to be read at the chip level @@ -41,6 +51,5 @@ if (! readout->data_exists) { continue; } -# if (0) - // uncomment to generate matched psfs + // PSF matching if (!psphotStackMatchPSFs (config, view)) { psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); @@ -48,5 +57,4 @@ return false; } -# endif // XXX for now, we assume there is only a single chip in the PHU: @@ -69,7 +77,13 @@ } } + + GetAstrometryChip(config, view, bilevelAstrometry); + // save output which is saved at the chip level if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); } + + SetAstrometryFPA(config, view, bilevelAstrometry); + // save output which is saved at the fpa level if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); @@ -87,2 +101,151 @@ */ +# define UPDATE_HEADER 0 + +bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) { + + bool status = false; + + bool foundAstrometry = false; + bool bilevelAstrometry = false; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // find the currently selected readout + pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest + psAssert (output, "missing file?"); + + pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest + pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest + pmFPAfile *input = inputRaw ? inputRaw : inputCnv; + psAssert (input, "missing input file"); + + // find the FPA phu + pmHDU *phu = pmFPAviewThisPHU(view, input->fpa); + if (!phu) { + psWarning("Unable to read bilevel mosaic astrometry for input FPA entry %d", i); + continue; + } + + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (!ctype) { + psWarning("Error in WCS keywords for input FPA entry %d", i); + continue; + } + + if (!foundAstrometry) { + bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } else { + if (bilevelAstrometry != !strcmp (&ctype[4], "-DIS")) { + psAbort("astrometry type mis-match"); + } + } + + if (bilevelAstrometry) { + // update the output structures + if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) { + psWarning("Unable to read bilevel mosaic astrometry for input FPA."); + } + } + } + return bilevelAstrometry; +} + +bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) { + + bool status = false; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // find the currently selected readout + pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest + psAssert (output, "missing file?"); + + pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest + pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest + pmFPAfile *input = inputRaw ? inputRaw : inputCnv; + psAssert (input, "missing input file"); + + // Need to update the output for astrometry. Read WCS data from the corresponding + // header and save in the output fpa + pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa); + pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output + +# if (UPDATE_HEADER) + pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa); + if (!outHDU) { + pmFPAAddSourceFromView(output->fpa, "name", view, output->format); + outHDU = pmFPAviewThisHDU (view, output->fpa); + psAssert (outHDU, "failed to make HDU"); + } +# endif + + if (bilevelAstrometry) { + if (!pmAstromReadBilevelChip (outChip, inHDU->header)) { + psWarning("Unable to read bilevel chip astrometry for input FPA."); + continue; + } +# if (UPDATE_HEADER) + if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) { + psWarning("Unable to generate WCS header."); + continue; + } +# endif + } else { + // we use a default FPA pixel scale of 1.0 + if (!pmAstromReadWCS (output->fpa, outChip, inHDU->header, 1.0)) { + psWarning("Unable to read WCS astrometry for input FPA."); + continue; + } +# if (UPDATE_HEADER) + if (UPDATE_HEADER && !pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) { + psWarning("Unable to generate WCS header."); + continue; + } +# endif + } + } + + return true; +} + +bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) { + + bool status = false; + + if (!bilevelAstrometry) return true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // find the currently selected readout + pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest + psAssert (output, "missing file?"); + +# if (UPDATE_HEADER) + pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa); + if (!PHU) { + pmFPAAddSourceFromView(output->fpa, "name", view, output->format); + PHU = pmFPAviewThisPHU (view, output->fpa); + psAssert (PHU, "failed to make PHU"); + } + + if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) { + psWarning("Unable to generate WCS header."); + } +# endif + } + + return true; +} Index: trunk/psphot/src/psphotStackMatchPSFs.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFs.c (revision 28010) +++ trunk/psphot/src/psphotStackMatchPSFs.c (revision 28013) @@ -1,71 +1,90 @@ # include "psphotInternal.h" -bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass) +bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view) { bool status = true; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT"); + psAssert(recipe, "We've thrown an error on this before."); int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - // loop over the available readouts + // 'options' carries info needed to perform the stack matching + psphotStackOptions *options = psphotStackOptionsAlloc(num); + + options->convolve = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.MATCH.PSF"); + psAssert (status, "PSPHOT.STACK.MATCH.PSF not in recipe"); + + if (options->convolve) { + char *convolveSource = psMetadataLookupStr (&status, recipe, "PSPHOT.STACK.MATCH.PSF.SOURCE"); + options->convolveSource = psphotStackConvolveSourceFromString (convolveSource); + if (options->convolveSource == PSPHOT_CNV_SRC_NONE) { + psError (PSPHOT_ERR_CONFIG, true, "stack convolution source not defined in recipe"); + return false; + } + } + + // loop over the available readouts (ignore chisq image) for (int i = 0; i < num; i++) { - if (!psphotStackMatchPSFsReadout (config, view, i)) { + if (!psphotStackMatchPSFsPrepare (config, view, options, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i); + return false; + } + } + + // Generate target PSF + if (options->convolve) { + options->psf = psphotStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask); + if (!options->psf) { + psError(psErrorCodeLast(), false, "Unable to determine output PSF."); + return false; + } + psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf); + options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target + psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing); + + // XXX is this needed to supply the psf to psphot?? + // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip + // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf); + // outChip->data_exists = true; + } + + // loop over the available readouts (ignore chisq image) + for (int i = 0; i < num; i++) { + if (!psphotStackMatchPSFsReadout (config, view, options, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); return false; } } + + psFree (options); return true; } // convolve the image to match desired PSF -bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) { +bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) { - bool status; - int pass; - float NSIGMA_PEAK = 25.0; - int NMAX = 0; + pmFPAfile *fileSrc = psphotStackGetConvolveSource(config, options, index); + if (!fileSrc) { + psError(PSPHOT_ERR_CONFIG, false, "desired convolution source is missing"); + return false; + } - // find the currently selected readout - pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest - psAssert (fileRaw, "missing file?"); + pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", index); // File of interest + psAssert (fileOut, "missing output file?"); - pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest - psAssert (fileCnv, "missing file?"); + pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa); + psAssert (readoutSrc, "missing readout?"); - pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa); - psAssert (readoutRaw, "missing readout?"); + pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa); + if (readoutOut == NULL) { + readoutOut = pmFPAGenerateReadout(config, view, "PSPHOT.STACK.OUTPUT.IMAGE", fileSrc->fpa, NULL, index); + psAssert (readoutOut, "missing readout?"); + } - pmReadout *readoutCnv = pmFPAviewThisReadout(view, fileCnv->fpa); - psAssert (readoutCnv, "missing readout?"); - - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - psAssert (maskVal, "missing mask value?"); - - /***** set up recipe options *****/ - - psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe - psAssert(stackRecipe, "We've thrown an error on this before."); - - // Look up appropriate values from the ppSub recipe - psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe - psAssert(subRecipe, "recipe missing"); - - int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size - - float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction - - psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in - psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch - psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor - psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels - psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - - bool mdok; // Status of MD lookup - float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness - int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads - - if (!pmReadoutMaskNonfinite(readout, maskVal)) { + // set NAN pixels to 'SAT' + psImageMaskType maskVal = pmConfigMaskGet("SAT", config); + if (!pmReadoutMaskNonfinite(readoutSrc, maskVal)) { psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout."); return false; @@ -74,79 +93,29 @@ // Image Matching (PSFs or just flux) if (options->convolve) { - // Full match of PSFs - pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily - - // Read previously produced kernel - if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { - loadKernel(); - } else { - matchKernel(); - } // !DEBUG.STACK - - saveMatchData(); - - saveChiSquare(); - - renormKernel(); - - // Reject image completely if the maximum deconvolution fraction exceeds the limit - float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction - if (deconv > deconvLimit) { - psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); - psFree(conv); - return NULL; - } - - readout->analysis = psMetadataCopy(readout->analysis, conv->analysis); - - psFree(conv); + matchKernel(config, readoutOut, readoutSrc, options, index); + saveMatchData(readoutOut, options, index); + // renormKernel(readoutCnv, options, index); } else { - // only match the flux - float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation - psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); + // only match the flux (NO! not for multi-filter, at least!) + // XXX do not generate readoutCnv in this case? + // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation + // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); } - // Ensure the background value is zero - psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psWarning("Can't measure background for image."); - psErrorClear(); - } else { - if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { - psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", - psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); - (void)psBinaryOp(readout->image, readout->image, "-", - psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); - } - } + rescaleData(readoutOut, config, options, index); - if (!stackRenormaliseReadout(config, readout)) { - psFree(rng); - psFree(bg); - return false; - } - - // Measure the variance level for the weighting - if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) { - if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { - psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image."); - psFree(rng); - psFree(bg); - return false; - } - options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); - } else { - options->weightings->data.F32[index] = 1.0; - } - psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", - index, options->weightings->data.F32[index]); - - psFree(rng); - psFree(bg); - - dumpImage3(); + dumpImage(readoutOut, readoutSrc, index, "convolved"); return true; } + + +# if (0) +// Read previously produced kernel +if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { + loadKernel(config, readoutCnv, options, index); +} else { + matchKernel(config, readoutCnv, readoutRaw, options, index); +} +# endif Index: trunk/psphot/src/psphotStackMatchPSFsPrepare.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFsPrepare.c (revision 28013) +++ trunk/psphot/src/psphotStackMatchPSFsPrepare.c (revision 28013) @@ -0,0 +1,114 @@ +# include "psphotInternal.h" + +pmPSF *selectPSF (pmFPAfile *fileSrc, const pmFPAview *view, psphotStackOptions *options, int index); +bool determineSeeing (pmPSF *psf, psphotStackOptions *options, int index); + +// set up the stacking parameters +bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) { + + if (!options->convolve) { + psLogMsg ("psphotStack", PS_LOG_INFO, "psf-matching NOT selected"); + return true; + } + + pmFPAfile *fileSrc = psphotStackGetConvolveSource(config, options, index); + if (!fileSrc) { + psError(PSPHOT_ERR_CONFIG, false, "desired convolution source is missing"); + return false; + } + + pmPSF *psf = selectPSF (fileSrc, view, options, index); + if (!psf) { + psError(PSPHOT_ERR_PSF, false, "Problem with PSF."); + return false; + } + + determineSeeing (psf, options, index); + + // load the sources (used to find reference sources for the kernel stamps) + { + pmFPAfile *inputSrc = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.SOURCES", index); // File of interest + pmReadout *readout = pmFPAviewThisReadout(view, inputSrc->fpa); // Readout with sources + pmDetections *detections = psMetadataLookupPtr(NULL, readout->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections || !detections->allSources) { + psWarning("No detections found for image %d --- rejecting.", index); + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0x01; + return true; // XXX not an error: we continue processing other images + } + psAssert (detections->allSources, "missing sources?"); + options->sourceLists->data[index] = psMemIncrRefCounter(detections->allSources); + } + + return true; +} + +// select the corresponding input psf +pmPSF *selectPSF (pmFPAfile *fileSrc, const pmFPAview *view, psphotStackOptions *options, int index) { + + bool status; + + pmChip *chip = pmFPAviewThisChip(view, fileSrc->fpa); // The chip holds the PSF + pmPSF *psf = psMetadataLookupPtr(&status, chip->analysis, "PSPHOT.PSF"); // PSF + if (!psf) { + // XXX if we were not supplied a PSF, we should be able to generate one by calling psphot + psError(PSPHOT_ERR_PROG, true, "Unable to find PSF."); + return NULL; + } + options->psfs->data[index] = psMemIncrRefCounter(psf); + + // find the image size + pmCell *cell = pmFPAviewThisCell(view, fileSrc->fpa); // Cell of interest + pmHDU *hdu = pmHDUFromCell(cell); + assert(hdu && hdu->header); + int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows + psAssert ((naxis1 > 0) && (naxis2 > 0), "Unable to determine size of image from PSF."); + if (!options->numCols) { + options->numCols = naxis1; + } else { + if (options->numCols != naxis1) { + psError (PSPHOT_ERR_CONFIG, true, "mismatched sizes (NAXIS1) for input images"); + return NULL; + } + } + if (!options->numRows) { + options->numRows = naxis2; + } else { + if (options->numRows != naxis2) { + psError (PSPHOT_ERR_CONFIG, true, "mismatched sizes (NAXIS2) for input images"); + return NULL; + } + } + return psf; +} + +// determine the input seeing +bool determineSeeing (pmPSF *psf, psphotStackOptions *options, int index) { + + // XXX set this based on the mode (+1 for poly, +0 for map) + float xNum = PS_MAX(psf->trendNx, 1); + float yNum = PS_MAX(psf->trendNy, 1); // Number of realisations + float sumFWHM = 0.0; // FWHM for image + int numFWHM = 0; // Number of FWHM measurements + for (float y = 0; y < yNum; y += 1.0) { + float yPos = options->numRows * ((y + 0.5) / yNum); + for (float x = 0; x < xNum; x++) { + float xPos = options->numCols * ((x + 0.5) / xNum); + float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image + if (isfinite(fwhm)) { + sumFWHM += fwhm; + numFWHM++; + } + } + } + if (numFWHM == 0) { + options->inputSeeing->data.F32[index] = NAN; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0x02; + psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index); + } else { + options->inputSeeing->data.F32[index] = sumFWHM / (float)numFWHM; + } + psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, options->inputSeeing->data.F32[index]); + return true; +} + Index: trunk/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 28010) +++ trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 28013) @@ -1,10 +1,4 @@ -/***** defines *****/ - -#define ARRAY_BUFFER 16 // Number to add to array at a time -#define MAG_IGNORE 50 // Ignore magnitudes fainter than this --- they're not real! -#define FAKE_SIZE 1 // Size of fake convolution kernel -#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources -#define NOISE_FRACTION 0.01 // Set minimum flux to this fraction of noise -#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix +# include "psphotInternal.h" +# define ARRAY_BUFFER 16 // Number to add to array at a time // XXX better name @@ -18,10 +12,10 @@ psFree(resolved); if (!fits) { - psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name); + psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name); return false; } psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest if (!image) { - psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name); + psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name); psFitsClose(fits); return false; @@ -51,6 +45,7 @@ } +# define SN_MIN 50.0 psArray *stackSourcesFilter(psArray *sources, // Source list to filter - int exclusion // Exclusion zone, pixels + int exclusion // Exclusion zone, pixels ) { @@ -68,4 +63,6 @@ continue; } + if (!source->peak) continue; + if (source->peak->SN < SN_MIN) continue; coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); numGood++; @@ -83,4 +80,6 @@ continue; } + if (!source->peak) continue; + if (source->peak->SN < SN_MIN) continue; float xSource, ySource; // Coordinates of source coordsFromSource(&xSource, &ySource, source); @@ -90,5 +89,5 @@ long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone - psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", + psTrace("psphotStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", coords->data.F64[0], coords->data.F64[1], numWithin); if (numWithin == 1) { @@ -104,5 +103,5 @@ psFree(y); - psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); + psLogMsg("psphotStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); return filtered; @@ -111,5 +110,5 @@ // Add background into the fake image // Based on ppSubBackground() -static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model +psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model const pmConfig *config // Configuration ) @@ -121,7 +120,7 @@ int numCols = image->numCols, numRows = image->numRows; // Size of image - psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); + psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK"); psAssert(ppStackRecipe, "Need PPSTACK recipe"); - psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, "PSPHOT"); psAssert(psphotRecipe, "Need PSPHOT recipe"); @@ -138,5 +137,5 @@ psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model if (!psImageUnbin(unbinned, binned, binning)) { - psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); + psError(PSPHOT_ERR_DATA, false, "Unable to unbin background model"); psFree(binned); psFree(unbinned); @@ -153,8 +152,7 @@ ) { -#if 1 bool mdok; // Status of metadata lookups - psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK"); // Recipe for ppStack psAssert(recipe, "Need PPSTACK recipe"); @@ -163,15 +161,15 @@ int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe"); + psError(PSPHOT_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe"); return false; } float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN"); if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe"); + psError(PSPHOT_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe"); return false; } float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); + psError(PSPHOT_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); return false; } @@ -181,7 +179,4 @@ psImageCovarianceTransfer(readout->variance, readout->covariance); return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); -#else - return true; -#endif } @@ -190,367 +185,424 @@ // It implicitly assumes the output root name is the same between invocations. -bool loadKernel () { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); - psAssert(file, "Require file"); - - pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest - view->chip = view->cell = view->readout = 0; - psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest - - // Read convolution kernel - psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename - psFree(filename); - psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel - psFree(resolved); - if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) { - psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel"); - psFitsClose(fits); - return false; - } - psFitsClose(fits); - - if (!readImage(&readout->image, options->convImages->data[index], config) || - !readImage(&readout->mask, options->convMasks->data[index], config) || - !readImage(&readout->variance, options->convVariances->data[index], config)) { - psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image."); - return false; - } - - psRegion *region = psMetadataLookupPtr(NULL, conv->analysis, - PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); - - pmSubtractionAnalysis(conv->analysis, NULL, kernels, region, - readout->image->numCols, readout->image->numRows); - - psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel - bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting - psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix - psImageCovarianceSetThreads(oldThreads); - psFree(readout->covariance); - readout->covariance = covar; - psFree(kernel); -} - -bool dumpImage() { - // XXX should be optional - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "fake_%03d.fits", index); - pmStackVisualPlotTestImage(fake->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); - psFitsClose(fits); - } - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "real_%03d.fits", index); - pmStackVisualPlotTestImage(readout->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); - psFitsClose(fits); - } -} - -bool dumpImage2() { - // XXX should be optional - - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "conv_%03d.fits", index); - pmStackVisualPlotTestImage(conv->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL); - psFitsClose(fits); - } - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "diff_%03d.fits", index); - pmStackVisualPlotTestImage(fake->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psBinaryOp(fake->image, conv->image, "-", fake->image); - psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); - psFitsClose(fits); - } -} - -bool dumpImage3() +# if (0) +bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) { + + // Read the convolution kernel from the saved file + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); + psAssert(file, "Require file"); + + pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest + view->chip = view->cell = view->readout = 0; + psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest + + // Read convolution kernel data + psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename + psFree(filename); + psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel + psFree(resolved); + if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) { + psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel"); + psFitsClose(fits); + return false; + } + psFitsClose(fits); + + // read the convolved pixels (image, mask, variance) -- names are pre-defined + if (!readImage(&readoutCnv->image, options->convImages->data[index], config) || + !readImage(&readoutCnv->mask, options->convMasks->data[index], config) || + !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) { + psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image."); + return false; + } + + // XXX ??? not sure what is happening here -- consult Paul + psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); + + pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows); + + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + + // update the covariance matrix + // XXX why is this needed if we have correctly read the saved data? + bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting + psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix + psImageCovarianceSetThreads(oldThreads); + psFree(readoutCnv->covariance); + readoutCnv->covariance = covar; + psFree(kernel); + return true; +} +# endif + +bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) { + + pmHDU *hdu = pmHDUFromCell(readoutRef->parent); + psString name = NULL; + psStringAppend(&name, "%s_%03d.fits", rootname, index); + pmStackVisualPlotTestImage(readoutOut->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, readoutOut->image, 0, NULL); + psFitsClose(fits); + return true; +} + +bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname) { + + pmHDU *hdu = pmHDUFromCell(readoutRef->parent); + psString name = NULL; + psStringAppend(&name, "%s_%03d.fits", rootname, index); + pmStackVisualPlotTestImage(readoutFake->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psBinaryOp(readoutFake->image, readoutConv->image, "-", readoutFake->image); + psFitsWriteImage(fits, hdu->header, readoutFake->image, 0, NULL); + psFitsClose(fits); + return true; +} + +// perform the bulk of the PSF-matching +bool matchKernel(pmConfig *config, pmReadout *readoutOut, pmReadout *readoutSrc, psphotStackOptions *options, int index) { + + bool mdok; + + psAssert(options->psf, "Require target PSF"); + psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); + + psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe + psAssert(stackRecipe, "We've thrown an error on this before."); + + // Look up appropriate values from the ppSub recipe + psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(subRecipe, "recipe missing"); + + psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in + psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor + psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad + + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness + int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads + + int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order + float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing + + int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size + int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size + + float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps + int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches + int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations + float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold + float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation + + const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type + pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type + psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths + psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders + int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius + int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order + int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel + float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction + float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search + int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search + float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor" + + bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE"); // Scale kernel parameters? + float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PSPHOT_ERR_CONFIG, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + // These values are specified specifically for stacking + const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS");// Stamps filename + + psVector *widthsCopy = NULL; + psVector *optWidths = NULL; + pmReadout *fake = NULL; + psArray *stampSources = NULL; + + bool optimum = false; + optWidths = SetOptWidths(&optimum, subRecipe); // Vector with FWHMs for optimum search + + // For the sake of stamps, remove nearby sources + stampSources = stackSourcesFilter(options->sourceLists->data[index], footprint); // Filtered list of sources + + fake = makeFakeReadout(config, readoutSrc, stampSources, options->psf, maskVal | maskBad, footprint + size); + if (!fake) goto escape; + + dumpImage(fake, readoutSrc, index, "fake"); + dumpImage(readoutSrc, readoutSrc, index, "real"); + + if (threads) pmSubtractionThreadsInit(); + + // Do the image matching + pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readoutSrc->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel + if (kernel) { + if (!pmSubtractionMatchPrecalc(NULL, readoutOut, fake, readoutSrc, readoutSrc->analysis, stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + goto escape; + } + } else { + // Scale the input parameters + widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths + if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, options->inputSeeing->data.F32[index], options->targetSeeing, scaleRef, scaleMin, scaleMax)) { + psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); + goto escape; + } + + if (!pmSubtractionMatch(NULL, readoutOut, fake, readoutSrc, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { + psError(psErrorCodeLast(), false, "Unable to match images."); + goto escape; + } + } + + // Reject image completely if the maximum deconvolution fraction exceeds the limit + float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction + float deconv = psMetadataLookupF32(NULL, readoutOut->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction + if (deconv > deconvLimit) { + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); + goto escape; + } + + dumpImage(readoutOut, readoutSrc, index, "conv"); + dumpImageDiff(readoutOut, fake, readoutSrc, index, "diff"); + + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(widthsCopy); + pmSubtractionThreadsFinalize(); + return true; + +escape: + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(widthsCopy); + pmSubtractionThreadsFinalize(); + return false; +} + +// Extract the regions and solutions used in the image matching +// This stops them from being freed when we iterate back up the FPA +// Record the chi-square value +// XXX this function may not be needed for psphotStack +bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index) { + + psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "convolved_%03d.fits", index); - pmStackVisualPlotTestImage(readout->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); - psFitsClose(fits); - } - -bool matchKernel() { - // Normal operations here - psAssert(options->psf, "Require target PSF"); - psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); - - int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order - float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs - float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing - int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size - float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps - int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches - int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations - float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold - float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel - float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw - float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images - float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky - float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation - - const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type - pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type - psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths - psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders - int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius - int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order - int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel - float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction - bool optimum = psMetadataLookupBool(&mdok, subRecipe, "OPTIMUM"); // Derive optimum parameters? - float optMin = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MIN"); // Minimum width for search - float optMax = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MAX"); // Maximum width for search - float optStep = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.STEP"); // Step for search - float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search - int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search - float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor" - - bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE"); // Scale kernel parameters? - float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling - float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling - float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling - if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { - psError(PPSTACK_ERR_CONFIG, false, - "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", - scaleRef, scaleMin, scaleMax); - return false; - } - - - // These values are specified specifically for stacking - const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename - - psVector *optWidths = NULL; // Vector with FWHMs for optimum search - if (optimum) { - optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); - } - - pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF - - psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); - psFree(fake); - psFree(optWidths); - psFree(conv); - psFree(bg); - psFree(rng); - return false; - } - float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_REGION); + regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); + } + psFree(iter); + } + + psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction + kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); + } + psFree(iter); + } + psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match"); + + // Record chi^2 + { + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernels = item->data.V; // Convolution kernels + sum += kernels->mean; + num++; + } + psFree(iter); + options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); + } + + return true; +} + +// Kernel normalisation for convolved readout +bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) { + + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_TYPE_F32); + float norm = item->data.F32; // Normalisation + sum += norm; + num++; + } + psFree(iter); + float conv = sum/num; // Mean normalisation from convolution + float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars + float renorm = stars / conv; // Renormalisation to apply + psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars); + + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); + return true; +} + +// adjust scaling for readout (remove background, ..., determine weighting) +bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) { + + psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe + psAssert(stackRecipe, "We've thrown an error on this before."); + + // Look up appropriate values from the ppSub recipe + psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(subRecipe, "recipe missing"); + + psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in + psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad + + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + // Ensure the background value is zero + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + + // XXX why is this in config->arguments and not recipe? + if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { + if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { + psAbort("Can't measure background for image."); + // XXX we used to clear error: why is this acceptable? psErrorClear(); + } + + float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN); + float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV); + + psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev); + psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32)); + } + + if (!stackRenormaliseReadout(config, readout)) { + psFree(rng); + psFree(bg); + return false; + } + + // Measure the variance level for the weighting + if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) { + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image."); psFree(rng); psFree(bg); - - // For the sake of stamps, remove nearby sources - psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index], - footprint); // Filtered list of sources - - bool oldThreads = pmReadoutFakeThreads(true); // Old threading state - if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, - stampSources, SOURCE_MASK, NULL, NULL, options->psf, - minFlux, footprint + size, false, true)) { - psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); - psFree(fake); - psFree(optWidths); - psFree(conv); - return false; - } - pmReadoutFakeThreads(oldThreads); - - fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK); - - // Add the background into the target image - psImage *bgImage = stackBackgroundModel(readout, config); // Image of background - psBinaryOp(fake->image, fake->image, "+", bgImage); - psFree(bgImage); - - dumpImage(); - - if (threads > 0) { - pmSubtractionThreadsInit(); - } - - // Do the image matching - pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel - if (kernel) { - if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis, - stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, - poorFrac, badFrac)) { - psError(psErrorCodeLast(), false, "Unable to convolve images."); - psFree(fake); - psFree(optWidths); - psFree(stampSources); - psFree(conv); - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - return false; - } - } else { - // Scale the input parameters - psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths - if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, - options->inputSeeing->data.F32[index], - options->targetSeeing, scaleRef, scaleMin, scaleMax)) { - psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); - psFree(fake); - psFree(optWidths); - psFree(stampSources); - psFree(conv); - psFree(widthsCopy); - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - return false; - } - - if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, - threshold, stampSources, stampsName, type, size, order, widthsCopy, - orders, inner, ringsOrder, binning, penalty, - optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, - sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, - poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { - psError(psErrorCodeLast(), false, "Unable to match images."); - psFree(fake); - psFree(optWidths); - psFree(stampSources); - psFree(conv); - psFree(widthsCopy); - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - return false; - } - psFree(widthsCopy); - } - - dumpImage2(); - - psFree(fake); - psFree(optWidths); - psFree(stampSources); - - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - - // Replace original images with convolved - psFree(readout->image); - psFree(readout->mask); - psFree(readout->variance); - psFree(readout->covariance); - readout->image = psMemIncrRefCounter(conv->image); - readout->mask = psMemIncrRefCounter(conv->mask); - readout->variance = psMemIncrRefCounter(conv->variance); - readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC); - -} - -bool saveMatchData () { - // Extract the regions and solutions used in the image matching - // This stops them from being freed when we iterate back up the FPA - psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions - { - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_REGION); - regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); - } - psFree(iter); + return false; } - psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels - { - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_UNKNOWN); - pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction - kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); - } - psFree(iter); - } - psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match"); -} - -bool saveChiSquare() { - // Record chi^2 - { - double sum = 0.0; // Sum of chi^2 - int num = 0; // Number of measurements of chi^2 - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_UNKNOWN); - pmSubtractionKernels *kernels = item->data.V; // Convolution kernels - sum += kernels->mean; - num++; - } - psFree(iter); - options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); - } - -} - -bool renormKernel() { - // Kernel normalisation - { - double sum = 0.0; // Sum of chi^2 - int num = 0; // Number of measurements of chi^2 - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_TYPE_F32); - float norm = item->data.F32; // Normalisation - sum += norm; - num++; - } - psFree(iter); - float conv = sum/num; // Mean normalisation from convolution - float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars - float renorm = stars / conv; // Renormalisation to apply - psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", - index, renorm, conv, stars); - psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); - psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); - } - -} + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); + } else { + options->weightings->data.F32[index] = 1.0; + } + psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]); + + psFree(rng); + psFree(bg); + return true; +} + +# define NOISE_FRACTION 0.01 // Set minimum flux to this fraction of noise +# define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources + +// generate a fake readout against which to PSF match +pmReadout *makeFakeReadout(pmConfig *config, pmReadout *readoutSrc, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize) { + + pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF + + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(bg, NULL, readoutSrc->image, readoutSrc->mask, maskVal, rng)) { + psError(PSPHOT_ERR_DATA, false, "Can't measure background for image."); + psFree(fake); + psFree(bg); + psFree(rng); + return NULL; + } + float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image + psFree(rng); + psFree(bg); + + bool oldThreads = pmReadoutFakeThreads(true); // Old threading state + if (!pmReadoutFakeFromSources(fake, readoutSrc->image->numCols, readoutSrc->image->numRows, sources, SOURCE_MASK, NULL, NULL, psf, minFlux, fullSize, false, true)) { + psError(PSPHOT_ERR_DATA, false, "Unable to generate fake image with target PSF."); + psFree(fake); + return NULL; + } + pmReadoutFakeThreads(oldThreads); + + fake->mask = psImageCopy(NULL, readoutSrc->mask, PS_TYPE_IMAGE_MASK); + + // Add the background into the target image + psImage *bgImage = stackBackgroundModel(readoutSrc, config); // Image of background + psBinaryOp(fake->image, fake->image, "+", bgImage); + psFree(bgImage); + + return fake; +} + +// set the widths +psVector *SetOptWidths (bool *optimum, psMetadata *recipe) { + + bool status; + + *optimum = psMetadataLookupBool(&status, recipe, "OPTIMUM"); // Derive optimum parameters? + psAssert (status, "missing recipe value %s", "OPTIMUM"); + + psVector *optWidths = NULL; // Vector with FWHMs for optimum search + + if (*optimum) { + float optMin = psMetadataLookupF32(&status, recipe, "OPTIMUM.MIN"); // Minimum width for search + psAssert (status, "missing recipe value %s", "OPTIMUM.MIN"); + + float optMax = psMetadataLookupF32(&status, recipe, "OPTIMUM.MAX"); // Maximum width for search + psAssert (status, "missing recipe value %s", "OPTIMUM.MAX"); + + float optStep = psMetadataLookupF32(&status, recipe, "OPTIMUM.STEP"); // Step for search + psAssert (status, "missing recipe value %s", "OPTIMUM.STEP"); + + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + + return optWidths; +} Index: trunk/psphot/src/psphotStackObjects.c =================================================================== --- trunk/psphot/src/psphotStackObjects.c (revision 28013) +++ trunk/psphot/src/psphotStackObjects.c (revision 28013) @@ -0,0 +1,48 @@ +# include "psphotInternal.h" + +// set a consistent position +bool psphotStackObjectsUnifyPosition (psArray *objects) { + + // consistent position: AVERAGE or CHISQ? + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + int npts = 0; + float x = 0.0; + float y = 0.0; + // measure the average position (weighted?) + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + if (!source) continue; + if (!source->peak) continue; + + x += source->peak->xf; + y += source->peak->yf; + npts ++; + } + if (npts == 0) continue; + + x /= (float) npts; + y /= (float) npts; + + // set the positions + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + if (!source) continue; + if (!source->peak) continue; + + source->peak->xf = x; + source->peak->yf = y; + npts ++; + } + object->x = x; + object->y = y; + } + + psLogMsg ("psphot", PS_LOG_INFO, "updated positions\n"); + return true; +} Index: trunk/psphot/src/psphotStackOptions.c =================================================================== --- trunk/psphot/src/psphotStackOptions.c (revision 28013) +++ trunk/psphot/src/psphotStackOptions.c (revision 28013) @@ -0,0 +1,98 @@ +# include "psphotInternal.h" + +static void psphotStackOptionsFree (psphotStackOptions *options) { + + if (options == NULL) return; + + // free the psf + psFree (options->psf); + + // free the array elements + psFree (options->psfs); + psFree (options->sourceLists); + psFree (options->kernels); + psFree (options->regions); + + // free the vector elements + psFree (options->inputMask); + psFree (options->inputSeeing); + psFree (options->norm); + psFree (options->matchChi2); + psFree (options->weightings); + + return; +} + +psphotStackOptions *psphotStackOptionsAlloc (int num) { + + psphotStackOptions *options = (psphotStackOptions *) psAlloc(sizeof(psphotStackOptions)); + psMemSetDeallocator(options, (psFreeFunc) psphotStackOptionsFree); + + options->numCols = 0; + options->numRows = 0; + + options->num = num; + options->psf = NULL; + options->convolve = false; + options->convolveSource = PSPHOT_CNV_SRC_NONE; + options->targetSeeing = NAN; + + options->psfs = psArrayAlloc(num); + options->sourceLists = psArrayAlloc(num); // Individual lists of sources for matching + options->kernels = psArrayAlloc(num); + options->regions = psArrayAlloc(num); + + options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs + options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32); + options->norm = psVectorAlloc(num, PS_TYPE_F32); + options->matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching + options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2) + + psVectorInit(options->inputMask, 0); + psVectorInit(options->inputSeeing, NAN); + psVectorInit(options->norm, NAN); + psVectorInit(options->matchChi2, NAN); + psVectorInit(options->weightings, NAN); + + return options; +} + +psphotStackConvolveSource psphotStackConvolveSourceFromString (const char *string) { + + if (!strcasecmp(string, "AUTO")) return PSPHOT_CNV_SRC_AUTO; + if (!strcasecmp(string, "CNV")) return PSPHOT_CNV_SRC_CNV; + if (!strcasecmp(string, "RAW")) return PSPHOT_CNV_SRC_RAW; + return PSPHOT_CNV_SRC_NONE; +} + +pmFPAfile *psphotStackGetConvolveSource (pmConfig *config, psphotStackOptions *options, int index) { + + // which image do we want to convolve? RAW, CNV, AUTO? + // find the currently selected readout + pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", index); // File of interest + pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", index); // File of interest + + pmFPAfile *fileSrc = NULL; + + switch (options->convolveSource) { + case PSPHOT_CNV_SRC_AUTO: + fileSrc = fileCnv ? fileCnv : fileRaw; + break; + + case PSPHOT_CNV_SRC_RAW: + fileSrc = fileRaw; + break; + + case PSPHOT_CNV_SRC_CNV: + fileSrc = fileCnv; + break; + + default: + psAbort("impossible case"); + } + if (!fileSrc) { + psError(PSPHOT_ERR_CONFIG, true, "desired convolution source is missing (cnv : %llx, raw : %llx)", (long long) fileCnv, (long long) fileRaw); + } + + return fileSrc; +} Index: trunk/psphot/src/psphotStackPSF.c =================================================================== --- trunk/psphot/src/psphotStackPSF.c (revision 28013) +++ trunk/psphot/src/psphotStackPSF.c (revision 28013) @@ -0,0 +1,64 @@ +# include "psphotInternal.h" + +pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask) +{ + bool mdok = false; + pmPSF *psf = NULL; + + // Get the recipe values + psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT"); // psphot recipe + psAssert(psphotRecipe, "We've thrown an error on this before."); + + bool autoPSF = psMetadataLookupBool (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.AUTO"); + + if (autoPSF) { + // Get the recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + int psfInstances = psMetadataLookupS32(NULL, recipe, "PSF.INSTANCES"); // Number of instances for PSF + float psfRadius = psMetadataLookupF32(NULL, recipe, "PSF.RADIUS"); // Radius for PSF + const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF + int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF + + psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in + if (!mdok || !maskValStr) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find MASK.VAL in recipe"); + return false; + } + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask + + for (int i = 0; i < psfs->n; i++) { + if (inputMask->data.U8[i]) { + psFree(psfs->data[i]); + psfs->data[i] = NULL; + } + } + + // Solve for the target PSF + psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal); + if (!psf) { + psError(PSPHOT_ERR_PSF, false, "Unable to determine output PSF."); + return NULL; + } + + } else { + + // externally-defined PSF + // XXX need to test for compatibility of target with inputs + + float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); + psAssert (isfinite(targetFWHM), "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM"); + + float Sxx = sqrt(2.0)*targetFWHM / 2.35; + + // XXX probably should make the model type (and par 7) optional from recipe + psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", Sxx, Sxx, 0.0, 1.0); + if (!psf) { + psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF."); + return NULL; + } + } + + return psf; +} Index: trunk/psphot/src/psphotStackParseCamera.c =================================================================== --- trunk/psphot/src/psphotStackParseCamera.c (revision 28010) +++ trunk/psphot/src/psphotStackParseCamera.c (revision 28013) @@ -25,40 +25,127 @@ psMetadata *input = item->data.md; // The input metadata of interest - // look for 'IMAGE', 'MASK', 'VARIANCE' in folder (only 'IMAGE' is required) - - psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image - if (!image || strlen(image) == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name); + pmFPAfile *rawInputFile = NULL; + pmFPAfile *cnvInputFile = NULL; + + // RAW (unconvolved) input data (RAW:IMAGE, RAW:MASK, RAW:VARIANCE, RAW:PSF) + psString rawImage = psMetadataLookupStr(&status, input, "RAW:IMAGE"); // Name of image + if (rawImage && strlen(rawImage) > 0) { + rawInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.RAW", rawImage, PM_FPA_FILE_IMAGE); // File for image + if (!rawInputFile) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, rawImage); + return false; + } + psString mask = psMetadataLookupStr(&status, input, "RAW:MASK"); // Name of mask + if (mask && strlen(mask) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.MASK.RAW", mask, PM_FPA_FILE_MASK)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); + return false; + } + } + psString variance = psMetadataLookupStr(&status, input, "RAW:VARIANCE"); // Name of variance map + if (variance && strlen(variance) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.VARIANCE.RAW", variance, PM_FPA_FILE_VARIANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); + return false; + } + } + psString psf = psMetadataLookupStr(&status, input, "RAW:PSF"); // Name of mask + if (psf && strlen(psf) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.PSF.RAW", psf, PM_FPA_FILE_PSF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + return false; + } + } + } + + // CNV (convolved) input data (CNV:IMAGE, CNV:MASK, CNV:VARIANCE, CNV:PSF) + psString cnvImage = psMetadataLookupStr(&status, input, "CNV:IMAGE"); // Name of image + if (cnvImage && strlen(cnvImage) > 0) { + cnvInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.CNV", cnvImage, PM_FPA_FILE_IMAGE); // File for image + if (!cnvInputFile) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, cnvImage); + return false; + } + psString mask = psMetadataLookupStr(&status, input, "CNV:MASK"); // Name of mask + if (mask && strlen(mask) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.MASK.CNV", mask, PM_FPA_FILE_MASK)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); + return false; + } + } + psString variance = psMetadataLookupStr(&status, input, "CNV:VARIANCE"); // Name of variance map + if (variance && strlen(variance) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.VARIANCE.CNV", variance, PM_FPA_FILE_VARIANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); + return false; + } + } + psString psf = psMetadataLookupStr(&status, input, "CNV:PSF"); // Name of mask + if (psf && strlen(psf) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.PSF.CNV", psf, PM_FPA_FILE_PSF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + return false; + } + } + } + + if (!rawInputFile && !cnvInputFile) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s (%d) lacks IMAGE of type STR", item->name, i); return false; } - pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image - if (!imageFile) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image); - return false; - } - - psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask - if (mask && strlen(mask) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); - return false; - } - } - - psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map - if (variance && strlen(variance) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); - return false; - } - } - // the output sources are carried on the input->fpa structures - pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT"); - if (!outsources) { - psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT"); - return false; - } - outsources->save = true; - outsources->fileID = i; // this is used to generate output names + + psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask + // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile; + if (sources && strlen(sources) > 0) { + if (!defineFile(config, NULL, "PSPHOT.STACK.SOURCES", sources, PM_FPA_FILE_CMF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)", i, sources); + return false; + } + } + + // generate an pmFPAimage for the output convolved image + // XXX output of these files should be optional + { + pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE"); + if (!outputImage) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.STACK.OUTPUT.IMAGE"); + return false; + } + outputImage->save = true; + outputImage->fileID = i; // this is used to generate output names + + pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK"); + if (!outputMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.MASK")); + return NULL; + } + if (outputMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.MASK is not of type MASK"); + return NULL; + } + outputMask->save = true; + outputMask->fileID = i; // this is used to generate output names + + pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE"); + if (!outputVariance) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.VARIANCE")); + return NULL; + } + if (outputVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.VARIANCE is not of type VARIANCE"); + return NULL; + } + outputVariance->save = true; + outputVariance->fileID = i; // this is used to generate output names + + // the output sources are carried on the outputImage->fpa structures + pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, outputImage, "PSPHOT.STACK.OUTPUT"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT"); + return false; + } + outsources->save = true; + outsources->fileID = i; // this is used to generate output names + } } psMetadataRemoveKey(config->arguments, "FILENAMES"); @@ -71,41 +158,35 @@ // generate an pmFPAimage for the chisqImage - pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE"); - if (!chisqImage) { - psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE"); - return false; - } - chisqImage->save = true; - - pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); - if (!chisqMask) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK")); - return NULL; - } - if (chisqMask->type != PM_FPA_FILE_MASK) { - psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK"); - return NULL; - } - chisqMask->save = true; - - pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); - if (!chisqVariance) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE")); - return NULL; - } - if (chisqVariance->type != PM_FPA_FILE_VARIANCE) { - psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE"); - return NULL; - } - chisqVariance->save = true; - -# if (0) - // define the additional input/output files associated with psphot - // XXX figure out which files are needed by psphotStack - if (false && !psphotDefineFiles (config, input)) { - psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files"); - return false; - } -# endif + // XXX output of these files should be optional + { + pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE"); + if (!chisqImage) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE"); + return false; + } + chisqImage->save = true; + + pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); + if (!chisqMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK")); + return NULL; + } + if (chisqMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK"); + return NULL; + } + chisqMask->save = true; + + pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); + if (!chisqVariance) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE")); + return NULL; + } + if (chisqVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE"); + return NULL; + } + chisqVariance->save = true; + } psTrace("psphot", 1, "Done with psphotStackParseCamera...\n"); @@ -145,2 +226,29 @@ } + +/*** + * + * psphotStack : + + * * inputs: + * * unconvolved images + * * raw convolved images + * * psfs (unconvolved or convolved?) + * * sources + + * optionally convolve the unconvolved or the raw inputs + * optionally perform no convolutions + * optionally save the psf-matched images + + */ + + +# if (0) + // define the additional input/output files associated with psphot + // XXX figure out which files are needed by psphotStack + if (false && !psphotDefineFiles (config, input)) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files"); + return false; + } +# endif + Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 28010) +++ trunk/psphot/src/psphotStackReadout.c (revision 28013) @@ -1,3 +1,6 @@ # include "psphotInternal.h" + +# define STACK_RAW "PSPHOT.STACK.INPUT.RAW" +# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE" bool psphotStackReadout (pmConfig *config, const pmFPAview *view) { @@ -17,6 +20,8 @@ PS_ASSERT_PTR_NON_NULL (breakPt, false); + // we have 3 relevant files: RAW, CNV, OUT + // set the photcode for each image - if (!psphotAddPhotcode (config, view)) { + if (!psphotAddPhotcode (config, view, STACK_OUT)) { psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -24,21 +29,23 @@ // Generate the mask and weight images - if (!psphotSetMaskAndVariance (config, view)) { - return psphotReadoutCleanup (config, view); + // XXX this should be done before we perform the convolutions + if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) { + return psphotReadoutCleanup (config, view, STACK_OUT); } if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } // generate a background model (median, smoothed image) // XXX I think this is not defined correctly for an array of images. - if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, view); + // XXX probably need to subtract the model (same model?) for both RAW and OUT + if (!psphotModelBackground (config, view, STACK_RAW)) { + return psphotReadoutCleanup (config, view, STACK_OUT); } - if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, view); + if (!psphotSubtractBackground (config, view, STACK_RAW)) { + return psphotReadoutCleanup (config, view, STACK_OUT); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } @@ -47,37 +54,43 @@ if (!psphotLoadPSF (config, view)) { psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } - if (!psphotStackChisqImage(config, view)) { + if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) { psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } if (!strcasecmp (breakPt, "CHISQ")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } // find the detections (by peak and/or footprint) in the image. // This finds the detections on Chisq image as well as the individuals - if (!psphotFindDetections (config, view, true)) { // pass 1 + if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1 // this only happens if we had an error in psphotFindDetections psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); + } + + // copy the detections from RAW to OUT + if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); + return psphotReadoutCleanup (config, view, STACK_OUT); } // construct sources and measure basic stats (saved on detections->newSources) // only run this on detections from the input images, not chisq image - if (!psphotSourceStats (config, view, true)) { // pass 1 + if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } - // *** generate the objects (which unify the sources from the different images) - psArray *objects = psphotMatchSources (config, view); + // generate the objects (object unify the sources from the different images) + psArray *objects = psphotMatchSources (config, view, STACK_OUT); // construct sources for the newly-generated sources (from other images) - if (!psphotSourceStats (config, view, false)) { // pass 1 + if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } @@ -85,5 +98,5 @@ // if (!psphotDeblendSatstars (config, view)) { // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - // return psphotReadoutCleanup (config, view); + // return psphotReadoutCleanup (config, view, STACK_OUT); // } @@ -91,61 +104,63 @@ // if (!psphotBasicDeblend (config, view)) { // psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); - // return psphotReadoutCleanup (config, view); + // return psphotReadoutCleanup (config, view, STACK_OUT); // } // classify sources based on moments, brightness // only run this on detections from the input images, not chisq image - if (!psphotRoughClass (config, view)) { + if (!psphotRoughClass (config, view, STACK_OUT)) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) // only run this on detections from the input images, not chisq image - if (!psphotImageQuality (config, view)) { // pass 1 + if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1 psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); - return psphotReadoutCleanup(config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, // this step is skipped - if (!psphotChoosePSF (config, view)) { // pass 1 + if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1 psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } - // include externally-supplied sources - // XXX fix this in the new multi-input context - // psphotLoadExtSources (config, view); // pass 1 - // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, view); + psphotGuessModels (config, view, STACK_OUT); // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - psphotMergeSources (config, view); + psphotMergeSources (config, view, STACK_OUT); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinearStack (config, objects, FALSE); - psFree (objects); // identify CRs and extended sources - psphotSourceSize (config, view, TRUE); + psphotSourceSize (config, view, STACK_OUT, TRUE); // measure aperture photometry corrections - if (!psphotApResid (config, view)) { + if (!psphotApResid (config, view, STACK_OUT)) { + psFree (objects); psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } + psphotStackObjectsUnifyPosition (objects); + psphotRadialAperturesByObject (config, objects, view, STACK_OUT); + + psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources) + psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources) + // calculate source magnitudes - psphotMagnitudes(config, view); + psphotMagnitudes(config, view, STACK_OUT); - if (!psphotEfficiency(config, view)) { + if (!psphotEfficiency(config, view, STACK_OUT)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -156,14 +171,16 @@ // replace background in residual image - psphotSkyReplace (config, view); + psphotSkyReplace (config, view, STACK_RAW); // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view); + psphotSourceFreePixels (config, view, STACK_OUT); // remove chisq image from config->file:PSPHOT.INPUT (why?) - psphotStackRemoveChisqFromInputs(config); + psphotStackRemoveChisqFromInputs(config, STACK_RAW); + + psFree (objects); // create the exported-metadata and free local data - return psphotReadoutCleanup (config, view); + return psphotReadoutCleanup (config, view, STACK_OUT); } Index: trunk/psphot/src/psphotSubtractBackground.c =================================================================== --- trunk/psphot/src/psphotSubtractBackground.c (revision 28010) +++ trunk/psphot/src/psphotSubtractBackground.c (revision 28013) @@ -4,5 +4,5 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) +bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { bool status = true; @@ -13,5 +13,5 @@ // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest pmFPA *inFPA = file->fpa; @@ -124,5 +124,5 @@ } -bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view) +bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = false; @@ -137,5 +137,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + if (!psphotSubtractBackgroundReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i); return false;