Index: trunk/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36658) +++ trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36707) @@ -53,5 +53,5 @@ #include "pmSourceOutputs.h" -// read in a readout from the fits file +// read in sources readout from a cff fits file psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header) { @@ -78,9 +78,11 @@ // 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++) { + long numRows = psFitsTableSize(fits); // Number of rows in table + psArray *sources = psArrayAllocEmpty(numRows); // Array of sources, to return + + // convert the table to the pmSource entries + unsigned int ID_last = 0; + pmSource *source = NULL; + for (int i = 0; i < numRows; i++) { psMetadata *row = psFitsReadTableRow(fits, i); // Table row if (!row) { @@ -123,79 +125,89 @@ float Sindex = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; -// RoughClass wants source type to be unknown -// source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - source->type = PM_SOURCE_TYPE_UNKNOWN; - - // 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; - } - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - source->seq = ID; - PAR[PM_PAR_XPOS] = X; - PAR[PM_PAR_YPOS] = Y; - - dPAR[PM_PAR_XPOS] = 0.0; - dPAR[PM_PAR_YPOS] = 0.0; - - PAR[PM_PAR_SKY] = 0.0; - dPAR[PM_PAR_SKY] = 0.0; - - PAR[PM_PAR_I0] = 1.0; - dPAR[PM_PAR_I0] = 0.0; - - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - source->psfMag = 0.0; - source->psfMagErr = 0.0; - source->apMag = 0.0; - source->apRadius = apRadius; - - // we generate a somewhat fake PSF model here -- - // in most (all?) contexts, we will replace this with a measured psf model - // elsewhere - axes.major = 1.0; - axes.minor = 1.0; - axes.theta = 0.0; - pmPSF_AxesToModel (PAR, axes, modelType); - - // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second - source->peak = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE); - source->peak->xf = X; // pmPeakAlloc converts X,Y to int, so reset here - source->peak->yf = Y; // pmPeakAlloc converts X,Y to int, so reset here - source->peak->dx = 0.0; - source->peak->dy = 0.0; - - source->moments = pmMomentsAlloc (); - source->moments->Mx = X; - source->moments->My = Y; - source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment - - // Don't mark the moments as measured because that causes many fields to be left blank. - // The moments code knows not to change the position or the Mrf for external sources - // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; - - if (isfinite(petRadius)) { - source->extpars = pmSourceExtendedParsAlloc (); - source->extpars->petrosianRadius = petRadius; - } - + if (!source || ID != source->seq) { + if (source) { + psArrayAdd (sources, 1, source); + psFree(source); + } + source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + // RoughClass wants source type to be unknown + // source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + source->type = PM_SOURCE_TYPE_UNKNOWN; + + // 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; + } + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + source->seq = ID; + ID_last = ID; + + PAR[PM_PAR_XPOS] = X; + PAR[PM_PAR_YPOS] = Y; + + dPAR[PM_PAR_XPOS] = 0.0; + dPAR[PM_PAR_YPOS] = 0.0; + + PAR[PM_PAR_SKY] = 0.0; + dPAR[PM_PAR_SKY] = 0.0; + + PAR[PM_PAR_I0] = 1.0; + dPAR[PM_PAR_I0] = 0.0; + + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + source->psfMag = 0.0; + source->psfMagErr = 0.0; + source->apMag = 0.0; + source->apRadius = apRadius; + + // we generate a somewhat fake PSF model here -- + // in most (all?) contexts, we will replace this with a measured psf model + // elsewhere + axes.major = 1.0; + axes.minor = 1.0; + axes.theta = 0.0; + pmPSF_AxesToModel (PAR, axes, modelType); + + // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second + source->peak = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE); + source->peak->xf = X; // pmPeakAlloc converts X,Y to int, so reset here + source->peak->yf = Y; // pmPeakAlloc converts X,Y to int, so reset here + source->peak->dx = 0.0; + source->peak->dy = 0.0; + + source->moments = pmMomentsAlloc (); + source->moments->Mx = X; + source->moments->My = Y; + source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment + + // Don't mark the moments as measured because that causes many fields to be left blank. + // The moments code knows not to change the position or the Mrf for external sources + // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + + if (isfinite(petRadius)) { + source->extpars = pmSourceExtendedParsAlloc (); + source->extpars->petrosianRadius = petRadius; + } + + } if (fitGalaxy && galaxyModelType >= 0) { - source->modelFits = psArrayAllocEmpty (1); - pmModel *model = pmModelAlloc(galaxyModelType); - psF32 *xPAR = model->params->data.F32; - - xPAR[PM_PAR_SKY] = 0.0; - xPAR[PM_PAR_I0] = 1.0; - xPAR[PM_PAR_XPOS] = X; + if (!source->modelFits) { + source->modelFits = psArrayAllocEmpty (1); + } + pmModel *model = pmModelAlloc(galaxyModelType); + psF32 *xPAR = model->params->data.F32; + + xPAR[PM_PAR_SKY] = 0.0; + xPAR[PM_PAR_I0] = 1.0; + xPAR[PM_PAR_XPOS] = X; xPAR[PM_PAR_YPOS] = Y; @@ -212,22 +224,13 @@ psArrayAdd (source->modelFits, 1, model); -#ifdef notyet - // XXX: set source->modelEXT to this model and flag as extended so that we have a better - // shot at subtracting extended sources? - - // This doesn't work right. Need to do some more work to flesh out the model before it can be - // used in subtaction. Better idea might be to do 2 passes in psphotFullForceReadout to - // First find the best extended models and then do everything else. - source->modelEXT = psMemIncrRefCounter(model); - // is this safe? (no see above) - source->type = PM_SOURCE_TYPE_EXTENDED; - source->mode |= PM_SOURCE_MODE_EXTMODEL; -#endif - psFree (model); } - sources->data[i] = source; psFree(row); + } + if (source) { + // close out last source + psArrayAdd (sources, 1, source); + psFree(source); } @@ -250,12 +253,30 @@ pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModelType devModelType = pmModelClassGetType("PS_MODEL_DEV"); + pmModelType selectedModelType = -1; + bool chooseBest = false; + bool chooseAll = false; psString modelToChoose = psMetadataLookupStr(&mdok, recipe, "EXT_MODEL_TYPE_FOR_CFF"); - pmModelType selectedModelType = -1; + if (mdok && modelToChoose != NULL) { - if (strcmp(modelToChoose, "BEST")) { + if (!strcmp(modelToChoose, "BEST")) { + chooseBest = true; + } else if (!strcmp(modelToChoose, "ALL")) { + chooseAll = true; + } else if (strcmp(modelToChoose, "PS_MODEL_SERSIC")) { + // We have selected a model type other than Sersic. + // Save it's type for use below. Sersic is handled specially selectedModelType = pmModelClassGetType(modelToChoose); } } + + // minimum sersic index to force devModel + psF32 sersicMinDev = psMetadataLookupF32(&mdok, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN"); + if (!mdok) { + sersicMinDev = NAN; + } + + sources = psArraySort (sources, pmSourceSortBySeq); for (int i = 0; i < sources->n; i++) { @@ -263,88 +284,169 @@ pmSource *source = thisSource->parent ? thisSource->parent : thisSource; - psF32 xPos, yPos, flux, rMajor, rMinor, theta; - psS32 modelType = 0; + #define MAX_ROWS_PER_SRC 10 + psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC]; + psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC]; + psS32 modelType[MAX_ROWS_PER_SRC]; + psF32 sersicIndex[MAX_ROWS_PER_SRC]; bool fitGalaxy = false; bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; - psF32 sersicIndex = 0; - // For now only perform galaxy fits on extended objects - if (source->modelEXT == NULL) { - pmModel *model = source->modelPSF; - if (model == NULL) continue; - psF32 *PAR = model->params->data.F32; - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || - !isfinite(source->psfFlux)) { - continue; - } - - xPos = model->params->data.F32[PM_PAR_XPOS]; - yPos = model->params->data.F32[PM_PAR_YPOS]; - flux = source->psfFlux; - rMajor = 0; - rMinor = 0; - theta = 0; - } else { - // Find the model with the selected type. If selected type is -1 choose the one selected ad - // modelEXT which was the best - int iModel = -1; - if (source->modelEXT) { - pmModelType ext_model_type = selectedModelType != -1 ? selectedModelType : source->modelEXT->type; + int n_rows = 0; + + // start with psf model + pmModel *model = source->modelPSF; + if (model == NULL) continue; + psF32 *PAR = model->params->data.F32; + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || + !isfinite(source->psfFlux)) { + continue; + } + + xPos[0] = model->params->data.F32[PM_PAR_XPOS]; + yPos[0] = model->params->data.F32[PM_PAR_YPOS]; + flux[0] = source->psfFlux; + rMajor[0] = 0; + rMinor[0] = 0; + theta[0] = 0; + modelType[0] = -1; + sersicIndex[0] = NAN; + + if (source->modelFits != NULL) { + // figure out which models to use based on recipe paramters + if (chooseAll) { + // Save parameters for all valid extended models + + // but make sure we aren't going to overflow our arrays + assert (source->modelFits->n < MAX_ROWS_PER_SRC); + for (int j=0; jmodelFits->n; j++) { - pmModel *aModel = source->modelFits->data[j]; - if (aModel->type == ext_model_type) { - iModel = j; - break; + pmModel *model = source->modelFits->data[j]; + psF32 *PAR = model->params->data.F32; + + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[n_rows] = PAR[PM_PAR_XPOS]; + yPos[n_rows] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor[n_rows] = axes.major; + rMinor[n_rows] = axes.minor; + theta[n_rows] = axes.theta*PS_DEG_RAD; + flux[n_rows] = pow(10.0, -0.4*model->mag); + modelType[n_rows] = model->type; + fitGalaxy = true; + if (model->params->n == 8) { + // this will save the sersic index in all models but that might + // be useful + sersicIndex[n_rows] = 0.5 / PAR[PM_PAR_7]; + } else { + sersicIndex[n_rows] = NAN; + } + + n_rows++; + } + } + } else { + int jModelSersic = -1; + int jModelDev = -1; + int jModelSelected = -1; + psF32 minChisq = NAN; + for (int j=0; jmodelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + if (chooseBest) { + // choose the model with lowest chisq + if (isfinite(model->chisq) && (!isfinite(minChisq) || model->chisq < minChisq)) { + jModelSelected = j; + minChisq = model->chisq; + } + } else { + // find the index of models of interest + if (model->type == selectedModelType) { + jModelSelected = j; + } else if (model->type == sersicModelType) { + jModelSersic = j; + } else if (model->type == devModelType) { + jModelDev = j; + } + } + } + if (jModelSelected >= 0 || jModelSersic >= 0) { + // If a specific non-sersic model we take paramers from that one. + // Otherwise we do the sersic model. + pmModel *model = jModelSelected >= 0 ? source->modelFits->data[jModelSelected] : + source->modelFits->data[jModelSersic]; + psF32 *PAR = model->params->data.F32; + + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[0] = PAR[PM_PAR_XPOS]; + yPos[0] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor[0] = axes.major; + rMinor[0] = axes.minor; + theta[0] = axes.theta*PS_DEG_RAD; + flux[0] = pow(10.0, -0.4*model->mag); + modelType[0] = model->type; + fitGalaxy = true; + if (model->type == sersicModelType) { + PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); + sersicIndex[0] = 0.5 / PAR[PM_PAR_7]; + } + + n_rows = 1; + + // Unless a specific non-sersic model type was selected do dev model for sources with + // sersic index above the recipe limit. + if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) && + sersicIndex[0] > sersicMinDev) { + + model = source->modelFits->data[jModelDev]; + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[1] = PAR[PM_PAR_XPOS]; + yPos[1] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor[1] = axes.major; + rMinor[1] = axes.minor; + theta[1] = axes.theta*PS_DEG_RAD; + flux[1] = pow(10.0, -0.4*model->mag); + modelType[1] = model->type; + sersicIndex[1] = NAN; + n_rows = 2; + } + } } } } - if (iModel == -1) { - // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits? - // XXX: Should this be an assert? - continue; - } - pmModel *model = source->modelFits->data[iModel]; - psF32 *PAR = model->params->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { - // XXX: do we need to handle this type - continue; - } else { - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || - !isfinite(model->mag)) { - // bad model - continue; - } - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); - rMajor = axes.major; - rMinor = axes.minor; - theta = axes.theta*PS_DEG_RAD; - flux = pow(10.0, -0.4*model->mag); - fitGalaxy = true; - } - modelType = model->type; - if (modelType == sersicModelType) { - PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); - sersicIndex = 0.5 / PAR[PM_PAR_7]; - } - } - psMetadata *row = psMetadataAlloc(); - psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux/exptime); - psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); - psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); - psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); - psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); - psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); - psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor); - psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor); - psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta); - psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType); - psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); - - psArrayAdd(table, 100, row); - psFree(row); + } + if (n_rows == 0) { + // didn't find a suitable extended model we have one row with the psf parameters + n_rows = 1; + } + + for (int j = 0; j < n_rows; j++) { + psMetadata *row = psMetadataAlloc(); + psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux[j]/exptime); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); + psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); + psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); + psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); + psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta[j]); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex[j]); + + psArrayAdd(table, 100, row); + psFree(row); + } }