Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 32970) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 32971) @@ -57,4 +57,7 @@ #define BLANK_HEADERS "BLANK.HEADERS" // Name of metadata in camera configuration containing header names // for putting values into a blank PHU +static bool pmReadoutReadXSRC(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXFIT(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXRAD(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); // lookup the EXTNAME values used for table data and image header segments @@ -961,5 +964,27 @@ psString dataname = NULL; psString deteffname = NULL; - if (!pmSourceIOextnames(&headname, &dataname, &deteffname, NULL, NULL, NULL, file, view)) { + psString xsrcname = NULL; + psString xfitname = NULL; + psString xradname = NULL; + + // determine the output table format. Assume if we need to output extendend source + // parameters that they may exist in the input. + // XXX: Perhaps we should use different recipe values. + // I.E. EXTENDED_SOURCE_ANALYSIS_READ or something like that + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data"); + return false; + } + // if this is not TRUE, the output files only contain the psf measurements. + bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); + bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); + bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES"); + + if (!pmSourceIOextnames(&headname, &dataname, &deteffname, + XSRC_OUTPUT ? &xsrcname : NULL, + XFIT_OUTPUT ? &xfitname : NULL, + XRAD_OUTPUT ? &xradname : NULL, + file, view)) { return false; } @@ -1039,4 +1064,45 @@ } + long *sourceIndex = NULL; + if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT) { + long seq_max = -1; + for (long i = sources->n -1; i >= 0; i--) { + pmSource *source = sources->data[i]; + if (source->seq > seq_max) { + seq_max = source->seq; + } + } + sourceIndex = psAlloc((seq_max + 1) * sizeof(long)); + for (long i = 0; i < seq_max; i++) { + sourceIndex[i] = -1; + } + for (long i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + sourceIndex[source->seq] = i; + } + } + if (XSRC_OUTPUT && xsrcname) { + if (!pmReadoutReadXSRC(file, exttype, hdu->header, xsrcname, sources, sourceIndex)) { + // XXX: is this an error? + psErrorClear(); + } + psFree(xsrcname); + } + if (XFIT_OUTPUT && xfitname) { + if (!pmReadoutReadXFIT(file, exttype, hdu->header, xfitname, sources, sourceIndex)) { + // XXX: is this an error? + psErrorClear(); + } + psFree(xfitname); + } + if (XRAD_OUTPUT && xradname) { + if (!pmReadoutReadXRAD(file, exttype, hdu->header, xradname, sources, sourceIndex)) { + // XXX: is this an error? + psErrorClear(); + } + psFree(xradname); + } + psFree(sourceIndex); + if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) { #if 0 @@ -1165,3 +1231,179 @@ } - +// XXX: We might be able to macroize this and reuse for the other types + +static bool pmReadoutReadXSRC(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtName (file->fits, xsrcname)) { + psError(PS_ERR_UNKNOWN, false, "cannot find xsrc extension %s in %s", xsrcname, file->filename); + return false; + } + + psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", xsrcname); + return false; + } + + // XXX these are case-sensitive since the EXTYPE is case-sensitive + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { +#ifdef notyet + if (!strcmp (exttype, "SMPDATA")) { + status = pmSourcesRead_SMPDATA_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DEV_0")) { + status = pmSourcesRead_PS1_DEV_0_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DEV_1")) { + status = pmSourcesRead_PS1_DEV_1_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V1")) { + status = pmSourcesRead_CMF_PS1_V1_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesRead_CMF_PS1_V2_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V3")) { + status = pmSourcesRead_CMF_PS1_V3_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V4")) { + status = pmSourcesRead_CMF_PS1_V4_XSRC (file->fits, hduHeader, sources); + } +#endif // notyet + if (!strcmp (exttype, "PS1_SV1")) { + status = pmSourcesRead_CMF_PS1_SV1_XSRC (file->fits, hduHeader, sources, sourceIndex); + } +#ifdef notyet + if (!strcmp (exttype, "PS1_DV1")) { + status = pmSourcesRead_CMF_PS1_DV1_XSRC (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DV2")) { + status = pmSourcesRead_CMF_PS1_DV2_XSRC (file->fits, hduHeader, sources); + } +#endif // notyet + } + psFree(tableHeader); + return status; +} + +static bool pmReadoutReadXFIT(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtName (file->fits, extname)) { + psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename); + return false; + } + + psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", extname); + return false; + } + + // XXX these are case-sensitive since the EXTYPE is case-sensitive + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { +#ifdef notyet + if (!strcmp (exttype, "SMPDATA")) { + status = pmSourcesRead_SMPDATA_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DEV_0")) { + status = pmSourcesRead_PS1_DEV_0_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DEV_1")) { + status = pmSourcesRead_PS1_DEV_1_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V1")) { + status = pmSourcesRead_CMF_PS1_V1_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesRead_CMF_PS1_V2_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V3")) { + status = pmSourcesRead_CMF_PS1_V3_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V4")) { + status = pmSourcesRead_CMF_PS1_V4_XFIT (file->fits, hduHeader, sources); + } +#endif // notyet + if (!strcmp (exttype, "PS1_SV1")) { + status = pmSourcesRead_CMF_PS1_SV1_XFIT (file->fits, hduHeader, sources, sourceIndex); + } +#ifdef notyet + if (!strcmp (exttype, "PS1_DV1")) { + status = pmSourcesRead_CMF_PS1_DV1_XFIT (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DV2")) { + status = pmSourcesRead_CMF_PS1_DV2_XFIT (file->fits, hduHeader, sources); + } +#endif // notyet + } + psFree(tableHeader); + return status; +} +static bool pmReadoutReadXRAD(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtName (file->fits, extname)) { + psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename); + return false; + } + + psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", extname); + return false; + } + + // XXX these are case-sensitive since the EXTYPE is case-sensitive + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { +#ifdef notyet + if (!strcmp (exttype, "SMPDATA")) { + status = pmSourcesRead_SMPDATA_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DEV_0")) { + status = pmSourcesRead_PS1_DEV_0_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DEV_1")) { + status = pmSourcesRead_PS1_DEV_1_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V1")) { + status = pmSourcesRead_CMF_PS1_V1_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesRead_CMF_PS1_V2_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V3")) { + status = pmSourcesRead_CMF_PS1_V3_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_V4")) { + status = pmSourcesRead_CMF_PS1_V4_XRAD (file->fits, hduHeader, sources); + } +#endif // notyet + if (!strcmp (exttype, "PS1_SV1")) { + status = pmSourcesRead_CMF_PS1_SV1_XRAD (file->fits, hduHeader, sources, sourceIndex); + } +#ifdef notyet + if (!strcmp (exttype, "PS1_DV1")) { + status = pmSourcesRead_CMF_PS1_DV1_XRAD (file->fits, hduHeader, sources); + } + if (!strcmp (exttype, "PS1_DV2")) { + status = pmSourcesRead_CMF_PS1_DV2_XRAD (file->fits, hduHeader, sources); + } +#endif // notyet + } + psFree(tableHeader); + return status; +} Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 32970) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 32971) @@ -93,4 +93,7 @@ psArray *pmSourcesRead_CMF_PS1_V4 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header); +bool pmSourcesRead_CMF_PS1_SV1_XSRC (psFits *fits, psMetadata *header, psArray *sources, long *); +bool pmSourcesRead_CMF_PS1_SV1_XFIT (psFits *fits, psMetadata *header, psArray *sources, long *); +bool pmSourcesRead_CMF_PS1_SV1_XRAD (psFits *fits, psMetadata *header, psArray *sources, long *); psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 32970) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 32971) @@ -539,4 +539,98 @@ } +bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n"); + return false; + } + + // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point + float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP"); + float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); + float magOffset = zeropt + 2.5*log10(exptime); + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + pmSource *source = NULL; +#ifndef ASSUME_SORTED + long j = seq < sources->n ? seq : sources->n - 1; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } +#else + long j = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + source = sources->data[j]; +#endif + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + + if (!source->extpars) { + source->extpars = pmSourceExtendedParsAlloc (); + } + pmSourceExtendedPars *extpars = source->extpars; + + // Assume that X_EXT Y_EXT and sigmas match the psf src so skip + + // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by + // psf size or something? + extpars->axes.major = 1.0; + extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO"); + extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); + + float mag = psMetadataLookupF32(&status, row, "PETRO_MAG"); + float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR"); + if (isfinite(mag)) { + extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); + if (isfinite(magErr)) { + extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; + } + } + + extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); + extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR"); + extpars->petrosianR50 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50"); + extpars->petrosianR50Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR"); + extpars->petrosianR90 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90"); + extpars->petrosianR90Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR"); + extpars->petrosianFill = psMetadataLookupF32(&status, row, "PETRO_FILL"); + + psVector *radSB = psMetadataLookupVector(&status, row, "PROF_SB"); + psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX"); + psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL"); + + if (radSB && radSB->n > 0) { + extpars->radProfile = pmSourceRadialProfileAlloc(); + extpars->radProfile->binSB = psMemIncrRefCounter(radSB); + extpars->radProfile->binSum = psMemIncrRefCounter(radFlux); + extpars->radProfile->binFill = psMemIncrRefCounter(radFill); + } + + psFree(row); + } + + 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, psMetadata *imageHeader, char *extname) @@ -684,4 +778,104 @@ psFree (outhead); psFree (table); + return true; +} + +bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n"); + return false; + } + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (!source->modelFits) { + // XXX: where to find the number of models to expect? + source->modelFits = psArrayAllocEmpty(5); + } + psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE"); + if (!modelName) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i); + psFree(row); + return false; + } + pmModelType modelType = pmModelClassGetType(modelName); + if (modelType < 0) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName); + psFree(row); + return false; + } + pmModel *model = pmModelAlloc(modelType); + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG"); + + model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG"); + model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); + + psEllipseAxes axes; + axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); + axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); + axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); + if (!pmPSF_AxesToModel(PAR, axes, modelType)) { + // Do we need to fail here or can this happen? + psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); + psFree(model); + psFree(row); + return false; + } + // XXX: clean this up + if (model->params->n > 7) { + PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); + } + // read the covariance matrix + int nparams = model->params->n; + psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); + for (int y = 0; y < nparams; y++) { + for (int x = 0; x < nparams; x++) { + char name[64]; + snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); + covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + } + } + model->covar = covar; + + psArrayAdd(source->modelFits, 1, model); + psFree(model); + + psFree(row); + } + return true; } @@ -831,2 +1025,57 @@ return true; } + +bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); + return false; + } + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (!source->radialAper) { + // XXX: where to find the number of models to expect? + source->radialAper = psArrayAllocEmpty(5); + } + pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc(); + + radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX")); + radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV")); + radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR")); + radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL")); + + psArrayAdd(source->radialAper, 1, radialAper); + + psFree(radialAper); + psFree(row); + } + + return true; +}