Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 13375) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 13900) @@ -5,12 +5,12 @@ // model with selected pixels, and the fit radius must be defined -// given the set of sources, each of which points to the pixels in the -// science image, we construct a set of simulated sources with their own pixels. -// these are used to determine the simultaneous linear fit of fluxes. +// given the set of sources, each of which points to the pixels in the +// science image, we construct a set of simulated sources with their own pixels. +// these are used to determine the simultaneous linear fit of fluxes. // the analysis is performed wrt the simulated pixel values static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal) { bool status; @@ -41,9 +41,9 @@ int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); if (!status) { - SKY_FIT_ORDER = 0; + SKY_FIT_ORDER = 0; } bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR"); if (!status) { - SKY_FIT_LINEAR = false; + SKY_FIT_LINEAR = false; } @@ -55,5 +55,5 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - if (source->type == PM_SOURCE_TYPE_STAR && + if (source->type == PM_SOURCE_TYPE_STAR && source->mode & PM_SOURCE_MODE_SATSTAR) continue; if (final) { @@ -67,5 +67,5 @@ y = source->peak->yf; - // is the source in the region of interest? + // is the source in the region of interest? if (x < AnalysisRegion.x0) continue; if (y < AnalysisRegion.y0) continue; @@ -97,11 +97,11 @@ psSparseMatrixElement (sparse, i, i, f); - // the formal error depends on the weighting scheme - if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false); - errors->data.F32[i] = 1.0 / sqrt(var); - } else { - errors->data.F32[i] = 1.0 / sqrt(f); - } + // the formal error depends on the weighting scheme + if (CONSTANT_PHOTOMETRIC_WEIGHTS) { + float var = pmSourceModelDotModel (SRCi, SRCi, false); + errors->data.F32[i] = 1.0 / sqrt(var); + } else { + errors->data.F32[i] = 1.0 / sqrt(f); + } @@ -110,21 +110,21 @@ psSparseVectorElement (sparse, i, f); - // add the per-source weights (border region) - switch (SKY_FIT_ORDER) { - case 1: - f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 1, f); - f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 2, f); - - case 0: - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 0, f); - break; - - default: - psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER); - break; - } + // add the per-source weights (border region) + switch (SKY_FIT_ORDER) { + case 1: + f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 1, f); + f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 2, f); + + case 0: + f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 0, f); + break; + + default: + psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER); + break; + } // loop over all other stars following this one @@ -159,9 +159,9 @@ psVector *skyfit = NULL; if (SKY_FIT_LINEAR) { - psSparseBorderSolve (&norm, &skyfit, constraint, border, 5); - fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); + psSparseBorderSolve (&norm, &skyfit, constraint, border, 5); + fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); } else { - norm = psSparseSolve (NULL, constraint, sparse, 5); - skyfit = NULL; + norm = psSparseSolve (NULL, constraint, sparse, 5); + skyfit = NULL; } @@ -169,5 +169,5 @@ for (int i = 0; i < fitSources->n; i++) { pmSource *source = fitSources->data[i]; - pmModel *model = pmSourceGetModel (NULL, source); + pmModel *model = pmSourceGetModel (NULL, source); // assign linearly-fitted normalization @@ -177,11 +177,11 @@ model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; - // XXX is the value of 'errors' modified by the sky fit? + // XXX is the value of 'errors' modified by the sky fit? // subtract object - pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; if (!final) source->mode |= PM_SOURCE_MODE_TEMPSUB; - // XXX not sure about the use of TEMPSUB + // XXX not sure about the use of TEMPSUB } @@ -190,5 +190,5 @@ pmSource *source = fitSources->data[i]; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->weight); + pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal); } @@ -219,10 +219,10 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; + pmModel *model = pmSourceGetModel (NULL, source); + if (model == NULL) continue; float x = model->params->data.F32[PM_PAR_XPOS]; float y = model->params->data.F32[PM_PAR_YPOS]; - psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); - } + psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + } // accumulate the image statistics from the masked regions @@ -233,31 +233,31 @@ double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; w = x = y = x2 = xy = y2 = fo = fx = fy = 0; - + int col0 = readout->image->col0; int row0 = readout->image->row0; for (int j = 0; j < readout->image->numRows; j++) { - for (int i = 0; i < readout->image->numCols; i++) { - if (mask[j][i]) continue; - if (constant_weights) { - wt = 1.0; - } else { - wt = weight[j][i]; - } - f = image[j][i]; - w += 1/wt; - fo += f/wt; - if (SKY_FIT_ORDER == 0) continue; - - xc = i + col0; - yc = j + row0; - x += xc/wt; - y += yc/wt; - x2 += xc*xc/wt; - xy += xc*yc/wt; - y2 += yc*yc/wt; - fx += f*xc/wt; - fy += f*yc/wt; - } + for (int i = 0; i < readout->image->numCols; i++) { + if (mask[j][i]) continue; + if (constant_weights) { + wt = 1.0; + } else { + wt = weight[j][i]; + } + f = image[j][i]; + w += 1/wt; + fo += f/wt; + if (SKY_FIT_ORDER == 0) continue; + + xc = i + col0; + yc = j + row0; + x += xc/wt; + y += yc/wt; + x2 += xc*xc/wt; + xy += xc*yc/wt; + y2 += yc*yc/wt; + fx += f*xc/wt; + fy += f*yc/wt; + } } @@ -269,15 +269,15 @@ psSparseBorderElementT (border, 0, 0, w); if (SKY_FIT_ORDER > 0) { - psSparseBorderElementG (border, 0, fx); - psSparseBorderElementG (border, 0, fy); - psSparseBorderElementT (border, 1, 0, x); - psSparseBorderElementT (border, 2, 0, y); - psSparseBorderElementT (border, 0, 1, x); - psSparseBorderElementT (border, 1, 1, x2); - psSparseBorderElementT (border, 2, 1, xy); - psSparseBorderElementT (border, 0, 2, y); - psSparseBorderElementT (border, 1, 2, xy); - psSparseBorderElementT (border, 2, 2, y2); - } + psSparseBorderElementG (border, 0, fx); + psSparseBorderElementG (border, 0, fy); + psSparseBorderElementT (border, 1, 0, x); + psSparseBorderElementT (border, 2, 0, y); + psSparseBorderElementT (border, 0, 1, x); + psSparseBorderElementT (border, 1, 1, x2); + psSparseBorderElementT (border, 2, 1, xy); + psSparseBorderElementT (border, 0, 2, y); + psSparseBorderElementT (border, 1, 2, xy); + psSparseBorderElementT (border, 2, 2, y2); + } return true; }