Changeset 9595
- Timestamp:
- Oct 16, 2006, 4:24:03 PM (20 years ago)
- Location:
- trunk/psphot/src
- Files:
-
- 4 edited
-
psphotEnsemblePSF.c (modified) (5 diffs)
-
psphotGrowthCurve.c (modified) (4 diffs)
-
psphotSourceFits.c (modified) (18 diffs)
-
psphotWeightBias.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/psphotEnsemblePSF.c
r9576 r9595 4 4 // 2006.02.07 : no leaks! 5 5 // fit all reasonable sources with the linear PSF model 6 bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { 6 bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { 7 7 8 8 bool status; … … 39 39 40 40 const bool CONSTANT_PHOTOMETRIC_WEIGHTS = 41 41 psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); 42 42 if (!status) { 43 43 psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); 44 44 } 45 45 46 46 for (int i = 0; i < sources->n; i++) { 47 48 49 50 51 if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 75 76 77 78 79 80 81 // ix,iy must land on inSource->pixels 82 83 84 85 86 87 88 psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 // build the model image 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 47 pmSource *inSource = sources->data[i]; 48 49 // skip non-astronomical objects (very likely defects) 50 // XXX EAM : should we try these anyway? 51 if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; 52 if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue; 53 if (final) { 54 if (inSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; 55 } else { 56 if (inSource->mode & PM_SOURCE_MODE_BLEND) continue; 57 } 58 59 if (inSource->moments->x < AnalysisRegion.x0) continue; 60 if (inSource->moments->y < AnalysisRegion.y0) continue; 61 if (inSource->moments->x > AnalysisRegion.x1) continue; 62 if (inSource->moments->y > AnalysisRegion.y1) continue; 63 64 pmSource *otSource = pmSourceAlloc (); 65 66 // really saturated stars should be re-measured for a better centroid 67 // XXX EAM : move this to a 'clear satstar function' 68 if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { 69 status = pmSourceMoments (inSource, INNER_RADIUS); 70 } 71 72 // XXX EAM : add option to use EXT or PSF form 73 // use the source moments, etc to guess basic model parameters 74 pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 75 if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { 76 modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x; 77 modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y; 78 } else { 79 // peak-up on peak (for non-sat objects) 80 81 // ix,iy must land on inSource->pixels 82 int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x)); 83 int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y)); 84 85 psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy); 86 psPlane min = psImageBicubeMin (bicube); 87 88 psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", 89 modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy); 90 91 // if min point is too deviant, keep the old value 92 if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { 93 modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix; 94 modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy; 95 } 96 psFree (bicube); 97 } 98 99 // set PSF parameters for this model 100 pmModel *model = pmModelFromPSF (modelEXT, psf); 101 psFree (modelEXT); 102 103 // save the original coords 104 x = model->params->data.F32[PM_PAR_XPOS]; 105 y = model->params->data.F32[PM_PAR_YPOS]; 106 107 // set the fit radius based on the object flux limit and the model 108 psphotCheckRadiusPSF (readout, inSource, model); 109 110 // make temporary copies of the image pixels and mask 111 otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); 112 otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); 113 otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32); 114 115 // build the model image 116 psImage *flux = otSource->pixels; 117 psImage *mask = otSource->mask; 118 119 // set model to unit peak, zero sky (we assume sky is subtracted) 120 model->params->data.F32[PM_PAR_SKY] = 0.0; 121 model->params->data.F32[PM_PAR_I0] = 1.0; 122 123 // fill in the model pixel values 124 psImageInit (flux, 0.0); 125 psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); 126 pmModelAdd (flux, mask, model, false, false); 127 128 // calculate nDOF (nPix - 1) 129 // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT); 130 // model->nDOF = mask->numCols*mask->numRows - Nmaskpix - 1; 131 132 // save source in list 133 otSource->modelPSF = model; 134 index->data.U32[models->n] = i; 135 psArrayAdd (models, 100, otSource); 136 psFree (otSource); 137 137 } 138 138 psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n); 139 139 140 140 // fill out the sparse matrix 141 141 psSparse *sparse = psSparseAlloc (models->n, 100); … … 144 144 145 145 for (int i = 0; i < models->n; i++) { 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 psSparseMatrixElement (sparse, j, i, f / r); 173 146 int N = index->data.U32[i]; 147 pmSource *Fi = sources->data[N]; 148 pmSource *Mi = models->data[i]; 149 150 // scale by diagonal element (auto-cross-product) 151 r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); 152 weight->data.F32[i] = r; 153 154 psSparseMatrixElement (sparse, i, i, 1.0); 155 156 // find the image x model value 157 f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); 158 psSparseVectorElement (sparse, i, f / r); 159 160 // loop over all other stars following this one 161 for (int j = i + 1; j < models->n; j++) { 162 pmSource *Mj = models->data[j]; 163 164 // skip over disjoint source images, break after last possible overlap 165 if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; 166 if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; 167 if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; 168 if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; 169 170 // got an overlap; calculate cross-product and add to output array 171 f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS); 172 psSparseMatrixElement (sparse, j, i, f / r); 173 } 174 174 } 175 175 psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); … … 186 186 // adjust models, set sources and subtract 187 187 for (int i = 0; i < models->n; i++) { 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 188 int N = index->data.U32[i]; 189 pmSource *Fi = sources->data[N]; 190 pmSource *Mi = models->data[i]; 191 192 // if we already have a PSF model, free it. 193 psFree (Fi->modelPSF); 194 195 // need to increment counter so we can free models here and sources above 196 Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF); 197 198 // assign linearly-fitted normalization 199 if (isnan(norm->data.F32[i])) { 200 psAbort ("psphot", "ensemble source is nan"); 201 } 202 Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; 203 Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i])); 204 // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it 205 206 // subtract object 207 pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); 208 Fi->mode |= PM_SOURCE_MODE_SUBTRACTED; 209 if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB; 210 210 } 211 211 212 212 // measure chisq for each source 213 213 for (int i = 0; final && (i < models->n); i++) { 214 215 216 217 218 219 220 221 222 223 psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND",NOT_U8(PM_MASK_MARK));214 int N = index->data.U32[i]; 215 pmSource *Fi = sources->data[N]; 216 pmModel *model = Fi->modelPSF; 217 218 x = model->params->data.F32[PM_PAR_XPOS]; 219 y = model->params->data.F32[PM_PAR_YPOS]; 220 221 psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); 222 pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight); 223 psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 224 224 } 225 225 … … 239 239 int Npix = 0; 240 240 for (int j = 0; j < image->numRows; j++) { 241 242 243 244 245 246 241 for (int i = 0; i < image->numCols; i++) { 242 if (mask->data.U8[j][i]) continue; 243 if (weight->data.F32[j][i] <= 0) continue; 244 dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i]; 245 Npix ++; 246 } 247 247 } 248 248 model->nDOF = Npix - 1; -
trunk/psphot/src/psphotGrowthCurve.c
r9529 r9595 5 5 // the 'center' option 6 6 7 bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf) { 7 bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf) { 8 8 9 9 // bool status; … … 15 15 pmModel *modelRef = pmModelAlloc(psf->type); 16 16 17 // use center of the image 17 // use center of the image 18 18 xc = 0.5*readout->image->numCols + readout->image->col0; 19 19 yc = 0.5*readout->image->numRows + readout->image->row0; … … 41 41 for (int i = 0; i < psf->growth->radius->n; i++) { 42 42 43 43 psImageInit (image, 0.0); 44 44 45 45 radius = psf->growth->radius->data.F32[i]; 46 46 47 47 psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK); 48 48 49 49 pmModelAdd (image, mask, model, false, false); 50 50 51 51 pmSourcePhotometryAper (&apMag, model, image, mask); 52 52 53 psImageKeepCircle (mask, xc, yc, radius, "AND",NOT_U8(PM_MASK_MARK));54 53 psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); 54 psf->growth->apMag->data.F32[i] = apMag; 55 55 } 56 56 … … 62 62 63 63 psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f\n", psf->growth->apLoss); 64 64 65 65 psFree (view); 66 66 psFree (image); -
trunk/psphot/src/psphotSourceFits.c
r9576 r9595 3 3 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful 4 4 5 bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 5 bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 6 6 7 7 float x, y, dR; … … 9 9 // if this source is not a possible blend, just fit as PSF 10 10 if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { 11 12 11 bool status = psphotFitPSF (readout, source, psf); 12 return status; 13 13 } 14 14 psTrace ("psphot", 5, "trying blend...\n"); … … 23 23 psArray *modelSet = psArrayAlloc (source->blends->n + 1); 24 24 psArrayAdd (modelSet, 16, PSF); 25 25 26 26 psArray *sourceSet = psArrayAlloc (source->blends->n + 1); 27 27 psArrayAdd (sourceSet, 16, source); … … 30 30 dR = 0; 31 31 for (int i = 0; i < source->blends->n; i++) { 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 // free to avoid double counting model 55 32 pmSource *blend = source->blends->data[i]; 33 34 // find the blend which is furthest from source 35 dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y)); 36 37 // create the model and guess parameters for this blend 38 pmModel *model = pmModelAlloc (PSF->type); 39 for (int j = 0; j < model->params->n; j++) { 40 model->params->data.F32[j] = PSF->params->data.F32[j]; 41 model->dparams->data.F32[j] = PSF->dparams->data.F32[j]; 42 } 43 44 // XXX assume local sky is 0.0? 45 model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky; 46 if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit"); 47 model->params->data.F32[2] = blend->peak->x; 48 model->params->data.F32[3] = blend->peak->y; 49 50 // add this blend to the list 51 psArrayAdd (modelSet, 16, model); 52 psArrayAdd (sourceSet, 16, blend); 53 54 // free to avoid double counting model 55 psFree (model); 56 56 } 57 57 … … 62 62 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); 63 63 pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); 64 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));64 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 65 65 66 66 // correct model chisq for flux trend … … 70 70 // evaluate the blend objects, subtract if good, free otherwise 71 71 for (int i = 1; i < modelSet->n; i++) { 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 72 pmSource *blend = sourceSet->data[i]; 73 pmModel *model = modelSet->data[i]; 74 75 // correct model chisq for flux trend 76 chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]); 77 model->chisqNorm = model->chisq / chiTrend; 78 79 // if this one failed, skip it 80 if (!psphotEvalPSF (blend, model)) { 81 psTrace ("psphot", 5, "failed on blend %d of %ld\n", i, modelSet->n); 82 continue; 83 } 84 85 // otherwise, supply the resulting model to the corresponding blend 86 psFree(blend->modelPSF); 87 blend->modelPSF = psMemIncrRefCounter (model); 88 psTrace ("psphot", 5, "fitted blend as PSF\n"); 89 pmModelSub (source->pixels, source->mask, model, false, false); 90 blend->mode |= PM_SOURCE_MODE_SUBTRACTED; 91 blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; 92 92 } 93 93 psFree (modelSet); … … 96 96 // evaluate the primary object 97 97 if (!psphotEvalPSF (source, PSF)) { 98 99 100 98 psTrace ("psphot", 5, "failed on blend 0 of %ld\n", modelSet->n); 99 psFree (PSF); 100 return false; 101 101 } 102 102 … … 110 110 } 111 111 112 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { 112 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { 113 113 114 114 float x, y; … … 128 128 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); 129 129 pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF); 130 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));131 130 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 131 132 132 // correct model chisq for flux trend 133 133 chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]); … … 136 136 // does the PSF model succeed? 137 137 if (!psphotEvalPSF (source, PSF)) { 138 139 138 psFree (PSF); 139 return false; 140 140 } 141 141 … … 143 143 pmModelSub (source->pixels, source->mask, PSF, false, false); 144 144 145 // free old model, save new model 145 // free old model, save new model 146 146 psFree (source->modelPSF); 147 147 source->modelPSF = PSF; … … 173 173 } 174 174 175 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { 175 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { 176 176 177 177 bool okEXT, okDBL; … … 216 216 ONE->chisqNorm = ONE->chisq / chiTrend; 217 217 218 psFree (tmpSrc); 218 psFree (tmpSrc); 219 219 220 220 if (okEXT && okDBL) { 221 222 // XXX EAM : a bogus bias: need to examine this better 223 224 221 psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL); 222 // XXX EAM : a bogus bias: need to examine this better 223 if (3*chiEXT > chiDBL) goto keepDBL; 224 goto keepEXT; 225 225 } 226 226 … … 276 276 277 277 // fit a double PSF source to an extended blob 278 psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { 278 psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { 279 279 280 280 float x, y, dx, dy; … … 284 284 psEllipseMoments moments; 285 285 psArray *modelSet; 286 286 287 287 // make a guess at the position of the two sources 288 288 moments.x2 = source->moments->Sx; … … 320 320 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); 321 321 pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); 322 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));322 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 323 323 324 324 return (modelSet); 325 325 } 326 326 327 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { 327 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { 328 328 329 329 float x, y; 330 330 331 331 // use the source moments, etc to guess basic model parameters 332 pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); 332 pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); 333 333 // if (isnan(EXT->params->data.F32[1])) psAbort ("psphot", "nan in ext fit"); 334 334 … … 339 339 340 340 if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) { 341 341 psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy); 342 342 } 343 343 … … 345 345 psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "OR", PM_MASK_MARK); 346 346 pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT); 347 psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));347 psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 348 348 349 349 return (EXT); -
trunk/psphot/src/psphotWeightBias.c
r9270 r9595 3 3 // select objects fitted with PSF model 4 4 // re-fit all of them with the non-poisson errors 5 // only allow the normalization to vary 5 // only allow the normalization to vary 6 6 // XXX eventually, we should be able to do this with linear fitting... 7 bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 7 bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 8 8 9 9 int Nfit = 0; … … 16 16 // source analysis is done in S/N order (brightest first) 17 17 sources = psArraySort (sources, psphotSortBySN); 18 18 19 19 // set fitting method to use non-poisson errors (local sky error) 20 20 pmSourceFitModelInit (15, 0.1, false); 21 21 22 22 // option to limit analysis to a specific region 23 23 char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); … … 29 29 30 30 for (int i = 0; i < sources->n; i++) { 31 31 pmSource *source = sources->data[i]; 32 32 33 34 if (source->type != PM_SOURCE_TYPE_STAR) continue; 35 36 37 38 39 33 // skip lower-quality objects 34 if (source->type != PM_SOURCE_TYPE_STAR) continue; 35 if (source->mode & PM_SOURCE_MODE_POOR) continue; 36 if (source->mode & PM_SOURCE_MODE_FAIL) continue; 37 if (source->mode & PM_SOURCE_MODE_PAIR) continue; 38 if (source->mode & PM_SOURCE_MODE_BLEND) continue; 39 if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; 40 40 41 42 41 // if model is NULL, we don't have a starting guess 42 if (source->modelPSF == NULL) continue; 43 43 44 45 46 47 44 if (source->moments->x < AnalysisRegion.x0) continue; 45 if (source->moments->y < AnalysisRegion.y0) continue; 46 if (source->moments->x > AnalysisRegion.x1) continue; 47 if (source->moments->y > AnalysisRegion.y1) continue; 48 48 49 50 49 // replace object in image 50 pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); 51 51 52 53 52 // make a temporary model (we don't keep the result of this analysis) 53 pmModel *PSF = pmModelCopy (source->modelPSF); 54 54 55 56 55 // extend source radius as needed 56 psphotCheckRadiusPSF (readout, source, PSF); 57 57 58 59 58 x = PSF->params->data.F32[2]; 59 y = PSF->params->data.F32[3]; 60 60 61 // fit PSF model (set/unset the pixel mask) 62 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); 63 pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM); 64 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK)); 65 66 // re-subtract PSF for object, leave local sky 67 pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); 61 // fit PSF model (set/unset the pixel mask) 62 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK); 63 pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM); 64 psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 68 65 69 PARp = source->modelPSF->params->data.F32; 70 PARc = PSF->params->data.F32; 71 fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky); 72 Nfit ++; 66 // re-subtract PSF for object, leave local sky 67 pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); 68 69 PARp = source->modelPSF->params->data.F32; 70 PARc = PSF->params->data.F32; 71 fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky); 72 Nfit ++; 73 73 } 74 74 fclose (f);
Note:
See TracChangeset
for help on using the changeset viewer.
