IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 10, 2013, 2:55:11 PM (12 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20130904

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psModules

  • trunk/psModules/src/objects/pmSource.c

    r35768 r36375  
    6666    psFree(tmp->extpars);
    6767    psFree(tmp->diffStats);
     68    psFree(tmp->galaxyFits);
    6869    psFree(tmp->radialAper);
    6970    psTrace("psModules.objects", 10, "---- end ----\n");
    … …  
    164165    source->extpars = NULL;
    165166    source->diffStats = NULL;
     167    source->galaxyFits = NULL;
    166168    source->radialAper = NULL;
    167169    source->parent = NULL;
    … …  
    690692            // why do we recalculate moments here?
    691693            // we already attempt to do this in psphotSourceStats
    692             // pmSourceMoments (source, INNER_RADIUS);
    693694            Nsatstar ++;
    694695            continue;
    … …  
    804805    return true;
    805806}
    806 
    807 /******************************************************************************
    808 pmSourceMoments(source, radius): this function takes a subImage defined in the
    809 pmSource data structure, along with the peak location, and determines the
    810 various moments associated with that peak.
    811 
    812 Requires the following to have been created:
    813     pmSource
    814     pmSource->peak
    815     pmSource->pixels
    816     pmSource->variance
    817     pmSource->mask
    818 
    819 XXX: The peak calculations are done in image coords, not subImage coords.
    820 
    821 XXX EAM : this version clips input pixels on S/N
    822 XXX EAM : this version returns false for several reasons
    823 *****************************************************************************/
    824 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    825 
    826 /*** this been moved to pmSourceMoments.c ***/
    827 # if (0)
    828 bool pmSourceMoments(pmSource *source,
    829                      psF32 radius)
    830 {
    831     psTrace("psModules.objects", 10, "---- begin ----\n");
    832     PS_ASSERT_PTR_NON_NULL(source, false);
    833     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    834     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    835     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    836 
    837     //
    838     // XXX: Verify the setting for sky if source->moments == NULL.
    839     //
    840     psF32 sky = 0.0;
    841     if (source->moments == NULL) {
    842         source->moments = pmMomentsAlloc();
    843     } else {
    844         sky = source->moments->Sky;
    845     }
    846 
    847     //
    848     // Sum = SUM (z - sky)
    849     // X1  = SUM (x - xc)*(z - sky)
    850     // X2  = SUM (x - xc)^2 * (z - sky)
    851     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    852     //
    853     psF32 peakPixel = -PS_MAX_F32;
    854     psS32 numPixels = 0;
    855     psF32 Sum = 0.0;
    856     psF32 Var = 0.0;
    857     psF32 X1 = 0.0;
    858     psF32 Y1 = 0.0;
    859     psF32 X2 = 0.0;
    860     psF32 Y2 = 0.0;
    861     psF32 XY = 0.0;
    862     psF32 x  = 0;
    863     psF32 y  = 0;
    864     psF32 R2 = PS_SQR(radius);
    865 
    866     psF32 xPeak = source->peak->x;
    867     psF32 yPeak = source->peak->y;
    868     psF32 xOff = source->pixels->col0 - source->peak->x;
    869     psF32 yOff = source->pixels->row0 - source->peak->y;
    870 
    871     // XXX why do I get different results for these two methods of finding Sx?
    872     // XXX Sx, Sy would be better measured if we clip pixels close to sky
    873     // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
    874     // We loop through all pixels in this subimage (source->pixels), and for each
    875     // pixel that is not masked, AND within the radius of the peak pixel, we
    876     // proceed with the moments calculation.  need to do two loops for a
    877     // numerically stable result.  first loop: get the sums.
    878     // XXX EAM : mask == 0 is valid
    879 
    880     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    881 
    882         psF32 *vPix = source->pixels->data.F32[row];
    883         psF32 *vWgt = source->variance->data.F32[row];
    884         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    885 
    886         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    887             if (vMsk) {
    888                 if (*vMsk) {
    889                     vMsk++;
    890                     psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to mask: %d\n",
    891                             col, row, (int)*vMsk);
    892                     continue;
    893                 }
    894                 vMsk++;
    895             }
    896             if (isnan(*vPix)) continue;
    897 
    898             psF32 xDiff = col + xOff;
    899             psF32 yDiff = row + yOff;
    900 
    901             // radius is just a function of (xDiff, yDiff)
    902             if (!VALID_RADIUS(xDiff, yDiff, R2)) {
    903 #if 1
    904                 psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to position: %f %f\n",
    905                         col, row, xDiff, yDiff);
    906 #endif
    907                 continue;
    908             }
    909 
    910             psF32 pDiff = *vPix - sky;
    911             psF32 wDiff = *vWgt;
    912 
    913             // XXX EAM : check for valid S/N in pixel
    914             // XXX EAM : should this limit be user-defined?
    915 #if 1
    916             if (PS_SQR(pDiff) < wDiff) {
    917                 psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to insignificance: %f, %f\n",
    918                         col, row, pDiff, wDiff);
    919                 continue;
    920             }
    921 #endif
    922 
    923             Var += wDiff;
    924             Sum += pDiff;
    925 
    926             psF32 xWght = xDiff * pDiff;
    927             psF32 yWght = yDiff * pDiff;
    928 
    929             X1  += xWght;
    930             Y1  += yWght;
    931 
    932             XY  += xDiff * yWght;
    933             X2  += xDiff * xWght;
    934             Y2  += yDiff * yWght;
    935 
    936             peakPixel = PS_MAX (*vPix, peakPixel);
    937             numPixels++;
    938         }
    939     }
    940 
    941     // if we have less than (1/4) of the possible pixels, force a retry
    942     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    943     if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    944         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n",
    945                  numPixels, (int)(0.75*R2), Sum);
    946         psTrace("psModules.objects", 10, "---- end (false) ----\n");
    947         return (false);
    948     }
    949 
    950     psTrace ("psModules.objects", 4, "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
    951              sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    952 
    953     //
    954     // first moment X  = X1/Sum + xc
    955     // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    956     // Sxy             = XY / Sum
    957     //
    958     x = X1/Sum;
    959     y = Y1/Sum;
    960     if ((fabs(x) > radius) || (fabs(y) > radius)) {
    961         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n",
    962                  source->peak->x, source->peak->y);
    963         psTrace("psModules.objects", 10, "---- end(false)  ----\n");
    964         return (false);
    965     }
    966 
    967     source->moments->Mx = x + xPeak;
    968     source->moments->My = y + yPeak;
    969 
    970     // XXX EAM : Sxy needs to have x*y subtracted
    971     source->moments->Mxy = XY/Sum - x*y;
    972     source->moments->Sum = Sum;
    973     source->moments->SN  = Sum / sqrt(Var);
    974     source->moments->Peak = peakPixel;
    975     source->moments->nPixels = numPixels;
    976 
    977     // XXX EAM : these values can be negative, so we need to limit the range
    978     // XXX EAM : make the use of this consistent: should this be the second moment or sqrt?
    979     // source->moments->Mxx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    980     // source->moments->Myy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    981     source->moments->Mxx = PS_MAX(X2/Sum - PS_SQR(x), 0);
    982     source->moments->Myy = PS_MAX(Y2/Sum - PS_SQR(y), 0);
    983 
    984     psTrace ("psModules.objects", 4,
    985              "sky: %f  Sum: %f  Mx: %f  My: %f  Mxx: %f  Myy: %f  Mxy: %f\n",
    986              sky, Sum, source->moments->Mx, source->moments->My,
    987              source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
    988 
    989     psTrace("psModules.objects", 10, "---- end ----\n");
    990     return(true);
    991 }
    992 # endif
    993807
    994808// construct a realization of the source model
Note: See TracChangeset for help on using the changeset viewer.