Index: trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- trunk/psModules/src/camera/pmFPAfile.c (revision 36375) +++ trunk/psModules/src/camera/pmFPAfile.c (revision 36441) @@ -565,4 +565,6 @@ case PM_FPA_FILE_CMF: return ("CMF"); + case PM_FPA_FILE_CFF: + return ("CFF"); case PM_FPA_FILE_WCS: return ("WCS"); Index: trunk/psModules/src/objects/pmModelClass.c =================================================================== --- trunk/psModules/src/objects/pmModelClass.c (revision 36375) +++ trunk/psModules/src/objects/pmModelClass.c (revision 36441) @@ -66,4 +66,5 @@ static pmModelClass *models = NULL; +static psVector *modelClassLookupTable = NULL; // translation between model types in header and here static int Nmodels = 0; @@ -135,4 +136,6 @@ models = NULL; Nmodels = 0; + psFree(modelClassLookupTable); + modelClassLookupTable = NULL; return; } @@ -193,2 +196,80 @@ } + +bool pmModelClassWriteHeader(psMetadata *header) +{ + psMetadataAddS32(header, PS_LIST_TAIL, "MTNUM", PS_META_REPLACE, "number of model types", Nmodels); + for (int i = 0; i < Nmodels; i++) { + char modelNameKey[16]; + char modelValKey[16]; + sprintf(modelNameKey, "MTNAM%02d", i); + sprintf(modelValKey, "MTVAL%02d", i); + psMetadataAddStr(header, PS_LIST_TAIL, modelNameKey, PS_META_REPLACE, "", models[i].name); + psMetadataAddS32(header, PS_LIST_TAIL, modelValKey, PS_META_REPLACE, "", i); + } + + return true; +} + +bool pmModelClassReadHeader(psMetadata *header) { + psFree(modelClassLookupTable); + + bool status; + int numHeaderModels = psMetadataLookupS32(&status, header, "MTNUM"); + if (!status) { + return false; + } + + psVector *inputTypes = psVectorAlloc(numHeaderModels, PS_TYPE_S32); + psVector *localTypes = psVectorAlloc(numHeaderModels, PS_TYPE_S32); + int max_val = -1; + for (int i = 0; i < numHeaderModels; i++) { + char modelNameKey[16]; + char modelValKey[16]; + sprintf(modelNameKey, "MTNAM%02d", i); + sprintf(modelValKey, "MTVAL%02d", i); + psString thisName = psMetadataLookupStr(&status, header, modelNameKey); + int thisVal = psMetadataLookupS32(&status, header, modelValKey); + if (thisVal > max_val) { + max_val = thisVal; + } + inputTypes->data.S32[i] = thisVal; + localTypes->data.S32[i] = pmModelClassGetType(thisName); + } + if (max_val < 0) { + psFree(inputTypes); + psFree(localTypes); + return false; + } + + modelClassLookupTable = psVectorAlloc(max_val + 1, PS_TYPE_S32); + psVectorInit(modelClassLookupTable, -1); + + for (int i = 0; i < numHeaderModels; i++) { + int thisVal = inputTypes->data.S32[i]; + int localVal = localTypes->data.S32[i]; + modelClassLookupTable->data.S32[thisVal] = localVal; + } + psFree(inputTypes); + psFree(localTypes); + + return true; +} + +pmModelType pmModelClassGetLocalType(pmModelType inputType) { + pmModelType localType = -1; + + if (modelClassLookupTable) { + if (inputType >= 0 && inputType < modelClassLookupTable->n) { + localType = modelClassLookupTable->data.S32[inputType]; + } + } else { + // no lookup table defined + // for backwards compatability if inputType refers to a defined model, return it + if (inputType >= 0 && pmModelClassGetName(inputType)) { + localType = inputType; + } + } + + return localType; +} Index: trunk/psModules/src/objects/pmModelClass.h =================================================================== --- trunk/psModules/src/objects/pmModelClass.h (revision 36375) +++ trunk/psModules/src/objects/pmModelClass.h (revision 36441) @@ -76,4 +76,10 @@ void pmModelClassSetLimits(pmModelLimitsType type); +// write keywords to header definining the model type values used by this program +bool pmModelClassWriteHeader(psMetadata *header); +// create a lookup table for translating input model type values to local model type values +bool pmModelClassReadHeader(psMetadata *header); +// translate input model type value to local value +pmModelType pmModelClassGetLocalType(pmModelType inputType); /// @} Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 36441) @@ -61,4 +61,5 @@ static bool pmReadoutReadXFIT(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXGAL(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); // lookup the EXTNAME values used for table data and image header segments @@ -374,5 +375,5 @@ } \ if (xgalname) { \ - status &= pmSourcesWrite_##TYPE##_XGAL (file->fits, sources, xgalname, recipe); \ + status &= pmSourcesWrite_##TYPE##_XGAL (file->fits, readout, sources, xgalname, recipe); \ } \ } @@ -1036,5 +1037,5 @@ bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES"); - bool XGAL_OUTPUT = false; // psMetadataLookupBool(&status, recipe, "GALAXY_SHAPES"); + bool XGAL_OUTPUT = psMetadataLookupBool(&status, recipe, "GALAXY_SHAPES"); if (!pmSourceIOextnames(&headname, &dataname, &deteffname, @@ -1122,5 +1123,9 @@ long *sourceIndex = NULL; - if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT) { + if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT || XGAL_OUTPUT) { + // Build sourceIndex. Lookup table from source->seq to index in sources array. + // Consists of an array of length max(source->seq) + 1. + + // find maximum sequence number long seq_max = -1; for (long i = sources->n -1; i >= 0; i--) { @@ -1135,8 +1140,10 @@ } } + // allocate and initialize the index sourceIndex = psAlloc((seq_max + 1) * sizeof(long)); for (long i = 0; i < seq_max; i++) { sourceIndex[i] = -1; } + // populate the index for (long i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; @@ -1165,4 +1172,11 @@ psFree(xradname); } + if (XGAL_OUTPUT && xgalname) { + // a cmf file may have an XGAL extension, but it is not required + if (!pmReadoutReadXGAL(file, readout, exttype, hdu->header, xgalname, sources, sourceIndex)) { + // do anything? + } + psFree(xgalname); + } psFree(sourceIndex); @@ -1461,2 +1475,40 @@ return status; } +static bool pmReadoutReadXGAL(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString xgalname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtNameClean (file->fits, xgalname)) { + psTrace ("pmFPAfile", 1, "cannot find xgal extension %s in %s, skipping", xgalname, 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")) { + psFree(tableHeader); + psWarning ("no binary table in extension %s, skipping\n", xgalname); + return false; + } + +# define PM_SOURCES_READ_XGAL(NAME,TYPE) \ + if (!strcmp (exttype, NAME)) { \ + status = pmSourcesRead_##TYPE##_XGAL(file->fits, readout, hduHeader, tableHeader, sources, sourceIndex); \ + } + + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { + PM_SOURCES_READ_XGAL("PS1_V1", CMF_PS1_V1); + PM_SOURCES_READ_XGAL("PS1_V2", CMF_PS1_V2); + PM_SOURCES_READ_XGAL("PS1_V3", CMF_PS1_V3); + PM_SOURCES_READ_XGAL("PS1_V4", CMF_PS1_V4); + PM_SOURCES_READ_XGAL("PS1_SV1", CMF_PS1_SV1); + PM_SOURCES_READ_XGAL("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_READ_XGAL("PS1_DV1", CMF_PS1_DV1); + PM_SOURCES_READ_XGAL("PS1_DV2", CMF_PS1_DV2); + PM_SOURCES_READ_XGAL("PS1_DV3", CMF_PS1_DV3); + } + psFree(tableHeader); + return status; +} Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 36441) @@ -21,9 +21,10 @@ bool pmSourcesWrite_##TYPE##_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); \ bool pmSourcesWrite_##TYPE##_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); \ - bool pmSourcesWrite_##TYPE##_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe); \ + bool pmSourcesWrite_##TYPE##_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe); \ psArray *pmSourcesRead_##TYPE (psFits *fits, psMetadata *header); \ bool pmSourcesRead_##TYPE##_XSRC (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \ bool pmSourcesRead_##TYPE##_XFIT (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \ bool pmSourcesRead_##TYPE##_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index);\ + bool pmSourcesRead_##TYPE##_XGAL (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index);\ // All of these functions need to use the same API, even if not all elements are used in a specific case Index: trunk/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36441) @@ -65,4 +65,7 @@ int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + // Read lookup table for model classes (if defined) + pmModelClassReadHeader(header); + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); if (PSF_NAME != NULL) { @@ -111,9 +114,10 @@ float theta = psMetadataLookupF32 (&status, row, "THETA"); - // XXX: we need to put a lookup table in the cff header to define the correspondence of the - // model type values in the cff with our models. (We want to use an interger for efficiency - // but the value for each model type is set on the organization of the the array in pmModelClass.c - // For now use the input values verbatim and trust the user that this is valid value int galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); + if (status) { + galaxyModelType = pmModelClassGetLocalType(galaxyModelType); + } else { + galaxyModelType = -1; + } float Sindex = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index @@ -123,5 +127,6 @@ source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS + // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and + // stars that are masked if (psfStar) { source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; @@ -180,7 +185,5 @@ } - // XXX: should use < 0 as invalid galaxyModelType - - if (fitGalaxy && galaxyModelType > 0) { + if (fitGalaxy && galaxyModelType >= 0) { source->modelFits = psArrayAllocEmpty (1); pmModel *model = pmModelAlloc(galaxyModelType); @@ -235,4 +238,8 @@ PS_ASSERT(mdok, false); + // write the definition of the model class type values to the header + psMetadata *outputHeader = psMetadataAlloc(); + pmModelClassWriteHeader(outputHeader); + psArray *table = psArrayAllocEmpty(sources->n); @@ -246,7 +253,8 @@ psS32 modelType = 0; bool fitGalaxy = false; - bool psfStar = false; + bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; psF32 sersicIndex = 0; - if (source->modelFits == NULL) { + // For now only perform galaxy fits on extended objects + if (source->modelEXT == NULL) { pmModel *model = source->modelPSF; if (model == NULL) continue; @@ -260,5 +268,4 @@ yPos = model->params->data.F32[PM_PAR_YPOS]; flux = source->psfFlux; - psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; rMajor = 0; rMinor = 0; @@ -329,14 +336,15 @@ psArrayAdd(table, 100, row); + psFree(row); } - if (!psFitsWriteTable(fits, NULL, table, extname)) { + if (!psFitsWriteTable(fits, outputHeader, table, extname)) { psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); psFree(table); - psFree(header); + psFree(outputHeader); return false; } psFree(table); - // psFree(header); + psFree(outputHeader); return true; Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36441) @@ -708,21 +708,9 @@ return false; } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq + // Find the source with this sequence number using the sourceIndex. 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 + pmSource *source = sources->data[j]; if (!source) { psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); @@ -793,4 +781,6 @@ // create a header to hold the output data psMetadata *outhead = psMetadataAlloc (); + + pmModelClassWriteHeader(outhead); // write the links to the image header @@ -1034,4 +1024,7 @@ return false; } + // set up the lookup table to translate between input model types and output model types + // if not defined it is assumed that the tables are the same + pmModelClassReadHeader(tableHeader); for (long i = 0; i < numSources; i++) { @@ -1042,15 +1035,8 @@ 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; - } - } + long j = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + pmSource *source = sources->data[j]; if (!source) { psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); @@ -1095,8 +1081,11 @@ // in the psf table. psS32 extModelType = psMetadataLookupS32(&status, row, "EXT_MODEL_TYPE"); - if (!status) { + if (status) { + // translate between the type value in xfit and values used by this program + extModelType = pmModelClassGetLocalType(extModelType); + } else { // older cmfs don't have this column extModelType = -1; - } + } psEllipseAxes axes; @@ -1339,15 +1328,8 @@ 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; - } - } + long j = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + pmSource *source = sources->data[j]; if (!source) { psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); @@ -1407,5 +1389,5 @@ // XXX where should I record the number of columns?? -bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe) { bool status = false; @@ -1422,4 +1404,6 @@ // write the links to the image header psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname); + + psMetadataAddStr (outhead, PS_LIST_TAIL, "HI", PS_META_REPLACE, "does this get through?", "THERE"); // let's write these out in S/N order @@ -1491,2 +1475,56 @@ } +bool pmSourcesRead_CMF_@CMFMODE@_XGAL(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psMetadata *tableHeader, 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, "XGAL 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 = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + + pmSource *source = sources->data[j]; + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + + psVector *Flux = psMetadataLookupVector(&status, row, "GAL_FLUX"); + psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR"); + psVector *chisq = psMetadataLookupVector(&status, row, "GAL_CHISQ"); + + if (Flux && Flux->n > 0) { + psFree(source->galaxyFits); + source->galaxyFits = pmSourceGalaxyFitsAlloc(); + source->galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX"); + + psFree(source->galaxyFits->Flux); + source->galaxyFits->Flux = psMemIncrRefCounter(Flux); + psFree(source->galaxyFits->dFlux); + source->galaxyFits->dFlux = psMemIncrRefCounter(dFlux); + psFree(source->galaxyFits->chisq); + source->galaxyFits->chisq = psMemIncrRefCounter(chisq); + } + + psFree(row); + } + + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 36441) @@ -714,5 +714,5 @@ } -bool pmSourcesWrite_PS1_CAL_0_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_PS1_CAL_0_XGAL (psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe) { return true; Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 36441) @@ -256,6 +256,6 @@ } -bool pmSourcesWrite_PS1_DEV_0_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe) -{ - return true; -} +bool pmSourcesWrite_PS1_DEV_0_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe) +{ + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 36441) @@ -596,5 +596,5 @@ } -bool pmSourcesWrite_PS1_DEV_1_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_PS1_DEV_1_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe) { return true; Index: trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 36441) @@ -226,5 +226,5 @@ } -bool pmSourcesWrite_SMPDATA_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_SMPDATA_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe) { return true;