Index: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 33837) +++ trunk/psphot/src/psphotKronIterate.c (revision 33877) @@ -1,6 +1,10 @@ # include "psphotInternal.h" +#define NO_USE_SE_KR 1 +#define NO_SAVE_IMAGES 1 + bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow); -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, bool smooth); +bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels); + bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule) @@ -41,5 +45,5 @@ // psAssert (psf, "missing psf?"); - if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) { + if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -53,5 +57,102 @@ bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max); -bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool pass1 = true; +// This is defined in pmSource.c +extern psString psphot_outroot; + +#ifdef USE_SE_KR +static int nKR = 0; +// input list indexed by list file's source->id +static float kr_sextractor[10000]; +static int used[10000]; + +int numMatched = 0; + +#define SE_KR_USED 0x800 + +static void find_se_kr(pmSource *source) { + int sid = source->id; + if (isfinite(kr_sextractor[sid])) { + source->moments->Mrf = kr_sextractor[sid] / 2.5; + source->tmpFlags |= SE_KR_USED; + numMatched++; + } +} + +static void getxy(pmSource *source, float *pX, float *pY) { + float x, y; + if (source->modelPSF) { + psF32 *PAR = source->modelPSF->params->data.F32; + x = PAR[PM_PAR_XPOS]; + y = PAR[PM_PAR_YPOS]; + } else if (pmSourcePositionUseMoments(source)) { + x = source->moments->Mx; + y = source->moments->My; + } else { + x = source->peak->x; + y = source->peak->y; + } + *pX = x; + *pY = y; +} + +static void build_kr_list(psArray *sources) { + // write list of source->id X Y to a file + pid_t pid = getpid(); + char filename[40]; + sprintf(filename, "sources.%d.list", pid); + FILE *out = fopen(filename, "w"); + for (int i=0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + float x, y; + getxy(source, &x, &y); + fprintf(out, "%6.3f %6.3f %5d\n", x, y, source->id); + } + fclose(out); + // run gcompare on that file + + char matchedfilename[20]; + sprintf(matchedfilename, "sources.%d.matched", pid); + char command[80]; + sprintf(command, "rungcompare %s %s", filename, matchedfilename); + int rc = system(command); + if (rc) { + fprintf(stderr, "failed to read matched list %d %d\n", rc, rc >> 8); + return; + } + + unlink(filename); + + for (int i = 0; i < 10000; i++) { + used[i] = -1; + kr_sextractor[i] = NAN; + } + int sid; + int ipp_idet; + float kr; + // read the results to build the kr_sextractor + FILE *in = fopen(matchedfilename, "r"); + while (fscanf(in, "%d %f %d", &sid, &kr, &ipp_idet) > 0) { + // take first match for each sid ... + if (!isfinite(kr_sextractor[sid])) { + // .. and for each ipp_idet Since the lists are sorted in + // order of SN this makes it more likely to get the right match + if (used[ipp_idet] == -1) { + kr_sextractor[sid] = kr; + used[ipp_idet] = sid; + nKR++; + } else { + fprintf(stderr, "match for %d %d is already used for %d\n", + sid, ipp_idet, used[ipp_idet]); + } + } + } + fclose(in); + unlink(matchedfilename); +} + +#endif // USE_SE_KR + +bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) { bool status = false; @@ -64,4 +165,10 @@ psTimerStart ("psphot.kron"); +#ifdef USE_SE_KR + if (!pass1) { + build_kr_list(sources); + } +#endif + // determine the number of allowed threads int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads @@ -97,4 +204,12 @@ if (!status) { KRON_SMOOTH = false; + } + float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA"); + if (!status) { + KRON_SMOOTH_SIGMA = 1.7; + } + float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA"); + if (!status) { + KRON_SMOOTH_NSIGMA = 2; } @@ -122,4 +237,7 @@ // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image +#ifdef USE_SE_KR + int previous = numMatched; +#endif for (int i = 0; i < sources->n; i++) { @@ -129,5 +247,70 @@ // (this skips really bad sources (no peak, no moments, DEFECT) psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW); - } +#ifdef USE_SE_KR + if (!pass1) { + find_se_kr(source); + } +#endif + } +#ifdef USE_SE_KR + if (!pass1) { + fprintf(stdout, "Matched %d sources previous %d\n", numMatched, previous); + } +#endif + + // We measure the Kron Radius on a smoothed copy of the readout image + psImage *smoothedImage = NULL; + if (KRON_SMOOTH) { + // Build the smoothed source image + // Replace the subtracted sources + psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false); + // Copy the image and smooth + psTimerStart ("psphot.kron.smooth"); + smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32); + psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA); + psLogMsg ("psphot.kron", PS_LOG_WARN, "Smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth")); + + // remove the sources + psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false ); + // Now subtract smooth versions of the sources from the smoothed image + psTimerStart ("psphot.kron.smooth.sources"); + for (int i=0; i< sources->n; i++) { + pmSource *source = sources->data[i]; + // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + // cache copy of smoothedPixels in the source + // tmpPtr is for use by a single "module" and must be null otherwise + psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!"); + + psImage *smoothedPixels = psImageSubset(smoothedImage, source->region); + source->tmpPtr = (psPtr) smoothedPixels; + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); + } + } + psLogMsg ("psphot.kron", PS_LOG_WARN, " removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources")); + + } + +#ifdef SAVE_IMAGES + psphot_outroot = psMetadataLookupStr(&status, config->arguments, "OUTPUT"); + { + // Save the background subtracted image + psString fn = NULL; + psStringAppend(&fn, "%s.p%d.src.sub.fits", psphot_outroot, pass1 ? 1 : 2); + psphotSaveImage(0, readout->image, fn); + psFree(fn); + fn = NULL; + if (KRON_SMOOTH) { + psStringAppend(&fn, "%s.p%d.src.sub.sm.fits", psphot_outroot, pass1 ? 1 : 2); + psphotSaveImage(0, smoothedImage, fn); + psFree(fn); + } + if (KRON_APPLY_WINDOW) { + psStringAppend(&fn, "%s.p%d.window.fits", psphot_outroot, pass1 ? 1 : 2); + psphotSaveImage(0, kronWindow, fn); + psFree(fn); + } + } +#endif // threaded measurement of the source magnitudes @@ -157,8 +340,9 @@ PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32); PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_SMOOTH, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA, PS_TYPE_F32); + psArrayAdd(job->args, 1, smoothedImage); // set this to 0 to run without threading -# if (0) +# if (1) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); @@ -192,4 +376,8 @@ psFree (cellGroups); psFree (kronWindow); + psFree (smoothedImage); + if (pass1) { + pass1 = false; + } psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); @@ -209,6 +397,8 @@ bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[8],S32); bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[9],S32); - bool KRON_SMOOTH = PS_SCALAR_VALUE(job->args->data[10],S32); - + float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[10],F32); + psImage *smoothedImage = job->args->data[11]; + + // psImage *smoothedPixels = NULL; for (int j = 0; j < KRON_ITERATIONS; j++) { for (int i = 0; i < sources->n; i++) { @@ -224,6 +414,14 @@ // replace object in image bool reSubtract = false; + psImage *smoothedPixels = NULL; if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + smoothedPixels = (psImage *) source->tmpPtr; + if (smoothedPixels) { + // psFree(source->tmpPtr); + // smoothedPixels = psImageSubset(smoothedImage, source->region); + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0); + // source->tmpPtr = smoothedPixels; + } reSubtract = true; } @@ -251,12 +449,36 @@ // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); psAssert (source->pixels, "WTF?"); + if (extend && smoothedPixels) { + psFree(source->tmpPtr); + smoothedPixels = psImageSubset(smoothedImage, source->region); + psAssert (smoothedPixels, "WTF?"); + source->tmpPtr = (psPtr) smoothedPixels ; + } + // clear the window function for this source based on the moments + // Note this function also applies cuts on the source and returns false if it + // does not meet the requirements for measuring the Kron Radius or Magnitude. + // Not that it performs that function even if KRON_APPLY_WINDOW is false if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) { +#ifdef SAVE_IMAGES +#ifdef SAVE_SOURCE_IMAGES + if (j == 1 && !pass1) { + char fn[80]; + if (smoothedPixels) { + sprintf(fn, "%s.s.%d.p%d.sm.fits", psphot_outroot, source->id, pass1 ? 1 : 2); + psphotSaveImage(0, smoothedPixels, fn); + } + sprintf(fn, "%s.s.%d.p%d.fits", psphot_outroot, source->id, pass1 ? 1 : 2); + psphotSaveImage(0, source->pixels, fn); + } +#endif +#endif + // this function populates moments->Mrf,KronFlux,KronFluxErr - psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, KRON_SMOOTH); + psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels); psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); @@ -267,13 +489,23 @@ // if we subtracted it above, re-subtract the object, leave local sky if (reSubtract) { - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (smoothedPixels) { + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); + if (j + 1 == KRON_ITERATIONS) { + // We're done with the smoothedPixels + psFree(source->tmpPtr); + source->tmpPtr = NULL; + smoothedPixels = NULL; + } + } } } } + // psFree(smoothedPixels); return true; } bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, -bool applyWeight, bool smooth) { + bool applyWeight, psImage *smoothedPixels) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -296,4 +528,5 @@ // Xn = SUM (x - xc)^n * (z - sky) + psF32 RF = 0.0; psF32 RS = 0.0; @@ -315,9 +548,6 @@ psF32 **vPix; - psImage *smoothedImage = NULL;; - if (smooth) { - smoothedImage = psImageCopy(NULL, source->pixels, PS_TYPE_F32); - psImageSmooth(smoothedImage, 1.7, 2); - vPix = smoothedImage->data.F32; + if (smoothedPixels) { + vPix = smoothedPixels->data.F32; } else { vPix = source->pixels->data.F32; @@ -376,4 +606,9 @@ psF32 rs = 0.5 * (fDiff1 + fDiff2); +#ifdef BE_SIMPLE_MINDED + rf = vPix[row][col] * sqrt(r2); + rs = vPix[row][col]; +#endif + RF += rf; RS += rs; @@ -386,4 +621,10 @@ Mrf = MIN (radius, Mrf); } +#ifdef USE_SE_KR + if (source->tmpFlags & SE_KR_USED) { + // we are using the sextratcor KR for this source set it back + Mrf = source->moments->Mrf; + } +#endif // Calculate the Kron magnitude (make this block optional?) @@ -399,4 +640,5 @@ // smoothed image above) vPix = source->pixels->data.F32; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -433,7 +675,4 @@ source->moments->KronFlux = Sum; source->moments->KronFluxErr = sqrt(Var); - if (smoothedImage) { - psFree(smoothedImage); - } return true; @@ -451,5 +690,5 @@ psAssert(kronWindow, "need a window"); - if (source->moments->Mrf < 0) return false; + if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false; if (!applyWindow) { @@ -478,5 +717,6 @@ } - float scale = fabs(0.5 * source->moments->Mrf / axes.major); + // Why this factor of 0.5 ? + float scale = 0.5 * source->moments->Mrf / axes.major; axes.major *= scale; axes.minor *= scale; @@ -489,7 +729,4 @@ } - float Sxx = PS_SQR(shape.sx); - float Syy = PS_SQR(shape.sy); - float Sxy = -1. * shape.sxy; float Smajor = axes.major; @@ -499,6 +736,8 @@ int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1); - float rMxx = 0.5 / Sxx; - float rMyy = 0.5 / Syy; + float rMxx = 0.5 / PS_SQR(shape.sx); + float rMyy = 0.5 / PS_SQR(shape.sy); + float Sxy = -1. * shape.sxy; // factor of -1 is included to match the previous window function + // implementation. XXX: Is this correct? for (int iy = minY; iy < maxY; iy++) {