Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 15016) +++ trunk/psModules/src/objects/Makefile.am (revision 15039) @@ -23,4 +23,5 @@ pmSourceIO_SMPDATA.c \ pmSourceIO_PS1_DEV_0.c \ + pmSourceIO_PS1_DEV_1.c \ pmSourcePlots.c \ pmSourcePlotPSFModel.c \ Index: trunk/psModules/src/objects/pmPeaks.c =================================================================== --- trunk/psModules/src/objects/pmPeaks.c (revision 15016) +++ trunk/psModules/src/objects/pmPeaks.c (revision 15039) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-24 00:11:02 $ + * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -51,4 +51,5 @@ // calculate peak position relative to ix,iy + // XXX these functions need to take a mask, weight, and calculate the errors psPolynomial2D *bicube = psImageBicubeFit (image, ix + image->col0, iy + image->row0); psPlane min = psImageBicubeMin (bicube); @@ -59,7 +60,11 @@ peak->xf = min.x + ix + image->col0; peak->yf = min.y + iy + image->row0; + peak->dx = 0.0; + peak->dy = 0.0; } else { peak->xf = ix; peak->yf = iy; + peak->dx = 1.0; + peak->dy = 1.0; } @@ -189,4 +194,36 @@ } psTrace("psModules.objects", 3, "---- %s(0) end ----\n", __func__); + return (0); +} + +// sort by SN (descending) +int pmPeakSortBySN (const void **a, const void **b) +{ + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 fA = A->SN; + psF32 fB = B->SN; + if (isnan (fA)) fA = 0; + if (isnan (fB)) fB = 0; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} + +// sort by Y (ascending) +int pmPeakSortByY (const void **a, const void **b) +{ + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 fA = A->y; + psF32 fB = B->y; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); return (0); } Index: trunk/psModules/src/objects/pmPeaks.h =================================================================== --- trunk/psModules/src/objects/pmPeaks.h (revision 15016) +++ trunk/psModules/src/objects/pmPeaks.h (revision 15039) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-24 00:11:02 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -53,4 +53,6 @@ float xf; ///< bicube fit to peak coord (x) float yf; ///< bicube fit to peak coord (y) + float dx; ///< bicube fit error on peak coord (x) + float dy; ///< bicube fit error on peak coord (y) float value; ///< level in detection image float flux; ///< level in unsmoothed sci image @@ -148,4 +150,7 @@ int pmPeaksCompareDescend (const void **a, const void **b); +int pmPeakSortBySN (const void **a, const void **b); +int pmPeakSortByY (const void **a, const void **b); + /// @} # endif /* PM_PEAKS_H */ Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 15016) +++ trunk/psModules/src/objects/pmSource.c (revision 15039) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:21:57 $ + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -87,4 +87,5 @@ pmSource *source = (pmSource *) psAlloc(sizeof(pmSource)); *(int *)&source->id = id++; + source->seq = -1; source->peak = NULL; source->pixels = NULL; @@ -111,4 +112,8 @@ source->skyErr = NAN; source->pixWeight = NAN; + + source->psfProb = NAN; + source->crNsigma = NAN; + source->extNsigma = NAN; psTrace("psModules.objects", 5, "---- end ----\n"); @@ -930,2 +935,35 @@ return NULL; } + +// sort by SN (descending) +int pmSourceSortBySN (const void **a, const void **b) +{ + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; + + psF32 fA = (A->peak == NULL) ? 0 : A->peak->SN; + psF32 fB = (B->peak == NULL) ? 0 : B->peak->SN; + if (isnan (fA)) fA = 0; + if (isnan (fB)) fB = 0; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} + +// sort by Y (ascending) +int pmSourceSortByY (const void **a, const void **b) +{ + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; + + psF32 fA = (A->peak == NULL) ? 0 : A->peak->y; + psF32 fB = (B->peak == NULL) ? 0 : B->peak->y; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); +} + Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 15016) +++ trunk/psModules/src/objects/pmSource.h (revision 15039) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-24 00:11:02 $ + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -59,4 +59,5 @@ struct pmSource { const int id; ///< Unique ID for object + int seq; ///< ID for output (generated on write) pmPeak *peak; ///< Description of peak pixel. psImage *pixels; ///< Rectangular region including object pixels. @@ -78,4 +79,7 @@ float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) float pixWeight; ///< model-weighted coverage of valid pixels + float psfProb; ///< probability of PSF + float crNsigma; ///< Nsigma deviation from PSF to CR + float extNsigma; ///< Nsigma deviation from PSF to EXT psRegion region; ///< area on image covered by selected pixels float sky, skyErr; ///< The sky and its error at the center of the object @@ -223,4 +227,7 @@ bool pmSourceCachePSF (pmSource *source, psMaskType maskVal); +int pmSourceSortBySN (const void **a, const void **b); +int pmSourceSortByY (const void **a, const void **b); + /// @} # endif /* PM_SOURCE_H */ Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 15016) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 15039) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:09:05 $ + * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -388,4 +388,7 @@ status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname); } + if (!strcmp (exttype, "PS1_DEV_1")) { + status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); + } if (!status) { @@ -708,4 +711,7 @@ sources = pmSourcesRead_PS1_DEV_0 (file->fits, hdu->header); } + if (!strcmp (exttype, "PS1_DEV_1")) { + sources = pmSourcesRead_PS1_DEV_1 (file->fits, hdu->header); + } psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname); Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 15016) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 15039) @@ -4,6 +4,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-14 03:20:44 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * @@ -26,4 +26,5 @@ bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); +bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); @@ -33,4 +34,5 @@ psArray *pmSourcesRead_SMPDATA (psFits *fits, psMetadata *header); psArray *pmSourcesRead_PS1_DEV_0 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_PS1_DEV_1 (psFits *fits, psMetadata *header); bool pmSourcesWritePSFs (psArray *sources, char *filename); Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 15039) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 15039) @@ -0,0 +1,224 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-27 03:35:29 $ + * + * 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) +{ + + psArray *table; + psMetadata *row; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // 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_PROBABILITY", PS_DATA_F32, "Probability of PSF-ness", source->psfProb); + 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); + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_PS1_DEV_1 (psFits *fits, psMetadata *header) +{ + + 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); + } + + psArray *table = psFitsReadTable (fits); + // validate a single row of the table (must match SMP) + + psArray *sources = psArrayAlloc (table->n); + + // convert the table to the pmSource entries + // XXX need to chooose PSF vs EXT, based on type? + for (int i = 0; i < table->n; i++) { + 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; + + psMetadata *row = table->data[i]; + + 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"); + + source->psfProb = psMetadataLookupF32 (&status, row, "PSF_PROBABILITY"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + + // XXX other values saved but not loaded? + // psMetadataLookupS64 (&status, row, "IPP_IDET"); + // psMetadataLookupF32 (&status, row, "N_FRAMES"); + + sources->data[i] = source; + } + psFree (table); + return (sources); +}