Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 21359) +++ trunk/psphot/src/psphot.h (revision 21366) @@ -90,5 +90,5 @@ void psphotModelClassInit (void); bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal); -bool psphotSetMaskAndWeight (pmConfig *config, pmReadout *readout, psMetadata *recipe); +bool psphotSetMaskAndVariance (pmConfig *config, pmReadout *readout, psMetadata *recipe); void psphotSourceFreePixels (psArray *sources); Index: trunk/psphot/src/psphotAnnuli.c =================================================================== --- trunk/psphot/src/psphotAnnuli.c (revision 21359) +++ trunk/psphot/src/psphotAnnuli.c (revision 21366) @@ -11,5 +11,5 @@ psVector *radius = source->extpars->profile->radius; - psVector *weight = source->extpars->profile->weight; + psVector *variance = source->extpars->profile->variance; psVector *flux = source->extpars->profile->flux; @@ -22,14 +22,14 @@ psVector *fluxValues = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); psVector *fluxSquare = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); - psVector *fluxWeight = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); + psVector *fluxVariance = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); psVector *pixelCount = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); psVectorInit (fluxValues, 0.0); psVectorInit (fluxSquare, 0.0); - psVectorInit (fluxWeight, 0.0); + psVectorInit (fluxVariance, 0.0); psVectorInit (pixelCount, 0.0); // XXX this code assumes the radii are in pixels. convert from arcsec with plate scale // XXX assume the annulii above are not overlapping? much faster... - // XXX this might be must faster in the reverse order: loop over annulii and use disection to + // XXX this might be must faster in the reverse order: loop over annulii and use disection to // skip to the start of the annulus. for (int i = 0; i < flux->n; i++) { @@ -39,5 +39,5 @@ fluxValues->data.F32[j] += flux->data.F32[i]; fluxSquare->data.F32[j] += PS_SQR(flux->data.F32[i]); - fluxWeight->data.F32[j] += weight->data.F32[i]; + fluxVariance->data.F32[j] += variance->data.F32[i]; pixelCount->data.F32[j] += 1.0; } @@ -49,8 +49,8 @@ fluxSquare->data.F32[j] -= PS_SQR(fluxValues->data.F32[j]); } - + source->extpars->annuli = pmSourceAnnuliAlloc (); source->extpars->annuli->flux = fluxValues; - source->extpars->annuli->fluxErr = fluxWeight; + source->extpars->annuli->fluxErr = fluxVariance; source->extpars->annuli->fluxVar = fluxSquare; Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 21359) +++ trunk/psphot/src/psphotBlendFit.c (revision 21366) @@ -21,5 +21,5 @@ int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads if (!status) { - nThreads = 0; + nThreads = 0; } @@ -44,5 +44,5 @@ psphotInitRadiusPSF (recipe, psf->type); - // starts the timer, sets up the array of fitSets + // starts the timer, sets up the array of fitSets psphotFitInit (nThreads); @@ -56,91 +56,89 @@ psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); - fprintf (stderr, "starting with %ld sources\n", sources->n); - for (int i = 0; i < cellGroups->n; i++) { - psArray *cells = cellGroups->data[i]; - - for (int j = 0; j < cells->n; j++) { - - if (nThreads) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT"); - psArray *newSources = psArrayAllocEmpty(16); - - psArrayAdd(job->args, 1, readout); - psArrayAdd(job->args, 1, recipe); - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - psArrayAdd(job->args, 1, newSources); // return for new sources - psFree (newSources); - - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return NULL; - } - psFree(job); - } else { - int nfit = 0; - int npsf = 0; - int next = 0; - int nfail = 0; - psArray *newSources = psArrayAllocEmpty(16); - - if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - Nfit += nfit; - Npsf += npsf; - Next += next; - Nfail += nfail; - - // add these back onto sources - for (int k = 0; k < newSources->n; k++) { - psArrayAdd (sources, 16, newSources->data[k]); - } - psFree (newSources); - } - } - - if (nThreads) { - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[5]; - Nfit += scalar->data.S32; - scalar = job->args->data[6]; - Npsf += scalar->data.S32; - scalar = job->args->data[7]; - Next += scalar->data.S32; - scalar = job->args->data[8]; - Nfail += scalar->data.S32; - - // add these back onto sources - psArray *newSources = job->args->data[4]; - for (int j = 0; j < newSources->n; j++) { - psArrayAdd (sources, 16, newSources->data[j]); - } - } - psFree(job); - } - } + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + if (nThreads) { + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT"); + psArray *newSources = psArrayAllocEmpty(16); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, recipe); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + psArrayAdd(job->args, 1, newSources); // return for new sources + psFree (newSources); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return NULL; + } + psFree(job); + } else { + int nfit = 0; + int npsf = 0; + int next = 0; + int nfail = 0; + psArray *newSources = psArrayAllocEmpty(16); + + if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + Nfit += nfit; + Npsf += npsf; + Next += next; + Nfail += nfail; + + // add these back onto sources + for (int k = 0; k < newSources->n; k++) { + psArrayAdd (sources, 16, newSources->data[k]); + } + psFree (newSources); + } + } + + if (nThreads) { + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[5]; + Nfit += scalar->data.S32; + scalar = job->args->data[6]; + Npsf += scalar->data.S32; + scalar = job->args->data[7]; + Next += scalar->data.S32; + scalar = job->args->data[8]; + Nfail += scalar->data.S32; + + // add these back onto sources + psArray *newSources = job->args->data[4]; + for (int j = 0; j < newSources->n; j++) { + psArrayAdd (sources, 16, newSources->data[j]); + } + } + psFree(job); + } + } } psFree (cellGroups); @@ -207,5 +205,5 @@ if (source->peak->SN < FIT_SN_LIM) continue; - // exclude sources outside optional analysis region + // exclude sources outside optional analysis region if (source->peak->xf < AnalysisRegion.x0) continue; if (source->peak->yf < AnalysisRegion.y0) continue; @@ -216,6 +214,6 @@ if (source->modelPSF == NULL) continue; - // skip sources which are insignificant flux? - // XXX this is somewhat ad-hoc + // skip sources which are insignificant flux? + // XXX this is somewhat ad-hoc if (source->modelPSF->params->data.F32[1] < 0.1) { psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", @@ -234,19 +232,19 @@ // try fitting PSFs or extended sources depending on source->mode // these functions subtract the resulting fitted source - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_EXTENDED; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); - Next ++; - continue; - } - } else { - if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_STAR; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); - Npsf ++; - continue; - } - } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); + Next ++; + continue; + } + } else { + if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { + source->type = PM_SOURCE_TYPE_STAR; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); + Npsf ++; + continue; + } + } psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); @@ -317,5 +315,5 @@ if (source->peak->SN < FIT_SN_LIM) continue; - // exclude sources outside optional analysis region + // exclude sources outside optional analysis region if (source->peak->xf < AnalysisRegion.x0) continue; if (source->peak->yf < AnalysisRegion.y0) continue; @@ -326,6 +324,6 @@ if (source->modelPSF == NULL) continue; - // skip sources which are insignificant flux? - // XXX this is somewhat ad-hoc + // skip sources which are insignificant flux? + // XXX this is somewhat ad-hoc if (source->modelPSF->params->data.F32[1] < 0.1) { psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", @@ -344,19 +342,19 @@ // try fitting PSFs or extended sources depending on source->mode // these functions subtract the resulting fitted source - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_EXTENDED; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); - Next ++; - continue; - } - } else { - if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_STAR; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); - Npsf ++; - continue; - } - } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); + Next ++; + continue; + } + } else { + if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { + source->type = PM_SOURCE_TYPE_STAR; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); + Npsf ++; + continue; + } + } psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); Index: trunk/psphot/src/psphotFindFootprints.c =================================================================== --- trunk/psphot/src/psphotFindFootprints.c (revision 21359) +++ trunk/psphot/src/psphotFindFootprints.c (revision 21366) @@ -68,5 +68,5 @@ } - psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); + psphotCullPeaks(readout->image, readout->variance, recipe, detections->footprints); detections->peaks = pmFootprintArrayToPeaks(detections->footprints); psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints")); Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 21359) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 21366) @@ -68,5 +68,5 @@ // if (source->type == PM_SOURCE_TYPE_STAR && - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; if (final) { @@ -76,8 +76,8 @@ } - // generate model for sources without, or skip if we can't - if (!source->modelFlux) { - if (!pmSourceCacheModel (source, maskVal)) continue; - } + // generate model for sources without, or skip if we can't + if (!source->modelFlux) { + if (!pmSourceCacheModel (source, maskVal)) continue; + } // save the original coords @@ -96,5 +96,5 @@ // vectors to store stats for each object - // psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32); + // psVector *variance = psVectorAlloc (fitSources->n, PS_TYPE_F32); psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32); @@ -128,5 +128,5 @@ psSparseVectorElement (sparse, i, f); - // add the per-source weights (border region) + // add the per-source variances (border region) switch (SKY_FIT_ORDER) { case 1: @@ -208,5 +208,5 @@ pmSource *source = fitSources->data[i]; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal); + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -253,5 +253,5 @@ // accumulate the image statistics from the masked regions psF32 **image = readout->image->data.F32; - psF32 **weight = readout->weight->data.F32; + psF32 **variance = readout->variance->data.F32; psImageMaskType **mask = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; @@ -268,5 +268,5 @@ wt = 1.0; } else { - wt = weight[j][i]; + wt = variance[j][i]; } f = image[j][i]; Index: trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- trunk/psphot/src/psphotMakeResiduals.c (revision 21359) +++ trunk/psphot/src/psphotMakeResiduals.c (revision 21366) @@ -58,14 +58,14 @@ // for each input source: // - construct a residual image, renormalized - // - construct a renormalized weight image + // - construct a renormalized variance image // - construct a new mask image // construct the output residual table (Nx*DX,Ny*DY) // for each output pixel: - // - construct a histogram of the values & weights (interpolate to the common pixel coordinate) + // - construct a histogram of the values & variances (interpolate to the common pixel coordinate) // - measure the robust median & sigma // - reject (mask) input pixels which are outliers // - re-measure the robust median & sigma - // - set output pixel, weight, and mask + // - set output pixel, variance, and mask // these mask values do not correspond to the recipe values: they @@ -76,7 +76,7 @@ // psVectorStats - const psImageMaskType badMask = 0x01; // mask bits - const psImageMaskType poorMask = 0x02; // from psImageInterpolate - const psImageMaskType clippedMask = 0x04; // mask bit set for clipped values + const psImageMaskType badMask = 0x01; // mask bits + const psImageMaskType poorMask = 0x02; // from psImageInterpolate + const psImageMaskType clippedMask = 0x04; // mask bit set for clipped values const psVectorMaskType fmaskVal = badMask | poorMask | clippedMask; @@ -101,17 +101,17 @@ psImage *image = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32); + psImage *variance = psImageCopy (NULL, source->variance, PS_TYPE_F32); psImage *mask = psImageCopy (NULL, source->maskView, PS_TYPE_IMAGE_MASK); pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal); - // re-normalize image and weight + // re-normalize image and variance float Io = model->params->data.F32[PM_PAR_I0]; psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32)); - psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32)); - - // we interpolate the image and weight - include the mask or not? - // XXX why not the mask? - // psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal, 0.0, 0.0, badMask, poorMask, 0.0, 0); - psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0); + psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32)); + + // we interpolate the image and variance - include the mask or not? + // XXX why not the mask? + // psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal, 0.0, 0.0, badMask, poorMask, 0.0, 0); + psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0); psArrayAdd (input, 100, interp); @@ -128,5 +128,5 @@ psFree (mask); psFree (image); - psFree (weight); + psFree (variance); psFree (interp); } @@ -176,8 +176,8 @@ fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask; } else { - fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = mflux; // XXX is mflux IMAGE or VECTOR type? - } - fluxes->data.F32[i] = flux; - dfluxes->data.F32[i] = dflux; + fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = mflux; // XXX is mflux IMAGE or VECTOR type? + } + fluxes->data.F32[i] = flux; + dfluxes->data.F32[i] = dflux; if (isnan(flux)) { fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask; @@ -204,5 +204,5 @@ // mark input pixels which are more than N sigma from the median - int nKeep = 0; + int nKeep = 0; for (int i = 0; i < fluxes->n; i++) { float delta = fluxes->data.F32[i] - fluxClip->robustMedian; @@ -210,9 +210,9 @@ float swing = fabs(delta) / sigma; - // mask pixels which are out of range + // mask pixels which are out of range if (swing > nSigma) { fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = clippedMask; } - if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++; + if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++; } @@ -225,5 +225,5 @@ resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption); resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; - //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev; + //resid->variance->data.F32[oy][ox] = fluxStats->sampleStdev; if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) { @@ -231,5 +231,5 @@ } - // fprintf (stderr, "res: %2d %2d : %6.4f %6.4f %6.4f %3d %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); + // fprintf (stderr, "res: %2d %2d : %6.4f %6.4f %6.4f %3d %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); } else { @@ -268,11 +268,11 @@ float dRo = sqrt(A->data.F32[0][0]); - // fprintf (stderr, "res: %2d %2d : %6.4f %6.4f %6.4f %3d %1d\n", - // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); + // fprintf (stderr, "res: %2d %2d : %6.4f %6.4f %6.4f %3d %1d\n", + // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) { resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; } - //resid->weight->data.F32[oy][ox] = XXX; + //resid->variance->data.F32[oy][ox] = XXX; } } @@ -301,5 +301,5 @@ psphotSaveImage (NULL, resid->Rx, "resid.rx.fits"); psphotSaveImage (NULL, resid->Ry, "resid.ry.fits"); - psphotSaveImage (NULL, resid->weight, "resid.wt.fits"); + psphotSaveImage (NULL, resid->variance, "resid.wt.fits"); psphotSaveImage (NULL, resid->mask, "resid.mk.fits"); } Index: trunk/psphot/src/psphotMaskReadout.c =================================================================== --- trunk/psphot/src/psphotMaskReadout.c (revision 21359) +++ trunk/psphot/src/psphotMaskReadout.c (revision 21366) @@ -1,6 +1,6 @@ # include "psphotInternal.h" -// generate mask and weight if not defined, additional mask for restricted subregion -bool psphotSetMaskAndWeight (pmConfig *config, pmReadout *readout, psMetadata *recipe) { +// generate mask and variance if not defined, additional mask for restricted subregion +bool psphotSetMaskAndVariance (pmConfig *config, pmReadout *readout, psMetadata *recipe) { bool status; @@ -15,5 +15,5 @@ psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.BAD", PS_META_REPLACE, "user-defined mask", maskBad); - // generate mask & weight images if they don't already exit + // generate mask & variance images if they don't already exit if (!readout->mask) { if (!pmReadoutGenerateMask(readout, maskSat, maskBad)) { @@ -22,7 +22,7 @@ } } - if (!readout->weight) { - if (!pmReadoutGenerateWeight(readout, true)) { - psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight"); + if (!readout->variance) { + if (!pmReadoutGenerateVariance(readout, true)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble creating variance"); return false; } @@ -49,5 +49,5 @@ psphotSaveImage (NULL, readout->image, "image.fits"); psphotSaveImage (NULL, readout->mask, "mask.fits"); - psphotSaveImage (NULL, readout->weight, "weight.fits"); + psphotSaveImage (NULL, readout->variance, "variance.fits"); } Index: trunk/psphot/src/psphotModelWithPSF.c =================================================================== --- trunk/psphot/src/psphotModelWithPSF.c (revision 21359) +++ trunk/psphot/src/psphotModelWithPSF.c (revision 21366) @@ -20,5 +20,5 @@ paramMask = constraint->paramMask; if (paramMask != NULL) { - PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false); + PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false); PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false); } @@ -42,7 +42,7 @@ // Alpha & Beta only contain elements to represent the unmasked parameters if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) { - psAbort ("programming error: no unmasked parameters to be fit\n"); - } - + psAbort ("programming error: no unmasked parameters to be fit\n"); + } + // allocate internal arrays (current vs Guess) psImage *alpha = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32); @@ -127,7 +127,7 @@ lambda *= 0.25; - // save the new convolved model image - psFree (source->modelFlux); - source->modelFlux = pmPCMDataSaveImage(pcm); + // save the new convolved model image + psFree (source->modelFlux); + source->modelFlux = pmPCMDataSaveImage(pcm); } else { lambda *= 10.0; @@ -142,18 +142,18 @@ psTrace ("psphot", 5, "failure to calculate covariance matrix\n"); } - // set covar values which are not masked - psImageInit (covar, 0.0); - for (int j = 0, J = 0; j < params->n; j++) { - if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) { - covar->data.F32[j][j] = 1.0; - continue; - } - for (int k = 0, K = 0; k < params->n; k++) { - if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue; - covar->data.F32[j][k] = Alpha->data.F32[J][K]; - K++; - } - J++; - } + // set covar values which are not masked + psImageInit (covar, 0.0); + for (int j = 0, J = 0; j < params->n; j++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) { + covar->data.F32[j][j] = 1.0; + continue; + } + for (int k = 0, K = 0; k < params->n; k++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue; + covar->data.F32[j][k] = Alpha->data.F32[J][K]; + K++; + } + J++; + } } @@ -192,5 +192,5 @@ PS_ASSERT_PTR_NON_NULL(source, NAN); PS_ASSERT_IMAGE_NON_NULL(source->pixels, NAN); - PS_ASSERT_IMAGE_NON_NULL(source->weight, NAN); + PS_ASSERT_IMAGE_NON_NULL(source->variance, NAN); PS_ASSERT_IMAGE_NON_NULL(source->maskObj, NAN); @@ -210,6 +210,6 @@ psImageInit (pcm->model, 0.0); for (int n = 0; n < params->n; n++) { - if (!pcm->dmodels->data[n]) continue; - psImageInit (pcm->dmodels->data[n], 0.0); + if (!pcm->dmodels->data[n]) continue; + psImageInit (pcm->dmodels->data[n], 0.0); } @@ -218,24 +218,24 @@ for (psS32 j = 0; j < source->pixels->numCols; j++) { - // XXX can we skip some of the data points where the model - // is not going to be fitted?? + // XXX can we skip some of the data points where the model + // is not going to be fitted?? // skip masked points - // XXX probably should not skipped masked points: - // XXX skip if convolution of unmasked pixels will not see this pixel + // XXX probably should not skipped masked points: + // XXX skip if convolution of unmasked pixels will not see this pixel // if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) { - // continue; - // } - - // skip zero-weight points - // XXX why is this not masked? - // if (source->weight->data.F32[i][j] == 0) { - // continue; - // } + // continue; + // } + + // skip zero-variance points + // XXX why is this not masked? + // if (source->variance->data.F32[i][j] == 0) { + // continue; + // } // skip nan value points - // XXX why is this not masked? + // XXX why is this not masked? // if (!isfinite(source->pixels->data.F32[i][j])) { - // continue; - // } + // continue; + // } // Convert i/j to image space: @@ -243,11 +243,11 @@ coord->data.F32[1] = (psF32) (i + source->pixels->row0); - pcm->model->data.F32[i][j] = func (deriv, params, coord); - - for (int n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } - psImage *dmodel = pcm->dmodels->data[n]; - dmodel->data.F32[i][j] = deriv->data.F32[n]; - } + pcm->model->data.F32[i][j] = func (deriv, params, coord); + + for (int n = 0; n < params->n; n++) { + if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } + psImage *dmodel = pcm->dmodels->data[n]; + dmodel->data.F32[i][j] = deriv->data.F32[n]; + } } } @@ -258,12 +258,12 @@ psImageConvolveDirect (pcm->modelConv, pcm->model, psf); for (int n = 0; n < pcm->dmodels->n; n++) { - if (pcm->dmodels->data[n] == NULL) continue; - psImage *dmodel = pcm->dmodels->data[n]; - psImage *dmodelConv = pcm->dmodelsConv->data[n]; - psImageConvolveDirect (dmodelConv, dmodel, psf); + if (pcm->dmodels->data[n] == NULL) continue; + psImage *dmodel = pcm->dmodels->data[n]; + psImage *dmodelConv = pcm->dmodelsConv->data[n]; + psImageConvolveDirect (dmodelConv, dmodel, psf); } // XXX TEST : SAVE IMAGES -# if (SAVE_IMAGES) +# if (SAVE_IMAGES) psphotSaveImage (NULL, psf->image, "psf.fits"); psphotSaveImage (NULL, pcm->model, "model.fits"); @@ -271,8 +271,8 @@ psphotSaveImage (NULL, source->pixels, "obj.fits"); psphotSaveImage (NULL, source->maskObj, "mask.fits"); - psphotSaveImage (NULL, source->weight, "weight.fits"); + psphotSaveImage (NULL, source->variance, "variance.fits"); # endif - // 2 *** accumulate alpha & beta + // 2 *** accumulate alpha & beta // zero alpha and beta for summing below @@ -283,11 +283,11 @@ for (psS32 i = 0; i < source->pixels->numRows; i++) { for (psS32 j = 0; j < source->pixels->numCols; j++) { - // XXX are we doing the right thing with the mask? + // XXX are we doing the right thing with the mask? // skip masked points if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) { continue; } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { + // skip zero-variance points + if (source->variance->data.F32[i][j] == 0) { continue; } @@ -297,28 +297,28 @@ } - float ymodel = pcm->modelConv->data.F32[i][j]; - float yweight = 1.0 / source->weight->data.F32[i][j]; - float delta = ymodel - source->pixels->data.F32[i][j]; - - chisq += PS_SQR(delta) * yweight; - - if (isnan(delta)) psAbort("nan in delta"); - if (isnan(chisq)) psAbort("nan in chisq"); - - // alpha & beta only contain unmasked elements - for (int n1 = 0, N1 = 0; n1 < params->n; n1++) { - if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue; - psImage *dmodel = pcm->dmodelsConv->data[n1]; - float weight = dmodel->data.F32[i][j] * yweight; - for (int n2 = 0, N2 = 0; n2 <= n1; n2++) { - if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue; - dmodel = pcm->dmodelsConv->data[n2]; - alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j]; - N2++; - } - beta->data.F32[N1] += weight * delta; - N1++; - } - } + float ymodel = pcm->modelConv->data.F32[i][j]; + float yweight = 1.0 / source->variance->data.F32[i][j]; + float delta = ymodel - source->pixels->data.F32[i][j]; + + chisq += PS_SQR(delta) * yweight; + + if (isnan(delta)) psAbort("nan in delta"); + if (isnan(chisq)) psAbort("nan in chisq"); + + // alpha & beta only contain unmasked elements + for (int n1 = 0, N1 = 0; n1 < params->n; n1++) { + if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue; + psImage *dmodel = pcm->dmodelsConv->data[n1]; + float weight = dmodel->data.F32[i][j] * yweight; + for (int n2 = 0, N2 = 0; n2 <= n1; n2++) { + if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue; + dmodel = pcm->dmodelsConv->data[n2]; + alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j]; + N2++; + } + beta->data.F32[N1] += weight * delta; + N1++; + } + } } @@ -356,7 +356,7 @@ pcm->dmodels = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - pcm->dmodels->data[n] = NULL; - if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } - pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodels->data[n] = NULL; + if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } + pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } @@ -365,7 +365,7 @@ pcm->dmodelsConv = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - pcm->dmodelsConv->data[n] = NULL; - if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } - pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodelsConv->data[n] = NULL; + if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } + pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } @@ -376,5 +376,5 @@ psImage *model = psImageCopy (NULL, pcm->modelConv, PS_TYPE_F32); - + return model; } @@ -383,11 +383,11 @@ * * we have a function func(param; value) - + * basic LMM: - + - fill in the data (x, y) - + chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func) - + while () { GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda) @@ -396,6 +396,6 @@ convergence tests... } - - + + ** GuessABP: Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 21359) +++ trunk/psphot/src/psphotOutput.c (revision 21366) @@ -18,5 +18,5 @@ pmReadout *psphotSelectBackgroundStdev (pmConfig *config, - const pmFPAview *view) { + const pmFPAview *view) { bool status; @@ -77,6 +77,6 @@ continue; } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { + // skip zero-variance points + if (source->variance->data.F32[i][j] == 0) { continue; } @@ -86,5 +86,5 @@ (i + source->pixels->row0), source->pixels->data.F32[i][j], - 1.0 / source->weight->data.F32[i][j], + 1.0 / source->variance->data.F32[i][j], source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]); } @@ -124,10 +124,10 @@ if (model == NULL) continue; - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - nEXT ++; - } - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { - nCR ++; - } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + nEXT ++; + } + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { + nCR ++; + } nSrc ++; } @@ -242,57 +242,57 @@ for (int i = 0; i < try->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) - continue; - - pmSource *source = try->sources->data[i]; - float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - - // set the mask and subtract the PSF model - // XXX should we be using maskObj? should we be unsetting the mask? - // use pmModelSub because modelFlux has not been generated - assert (source->maskObj); - psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal); - pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal); - psImageKeepCircle (source->maskObj, x, y, radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + // masked for: bad model fit, outlier in parameters + if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = try->sources->data[i]; + float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + + // set the mask and subtract the PSF model + // XXX should we be using maskObj? should we be unsetting the mask? + // use pmModelSub because modelFlux has not been generated + assert (source->maskObj); + psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal); + pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal); + psImageKeepCircle (source->maskObj, x, y, radius, "AND", PS_NOT_IMAGE_MASK(markVal)); } FILE *f = fopen ("shapes.dat", "w"); for (int i = 0; i < try->sources->n; i++) { - psF32 inPar[10]; // must be psF32 to pmPSF_FitToModel - - // masked for: bad model fit, outlier in parameters - if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; - - pmSource *source = try->sources->data[i]; - psF32 *outPar = source->modelEXT->params->data.F32; - - psEllipseShape shape; - - shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = outPar[PM_PAR_SXY]; - - psEllipsePol pol = pmPSF_ModelToFit (outPar); - inPar[PM_PAR_E0] = pol.e0; - inPar[PM_PAR_E1] = pol.e1; - inPar[PM_PAR_E2] = pol.e2; - pmPSF_FitToModel (inPar, 0.1); - - psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); - psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1); - - psEllipsePol pol2 = psEllipseAxesToPol (axes1); - - fprintf (f, "%3d %7.2f %7.2f %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n", - i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS], - outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], - pol.e0, pol.e1, pol.e2, - pol2.e0, pol2.e1, pol2.e2, - inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY], - axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD, - axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD - ); + psF32 inPar[10]; // must be psF32 to pmPSF_FitToModel + + // masked for: bad model fit, outlier in parameters + if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = try->sources->data[i]; + psF32 *outPar = source->modelEXT->params->data.F32; + + psEllipseShape shape; + + shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; + shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; + shape.sxy = outPar[PM_PAR_SXY]; + + psEllipsePol pol = pmPSF_ModelToFit (outPar); + inPar[PM_PAR_E0] = pol.e0; + inPar[PM_PAR_E1] = pol.e1; + inPar[PM_PAR_E2] = pol.e2; + pmPSF_FitToModel (inPar, 0.1); + + psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1); + + psEllipsePol pol2 = psEllipseAxesToPol (axes1); + + fprintf (f, "%3d %7.2f %7.2f %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n", + i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS], + outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], + pol.e0, pol.e1, pol.e2, + pol2.e0, pol2.e1, pol2.e2, + inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY], + axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD, + axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD + ); } fclose (f); Index: trunk/psphot/src/psphotRadialProfile.c =================================================================== --- trunk/psphot/src/psphotRadialProfile.c (revision 21359) +++ trunk/psphot/src/psphotRadialProfile.c (revision 21366) @@ -11,7 +11,7 @@ flux->data.F32[A] = flux->data.F32[B]; \ flux->data.F32[B] = tmp; \ - tmp = weight->data.F32[A]; \ - weight->data.F32[A] = weight->data.F32[B]; \ - weight->data.F32[B] = tmp; \ + tmp = variance->data.F32[A]; \ + variance->data.F32[A] = variance->data.F32[B]; \ + variance->data.F32[B] = tmp; \ } \ } @@ -31,9 +31,9 @@ source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32); source->extpars->profile->flux = psVectorAllocEmpty (nPts, PS_TYPE_F32); - source->extpars->profile->weight = psVectorAllocEmpty (nPts, PS_TYPE_F32); + source->extpars->profile->variance = psVectorAllocEmpty (nPts, PS_TYPE_F32); psVector *radius = source->extpars->profile->radius; psVector *flux = source->extpars->profile->flux; - psVector *weight = source->extpars->profile->weight; + psVector *variance = source->extpars->profile->variance; // XXX use the extended source model here for Xo, Yo? @@ -41,5 +41,5 @@ int n = 0; - + float Xo = 0.0; float Yo = 0.0; @@ -57,10 +57,10 @@ radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ; flux->data.F32[n] = source->pixels->data.F32[iy][ix]; - weight->data.F32[n] = source->weight->data.F32[iy][ix]; + variance->data.F32[n] = source->variance->data.F32[iy][ix]; n++; } } radius->n = n; - weight->n = n; + variance->n = n; flux->n = n; Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 21359) +++ trunk/psphot/src/psphotReadout.c (revision 21366) @@ -40,5 +40,5 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndWeight (config, readout, recipe); + psphotSetMaskAndVariance (config, readout, recipe); if (!strcasecmp (breakPt, "NOTHING")) { return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); Index: trunk/psphot/src/psphotReadoutFindPSF.c =================================================================== --- trunk/psphot/src/psphotReadoutFindPSF.c (revision 21359) +++ trunk/psphot/src/psphotReadoutFindPSF.c (revision 21366) @@ -18,6 +18,6 @@ PS_ASSERT_PTR_NON_NULL (readout, false); - // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndWeight (config, readout, recipe); + // Generate the mask and variance images, including the user-defined analysis region of interest + psphotSetMaskAndVariance (config, readout, recipe); // display the image, weight, mask (ch 1,2,3) Index: trunk/psphot/src/psphotReadoutKnownSources.c =================================================================== --- trunk/psphot/src/psphotReadoutKnownSources.c (revision 21359) +++ trunk/psphot/src/psphotReadoutKnownSources.c (revision 21366) @@ -19,5 +19,5 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndWeight (config, readout, recipe); + psphotSetMaskAndVariance (config, readout, recipe); // display the image, weight, mask (ch 1,2,3) @@ -41,8 +41,8 @@ // use the peak measured in the moments analysis: for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - source->peak->flux = source->moments->Peak; + pmSource *source = sources->data[i]; + source->peak->flux = source->moments->Peak; } - + // classify sources based on moments, brightness (psf is not known) if (!psphotRoughClass (readout, sources, recipe, false)) { @@ -53,6 +53,6 @@ pmPSF *psf = psphotChoosePSF (readout, sources, recipe); if (!psf) { - psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + psLogMsg ("psphot", 3, "failure to construct a psf model"); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } psphotVisualShowPSFModel (readout, psf); Index: trunk/psphot/src/psphotReadoutMinimal.c =================================================================== --- trunk/psphot/src/psphotReadoutMinimal.c (revision 21359) +++ trunk/psphot/src/psphotReadoutMinimal.c (revision 21366) @@ -5,5 +5,5 @@ // used by ppSub. -// NOTE: ppSub needs to perform extended source analysis for comets and trails. +// NOTE: ppSub needs to perform extended source analysis for comets and trails. bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view) { @@ -27,5 +27,5 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndWeight (config, readout, recipe); + psphotSetMaskAndVariance (config, readout, recipe); // display the image, weight, mask (ch 1,2,3) @@ -78,5 +78,5 @@ // XXX eventually, add the extended source fits here -# if (0) +# if (0) // measure source size for the remaining sources psphotSourceSize (config, readout, sources, recipe, 0); Index: trunk/psphot/src/psphotSignificanceImage.c =================================================================== --- trunk/psphot/src/psphotSignificanceImage.c (revision 21359) +++ trunk/psphot/src/psphotSignificanceImage.c (revision 21366) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -// In this function, we smooth the image and weight, then generate the significance image : +// In this function, we smooth the image and variance, then generate the significance image : // (S/N)^2. If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the // smoothing kernel. @@ -10,5 +10,5 @@ bool guess = false; - // smooth the image and weight map + // smooth the image and variance map psTimerStart ("psphot.smooth"); bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve @@ -48,5 +48,5 @@ psLogMsg("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("psphot.smooth")); - // Smooth the weight, applying the mask as we go. The variance is smoothed by the PSF^2, + // Smooth the variance, applying the mask as we go. The variance is smoothed by the PSF^2, // renomalized to maintain the input level of the variance. We achieve this by smoothing // with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity normalization. Note that @@ -55,8 +55,8 @@ // measurements based on apertures comparable to or larger than the smoothing kernel, the // effective per-pixel variance is maintained. - psImage *smooth_wt = psImageCopy(NULL, readout->weight, PS_TYPE_F32); + psImage *smooth_wt = psImageCopy(NULL, readout->variance, PS_TYPE_F32); psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss); - psLogMsg("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark("psphot.smooth")); + psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth")); psImage *mask = readout->mask; @@ -100,4 +100,12 @@ float factor = guess ? 4.0 * M_PI * PS_SQR(SIGMA_SMTH) : 4.0 * M_PI * PS_SQR(SIGMA_SMTH); + + // Correct the correction factor for the covariance produced by the (potentially multiple) smoothing + psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing + psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + psFree(kernel); + factor /= psImageCovarianceFactor(covar); + psFree(covar); + // record the effective area and significance scaling factor float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH); Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 21359) +++ trunk/psphot/src/psphotSourceSize.c (revision 21366) @@ -2,5 +2,5 @@ # include -static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask, +static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, psImageMaskType maskVal, const pmModel *model, float Ro); @@ -62,9 +62,9 @@ psF32 **resid = source->pixels->data.F32; - psF32 **weight = source->weight->data.F32; + psF32 **variance = source->variance->data.F32; psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; // check for extendedness: measure the delta flux significance at the 1 sigma contour - source->extNsigma = psphotModelContour(source->pixels, source->weight, source->maskObj, maskVal, + source->extNsigma = psphotModelContour(source->pixels, source->variance, source->maskObj, maskVal, source->modelPSF, 1.0); @@ -103,25 +103,25 @@ // Compare the central pixel with those on either side, for the four possible lines through it. - // Soften weights (add systematic error) - float softening = soft * PS_SQR(source->peak->flux); // Softening for weights + // Soften variances (add systematic error) + float softening = soft * PS_SQR(source->peak->flux); // Softening for variances // Across the middle: y = 0 float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; - float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1]; + float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; float nX = cX / sqrtf(dcX + softening); // Up the centre: x = 0 float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; - float dcY = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak+0] + weight[yPeak+1][xPeak+0]; + float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; float nY = cY / sqrtf(dcY + softening); // Diagonal: x = y float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; - float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1]; + float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; float nL = cL / sqrtf(dcL + softening); // Diagonal: x = - y float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; - float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1]; + float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; float nR = cR / sqrtf(dcR + softening); @@ -160,6 +160,6 @@ // this source is thought to be a cosmic ray. flag the detection and mask the pixels if (source->crNsigma > CR_NSIGMA_LIMIT) { - // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); - psphotMaskCosmicRay_Old (source, maskVal, crMask); + // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); + psphotMaskCosmicRay_Old (source, maskVal, crMask); } } @@ -190,5 +190,5 @@ // deviation in sigmas. This is measured on the residual image - should we ignore negative // deviations? -static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask, +static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, psImageMaskType maskVal, const pmModel *model, float Ro) { @@ -240,5 +240,5 @@ if (yPixM >= 0 && yPixM < image->numRows && !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixM][xPix] / sqrtf(weight->data.F32[yPixM][xPix]); + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); nSigma += dSigma; nPts++; @@ -251,5 +251,5 @@ if (yPixP >= 0 && yPixP < image->numRows && !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixP][xPix] / sqrtf(weight->data.F32[yPixP][xPix]); + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); nSigma += dSigma; nPts++; @@ -274,26 +274,26 @@ pmFootprint *footprint = peak->footprint; if (!footprint) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); - return true; + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRay_Old (source, maskVal, crMask); + return true; } if (!footprint->spans) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); - return true; + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRay_Old (source, maskVal, crMask); + return true; } // mask all of the pixels covered by the spans of the footprint for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int iy = span1->y; - int xs = span1->x0; - int xe = span1->x1; - - for (int ix = xs; ix < xe; ix++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + pmSpan *span1 = footprint->spans->data[j]; + + int iy = span1->y; + int xs = span1->x0; + int xe = span1->x1; + + for (int ix = xs; ix < xe; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } return true; @@ -308,5 +308,5 @@ psImage *mask = source->maskView; psImage *pixels = source->pixels; - psImage *weight = source->weight; + psImage *variance = source->variance; // XXX This should be a recipe variable @@ -318,14 +318,14 @@ // mark the pixels in this row to the left, then the right for (int ix = xo; ix >= 0; ix--) { - float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } } for (int ix = xo + 1; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } } @@ -333,33 +333,33 @@ // first go up: for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } // next go down: for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } return true; Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 21359) +++ trunk/psphot/src/psphotVisual.c (revision 21366) @@ -2,11 +2,11 @@ // this function displays representative images as the psphot analysis progresses: -// 0 : image, 1 : weight -// 0 : backsub, 1 : weight, 2 : backgnd -// 0 : backsub, 1 : weight, 2 : signif +// 0 : image, 1 : variance +// 0 : backsub, 1 : variance, 2 : backgnd +// 0 : backsub, 1 : variance, 2 : signif // (overlay peaks on images) // (overlay footprints on images) // (overlay moments on images) -// (overlay rough class on images) +// (overlay rough class on images) // 0 : backsub, 1 : psfpos, 2: psfsub // 0 : backsub, 1 : lin_resid, 2: psfsub @@ -40,6 +40,6 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) { - fprintf (stderr, "failed to get background values\n"); - return false; + fprintf (stderr, "failed to get background values\n"); + return false; } @@ -49,8 +49,8 @@ ALLOCATE (image.data2d, float *, image.Ny); for (int iy = 0; iy < image.Ny; iy++) { - ALLOCATE (image.data2d[iy], float, image.Nx); - for (int ix = 0; ix < image.Nx; ix++) { - image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]; - } + ALLOCATE (image.data2d[iy], float, image.Nx); + for (int ix = 0; ix < image.Nx; ix++) { + image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]; + } } @@ -60,10 +60,10 @@ data.range = 32; data.logflux = 0; - + KiiSetChannel (kapaFD, channel); KiiNewPicture2D (kapaFD, &image, &data, &coords); for (int iy = 0; iy < image.Ny; iy++) { - free (image.data2d[iy]); + free (image.data2d[iy]); } free (image.data2d); @@ -86,6 +86,6 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) { - fprintf (stderr, "failed to get background values\n"); - return false; + fprintf (stderr, "failed to get background values\n"); + return false; } @@ -99,5 +99,5 @@ data.range = 5*stats->robustStdev; data.logflux = 0; - + KiiSetChannel (kapaFD, channel); KiiNewPicture2D (kapaFD, &image, &data, &coords); @@ -126,5 +126,5 @@ data.range = max - min; data.logflux = 0; - + KiiSetChannel (kapaFD, channel); KiiNewPicture2D (kapaFD, &image, &data, &coords); @@ -139,13 +139,13 @@ if (kapa == -1) { kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } // psphotVisualShowMask (kapa, readout->mask, "mask", 2); - psphotVisualScaleImage (kapa, readout->weight, "weight", 1); + psphotVisualScaleImage (kapa, readout->variance, "variance", 1); psphotVisualScaleImage (kapa, readout->image, "image", 0); @@ -155,5 +155,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -168,10 +168,10 @@ if (kapa == -1) { kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } bool status = false; @@ -179,7 +179,7 @@ if (file->mode == PM_FPA_MODE_INTERNAL) { - backgnd = file->readout; + backgnd = file->readout; } else { - backgnd = pmFPAviewThisReadout (view, file->fpa); + backgnd = pmFPAviewThisReadout (view, file->fpa); } @@ -192,5 +192,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -203,10 +203,10 @@ if (kapa == -1) { kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } // XXX test: image->data.F32[10][10] = 10000; @@ -218,5 +218,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -227,10 +227,10 @@ int Noverlay; KiiOverlay *overlay; - + if (!isVisual) return true; if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; + fprintf (stderr, "kapa not opened, skipping\n"); + return false; } @@ -244,34 +244,34 @@ for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - if (peak == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_BOX; - overlay[Noverlay].x = peak->xf; - overlay[Noverlay].y = peak->yf; - overlay[Noverlay].dx = 2.0; - overlay[Noverlay].dy = 2.0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmPeak *peak = peaks->data[i]; + if (peak == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 2.0; + overlay[Noverlay].dy = 2.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; # if (0) - overlay[Noverlay].type = KII_OVERLAY_BOX; - overlay[Noverlay].x = peak->x; - overlay[Noverlay].y = peak->y; - overlay[Noverlay].dx = 1.0; - overlay[Noverlay].dy = 1.0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = peak->xf; - overlay[Noverlay].y = peak->yf; - overlay[Noverlay].dx = 2.0; - overlay[Noverlay].dy = 2.0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->x; + overlay[Noverlay].y = peak->y; + overlay[Noverlay].dx = 1.0; + overlay[Noverlay].dy = 1.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 2.0; + overlay[Noverlay].dy = 2.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; # endif } @@ -296,5 +296,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -305,10 +305,10 @@ int Noverlay; KiiOverlay *overlay; - + if (!isVisual) return true; if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; + fprintf (stderr, "kapa not opened, skipping\n"); + return false; } @@ -323,81 +323,81 @@ for (int i = 0; i < footprints->n; i++) { - pmSpan *span = NULL; - - pmFootprint *footprint = footprints->data[i]; - if (footprint == NULL) continue; - if (footprint->spans == NULL) continue; - if (footprint->spans->n < 1) continue; - - // draw the top - span = footprint->spans->data[0]; - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = span->x0; - overlay[Noverlay].y = span->y; - overlay[Noverlay].dx = span->x1 - span->x0; - overlay[Noverlay].dy = 0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - int ys = span->y; - int x0s = span->x0; - int x1s = span->x1; - - // draw the outer span edges - for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int ye = span1->y; - int x0e = span1->x0; - int x1e = span1->x1; - - // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) - // find all of the spans in this row and generate x0e, x01: - for (int k = j + 1; k < footprint->spans->n; k++) { - pmSpan *span2 = footprint->spans->data[k]; - if (span2->y > span1->y) break; - x0e = PS_MIN (x0e, span2->x0); - x1e = PS_MAX (x1e, span2->x1); - j++; - } - - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = x0s; - overlay[Noverlay].y = ys; - overlay[Noverlay].dx = x0e - x0s; - overlay[Noverlay].dy = ye - ys; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = x1s; - overlay[Noverlay].y = ys; - overlay[Noverlay].dx = x1e - x1s; - overlay[Noverlay].dy = ye - ys; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - ys = ye; - x0s = x0e; - x1s = x1e; - } - - // draw the bottom - span = footprint->spans->data[footprint->spans->n - 1]; - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = span->x0; - overlay[Noverlay].y = span->y; - overlay[Noverlay].dx = span->x1 - span->x0; - overlay[Noverlay].dy = 0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + pmSpan *span = NULL; + + pmFootprint *footprint = footprints->data[i]; + if (footprint == NULL) continue; + if (footprint->spans == NULL) continue; + if (footprint->spans->n < 1) continue; + + // draw the top + span = footprint->spans->data[0]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + int ys = span->y; + int x0s = span->x0; + int x1s = span->x1; + + // draw the outer span edges + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int ye = span1->y; + int x0e = span1->x0; + int x1e = span1->x1; + + // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) + // find all of the spans in this row and generate x0e, x01: + for (int k = j + 1; k < footprint->spans->n; k++) { + pmSpan *span2 = footprint->spans->data[k]; + if (span2->y > span1->y) break; + x0e = PS_MIN (x0e, span2->x0); + x1e = PS_MAX (x1e, span2->x1); + j++; + } + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x0s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x0e - x0s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x1s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x1e - x1s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + ys = ye; + x0s = x0e; + x1s = x1e; + } + + // draw the bottom + span = footprint->spans->data[footprint->spans->n - 1]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); } @@ -410,5 +410,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -419,5 +419,5 @@ int Noverlay; KiiOverlay *overlay; - + psEllipseMoments emoments; psEllipseAxes axes; @@ -426,6 +426,6 @@ if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; + fprintf (stderr, "kapa not opened, skipping\n"); + return false; } @@ -436,25 +436,25 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - pmMoments *moments = source->moments; - if (moments == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = moments->Mx; - overlay[Noverlay].y = moments->My; - - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - - overlay[Noverlay].dx = 2.0*axes.major; - overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; // XXXXXXXX the axes angle is negative to display of object on kapa - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->Mx; + overlay[Noverlay].y = moments->My; + + emoments.x2 = moments->Mxx; + emoments.xy = moments->Mxy; + emoments.y2 = moments->Myy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; // XXXXXXXX the axes angle is negative to display of object on kapa + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -467,5 +467,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } @@ -482,10 +482,10 @@ if (kapa3 == -1) { kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } KapaClearPlots (kapa3); @@ -494,5 +494,5 @@ // there are N regions: use the first (guaranteed to exist) to get the overall limits psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000"); - + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); @@ -558,42 +558,42 @@ // XXX draw N circles to outline the clumps - { - // draw a circle centered on psfX,Y with size of the psf limit - psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); - psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); - - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); - float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); - - graphdata.color = KapaColorByName ("blue"); - graphdata.style = 0; - - for (int n = 0; n < nRegions; n++) { - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); - float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); - float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - float Rx = psfdX * PSF_CLUMP_NSIGMA; - float Ry = psfdY * PSF_CLUMP_NSIGMA; - - for (int i = 0; i < xLimit->n; i++) { - xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; - yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; - } - KapaPrepPlot (kapa3, xLimit->n, &graphdata); - KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x"); - KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y"); - } - psFree (xLimit); - psFree (yLimit); + { + // draw a circle centered on psfX,Y with size of the psf limit + psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); + psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); + + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 0; + + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + float Rx = psfdX * PSF_CLUMP_NSIGMA; + float Ry = psfdY * PSF_CLUMP_NSIGMA; + + for (int i = 0; i < xLimit->n; i++) { + xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; + yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; + } + KapaPrepPlot (kapa3, xLimit->n, &graphdata); + KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y"); + } + psFree (xLimit); + psFree (yLimit); } # if (0) - // *** make a histogram of the source counts in the x and y directions + // *** make a histogram of the source counts in the x and y directions psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0); psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0); @@ -605,10 +605,10 @@ psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); for (int i = 0; i < nX->nums->n; i++) { - dX->data.F32[i] = nX->nums->data.S32[i]; - vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); + dX->data.F32[i] = nX->nums->data.S32[i]; + vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); } for (int i = 0; i < nY->nums->n; i++) { - dY->data.F32[i] = nY->nums->data.S32[i]; - vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); + dY->data.F32[i] = nY->nums->data.S32[i]; + vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); } @@ -640,5 +640,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -650,5 +650,5 @@ int Noverlay; KiiOverlay *overlay; - + psEllipseMoments emoments; psEllipseAxes axes; @@ -660,28 +660,28 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - if (source->type != type) continue; - if (mode && !(source->mode & mode)) continue; - - pmMoments *moments = source->moments; - if (moments == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = moments->Mx; - overlay[Noverlay].y = moments->My; - - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - - overlay[Noverlay].dx = 2.0*axes.major; - overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + if (source->type != type) continue; + if (mode && !(source->mode & mode)) continue; + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->Mx; + overlay[Noverlay].y = moments->My; + + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -697,6 +697,6 @@ if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; + fprintf (stderr, "kapa not opened, skipping\n"); + return false; } @@ -716,5 +716,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -727,10 +727,10 @@ if (kapa2 == -1) { kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } int DX = 64; @@ -739,5 +739,5 @@ psImage *psfMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32); psImageInit (psfMosaic, 0.0); - + psImage *funMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32); psImageInit (funMosaic, 0.0); @@ -750,28 +750,28 @@ // generate a fake model at each of the 3x3 image grid positions for (int x = -2; x <= +2; x ++) { - for (int y = -2; y <= +2; y ++) { - // use the center of the center pixel of the image - float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5; - float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5; - - // assign the x and y coords to the image center - // create an object with center intensity of 1000 - modelRef->params->data.F32[PM_PAR_SKY] = 0; - modelRef->params->data.F32[PM_PAR_I0] = 1000; - modelRef->params->data.F32[PM_PAR_XPOS] = xc; - modelRef->params->data.F32[PM_PAR_YPOS] = yc; - - // create modelPSF from this model - pmModel *model = pmModelFromPSF (modelRef, psf); - if (!model) continue; - - // place the reference object in the image center - // no need to mask the source here - // XXX should we measure this for the analytical model only or the full model? - pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - psFree (model); - } + for (int y = -2; y <= +2; y ++) { + // use the center of the center pixel of the image + float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5; + float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + modelRef->params->data.F32[PM_PAR_XPOS] = xc; + modelRef->params->data.F32[PM_PAR_YPOS] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + if (!model) continue; + + // place the reference object in the image center + // no need to mask the source here + // XXX should we measure this for the analytical model only or the full model? + pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + psFree (model); + } } @@ -792,5 +792,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -805,10 +805,10 @@ if (kapa2 == -1) { kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -827,8 +827,8 @@ // counters to track the size of the image and area used in a row - int dX = 0; // starting corner of next box - int dY = 0; // height of row so far - int NX = 20*DX; // full width of output image - int NY = 0; // total height of output image + int dX = 0; // starting corner of next box + int dY = 0; // height of row so far + int NX = 20*DX; // full width of output image + int NY = 0; // total height of output image // first, examine the PSF stars: @@ -838,30 +838,30 @@ pmSource *source = sources->data[i]; - bool keep = false; + bool keep = false; keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); - if (!keep) continue; - - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows - - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + if (!keep) continue; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } NY += DY; @@ -873,7 +873,7 @@ psImageInit (outsub, 0.0); - int Xo = 0; // starting corner of next box - int Yo = 0; // starting corner of next box - dY = 0; // height of row so far + int Xo = 0; // starting corner of next box + int Yo = 0; // starting corner of next box + dY = 0; // height of row so far int nPSF = 0; @@ -885,52 +885,52 @@ pmSource *source = sources->data[i]; - bool keep = false; + bool keep = false; if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - nPSF ++; - keep = true; - } - if (!keep) continue; - - if (Xo + DX > NX) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - - psphotSubWithTest (source, false, maskVal); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded - - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - - psphotAddWithTest (source, true, maskVal); // replace source if subtracted - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - - psphotSubWithTest (source, false, maskVal); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) - - Xo = DX; - dY = DY; - } - } else { - // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - - psphotSubWithTest (source, false, maskVal); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) - - Xo += DX; - dY = PS_MAX (dY, DY); - } + nPSF ++; + keep = true; + } + if (!keep) continue; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) + + Xo = DX; + dY = DY; + } + } else { + // extend this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) + + Xo += DX; + dY = PS_MAX (dY, DY); + } } @@ -943,5 +943,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } @@ -965,10 +965,10 @@ if (kapa2 == -1) { kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -987,8 +987,8 @@ // counters to track the size of the image and area used in a row - int dX = 0; // starting corner of next box - int dY = 0; // height of row so far - int NX = 10*DX; // full width of output image - int NY = 0; // total height of output image + int dX = 0; // starting corner of next box + int dY = 0; // height of row so far + int NX = 10*DX; // full width of output image + int NY = 0; // total height of output image // first, examine the PSF and SAT stars: @@ -998,30 +998,30 @@ pmSource *source = sources->data[i]; - bool keep = false; + bool keep = false; keep |= (source->mode & PM_SOURCE_MODE_SATSTAR); - if (!keep) continue; - - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows - - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + if (!keep) continue; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } NY += DY; @@ -1031,7 +1031,7 @@ psImageInit (outsat, 0.0); - int Xo = 0; // starting corner of next box - int Yo = 0; // starting corner of next box - dY = 0; // height of row so far + int Xo = 0; // starting corner of next box + int Yo = 0; // starting corner of next box + dY = 0; // height of row so far int nSAT = 0; @@ -1043,42 +1043,42 @@ pmSource *source = sources->data[i]; - bool keep = false; + bool keep = false; if (source->mode & PM_SOURCE_MODE_SATSTAR) { - nSAT ++; - keep = true; - } - if (!keep) continue; - - if (Xo + DX > NX) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted - psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded - - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - psphotAddWithTest (source, true, maskVal); // replace source if subtracted - psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // replace source (has been subtracted) - - Xo = DX; - dY = DY; - } - } else { - // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted - psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // replace source (has been subtracted) - - Xo += DX; - dY = PS_MAX (dY, DY); - } + nSAT ++; + keep = true; + } + if (!keep) continue; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + psphotSetState (source, true, maskVal); // replace source (has been subtracted) + + Xo = DX; + dY = DY; + } + } else { + // extend this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + psphotSetState (source, true, maskVal); // replace source (has been subtracted) + + Xo += DX; + dY = PS_MAX (dY, DY); + } } @@ -1090,5 +1090,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } @@ -1117,23 +1117,23 @@ float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; for (int iy = 0; iy < source->pixels->numRows; iy++) { - for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { - // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; - rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; - Rb->data.F32[nb] = log10(rb->data.F32[nb]); - fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); - nb++; - } else { - // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; - rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; - Rg->data.F32[ng] = log10(rg->data.F32[ng]); - fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]); - ng++; - } - } - } - + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; + Rb->data.F32[nb] = log10(rb->data.F32[nb]); + fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); + nb++; + } else { + // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; + Rg->data.F32[ng] = log10(rg->data.F32[ng]); + fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]); + ng++; + } + } + } + // reset source Add/Sub state to recorded - psphotSetState (source, state, maskVal); + psphotSetState (source, state, maskVal); KapaInitGraph (&graphdata); @@ -1148,10 +1148,10 @@ graphdata.ymax = +5.05; KapaSetLimits (myKapa, &graphdata); - + KapaSetFont (myKapa, "helvetica", 14); KapaBox (myKapa, &graphdata); KapaSendLabel (myKapa, "radius (pixels)", KAPA_LABEL_XM); KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM); - + graphdata.color = KapaColorByName ("black"); graphdata.ptype = 2; @@ -1169,5 +1169,5 @@ KapaPlotVector (myKapa, nb, rb->data.F32, "x"); KapaPlotVector (myKapa, nb, fb->data.F32, "y"); - + // ** loglog ** KapaSelectSection (myKapa, "loglog"); @@ -1179,10 +1179,10 @@ graphdata.ymax = +5.05; KapaSetLimits (myKapa, &graphdata); - + KapaSetFont (myKapa, "helvetica", 14); KapaBox (myKapa, &graphdata); KapaSendLabel (myKapa, "log radius (pixels)", KAPA_LABEL_XM); KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM); - + graphdata.color = KapaColorByName ("black"); graphdata.ptype = 2; @@ -1212,5 +1212,5 @@ bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) { - KapaSection section; // put the positive profile in one and the residuals in another? + KapaSection section; // put the positive profile in one and the residuals in another? if (!isVisual) return true; @@ -1218,10 +1218,10 @@ if (kapa3 == -1) { kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1257,20 +1257,20 @@ if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - psphotVisualPlotRadialProfile (kapa3, source, maskVal); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - if (key[0] == 'e') { - KapaClearPlots (kapa3); - } - if (key[0] == 's') { - break; - } - } + psphotVisualPlotRadialProfile (kapa3, source, maskVal); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + if (key[0] == 'e') { + KapaClearPlots (kapa3); + } + if (key[0] == 's') { + break; + } + } return true; @@ -1282,5 +1282,5 @@ int NoverlayO, NOVERLAYO; KiiOverlay *overlayE, *overlayO; - + psEllipseMoments emoments; psEllipseAxes axes; @@ -1289,6 +1289,6 @@ if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; + fprintf (stderr, "kapa not opened, skipping\n"); + return false; } @@ -1304,62 +1304,62 @@ for (int i = 0; i < sources->n; i++) { - float Xo, Yo, Rmaj, Rmin, cs, sn; - - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - pmMoments *moments = source->moments; - if (0) { - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - Xo = moments->Mx; - Yo = moments->My; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - Rmaj = 2.0*axes.major; - Rmin = 2.0*axes.minor; - cs = cos(axes.theta); - sn = sin(axes.theta); - } else { - Rmaj = Rmin = 5.0; - cs = 1.0; - sn = 0.0; - Xo = source->peak->xf; - Yo = source->peak->yf; - } - - unsigned short int flagMask = 0x01; - for (int j = 0; j < 8; j++) { - if (source->mode & flagMask) { - overlayE[NoverlayE].type = KII_OVERLAY_LINE; - overlayE[NoverlayE].x = Xo; - overlayE[NoverlayE].y = Yo; - - float phi = j*M_PI/4.0; - overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; - overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; - overlayE[NoverlayE].angle = 0; - overlayE[NoverlayE].text = NULL; - NoverlayE ++; - CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100); - } - flagMask <<= 1; - - if (source->mode & flagMask) { - overlayO[NoverlayO].type = KII_OVERLAY_LINE; - overlayO[NoverlayO].x = Xo + 1; - overlayO[NoverlayO].y = Yo; - - float phi = j*M_PI/4.0; - overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; - overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; - overlayO[NoverlayO].angle = 0; - overlayO[NoverlayO].text = NULL; - NoverlayO ++; - CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100); - } - flagMask <<= 1; - } + float Xo, Yo, Rmaj, Rmin, cs, sn; + + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + pmMoments *moments = source->moments; + if (0) { + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + Xo = moments->Mx; + Yo = moments->My; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + Rmaj = 2.0*axes.major; + Rmin = 2.0*axes.minor; + cs = cos(axes.theta); + sn = sin(axes.theta); + } else { + Rmaj = Rmin = 5.0; + cs = 1.0; + sn = 0.0; + Xo = source->peak->xf; + Yo = source->peak->yf; + } + + unsigned short int flagMask = 0x01; + for (int j = 0; j < 8; j++) { + if (source->mode & flagMask) { + overlayE[NoverlayE].type = KII_OVERLAY_LINE; + overlayE[NoverlayE].x = Xo; + overlayE[NoverlayE].y = Yo; + + float phi = j*M_PI/4.0; + overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; + overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; + overlayE[NoverlayE].angle = 0; + overlayE[NoverlayE].text = NULL; + NoverlayE ++; + CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100); + } + flagMask <<= 1; + + if (source->mode & flagMask) { + overlayO[NoverlayO].type = KII_OVERLAY_LINE; + overlayO[NoverlayO].x = Xo + 1; + overlayO[NoverlayO].y = Yo; + + float phi = j*M_PI/4.0; + overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; + overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; + overlayO[NoverlayO].angle = 0; + overlayO[NoverlayO].text = NULL; + NoverlayO ++; + CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100); + } + flagMask <<= 1; + } } @@ -1376,5 +1376,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } @@ -1390,6 +1390,6 @@ if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; + fprintf (stderr, "kapa not opened, skipping\n"); + return false; } @@ -1402,19 +1402,19 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue; - - overlay[Noverlay].type = KII_OVERLAY_BOX; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - - overlay[Noverlay].dx = 4; - overlay[Noverlay].dy = 4; - overlay[Noverlay].angle = 0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + + overlay[Noverlay].dx = 4; + overlay[Noverlay].dy = 4; + overlay[Noverlay].angle = 0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); } KiiLoadOverlay (kapa, overlay, Noverlay, "red"); @@ -1424,20 +1424,20 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - // mark EXTs with yellow circles - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - - overlay[Noverlay].dx = 10; - overlay[Noverlay].dy = 10; - overlay[Noverlay].angle = 0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + // mark EXTs with yellow circles + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + + overlay[Noverlay].dx = 10; + overlay[Noverlay].dy = 10; + overlay[Noverlay].angle = 0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); } @@ -1453,5 +1453,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } @@ -1468,10 +1468,10 @@ if (kapa3 == -1) { kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } KapaClearPlots (kapa3); @@ -1500,7 +1500,7 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - if (!source) continue; + if (!source) continue; if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!isfinite (source->crNsigma)) continue; + if (!isfinite (source->crNsigma)) continue; x->data.F32[n] = -2.5*log10(source->peak->flux); @@ -1555,7 +1555,7 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - if (!source) continue; + if (!source) continue; if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!isfinite (source->extNsigma)) continue; + if (!isfinite (source->extNsigma)) continue; x->data.F32[n] = -2.5*log10(source->peak->flux); @@ -1597,5 +1597,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -1608,10 +1608,10 @@ if (kapa == -1) { kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } psphotVisualScaleImage (kapa, readout->image, "resid", 1); @@ -1622,5 +1622,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -1635,10 +1635,10 @@ if (kapa3 == -1) { kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - isVisual = false; - return false; - } - } + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } KapaClearPlots (kapa3); @@ -1657,8 +1657,8 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - if (!source) continue; + if (!source) continue; if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!isfinite (source->apMag)) continue; - if (!isfinite (source->psfMag)) continue; + if (!isfinite (source->apMag)) continue; + if (!isfinite (source->psfMag)) continue; x->data.F32[n] = source->psfMag; @@ -1705,5 +1705,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true;