Index: /branches/unlabeled-1.8.4/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- /branches/unlabeled-1.8.4/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 16295) +++ /branches/unlabeled-1.8.4/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 16295) @@ -0,0 +1,402 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-03 22:04:47 $ + * + * 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 "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 + +// this output format is valid for psphot analysis of an image, and does not include calibrated +// values derived in the DVO database. +// XXX how do I generate the source tables which I need to send to PSPS? + +bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, + psMetadata *tableHeader, char *extname, char *xsrcname) +{ + 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; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // if we request XSRC output, add the XSRC name to this header + if (xsrcname) { + psMetadataAddStr (tableHeader, PS_LIST_TAIL, "XSRCNAME", PS_META_REPLACE, "name of XSRC table extension", xsrcname); + } + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + source->seq = i; + + // no difference between PSF and non-PSF model + pmModel *model = pmSourceGetModel (NULL, source); + + if (model != NULL) { + 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); + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + axes.major = 0.0; + axes.minor = 0.0; + axes.theta = 0.0; + } + + float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + psS16 nImageOverlap = 1; + + 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_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); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", PS_MIN (99.0, source->psfMag)); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", PS_MIN (99.0, source->errMag)); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", PS_MIN (99.0, peakMag)); + 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 to top 9 pixels", source->psfChisq); + 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); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_X", PS_DATA_F32, "PSF width in x coordinate", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_Y", PS_DATA_F32, "PSF width in y coordinate", 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); + + // 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, "FLAGS", PS_DATA_U16, "psphot analysis flags", source->mode); + + // XXX these calibrated values are not supplied by psphot analysis + // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "RA from PSF fit", RA); + // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "DEC from PSF fit", DEC); + // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF_SIG", PS_DATA_F32, "Sigma of PSF fit RA", dRA); + // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF_SIG", PS_DATA_F32, "Sigma of PSF fit DEC", dDEC); + // psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "Calibrated magnitude", calMag); + // psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_sIG", PS_DATA_F32, "Sigma of calibrated magnitude", calMagErr); + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, tableHeader, extname); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, tableHeader, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + return false; + } + psFree (table); + + if (xsrcname) { + pmSourcesWrite_PS1_DEV_1_XSRC (fits, sources, xsrcname); + } + + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_PS1_DEV_1 (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 + 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); + + // XXX need to look up the XSRCNAME entries + + // validate a single row of the table (must match SMP) + + // XXX test return values + + // XXX we have a memory problem, which is illustrated here: if I allocate the sources array + // (line 1) before freeing the table, the data is not really freed (but not a leak?) if I + // allocate the sources array *after* freeing the table, the data is actually freed + // psArray *fooSources = psArrayAllocEmpty (10000); + // psFree (table); + // return (fooSources); + + + // 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. 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 entries + // XXX need to chooose PSF vs EXT, based on type? + for (int i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; + + // 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; + + 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_WIDTH_X"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_WIDTH_Y"); + 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"); + + pmPSF_AxesToModel (PAR, axes); + + float lflux = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + float flux = (isfinite(lflux)) ? pow(10.0, -0.4*lflux) : NAN; + source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], flux, PM_PEAK_LONE); + source->peak->flux = flux; + + source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + source->psfChisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + + source->mode = psMetadataLookupU16 (&status, row, "FLAGS"); + assert (status); + + // XXX other values saved but not loaded? + // psMetadataLookupS64 (&status, row, "IPP_IDET"); + // psMetadataLookupF32 (&status, row, "N_FRAMES"); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname) +{ + + psArray *table; + psMetadata *row; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // 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); + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (i = 0; i < sources->n; i++) { + // source->seq is set during the PSF output step, which must be called + // XXX test here if seq is set? + + // 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]; + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelConv; + if (model == NULL) { + model = source->modelEXT; + } + if (model == NULL) continue; + + // XXX do we need to do this? + if (source->extpars == NULL) continue; + + 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) + 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); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); + + // XXX at the moment, this will be a fixed size table, but perhaps have multiple output formats depending on what is measured? + + // other values that I need to report: + // R, Mag Petrosian, ..... + if (source->extpars) { + + // Petrosian measurements + pmSourcePetrosianValues *petrosian = source->extpars->petrosian; + if (petrosian) { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + } + + // Kron measurements + 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); + } + + // Isophot measurements + // XXX insert header data: isophotal level + 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); + } + + pmSourceAnnuli *annuli = source->extpars->annuli; + if (annuli) { + psVector *fluxVal = annuli->flux; + psVector *fluxErr = annuli->fluxErr; + psVector *fluxVar = annuli->fluxVar; + + for (int j = 0; j < fluxVal->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); + } + } + } + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, outhead, extname); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + return false; + } + psFree (table); + + return true; +}