IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2011, 1:32:31 PM (15 years ago)
Author:
eugene
Message:
  • clarified the output to show the stages and their timing more clearly
  • psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
  • psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
  • for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
  • extended source analysis (petrosians) is threaded
  • iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
  • the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
  • dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
  • fake sources inserted for the efficiency analysis are removed at the end of the analysis
  • radial apertures is now threaded
  • edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
  • new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
  • define the fit region and window region based on the 1st radial moments
  • use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
  • remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis
Location:
trunk/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        2222psphotMakePSF
        2323psphotStack
         24psphotModelTest

  • trunk/psphot/src/psphotKronMasked.c

    r31673 r32348  
    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
    98         // replace object in image
    99         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    100             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    101         }
    102 
    103         psphotKronMag (source, RADIUS, MIN_KRON_RADIUS, maskVal);
    104 
    105         // re-subtract the object, leave local sky
    106         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    107 
    108         continue;
    109 
    110         // XXX skip this code
    111         // generate the pixel masks
    112         // int Xo = source->moments->Mx;
    113         // int Yo = source->moments->My;
    114         float dX = source->moments->Mx - source->peak->xf;
    115         float dY = source->moments->My - source->peak->yf;
    116         float dR = hypot(dX, dY);
    117        
    118         float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    119         float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    120        
    121         int Kr = 2.5*source->moments->Mrf;
    122         int Kr2 = Kr*Kr;
    123        
    124         for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
    125             if (iy < 0) continue;
    126             if (iy >= Ny) continue;
    127             for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
    128                 if (ix < 0) continue;
    129                 if (ix >= Nx) continue;
    130 
    131                 if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
    132                 if (kronMask->data.U8[iy][ix] < 0xff) {
    133                     kronMask->data.U8[iy][ix] ++;
    134                 }
    135             }
    136         }
    137     }           
     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
     124        float windowRadius = RADIUS;
     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");
    138166    // psphotSaveImage (NULL, kronMask, "kronmask.fits");
     167
    139168    psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    140169
    … …  
    142171    return true;
    143172}
     173
     174# define WEIGHTED 0
    144175
    145176bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
    … …  
    150181    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    151182
    152     psF32 Sum = 0.0;
    153     psF32 Var = 0.0;
    154183    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
    155188
    156189    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    168201    psF32 RS = 0.0;
    169202
    170 # if (1)
     203    // the peak position is less accurate but less subject to extreme deviations
    171204    float dX = source->moments->Mx - source->peak->xf;
    172205    float dY = source->moments->My - source->peak->yf;
    … …  
    174207    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    175208    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    176 # else
    177     float Xo = source->moments->Mx;
    178     float Yo = source->moments->My;
    179 # endif
    180209
    181210    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
    … …  
    210239            if (r2 > R2) continue;
    211240
     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
    212247            psF32 pDiff = *vPix;
    213 
    214             Sum += pDiff;
    215 
    216             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
     248# endif
     249
     250            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    217251            psF32 rf = pDiff * sqrt(r2);
    218252            psF32 rs = pDiff;
    … …  
    234268
    235269    int nKronPix = 0;
    236     Sum = Var = 0.0;
     270    float Sum = 0.0;
     271    float Var = 0.0;
    237272
    238273    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    263298            if (r2 > radKron2) continue;
    264299
     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
    265307            psF32 pDiff = *vPix;
    266308            psF32 wDiff = *vWgt;
     309# endif
    267310
    268311            Sum += pDiff;
    … …  
    272315    }
    273316
    274     source->moments->Mrh = Mrf;
    275 
    276     // XXX for a test, save the old values here:
    277     source->moments->KronFouter    = source->moments->KronFlux;
    278     source->moments->KronCoreErr = source->moments->KronFluxErr;
    279 
    280     source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
    281     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;
    282322
    283323    return true;
    284324}
     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}
Note: See TracChangeset for help on using the changeset viewer.