Index: /trunk/psModules/src/objects/pmSource.c =================================================================== --- /trunk/psModules/src/objects/pmSource.c (revision 33878) +++ /trunk/psModules/src/objects/pmSource.c (revision 33879) @@ -39,4 +39,5 @@ #include "pmSourceExtendedPars.h" #include "pmSourceDiffStats.h" +#include "pmSourcePhotometry.h" #include "pmSource.h" @@ -159,4 +160,5 @@ source->radialAper = NULL; source->parent = NULL; + source->tmpPtr = NULL; source->imageID = -1; source->nFrames = 0; @@ -1157,4 +1159,92 @@ PAR[PM_PAR_SYY] = oldshape.sy; PAR[PM_PAR_SXY] = oldshape.sxy; + + return true; +} + +bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy) +{ +// assert (mode & PM_MODEL_OP_NOISE); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(target, false); + + if (add) { + psTrace ("psphot", 3, "adding smoothed object at %f,%f\n", source->peak->xf, source->peak->yf); + } else { + psTrace ("psphot", 3, "removing smooted object at %f,%f\n", source->peak->xf, source->peak->yf); + } + + pmModel *model = pmSourceGetModel(NULL, source); + if (!model) { + return false; + } + pmSourceSmoothOpModel (model, source, mode, target, sigma, add, maskVal, dx, dy); + + return true; +} + +bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy) +{ + bool status; + psEllipseShape oldshape; + psEllipseShape newshape; + psEllipseAxes axes; + + if (add) { + psTrace ("psphot", 4, "adding smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + } else { + psTrace ("psphot", 4, "removing smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + } + + psF32 *PAR = model->params->data.F32; + + // Isn't this hanging around somewhere? + float oldFlux = NAN; + if (!pmSourcePhotometryModel (NULL, &oldFlux, model)) return false; + + // save original values + float oldI0 = PAR[PM_PAR_I0]; + float oldsxx = PAR[PM_PAR_SXX]; + float oldsyy = PAR[PM_PAR_SYY]; + float oldsxy = PAR[PM_PAR_SXY]; + + // Since we are going to scale the flux correctly we need to get our + // factors of sqrt(2) right + oldshape.sx = oldsxx / M_SQRT2; + oldshape.sy = oldsyy / M_SQRT2; + oldshape.sxy = oldsxy; + + // XXX can this be done more intelligently? + if (oldI0 == 0.0) return false; + if (!isfinite(oldI0)) return false; + + // increase size and height of source + axes = psEllipseShapeToAxes (oldshape, 20.0); + axes.major = sqrt(PS_SQR(axes.major) + PS_SQR(sigma)); + axes.minor = sqrt(PS_SQR(axes.minor) + PS_SQR(sigma)); + newshape = psEllipseAxesToShape (axes); + PAR[PM_PAR_SXX] = newshape.sx * M_SQRT2; + PAR[PM_PAR_SYY] = newshape.sy * M_SQRT2; + PAR[PM_PAR_SXY] = newshape.sxy; + + PAR[PM_PAR_I0] = 1.0; + + float newFlux; + if (!pmSourcePhotometryModel (NULL, &newFlux, model)) return false; + + PAR[PM_PAR_I0] = oldFlux / newFlux; + + if (add) { + status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy); + } else { + status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy); + } + + // restore original values + PAR[PM_PAR_I0] = oldI0; + PAR[PM_PAR_SXX] = oldsxx; + PAR[PM_PAR_SYY] = oldsyy; + PAR[PM_PAR_SXY] = oldsxy; return true; Index: /trunk/psModules/src/objects/pmSource.h =================================================================== --- /trunk/psModules/src/objects/pmSource.h (revision 33878) +++ /trunk/psModules/src/objects/pmSource.h (revision 33879) @@ -117,4 +117,5 @@ psArray *radialAper; ///< radial flux in circular apertures pmSource *parent; ///< reference to the master source from which this is derived + psPtr *tmpPtr; ///< pointer that may be used to store data in a particular module. e.g. psphotKronIterate. int imageID; psU16 nFrames; @@ -302,4 +303,7 @@ bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy); +bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy); +bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy); + bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy); bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal); Index: /trunk/psphot/src/psphotKronIterate.c =================================================================== --- /trunk/psphot/src/psphotKronIterate.c (revision 33878) +++ /trunk/psphot/src/psphotKronIterate.c (revision 33879) @@ -1,6 +1,3 @@ # 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); @@ -57,101 +54,4 @@ bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max); -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) { @@ -165,9 +65,4 @@ psTimerStart ("psphot.kron"); -#ifdef USE_SE_KR - if (!pass1) { - build_kr_list(sources); - } -#endif // determine the number of allowed threads @@ -237,7 +132,4 @@ // 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++) { @@ -247,15 +139,5 @@ // (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 @@ -291,26 +173,4 @@ } - -#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 @@ -377,7 +237,4 @@ 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); @@ -465,18 +322,4 @@ 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, smoothedPixels); @@ -606,9 +449,4 @@ 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; @@ -621,10 +459,4 @@ 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?)