Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 9594) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 9595) @@ -4,5 +4,5 @@ // 2006.02.07 : no leaks! // fit all reasonable sources with the linear PSF model -bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { +bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { bool status; @@ -39,103 +39,103 @@ const bool CONSTANT_PHOTOMETRIC_WEIGHTS = - psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); if (!status) { - psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); + psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); } for (int i = 0; i < sources->n; i++) { - pmSource *inSource = sources->data[i]; - - // skip non-astronomical objects (very likely defects) - // XXX EAM : should we try these anyway? - if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; - if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue; - if (final) { - if (inSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; - } else { - if (inSource->mode & PM_SOURCE_MODE_BLEND) continue; - } - - if (inSource->moments->x < AnalysisRegion.x0) continue; - if (inSource->moments->y < AnalysisRegion.y0) continue; - if (inSource->moments->x > AnalysisRegion.x1) continue; - if (inSource->moments->y > AnalysisRegion.y1) continue; - - pmSource *otSource = pmSourceAlloc (); - - // really saturated stars should be re-measured for a better centroid - // XXX EAM : move this to a 'clear satstar function' - if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { - status = pmSourceMoments (inSource, INNER_RADIUS); - } - - // XXX EAM : add option to use EXT or PSF form - // use the source moments, etc to guess basic model parameters - pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); - if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { - modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x; - modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y; - } else { - // peak-up on peak (for non-sat objects) - - // ix,iy must land on inSource->pixels - int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x)); - int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y)); - - psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy); - psPlane min = psImageBicubeMin (bicube); - - psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", - modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy); - - // if min point is too deviant, keep the old value - if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { - modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix; - modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy; - } - psFree (bicube); - } - - // set PSF parameters for this model - pmModel *model = pmModelFromPSF (modelEXT, psf); - psFree (modelEXT); - - // save the original coords - x = model->params->data.F32[PM_PAR_XPOS]; - y = model->params->data.F32[PM_PAR_YPOS]; - - // set the fit radius based on the object flux limit and the model - psphotCheckRadiusPSF (readout, inSource, model); - - // make temporary copies of the image pixels and mask - otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); - otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); - otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32); - - // build the model image - psImage *flux = otSource->pixels; - psImage *mask = otSource->mask; - - // set model to unit peak, zero sky (we assume sky is subtracted) - model->params->data.F32[PM_PAR_SKY] = 0.0; - model->params->data.F32[PM_PAR_I0] = 1.0; - - // fill in the model pixel values - psImageInit (flux, 0.0); - psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); - pmModelAdd (flux, mask, model, false, false); - - // calculate nDOF (nPix - 1) - // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT); - // model->nDOF = mask->numCols*mask->numRows - Nmaskpix - 1; - - // save source in list - otSource->modelPSF = model; - index->data.U32[models->n] = i; - psArrayAdd (models, 100, otSource); - psFree (otSource); + pmSource *inSource = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; + if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue; + if (final) { + if (inSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; + } else { + if (inSource->mode & PM_SOURCE_MODE_BLEND) continue; + } + + if (inSource->moments->x < AnalysisRegion.x0) continue; + if (inSource->moments->y < AnalysisRegion.y0) continue; + if (inSource->moments->x > AnalysisRegion.x1) continue; + if (inSource->moments->y > AnalysisRegion.y1) continue; + + pmSource *otSource = pmSourceAlloc (); + + // really saturated stars should be re-measured for a better centroid + // XXX EAM : move this to a 'clear satstar function' + if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { + status = pmSourceMoments (inSource, INNER_RADIUS); + } + + // XXX EAM : add option to use EXT or PSF form + // use the source moments, etc to guess basic model parameters + pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); + if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { + modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x; + modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y; + } else { + // peak-up on peak (for non-sat objects) + + // ix,iy must land on inSource->pixels + int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x)); + int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y)); + + psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy); + psPlane min = psImageBicubeMin (bicube); + + psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", + modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy); + + // if min point is too deviant, keep the old value + if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { + modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix; + modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy; + } + psFree (bicube); + } + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelEXT, psf); + psFree (modelEXT); + + // save the original coords + x = model->params->data.F32[PM_PAR_XPOS]; + y = model->params->data.F32[PM_PAR_YPOS]; + + // set the fit radius based on the object flux limit and the model + psphotCheckRadiusPSF (readout, inSource, model); + + // make temporary copies of the image pixels and mask + otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); + otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); + otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32); + + // build the model image + psImage *flux = otSource->pixels; + psImage *mask = otSource->mask; + + // set model to unit peak, zero sky (we assume sky is subtracted) + model->params->data.F32[PM_PAR_SKY] = 0.0; + model->params->data.F32[PM_PAR_I0] = 1.0; + + // fill in the model pixel values + psImageInit (flux, 0.0); + psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); + pmModelAdd (flux, mask, model, false, false); + + // calculate nDOF (nPix - 1) + // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT); + // model->nDOF = mask->numCols*mask->numRows - Nmaskpix - 1; + + // save source in list + otSource->modelPSF = model; + index->data.U32[models->n] = i; + psArrayAdd (models, 100, otSource); + psFree (otSource); } psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n); - + // fill out the sparse matrix psSparse *sparse = psSparseAlloc (models->n, 100); @@ -144,32 +144,32 @@ for (int i = 0; i < models->n; i++) { - int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmSource *Mi = models->data[i]; - - // scale by diagonal element (auto-cross-product) - r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); - weight->data.F32[i] = r; - - psSparseMatrixElement (sparse, i, i, 1.0); - - // find the image x model value - f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseVectorElement (sparse, i, f / r); - - // loop over all other stars following this one - for (int j = i + 1; j < models->n; j++) { - pmSource *Mj = models->data[j]; - - // skip over disjoint source images, break after last possible overlap - if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; - if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; - if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; - if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; - - // got an overlap; calculate cross-product and add to output array - f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseMatrixElement (sparse, j, i, f / r); - } + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; + pmSource *Mi = models->data[i]; + + // scale by diagonal element (auto-cross-product) + r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); + weight->data.F32[i] = r; + + psSparseMatrixElement (sparse, i, i, 1.0); + + // find the image x model value + f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseVectorElement (sparse, i, f / r); + + // loop over all other stars following this one + for (int j = i + 1; j < models->n; j++) { + pmSource *Mj = models->data[j]; + + // skip over disjoint source images, break after last possible overlap + if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; + if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; + if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; + if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; + + // got an overlap; calculate cross-product and add to output array + f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseMatrixElement (sparse, j, i, f / r); + } } psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); @@ -186,40 +186,40 @@ // adjust models, set sources and subtract for (int i = 0; i < models->n; i++) { - int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmSource *Mi = models->data[i]; - - // if we already have a PSF model, free it. - psFree (Fi->modelPSF); - - // need to increment counter so we can free models here and sources above - Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF); - - // assign linearly-fitted normalization - if (isnan(norm->data.F32[i])) { - psAbort ("psphot", "ensemble source is nan"); - } - Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; - Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i])); - // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it - - // subtract object - pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); - Fi->mode |= PM_SOURCE_MODE_SUBTRACTED; - if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB; + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; + pmSource *Mi = models->data[i]; + + // if we already have a PSF model, free it. + psFree (Fi->modelPSF); + + // need to increment counter so we can free models here and sources above + Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF); + + // assign linearly-fitted normalization + if (isnan(norm->data.F32[i])) { + psAbort ("psphot", "ensemble source is nan"); + } + Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; + Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i])); + // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it + + // subtract object + pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); + Fi->mode |= PM_SOURCE_MODE_SUBTRACTED; + if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB; } // measure chisq for each source for (int i = 0; final && (i < models->n); i++) { - int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmModel *model = Fi->modelPSF; - - x = model->params->data.F32[PM_PAR_XPOS]; - y = model->params->data.F32[PM_PAR_YPOS]; - - psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); - pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight); - psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; + pmModel *model = Fi->modelPSF; + + x = model->params->data.F32[PM_PAR_XPOS]; + y = model->params->data.F32[PM_PAR_YPOS]; + + psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); + pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight); + psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); } @@ -239,10 +239,10 @@ int Npix = 0; for (int j = 0; j < image->numRows; j++) { - for (int i = 0; i < image->numCols; i++) { - if (mask->data.U8[j][i]) continue; - if (weight->data.F32[j][i] <= 0) continue; - dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i]; - Npix ++; - } + for (int i = 0; i < image->numCols; i++) { + if (mask->data.U8[j][i]) continue; + if (weight->data.F32[j][i] <= 0) continue; + dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i]; + Npix ++; + } } model->nDOF = Npix - 1; Index: /trunk/psphot/src/psphotGrowthCurve.c =================================================================== --- /trunk/psphot/src/psphotGrowthCurve.c (revision 9594) +++ /trunk/psphot/src/psphotGrowthCurve.c (revision 9595) @@ -5,5 +5,5 @@ // the 'center' option -bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf) { +bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf) { // bool status; @@ -15,5 +15,5 @@ pmModel *modelRef = pmModelAlloc(psf->type); - // use center of the image + // use center of the image xc = 0.5*readout->image->numCols + readout->image->col0; yc = 0.5*readout->image->numRows + readout->image->row0; @@ -41,16 +41,16 @@ for (int i = 0; i < psf->growth->radius->n; i++) { - psImageInit (image, 0.0); + psImageInit (image, 0.0); - radius = psf->growth->radius->data.F32[i]; + radius = psf->growth->radius->data.F32[i]; - psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK); + psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK); - pmModelAdd (image, mask, model, false, false); + pmModelAdd (image, mask, model, false, false); - pmSourcePhotometryAper (&apMag, model, image, mask); + pmSourcePhotometryAper (&apMag, model, image, mask); - psImageKeepCircle (mask, xc, yc, radius, "AND", NOT_U8(PM_MASK_MARK)); - psf->growth->apMag->data.F32[i] = apMag; + psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); + psf->growth->apMag->data.F32[i] = apMag; } @@ -62,5 +62,5 @@ psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f\n", psf->growth->apLoss); - + psFree (view); psFree (image); Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 9594) +++ /trunk/psphot/src/psphotSourceFits.c (revision 9595) @@ -3,5 +3,5 @@ // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { float x, y, dR; @@ -9,6 +9,6 @@ // if this source is not a possible blend, just fit as PSF if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { - bool status = psphotFitPSF (readout, source, psf); - return status; + bool status = psphotFitPSF (readout, source, psf); + return status; } psTrace ("psphot", 5, "trying blend...\n"); @@ -23,5 +23,5 @@ psArray *modelSet = psArrayAlloc (source->blends->n + 1); psArrayAdd (modelSet, 16, PSF); - + psArray *sourceSet = psArrayAlloc (source->blends->n + 1); psArrayAdd (sourceSet, 16, source); @@ -30,28 +30,28 @@ dR = 0; for (int i = 0; i < source->blends->n; i++) { - pmSource *blend = source->blends->data[i]; - - // find the blend which is furthest from source - dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y)); - - // create the model and guess parameters for this blend - pmModel *model = pmModelAlloc (PSF->type); - for (int j = 0; j < model->params->n; j++) { - model->params->data.F32[j] = PSF->params->data.F32[j]; - model->dparams->data.F32[j] = PSF->dparams->data.F32[j]; - } - - // XXX assume local sky is 0.0? - model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky; - if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit"); - model->params->data.F32[2] = blend->peak->x; - model->params->data.F32[3] = blend->peak->y; - - // add this blend to the list - psArrayAdd (modelSet, 16, model); - psArrayAdd (sourceSet, 16, blend); - - // free to avoid double counting model - psFree (model); + pmSource *blend = source->blends->data[i]; + + // find the blend which is furthest from source + dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y)); + + // create the model and guess parameters for this blend + pmModel *model = pmModelAlloc (PSF->type); + for (int j = 0; j < model->params->n; j++) { + model->params->data.F32[j] = PSF->params->data.F32[j]; + model->dparams->data.F32[j] = PSF->dparams->data.F32[j]; + } + + // XXX assume local sky is 0.0? + model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky; + if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit"); + model->params->data.F32[2] = blend->peak->x; + model->params->data.F32[3] = blend->peak->y; + + // add this blend to the list + psArrayAdd (modelSet, 16, model); + psArrayAdd (sourceSet, 16, blend); + + // free to avoid double counting model + psFree (model); } @@ -62,5 +62,5 @@ psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); - psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); + psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); // correct model chisq for flux trend @@ -70,24 +70,24 @@ // evaluate the blend objects, subtract if good, free otherwise for (int i = 1; i < modelSet->n; i++) { - pmSource *blend = sourceSet->data[i]; - pmModel *model = modelSet->data[i]; - - // correct model chisq for flux trend - chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]); - model->chisqNorm = model->chisq / chiTrend; - - // if this one failed, skip it - if (!psphotEvalPSF (blend, model)) { - psTrace ("psphot", 5, "failed on blend %d of %ld\n", i, modelSet->n); - continue; - } - - // otherwise, supply the resulting model to the corresponding blend - psFree(blend->modelPSF); - blend->modelPSF = psMemIncrRefCounter (model); - psTrace ("psphot", 5, "fitted blend as PSF\n"); - pmModelSub (source->pixels, source->mask, model, false, false); - blend->mode |= PM_SOURCE_MODE_SUBTRACTED; - blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; + pmSource *blend = sourceSet->data[i]; + pmModel *model = modelSet->data[i]; + + // correct model chisq for flux trend + chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]); + model->chisqNorm = model->chisq / chiTrend; + + // if this one failed, skip it + if (!psphotEvalPSF (blend, model)) { + psTrace ("psphot", 5, "failed on blend %d of %ld\n", i, modelSet->n); + continue; + } + + // otherwise, supply the resulting model to the corresponding blend + psFree(blend->modelPSF); + blend->modelPSF = psMemIncrRefCounter (model); + psTrace ("psphot", 5, "fitted blend as PSF\n"); + pmModelSub (source->pixels, source->mask, model, false, false); + blend->mode |= PM_SOURCE_MODE_SUBTRACTED; + blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; } psFree (modelSet); @@ -96,7 +96,7 @@ // evaluate the primary object if (!psphotEvalPSF (source, PSF)) { - psTrace ("psphot", 5, "failed on blend 0 of %ld\n", modelSet->n); - psFree (PSF); - return false; + psTrace ("psphot", 5, "failed on blend 0 of %ld\n", modelSet->n); + psFree (PSF); + return false; } @@ -110,5 +110,5 @@ } -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { float x, y; @@ -128,6 +128,6 @@ psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF); - psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); - + psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); + // correct model chisq for flux trend chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]); @@ -136,6 +136,6 @@ // does the PSF model succeed? if (!psphotEvalPSF (source, PSF)) { - psFree (PSF); - return false; + psFree (PSF); + return false; } @@ -143,5 +143,5 @@ pmModelSub (source->pixels, source->mask, PSF, false, false); - // free old model, save new model + // free old model, save new model psFree (source->modelPSF); source->modelPSF = PSF; @@ -173,5 +173,5 @@ } -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { bool okEXT, okDBL; @@ -216,11 +216,11 @@ ONE->chisqNorm = ONE->chisq / chiTrend; - psFree (tmpSrc); + psFree (tmpSrc); if (okEXT && okDBL) { - psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL); - // XXX EAM : a bogus bias: need to examine this better - if (3*chiEXT > chiDBL) goto keepDBL; - goto keepEXT; + psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL); + // XXX EAM : a bogus bias: need to examine this better + if (3*chiEXT > chiDBL) goto keepDBL; + goto keepEXT; } @@ -276,5 +276,5 @@ // fit a double PSF source to an extended blob -psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { +psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { float x, y, dx, dy; @@ -284,5 +284,5 @@ psEllipseMoments moments; psArray *modelSet; - + // make a guess at the position of the two sources moments.x2 = source->moments->Sx; @@ -320,15 +320,15 @@ psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); - psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); + psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); return (modelSet); } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { float x, y; // use the source moments, etc to guess basic model parameters - pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); + pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); // if (isnan(EXT->params->data.F32[1])) psAbort ("psphot", "nan in ext fit"); @@ -339,5 +339,5 @@ if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) { - psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy); + psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy); } @@ -345,5 +345,5 @@ psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "OR", PM_MASK_MARK); pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT); - psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); + psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); return (EXT); Index: /trunk/psphot/src/psphotWeightBias.c =================================================================== --- /trunk/psphot/src/psphotWeightBias.c (revision 9594) +++ /trunk/psphot/src/psphotWeightBias.c (revision 9595) @@ -3,7 +3,7 @@ // select objects fitted with PSF model // re-fit all of them with the non-poisson errors -// only allow the normalization to vary +// only allow the normalization to vary // XXX eventually, we should be able to do this with linear fitting... -bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { +bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { int Nfit = 0; @@ -16,8 +16,8 @@ // source analysis is done in S/N order (brightest first) sources = psArraySort (sources, psphotSortBySN); - + // set fitting method to use non-poisson errors (local sky error) pmSourceFitModelInit (15, 0.1, false); - + // option to limit analysis to a specific region char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); @@ -29,46 +29,46 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + pmSource *source = sources->data[i]; - // skip lower-quality objects - if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (source->mode & PM_SOURCE_MODE_POOR) continue; - if (source->mode & PM_SOURCE_MODE_FAIL) continue; - if (source->mode & PM_SOURCE_MODE_PAIR) continue; - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + // skip lower-quality objects + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (source->mode & PM_SOURCE_MODE_POOR) continue; + if (source->mode & PM_SOURCE_MODE_FAIL) continue; + if (source->mode & PM_SOURCE_MODE_PAIR) continue; + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; - // if model is NULL, we don't have a starting guess - if (source->modelPSF == NULL) continue; + // if model is NULL, we don't have a starting guess + if (source->modelPSF == NULL) continue; - if (source->moments->x < AnalysisRegion.x0) continue; - if (source->moments->y < AnalysisRegion.y0) continue; - if (source->moments->x > AnalysisRegion.x1) continue; - if (source->moments->y > AnalysisRegion.y1) continue; + if (source->moments->x < AnalysisRegion.x0) continue; + if (source->moments->y < AnalysisRegion.y0) continue; + if (source->moments->x > AnalysisRegion.x1) continue; + if (source->moments->y > AnalysisRegion.y1) continue; - // replace object in image - pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + // replace object in image + pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); - // make a temporary model (we don't keep the result of this analysis) - pmModel *PSF = pmModelCopy (source->modelPSF); + // make a temporary model (we don't keep the result of this analysis) + pmModel *PSF = pmModelCopy (source->modelPSF); - // extend source radius as needed - psphotCheckRadiusPSF (readout, source, PSF); + // extend source radius as needed + psphotCheckRadiusPSF (readout, source, PSF); - x = PSF->params->data.F32[2]; - y = PSF->params->data.F32[3]; + x = PSF->params->data.F32[2]; + y = PSF->params->data.F32[3]; - // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); - pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM); - psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); - - // re-subtract PSF for object, leave local sky - pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); + pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM); + psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); - PARp = source->modelPSF->params->data.F32; - PARc = PSF->params->data.F32; - fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky); - Nfit ++; + // re-subtract PSF for object, leave local sky + pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + + PARp = source->modelPSF->params->data.F32; + PARc = PSF->params->data.F32; + fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky); + Nfit ++; } fclose (f);