IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32205


Ignore:
Timestamp:
Aug 29, 2011, 2:59:09 AM (15 years ago)
Author:
eugene
Message:

find the footprint saddles; iterative search for kron (in masked); add iterative search for kron (downweighting); add footprint based mask

Location:
branches/eam_branches/ipp-20110710/psphot/src
Files:
2 added
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110710/psphot/src/Makefile.am

    r32154 r32205  
    130130        psphotVisual.c                 \
    131131        psphotCullPeaks.c              \
     132        psphotFootprintSaddles.c       \
    132133        psphotVersion.c                \
    133134        psphotModelGroupInit.c         \
    … …  
    183184        psphotRadialPlot.c             \
    184185        psphotKronMasked.c             \
     186        psphotKronIterate.c            \
    185187        psphotDeblendSatstars.c        \
    186188        psphotMosaicSubimage.c         \

  • branches/eam_branches/ipp-20110710/psphot/src/psphot.h

    r32157 r32205  
    482482bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
    483483
     484bool psphotFootprintSaddles(pmReadout *readout, psArray *footprints);
     485bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal);
     486
     487bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
     488bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
     489
    484490#endif

  • branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c

    r32157 r32205  
    319319        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal);
    320320
    321         // Recalculate the source moments using the larger extended-source moments radius.  At
    322         // this stage, skip Gaussian windowing, and do not clip pixels by S/N.  Save the
    323         // psf-based moments for output
    324         psfMoments = *source->moments;
    325         if (!pmSourceMoments (source, windowRadius, 0.0, 0.0, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) {
    326             // subtract the best fit from the object, leave local sky
    327             fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
    328             *source->moments = psfMoments;
    329 
    330             // probably need to bump up the size for subtraction
    331             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    332             // XXX raise an error flag of some kind
    333             continue;
    334         }
     321        // UPDATE : we have changed the moments calculation.  There is now an iteration within
     322        // psphotKronMasked to determine moments appropriate for a larger object.  The values
     323        // Mrf, KronFlux, and KronFluxErr are calculated for the iterated radius.  The other
     324        // values are left at the psf-based values.
    335325
    336326        // allocate the array to store the model fits
    … …  
    338328            source->modelFits = psArrayAllocEmpty (models->list->n);
    339329        }
    340         // extFitPars are the non-parametric data measured for this model fit (moments, kron, etc)
    341         if (source->extFitPars == NULL) {
    342             source->extFitPars = psArrayAllocEmpty (models->list->n);
    343         }
     330
     331        // XXX deprecate?
     332        // XXX // extFitPars are the non-parametric data measured for this model fit (moments, kron, etc)
     333        // XXX if (source->extFitPars == NULL) {
     334        // XXX     source->extFitPars = psArrayAllocEmpty (models->list->n);
     335        // XXX }
    344336
    345337        // loop here over the models chosen for each source (exclude by S/N)
    … …  
    405397          }
    406398
    407           pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
    408           extFitPars->Mxx = source->moments->Mxx;
    409           extFitPars->Mxy = source->moments->Mxy;
    410           extFitPars->Myy = source->moments->Myy;
    411           extFitPars->Mrf = source->moments->Mrf;
    412           extFitPars->Mrh = source->moments->Mrh;
    413           extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
     399          // XXX deprecate?
     400          // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
     401          // XXX extFitPars->Mxx = source->moments->Mxx;
     402          // XXX extFitPars->Mxy = source->moments->Mxy;
     403          // XXX extFitPars->Myy = source->moments->Myy;
     404          // XXX extFitPars->Mrf = source->moments->Mrf;
     405          // XXX extFitPars->Mrh = source->moments->Mrh;
     406          // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
    414407         
    415408          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
    416           extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
    417 
    418           psArrayAdd (source->extFitPars, 4, extFitPars);
    419           psFree (extFitPars);
     409          // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
     410
     411          // psArrayAdd (source->extFitPars, 4, extFitPars);
     412          // psFree (extFitPars);
    420413
    421414          // test for fit quality / result

  • branches/eam_branches/ipp-20110710/psphot/src/psphotFindFootprints.c

    r31154 r32205  
    5858    psphotCullPeaks(readout, significance, recipe, detections->footprints);
    5959    detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
     60
     61    psphotFootprintSaddles (readout, detections->footprints);
     62
     63    int nSaddle = 0;
     64    for (int i = 0; i < detections->peaks->n; i++) {
     65        pmPeak *peak = detections->peaks->data[i];
     66       
     67        if (peak->saddlePoints) nSaddle += peak->saddlePoints->n;
     68    }
     69    fprintf (stderr, "%d saddle points for %ld peaks\n", nSaddle, detections->peaks->n);
     70
    6071    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
    6172

  • branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c

    r32157 r32205  
    11# include "psphotInternal.h"
    2 
     2# ifndef ROUND
     3# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
     4# endif
     5
     6bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal);
    37bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
    48
    … …  
    4549}
    4650
     51int psphotKapaChannel (int channel);
     52bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
     53bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
     54
    4755bool psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
    4856
    … …  
    7280    }
    7381
     82    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
     83
    7484    // bit-masks to test for good/bad pixels
    7585    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    … …  
    8393    maskVal |= markVal;
    8494
     95    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
     96
     97    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
     98
     99    // XXX tmp visualization
     100    // int kapa = psphotKapaChannel (1);
     101
    85102    // generate the mask image: increment counter for every source overlapping the pixel
    86     psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
     103    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
    87104    psImageInit (kronMask, 0);
    88     int Nx = kronMask->numCols;
    89     int Ny = kronMask->numRows;
    90105    for (int i = 0; i < sources->n; i++) {
    91106
    … …  
    96111        if (!source->moments) continue;
    97112
     113        // replace object in image
     114        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     115            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     116        }
     117
     118        // save the PSF-based values
     119        // source->moments->KronFluxPSF    = source->moments->KronFlux;
     120        // source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
     121        // source->moments->KronRadiusPSF  = source->moments->Mrf;
     122
     123        // iterate to the window radius
    98124        float windowRadius = RADIUS;
    99         if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
    100             windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
    101         }
    102 
    103         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    104         pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    105 
    106         pmMoments psfMoments = *source->moments;
    107 
    108         // replace object in image
    109         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    110             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    111         }
    112 
    113         psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
    114 
    115         // re-subtract the object, leave local sky
    116         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    117 
    118         continue;
    119 
    120         // XXX skip this code
    121         // generate the pixel masks
    122         // int Xo = source->moments->Mx;
    123         // int Yo = source->moments->My;
    124         float dX = source->moments->Mx - source->peak->xf;
    125         float dY = source->moments->My - source->peak->yf;
    126         float dR = hypot(dX, dY);
    127        
    128         float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    129         float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    130        
    131         int Kr = 2.5*source->moments->Mrf;
    132         int Kr2 = Kr*Kr;
    133        
    134         for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
    135             if (iy < 0) continue;
    136             if (iy >= Ny) continue;
    137             for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
    138                 if (ix < 0) continue;
    139                 if (ix >= Nx) continue;
    140 
    141                 if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
    142                 if (kronMask->data.U8[iy][ix] < 0xff) {
    143                     kronMask->data.U8[iy][ix] ++;
    144                 }
    145             }
    146         }
    147     }           
    148 
     125        for (int j = 0; j < 4; j++) {
     126            // XXX use some S/N criterion to limit?
     127            // if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
     128            // windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
     129            // }
     130
     131            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
     132            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     133
     134            // mask the pixels not contained by the footprint
     135            // if (psphotMaskFootprint (readout, source, markVal)) {
     136            //  // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
     137            //  // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
     138            //  // fprintf (stderr, "masked\n");
     139            // }           
     140
     141            // mask the pixels associated with near neighbors
     142            if (psphotKronMask (source, kronMask, markVal)) {
     143                // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
     144                // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
     145                // fprintf (stderr, "masked\n");
     146            }       
     147            // this function populates moments->Mrf,KronFlux,KronFluxErr
     148            // XXX what about KronInner, KronOuter, etc?
     149            psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
     150            windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
     151            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     152        }
     153        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
     154        float kmag = -2.5*log10(source->moments->KronFlux);
     155        if (source->psfMag - kmag > 0.25) {
     156            // psphotVisualShowMask (kapa, source->maskObj, "mask", 1);
     157            // psphotVisualRangeImage (kapa, source->pixels, "image", 0, -50, 300);
     158            // fprintf (stderr, "psf: %f, kron: %f, dmag: %f\n", source->psfMag, kmag, source->psfMag - kmag);
     159            // fprintf (stderr, "continue\n");
     160        }
     161        // re-subtract the object, leave local sky
     162        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     163    }
     164
     165    // psphotSaveImage (NULL, readout->mask, "kron.masked.fits");
    149166    // psphotSaveImage (NULL, kronMask, "kronmask.fits");
     167
    150168    psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    151169
    … …  
    153171    return true;
    154172}
     173
     174# define WEIGHTED 0
    155175
    156176bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
    … …  
    161181    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    162182
    163     psF32 Sum = 0.0;
    164     psF32 Var = 0.0;
    165183    psF32 R2 = PS_SQR(radius);
     184
     185# if (WEIGHTED)
     186    float rsigma2 = 0.5 / PS_SQR(radius / 4.0); // use a Gaussian window with sigma = R_window / 2
     187# endif
    166188
    167189    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    179201    psF32 RS = 0.0;
    180202
    181 # if (1)
     203    // the peak position is less accurate but less subject to extreme deviations
    182204    float dX = source->moments->Mx - source->peak->xf;
    183205    float dY = source->moments->My - source->peak->yf;
    … …  
    185207    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    186208    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    187 # else
    188     float Xo = source->moments->Mx;
    189     float Yo = source->moments->My;
    190 # endif
    191209
    192210    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
    … …  
    221239            if (r2 > R2) continue;
    222240
     241# if (WEIGHTED)
     242            float z = r2 * rsigma2;
     243            assert (z >= 0.0);
     244            float weight  = exp(-z);
     245            psF32 pDiff = *vPix * weight;
     246# else
    223247            psF32 pDiff = *vPix;
    224 
    225             Sum += pDiff;
    226 
    227             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
     248# endif
     249
     250            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    228251            psF32 rf = pDiff * sqrt(r2);
    229252            psF32 rs = pDiff;
    … …  
    245268
    246269    int nKronPix = 0;
    247     Sum = Var = 0.0;
     270    float Sum = 0.0;
     271    float Var = 0.0;
    248272
    249273    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    274298            if (r2 > radKron2) continue;
    275299
     300# if (WEIGHTED)
     301            float z = r2 * rsigma2;
     302            assert (z >= 0.0);
     303            float weight  = exp(-z);
     304            psF32 pDiff = *vPix * weight;
     305            psF32 wDiff = *vWgt * weight;
     306# else
    276307            psF32 pDiff = *vPix;
    277308            psF32 wDiff = *vWgt;
     309# endif
    278310
    279311            Sum += pDiff;
    … …  
    283315    }
    284316
    285     source->moments->Mrh = Mrf;
    286 
    287     // XXX for a test, save the old values here:
    288     source->moments->KronFouter    = source->moments->KronFlux;
    289     source->moments->KronCoreErr = source->moments->KronFluxErr;
    290 
    291     source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
    292     source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
     317    source->moments->Mrf = Mrf;
     318    source->moments->KronFlux    = Sum;
     319    source->moments->KronFluxErr = sqrt(Var);
     320    // source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
     321    // source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
    293322
    294323    return true;
    295324}
     325
     326enum {Y_P, Y_M, X_P, X_M};
     327
     328bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal) {
     329
     330    int X, Y, e, dXs, dYs;
     331
     332    if (!source->peak) return false;
     333    if (!source->peak->saddlePoints) return false;
     334    if (!source->peak->saddlePoints->n) return false;
     335
     336    // return true;
     337
     338    int xMax = source->maskObj->numCols;
     339    int yMax = source->maskObj->numRows;
     340    // int xOff = source->maskObj->col0;
     341    // int yOff = source->maskObj->row0;
     342
     343    int Xo = ROUND(source->peak->xf) - source->maskObj->col0;
     344    int Yo = ROUND(source->peak->yf) - source->maskObj->row0;
     345
     346    for (int i = 0; i < source->peak->saddlePoints->n; i++) {
     347
     348        psVector *saddlePoint = source->peak->saddlePoints->data[i];
     349        int Xs = ROUND(saddlePoint->data.S32[0]) - source->maskObj->col0;
     350        int Ys = ROUND(saddlePoint->data.S32[1]) - source->maskObj->row0;
     351       
     352        // fprintf (stderr, "mask %d,%d @ %d,%d (%d,%d)\n", Xo, Yo, Xs, Ys, saddlePoint->data.S32[0], saddlePoint->data.S32[1]);
     353
     354        // We want to mask the pixels between the edge of the image and the line perpendicular
     355        // to the connecting line.  Call dX,dY = (Xs-Xo,Ys-Yo).  There are 4 cases:
     356        // +y : dY > 0, fabs(dY) > fabs(dX)
     357        // -y : dY < 0, fabs(dY) > fabs(dX)
     358        // +x : dX > 0, fabs(dX) > fabs(dY)
     359        // -x : dX < 0, fabs(dX) > fabs(dY)
     360
     361        int dX = Xs - Xo;
     362        int dY = Ys - Yo;
     363
     364        // +y : dY > 0, fabs(dY) > fabs(dX)
     365        // -y : dY < 0, fabs(dY) > fabs(dX)
     366        if (fabs(dY) > fabs(dX)) {
     367            // points to right of saddle
     368            Y = Ys;
     369            e = 0;
     370            dXs = (dY > 0) ? +dY : -dY;  // this forces dXs > 0
     371            dYs = (dY > 0) ? -dX : +dX;
     372            for (X = Xs; X < xMax; X++) {
     373                if (dY > 0) {
     374                    for (int Ym = Y; Ym < yMax; Ym++) {
     375                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     376                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     377                    }
     378                } else {
     379                    for (int Ym = Y; Ym >= 0; Ym--) {
     380                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     381                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     382                    }
     383                }
     384                e += dYs;
     385                int e2 = 2 * e;
     386                if (e2 > dXs) {
     387                    Y ++;
     388                    e -= dXs;
     389                }
     390                if (e2 < -dXs) {
     391                    Y --;
     392                    e += dXs;
     393                }
     394                if (Y >= yMax) break;
     395                if (Y < 0) break;
     396            }
     397            // points to left of saddle
     398            Y = Ys;
     399            e = 0;
     400            for (X = Xs; X >= 0; X--) {
     401                if (dY > 0) {
     402                    for (int Ym = Y; Ym < yMax; Ym++) {
     403                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     404                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     405                    }
     406                } else {
     407                    for (int Ym = Y; Ym >= 0; Ym--) {
     408                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     409                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     410                    }
     411                }
     412                e -= dYs;
     413                int e2 = 2 * e;
     414                if (e2 > dXs) {
     415                    Y ++;
     416                    e -= dXs;
     417                }
     418                if (e2 < -dXs) {
     419                    Y --;
     420                    e += dXs;
     421                }
     422                if (Y >= yMax) break;
     423                if (Y < 0) break;
     424            }
     425        } else {
     426            // points to right of saddle
     427            X = Xs;
     428            e = 0;
     429            dXs = (dX > 0) ? -dY : +dY;
     430            dYs = (dX > 0) ? +dX : -dX;  // this forces dYs > 0
     431            for (Y = Ys; Y < yMax; Y++) {
     432                if (dX > 0) {
     433                    for (int Xm = X; Xm < xMax; Xm++) {
     434                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     435                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     436                    }
     437                } else {
     438                    for (int Xm = X; Xm >= 0; Xm--) {
     439                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     440                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     441                    }
     442                }
     443                e += dXs;
     444                int e2 = 2 * e;
     445                if (e2 > dYs) {
     446                    X ++;
     447                    e -= dYs;
     448                }
     449                if (e2 < -dYs) {
     450                    X --;
     451                    e += dYs;
     452                }
     453                if (X >= yMax) break;
     454                if (X < 0) break;
     455            }
     456            // points to left of saddle
     457            X = Xs;
     458            e = 0;
     459            for (Y = Ys; Y >= 0; Y--) {
     460                if (dX > 0) {
     461                    for (int Xm = X; Xm < xMax; Xm++) {
     462                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     463                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     464                    }
     465                } else {
     466                    for (int Xm = X; Xm >= 0; Xm--) {
     467                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     468                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     469                    }
     470                }
     471                e -= dXs;
     472                int e2 = 2 * e;
     473                if (e2 > dYs) {
     474                    X ++;
     475                    e -= dYs;
     476                }
     477                if (e2 < -dYs) {
     478                    X --;
     479                    e += dYs;
     480                }
     481                if (X >= yMax) break;
     482                if (X < 0) break;
     483            }
     484        }
     485    }
     486    return true;
     487}

  • branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c

    r32154 r32205  
    8484
    8585    // define the source of interest
    86     float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     86   float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
    8787    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
    8888    if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
    … …  
    102102    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    103103    if (!status) psAbort("pmSourceLocalSky error");
     104
     105    {
     106        // XXX I want to test an iterative aperture for brighter sources
     107        float radius = mRADIUS;
     108        for (int i = 0; i < 10; i++) {
     109
     110            // get the source moments
     111            status = pmSourceMoments (source, radius, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
     112            if (!status) psAbort("psSourceMoments error");
     113
     114            float oldRadius = radius;
     115            radius = source->moments->Mrf * RADIUS;
     116           
     117            fprintf (stderr, "%d %f  %f  %f\n", i, oldRadius, radius, source->moments->Mrf);
     118        }
     119        exit (0);
     120    }
    104121
    105122    // get the source moments

  • branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c

    r32157 r32205  
    241241}
    242242
     243// call this function whenever you (re)-define the EXT model
     244bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
     245
     246    psAssert (source, "source not defined??");
     247    psAssert (source->peak, "peak not defined??");
     248
     249    pmPeak *peak = source->peak;
     250
     251    // set the radius based on the footprint:
     252    if (!peak->footprint) return false;
     253    pmFootprint *footprint = peak->footprint;
     254    if (!footprint->spans) return false;
     255    if (footprint->spans->n < 1) return false;
     256
     257    int Xo = source->maskObj->col0;
     258    int Yo = source->maskObj->row0;
     259
     260    // mark all pixels
     261    for (int j = 0; j < source->maskObj->numRows; j++) {
     262        for (int i = 0; i < source->maskObj->numCols; i++) {
     263            source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= markVal;
     264        }
     265    }
     266
     267    psImageMaskType clearVal = PS_NOT_IMAGE_MASK(markVal);
     268
     269    for (int j = 0; j < footprint->spans->n; j++) {
     270        pmSpan *span = footprint->spans->data[j];
     271
     272        // mask the rows before and after each span
     273        int minX = span->x0 - Xo - 2;
     274        int maxX = span->x1 - Xo + 2;
     275        int myY = span->y - Yo;
     276
     277        // unmark pixels inside the footprint
     278        for (int i = minX; i <= maxX; i++) {
     279            source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[myY][i] &= clearVal;
     280        }
     281    }
     282    return true;
     283}
     284
    243285// alternative EXT radius based on model guess (for use without footprints)
    244286bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {

  • branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c

    r32157 r32205  
    205205   
    206206    // but this is chosen above to be appropriate for the PSF objects (not galaxies)
    207     psphotKronMasked(config, view, filerule);
     207    // psphotKronMasked(config, view, filerule);
     208    psphotKronIterate(config, view, filerule);
    208209
    209210    // identify CRs and extended sources (only unmeasured sources are measured)
    … …  
    321322    // XXX re-measure the kron mags with models subtracted
    322323    // Note that this uses the PSF_MOMENTS_RADIUS as a window
    323     psphotKronMasked(config, view, filerule);
     324    // psphotKronMasked(config, view, filerule);
    324325
    325326    // measure source size for the remaining sources

  • branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c

    r31990 r32205  
    542542   
    543543    // get the PSFratio and the ASPratio : we use these to choose the model class
    544     float PSFratio = momentAxes.major / (2.35 * psfAxes.major) ;
     544    float PSFratio = source->moments->Mrf / (2.35 * psfAxes.major) ;
    545545    float ASPratio = momentAxes.minor / momentAxes.major;
    546546

  • branches/eam_branches/ipp-20110710/psphot/src/psphotSourceSize.c

    r31452 r32205  
    209209        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    210210
    211         float kMag = -2.5*log10(source->moments->KronFlux);
     211        float kMag = -2.5*log10(source->moments->KronFluxPSF);
    212212        float dMag = source->psfMag - kMag;
    213213
    … …  
    334334        psF32 Mxy = source->moments->Mxy;
    335335
    336         float KronMag = -2.5*log10(source->moments->KronFlux);
     336        float KronMag = -2.5*log10(source->moments->KronFluxPSF);
    337337
    338338        float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    … …  
    509509        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    510510
    511         float kMag = -2.5*log10(source->moments->KronFlux);
     511        float kMag = -2.5*log10(source->moments->KronFluxPSF);
    512512        float dMag = source->psfMag - kMag;
    513513

  • branches/eam_branches/ipp-20110710/psphot/src/psphotSourceStats.c

    r31673 r32205  
    554554        // determine the PSF parameters from the source moment values
    555555        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
    556         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
     556        psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    557557
    558558        Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
Note: See TracChangeset for help on using the changeset viewer.