Index: /branches/eam_branches/20100225/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- /branches/eam_branches/20100225/psphot/src/psphotFitSourcesLinearStack.c (revision 27385) +++ /branches/eam_branches/20100225/psphot/src/psphotFitSourcesLinearStack.c (revision 27386) @@ -9,6 +9,4 @@ // 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, psImageMaskType markVal); bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { @@ -194,58 +192,2 @@ return true; } - -// Calculate the weight terms for the sky fit component of the matrix. This function operates -// on the pixels which correspond to all of the sources of interest. These elements fill in -// the border matrix components in the sparse matrix equation. -static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal) { - - // generate the image-wide weight terms - // turn on MARK for all image pixels - psRegion fullArray = psRegionSet (0, 0, 0, 0); - fullArray = psRegionForImage (readout->mask, fullArray); - psImageMaskRegion (readout->mask, fullArray, "OR", markVal); - - // turn off MARK for all object pixels - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - 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_IMAGE_MASK(markVal)); - } - - // accumulate the image statistics from the masked regions - psF32 **image = readout->image->data.F32; - psF32 **variance = readout->variance->data.F32; - psImageMaskType **mask = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; - - 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 = variance[j][i]; - } - f = image[j][i]; - w += 1/wt; - fo += f/wt; - } - } - - // turn off MARK for all image pixels - psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // set the Border T elements - psSparseBorderElementG (border, 0, fo); - psSparseBorderElementT (border, 0, 0, w); - - return true; -} Index: /branches/eam_branches/20100225/psphot/src/psphotSourceMerge.c =================================================================== --- /branches/eam_branches/20100225/psphot/src/psphotSourceMerge.c (revision 27386) +++ /branches/eam_branches/20100225/psphot/src/psphotSourceMerge.c (revision 27386) @@ -0,0 +1,49 @@ + +// we need a couple of functions to distinguish coincident sources: +// XXX these are similar (identical?) to the goals of pmSourceMatch.c + +# define NEXT1 { i++; continue; } +# define NEXT2 { j++; continue; } + +psphotSourceMerge () { + + float dx, dy; + + // sort the source list by X + pmSourceSortByX (sources1); + pmSourceSortByX (sources2); + + int i, j; + for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) { + + pmSource *src1 = sources1->data[i]; + pmSource *src2 = sources2->data[j]; + + if (!src1) NEXT1; + if (!src1->peak) NEXT1; + if (!finite(src1->peak->xf)) NEXT1; + if (!finite(src1->peak->yf)) NEXT1; + + if (!src2) NEXT2; + if (!src2->peak) NEXT2; + if (!finite(src2->peak->xf)) NEXT2; + if (!finite(src2->peak->yf)) NEXT2; + + dx = src1->peak->xf - src2->peak->xf; + if (dx < -1.02*RADIUS) NEXT1; + if (dx > +1.02*RADIUS) NEXT2; + + // we are within match range, look for matches: + for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { + + dx = src1->peak->xf - src2->peak->xf; + dy = src1->peak->yf - src2->peak->yf; + + dr = dx*dx + dy*dy; + if (dr > RADIUS2) continue; + + // add to group? + } + i++; + } +} Index: /branches/eam_branches/20100225/psphot/src/psphotStackChisqImage.c =================================================================== --- /branches/eam_branches/20100225/psphot/src/psphotStackChisqImage.c (revision 27386) +++ /branches/eam_branches/20100225/psphot/src/psphotStackChisqImage.c (revision 27386) @@ -0,0 +1,85 @@ +# include "psphotInternal.h" + +// Chi = sqrt(\sum_i im_i^2/var_i) + +static int npass = 0; + +// XXX supply filename or keep PSPHOT.INPUT fixed? +bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view) +{ + bool status = false; + + psTimerStart ("psphot.chisq.image"); + + // create a holder for the image + pmReadout *chiImage = pmReadoutAlloc(); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotStackChisqImageAddReadout(config, view, chiImage, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); + return false; + } + } + + // need to save the resulting image somewhere (PSPHOT.INPUT?) + + psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image")); + + return true; +} + +bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration + pmFPAview *view, + pmReadout *chiReadout, + char *filename, + int index) +{ + bool status = true; + + // find the currently selected readout + pmFPAfile *input = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (input, "missing file?"); + + pmReadout *inReadout = pmFPAviewThisReadout(view, input->fpa); + psAssert (inReadout, "missing readout?"); + + psImage *inImage = inReadout->image; + psImage *inVariance = inReadout->variance; + psImage *inMask = inReadout->mask; + + psImage *chiImage = chiReadout->image; + psImage *chiVariance = chiReadout->variance; + psImage *chiMask = chiReadout->mask; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // generate the chisq image: + for (int iy = 0; iy < inImage->numRows; iy++) { + for (int ix = 0; ix < inImage->numCols; ix++) { + if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + chiImage->data.F32[iy][ix] += PS_SQR(inImage->data.F32[iy][ix]) / inVariance->data.F32[iy][ix]; + chiVariance->data.F32[iy][ix] = 1.0; // ?? what is the right value? just init to this? + chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x00; // we have valid data so unmask this pixel + } + } + + if (psTraceGetLevel("psphot") > 5) { + char name[256]; + sprintf (name, "chiImage.%02d.fits", npass); + psphotSaveImage (NULL, chiImage, name); + } + npass ++; + + return true; +} +