IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psModules

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1212pmSourceIO_CMF_PS1_V1.v1.c
        1313pmSourceIO_CMF_PS1_V4.c
         14pmSourceIO_CMF_PS1_V5.c
        1415pmSourceIO_CMF_PS1_SV1.c
        1516pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
        1618pmSourceIO_CMF_PS1_DV1.c
        1719pmSourceIO_CMF_PS1_DV2.c
        1820pmSourceIO_CMF_PS1_DV3.c
        19 
         21pmSourceIO_CMF_PS1_DV4.c

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects/pmSource.c

    r35768 r37066  
    3333#include "pmMoments.h"
    3434#include "pmModelFuncs.h"
     35#include "pmModelClass.h"
    3536#include "pmModel.h"
    3637#include "pmModelUtils.h"
    37 #include "pmModelClass.h"
    3838#include "pmSourceMasks.h"
    3939#include "pmSourceExtendedPars.h"
    … …  
    4141#include "pmSourcePhotometry.h"
    4242#include "pmSourceSatstar.h"
     43#include "pmSourceLensing.h"
    4344#include "pmSource.h"
    4445
    … …  
    6667    psFree(tmp->extpars);
    6768    psFree(tmp->diffStats);
     69    psFree(tmp->galaxyFits);
    6870    psFree(tmp->radialAper);
     71    psFree(tmp->lensingOBJ);
     72    psFree(tmp->lensingPSF);
    6973    psTrace("psModules.objects", 10, "---- end ----\n");
    7074}
    … …  
    142146    source->apMagRaw         = NAN;
    143147    source->apRadius         = NAN;
     148    source->apNpixels        = 0;
    144149    source->apFlux           = NAN;
    145150    source->apFluxErr        = NAN;
    … …  
    158163    source->sky              = NAN;
    159164    source->skyErr           = NAN;   
     165    source->extSN            = NAN;   
    160166
    161167    source->region = psRegionSet(NAN, NAN, NAN, NAN);
    … …  
    164170    source->extpars = NULL;
    165171    source->diffStats = NULL;
     172    source->galaxyFits = NULL;
     173    source->lensingOBJ = NULL;
     174    source->lensingPSF = NULL;
    166175    source->radialAper = NULL;
    167176    source->parent = NULL;
    … …  
    229238    source->apMagRaw         = in->apMagRaw;
    230239    source->apRadius         = in->apRadius;
     240    source->apNpixels        = in->apNpixels;
    231241    source->apFlux           = in->apFlux;
    232242    source->apFluxErr        = in->apFluxErr;
    … …  
    247257
    248258    source->region           = in->region;
     259
     260    // XXX I am not copying the pointers to things like the blends, satstar profile, galaxyFits, etc
    249261
    250262    return(source);
    … …  
    690702            // why do we recalculate moments here?
    691703            // we already attempt to do this in psphotSourceStats
    692             // pmSourceMoments (source, INNER_RADIUS);
    693704            Nsatstar ++;
    694705            continue;
    … …  
    804815    return true;
    805816}
    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
    993817
    994818// construct a realization of the source model
    … …  
    1165989    if (!isfinite(oldI0)) return false;
    1166990
    1167     bool useReff = pmModelUseReff (model->type);
     991    bool useReff = model->class->useReff;
    1168992    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
    1169993
Note: See TracChangeset for help on using the changeset viewer.