IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 29, 2011, 1:26:12 PM (15 years ago)
Author:
eugene
Message:

improvements (?) to the kron flux analysis : attempt to window neighbors based on the symmetry of the object

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r32205 r32216  
    1010{
    1111    bool status = true;
     12
     13    // return true;
    1214
    1315    // select the appropriate recipe information
    … …  
    127129    psphotSaveImage (NULL, kronWindow, "kron.window.v0.fits");
    128130
    129     for (int i = 0; i < sources->n; i++) {
    130 
    131         pmSource *source = sources->data[i];
    132         if (!source->peak) continue; // XXX how can we have a peak-less source?
    133 
    134         // allocate space for moments
    135         if (!source->moments) continue;
    136 
    137         // replace object in image
    138         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    139             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     131    for (int j = 0; j < 5; j++) {
     132        for (int i = 0; i < sources->n; i++) {
     133
     134            pmSource *source = sources->data[i];
     135            if (!source->peak) continue; // XXX how can we have a peak-less source?
     136
     137            // allocate space for moments
     138            if (!source->moments) continue;
     139
     140            // replace object in image
     141            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     142                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     143            }
     144
     145            // iterate to the window radius
     146            float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
     147
     148            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
     149            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     150
     151            // clear the window function for this source based on the moments
     152            psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
     153            // psphotKronWindowSetSource (source, kronWindow, false, false);
     154            // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
     155
     156            // 165, 539;
     157            if ((fabs(source->peak->xf - 165) < 3) && (fabs(source->peak->yf - 539) < 3)) {
     158                fprintf (stderr, "test obj\n");
     159            }
     160
     161            // this function populates moments->Mrf,KronFlux,KronFluxErr
     162            psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
     163            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     164
     165            // set a window function for each source based on the moments
     166            psphotKronWindowSetSource (source, kronWindow, true, true);
     167            // psphotKronWindowSetSource (source, kronWindow, false, true);
     168
     169            // test source fluxes
     170            pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
     171            float kmag = -2.5*log10(source->moments->KronFlux);
     172# define TEST_X1 167
     173# define TEST_Y1 299
     174# define TEST_X2 180
     175# define TEST_Y2 300
     176            if ((fabs(source->peak->xf - TEST_X1) < 3) && (fabs(source->peak->yf - TEST_Y1) < 3)) {
     177                fprintf (stderr, "R1: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
     178            }
     179            if ((fabs(source->peak->xf - TEST_X2) < 3) && (fabs(source->peak->yf - TEST_Y2) < 3)) {
     180                fprintf (stderr, "R2: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
     181            }
     182
     183            // re-subtract the object, leave local sky
     184            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    140185        }
    141 
    142         // iterate to the window radius
    143         float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
    144 
    145         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    146         pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    147 
    148         // clear the window function for this source based on the moments
    149         psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
    150         // psphotKronWindowSetSource (source, kronWindow, false, false);
    151         // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
    152 
    153         // this function populates moments->Mrf,KronFlux,KronFluxErr
    154         psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    155         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    156 
    157         // set a window function for each source based on the moments
    158         psphotKronWindowSetSource (source, kronWindow, true, true);
    159         // psphotKronWindowSetSource (source, kronWindow, false, true);
    160 
    161         // test source fluxes
    162         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    163         float kmag = -2.5*log10(source->moments->KronFlux);
    164         if (source->psfMag - kmag > 0.25) {
    165             // fprintf (stderr, "continue\n");
    166         }
    167 
    168         // re-subtract the object, leave local sky
    169         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    170     }
    171 
    172     psphotSaveImage (NULL, kronWindow, "kron.window.v1.fits");
    173 
    174     for (int i = 0; i < sources->n; i++) {
    175 
    176         pmSource *source = sources->data[i];
    177         if (!source->peak) continue; // XXX how can we have a peak-less source?
    178 
    179         // allocate space for moments
    180         if (!source->moments) continue;
    181 
    182         // replace object in image
    183         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    184             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    185         }
    186 
    187         // iterate to the window radius
    188         float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
    189 
    190         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    191         pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    192 
    193         // clear the window function for this source based on the moments
    194         psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
    195         // psphotKronWindowSetSource (source, kronWindow, false, false);
    196         // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
    197 
    198         // this function populates moments->Mrf,KronFlux,KronFluxErr
    199         psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    200         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    201 
    202         // set a window function for each source based on the moments
    203         psphotKronWindowSetSource (source, kronWindow, true, true);
    204         // psphotKronWindowSetSource (source, kronWindow, false, true);
    205 
    206         // test source fluxes
    207         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    208         float kmag = -2.5*log10(source->moments->KronFlux);
    209         if (source->psfMag - kmag > 0.25) {
    210             // fprintf (stderr, "continue\n");
    211         }
    212 
    213         // re-subtract the object, leave local sky
    214         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    215     }
    216 
    217     psphotSaveImage (NULL, kronWindow, "kron.window.v2.fits");
     186        char name[64];
     187        sprintf (name, "kron.window.v%d.fits", j+1);
     188        psphotSaveImage (NULL, kronWindow, name);
     189    }
    218190    psFree (kronWindow);
    219191
    … …  
    231203
    232204    psF32 R2 = PS_SQR(radius);
     205    float rsigma2 = 0.5 / PS_SQR(radius/2.0);
    233206
    234207    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    261234    int Ywo = source->pixels->row0;
    262235
     236    psF32 **vPix = source->pixels->data.F32;
     237    psF32 **vWin = kronWindow->data.F32;
     238    psF32 **vWgt = source->variance->data.F32;
     239   
     240    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     241
    263242    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    264243
    … …  
    266245        if (fabs(yDiff) > radius) continue;
    267246
    268         psF32 *vPix = source->pixels->data.F32[row];
    269         psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
    270 
    271         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    272 
    273         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWin++) {
    274             if (vMsk) {
    275                 if (*vMsk & maskVal) {
    276                     vMsk++;
    277                     continue;
    278                 }
    279                 vMsk++;
    280             }
    281             if (isnan(*vPix)) continue;
     247        // coordinate of mirror pixel
     248        int yFlip = yCM - yDiff;
     249        if (yFlip < 0) continue;
     250        if (yFlip >= source->pixels->numRows) continue;
     251
     252        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     253            // check mask and value for this pixel
     254            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     255            if (isnan(vPix[row][col])) continue;
    282256
    283257            psF32 xDiff = col - xCM;
    284258            if (fabs(xDiff) > radius) continue;
     259
     260            // coordinate of mirror pixel
     261            int xFlip = xCM - xDiff;
     262            if (xFlip < 0) continue;
     263            if (xFlip >= source->pixels->numCols) continue;
     264
     265            // check mask and value for mirror pixel
     266            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     267            if (isnan(vPix[yFlip][xFlip])) continue;
    285268
    286269            // radius is just a function of (xDiff, yDiff)
    … …  
    289272
    290273            // flux * window
    291             float weight  = *vWin;
     274            float z = r2 * rsigma2;
     275            assert (z >= 0.0);
     276
    292277            // float weight  = 1.0;
    293             psF32 pDiff = *vPix * weight;
     278            float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
     279            float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
     280            // float weight1  = vWin[row+Ywo][col+Xwo];
     281            // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
     282
     283            float fDiff1 = vPix[row][col]*weight1;
     284            float fDiff2 = vPix[yFlip][xFlip]*weight2;
     285
     286            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    294287
    295288            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    296289            psF32 rf = pDiff * sqrt(r2);
    297             psF32 rs = pDiff;
     290            psF32 rs = 0.5 * (fDiff1 + fDiff2);
    298291
    299292            RF  += rf;
    … …  
    322315        if (fabs(yDiff) > radKron) continue;
    323316
    324         psF32 *vPix = source->pixels->data.F32[row];
    325         psF32 *vWgt = source->variance->data.F32[row];
    326         psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
    327 
    328         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    329 
    330         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++, vWin++) {
    331             if (vMsk) {
    332                 if (*vMsk & maskVal) {
    333                     vMsk++;
    334                     continue;
    335                 }
    336                 vMsk++;
    337             }
    338             if (isnan(*vPix)) continue;
     317        // coordinate of mirror pixel
     318        int yFlip = yCM - yDiff;
     319        if (yFlip < 0) continue;
     320        if (yFlip >= source->pixels->numRows) continue;
     321
     322        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     323            // check mask and value for this pixel
     324            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     325            if (isnan(vPix[row][col])) continue;
    339326
    340327            psF32 xDiff = col - xCM;
    341328            if (fabs(xDiff) > radKron) continue;
     329
     330            // coordinate of mirror pixel
     331            int xFlip = xCM - xDiff;
     332            if (xFlip < 0) continue;
     333            if (xFlip >= source->pixels->numCols) continue;
     334
     335            // check mask and value for mirror pixel
     336            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     337            if (isnan(vPix[yFlip][xFlip])) continue;
    342338
    343339            // radKron is just a function of (xDiff, yDiff)
    … …  
    345341            if (r2 > radKron2) continue;
    346342
    347             float weight  = *vWin;
     343            // float z = r2 * rsigma2;
     344            // assert (z >= 0.0);
     345
    348346            // float weight  = 1.0;
    349             psF32 pDiff = *vPix * weight;
    350             psF32 wDiff = *vWgt * weight;
     347            // float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
     348            // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
     349            float weight1  = vWin[row+Ywo][col+Xwo];
     350            float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
     351
     352            float fDiff1 = vPix[row][col]*weight1;
     353            float fDiff2 = vPix[yFlip][xFlip]*weight2;
     354
     355            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     356            psF32 wDiff = vWgt[row][col] * weight1;
    351357
    352358            Sum += pDiff;
    353359            Var += wDiff;
    354             Win += weight;
     360            Win += weight1;
    355361            nKronPix ++;
    356362        }
    … …  
    387393    float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    388394
     395    // float kratio = source->moments->KronFinner / source->moments->KronFlux;
     396
     397    float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
    389398    // float scale = useKronRadius ? 2.0 * source->moments->Mrf / Mmajor : 2.0;
    390     float scale = 3.0 * source->moments->Mrf / Mmajor;
    391     // float scale = 2.0;
     399    // float scale = (kratio > 0.4) ? 9.0 * source->moments->Mrf / Mmajor : 3.0 * source->moments->Mrf / Mmajor;
    392400
    393401    float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
Note: See TracChangeset for help on using the changeset viewer.