IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5048


Ignore:
Timestamp:
Sep 13, 2005, 3:35:20 PM (21 years ago)
Author:
eugene
Message:

substantial work to match with psLib v.7

Location:
trunk/psphot
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/doc/notes.txt

    r4983 r5048  
    5050  - fixed 'center' option in pmSourceAdd/Sub
    5151
     52  - pmPSF params should be F64
     53
    5254  lingering concerns:
    5355  - Polynomials and related functions are STILL messed up.  this is a

  • trunk/psphot/src/pmObjects_EAM.c

    r4977 r5048  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-09-08 04:32:40 $
     8 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-09-14 01:35:20 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    603603     members.
    604604*****************************************************************************/
    605 pmSource *pmSourceLocalSky(const psImage *image,
    606                            const pmPeak *peak,
    607                            psStatsOptions statsOptions,
    608                            psF32 innerRadius,
    609                            psF32 outerRadius)
    610 {
    611     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    612     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
    613     PS_ASSERT_PTR_NON_NULL(peak, NULL);
    614     PS_FLOAT_COMPARE(0.0, innerRadius, NULL);
    615     PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL);
    616     psS32 innerRadiusS32 = (psS32) innerRadius;
    617     psS32 outerRadiusS32 = (psS32) outerRadius;
    618 
    619     //
    620     // We define variables for code readability.
    621     //
    622     // XXX: Since the peak->xy coords are in image, not subImage coords,
    623     // these variables should be renamed for clarity (imageCenterRow, etc).
    624     //
    625     // peak->x,y is guaranteed to be on image
    626     psS32 SubImageCenterRow = peak->y;
    627     psS32 SubImageCenterCol = peak->x;
    628 
    629     // XXX EAM : I added this code to stay on the image. So did George
    630     // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
    631     // XXX EAM : dropped the Annulus width (not needed)
    632     // XXX EAM : dropped off-boundary tests (not needed)
    633     psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
    634     psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1);
    635     psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
    636     psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1);
    637     //
    638     // Grab a subimage of the original image of size (2 * outerRadius).
    639     //
    640     // XXX EAM : merged psImageSubset & psRegionSet
    641     // XXX EAM : cast for image is needed because of const (above)
    642     psImage *subImage = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
    643                                                          SubImageEndCol,
    644                                                          SubImageStartRow,
    645                                                          SubImageEndRow));
    646 
    647     psImage *subImageMask = psImageAlloc(subImage->numCols,
    648                                          subImage->numRows,
    649                                          PS_TYPE_U8);
    650     // XXX EAM : for consistency, mask needs col0,row0 set to match image
    651     // XXX EAM : CONFLICT between psLib and code requirement. FIX PSLIB!!!
    652     // subImageMask->col0 = subImage->col0;
    653     // subImageMask->row0 = subImage->row0;
    654     psAbort ("pmObjects", "must fix this error before psphot will work!!!\n");
    655 
    656     //
    657     // Loop through the subimage mask & initialize PSPHOT_MASK_INVALID
    658     //
    659     for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
    660         for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
    661             subImageMask->data.U8[row][col] = PSPHOT_MASK_INVALID;
    662         }
    663     }
    664 
    665     //
    666     // Loop through the subimage, mask off pixels in the inner square.
    667     // XXX this uses a static mask value of 0
    668     // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image
    669     // XXX EAM : this masks the pixels in the center region
    670     // XXX EAM : should we be using psImageMaskRegion ???
    671     psS32 StartRow  = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32);
    672     psS32 EndRow    = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1);
    673     psS32 StartCol  = PS_MAX (0, SubImageCenterCol - subImageMask->col0 - innerRadiusS32);
    674     psS32 EndCol    = PS_MIN (subImageMask->numCols, SubImageCenterCol - subImageMask->col0 + innerRadiusS32 + 1);
    675     for (psS32 row = StartRow; row < EndRow; row++) {
    676         for (psS32 col = StartCol; col < EndCol; col++) {
    677             subImageMask->data.U8[row][col] = PSPHOT_MASK_CLEAR;
    678         }
    679     }
    680 
    681     // XXX EAM : make this trace information??
    682     //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
    683     //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
    684     //            printf("(%d) ", subImageMask->data.U8[row][col]);
    685     //        }
    686     //        printf("\n");
    687     //    }
    688 
    689     //
    690     // Allocate the myStats structure, then call psImageStats(), which will
    691     // calculate the specified statistic.
    692     //
     605bool pmSourceLocalSky (pmSource *source,
     606                           psStatsOptions statsOptions,
     607                           psF32 Radius)
     608{
     609
     610    psImage *image = source->pixels;
     611    psImage *mask  = source->mask;
     612    pmPeak *peak  = source->peak;
     613    psRegion srcRegion;
     614
     615    // XXX EAM : psRegionXXX funcs should take value not ptr
     616    srcRegion = psRegionForSquare (peak->x, peak->y, Radius);
     617    srcRegion = psRegionForImage (mask, &srcRegion);
     618
     619    psImageMaskRegion (mask, &srcRegion, "OR", PSPHOT_MASK_MARKED);
    693620    psStats *myStats = psStatsAlloc(statsOptions);
    694     myStats = psImageStats(myStats, subImage, subImageMask, 1);
    695 
    696     //
    697     // Create the output mySource, and set appropriate members.
    698     //
    699     pmSource *mySource = pmSourceAlloc();
    700     mySource->peak = (pmPeak *) peak;
    701     mySource->moments = pmMomentsAlloc();
     621    myStats = psImageStats(myStats, image, mask, 0xff);
     622    psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED);
     623
    702624    psF64 tmpF64;
    703625    p_psGetStatValue(myStats, &tmpF64);
    704     mySource->moments->Sky = (psF32) tmpF64;
    705     mySource->pixels = subImage;
    706     mySource->mask = subImageMask;
    707 
    708     //
    709     // Free things.  XXX: This should be static memory.
    710     //
    711     psFree(myStats);
    712 
    713     return(mySource);
    714 }
    715 
    716 /******************************************************************************
    717 bool checkRadius(*peak, radius, x, y): private function which simply
    718 determines if the (x, y) point is within the radius of the specified peak.
    719  
    720 XXX: macro this for performance.
    721 *****************************************************************************/
    722 static bool checkRadius(pmPeak *peak,
    723                         psF32 radius,
    724                         psS32 x,
    725                         psS32 y)
    726 {
    727     if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) {
    728         return(true);
    729     }
    730 
    731     return(false);
     626    psFree (myStats);
     627
     628    if (isnan(tmpF64)) return (false);
     629    source->moments = pmMomentsAlloc();
     630    source->moments->Sky = (psF32) tmpF64;
     631    return (true);
    732632}
    733633
    … …  
    762662    pmSource->peak
    763663    pmSource->pixels
     664    pmSource->noise
     665    pmSource->mask
    764666 
    765667XXX: The peak calculations are done in image coords, not subImage coords.
    766668 
    767 XXX: mask values?
    768 *****************************************************************************/
     669XXX EAM : this version clips input pixels on S/N
     670XXX EAM : this version returns false for several reasons
     671*****************************************************************************/
     672# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
     673
    769674bool pmSourceMoments(pmSource *source,
    770675                     psF32 radius)
    771676{
    772     PS_ASSERT_PTR_NON_NULL(source, NULL);
    773     PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
    774     PS_ASSERT_PTR_NON_NULL(source->pixels, NULL);
     677    // PS_PTR_CHECK_NULL(source, NULL);
     678    // PS_PTR_CHECK_NULL(source->peak, NULL);
     679    // PS_PTR_CHECK_NULL(source->pixels, NULL);
     680    // PS_PTR_CHECK_NULL(source->noise, NULL);
    775681    PS_FLOAT_COMPARE(0.0, radius, NULL);
    776682
    … …  
    791697    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    792698    //
    793     psF32 Sum = 0.0;
    794699    psF32 peakPixel = -PS_MAX_F32;
    795700    psS32 numPixels = 0;
     701    psF32 Sum = 0.0;
    796702    psF32 X1 = 0.0;
    797703    psF32 Y1 = 0.0;
    … …  
    799705    psF32 Y2 = 0.0;
    800706    psF32 XY = 0.0;
    801     psF32 x = 0;
    802     psF32 y = 0;
    803     //
     707    psF32 x  = 0;
     708    psF32 y  = 0;
     709    psF32 R2 = PS_SQR(radius);
     710
     711    psF32 xPeak = source->peak->x;
     712    psF32 yPeak = source->peak->y;
     713
    804714    // XXX why do I get different results for these two methods of finding Sx?
    805715    // XXX Sx, Sy would be better measured if we clip pixels close to sky
    … …  
    809719    // proceed with the moments calculation.  need to do two loops for a
    810720    // numerically stable result.  first loop: get the sums.
    811     //
     721    // XXX EAM : mask == 0 is valid
     722
    812723    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    813724        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    814             // XXX EAM : mask == 0 is valid
    815725            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue;
    816             psS32 imgColCoord = col + source->pixels->col0;
    817             psS32 imgRowCoord = row + source->pixels->row0;
    818             if (checkRadius(source->peak,
    819                             radius,
    820                             imgColCoord,
    821                             imgRowCoord)) {
    822                 psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
    823                 psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
    824                 psF32 pDiff = source->pixels->data.F32[row][col] - sky;
    825 
    826                 Sum+= pDiff;
    827                 X1+= xDiff * pDiff;
    828                 Y1+= yDiff * pDiff;
    829                 XY+= xDiff * yDiff * pDiff;
    830 
    831                 X2+= PS_SQR(xDiff) * pDiff;
    832                 Y2+= PS_SQR(yDiff) * pDiff;
    833 
    834                 if (source->pixels->data.F32[row][col] > peakPixel) {
    835                     peakPixel = source->pixels->data.F32[row][col];
    836                 }
    837                 numPixels++;
    838             }
     726
     727            psF32 xDiff = col + source->pixels->col0 - xPeak;
     728            psF32 yDiff = row + source->pixels->row0 - yPeak;
     729
     730            // XXX EAM : calculate xDiff, yDiff up front;
     731            //           radius is just a function of (xDiff, yDiff)
     732            if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
     733
     734            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     735
     736            // XXX EAM : check for valid S/N in pixel
     737            // XXX EAM : should this limit be user-defined?
     738            if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue;
     739           
     740            Sum += pDiff;
     741            X1  += xDiff * pDiff;
     742            Y1  += yDiff * pDiff;
     743            XY  += xDiff * yDiff * pDiff;
     744           
     745            X2  += PS_SQR(xDiff) * pDiff;
     746            Y2  += PS_SQR(yDiff) * pDiff;
     747           
     748            peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
     749            numPixels++;
    839750        }
    840751    }
    841752    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    842753    if ((numPixels < 3) || (Sum <= 0)) {
    843         psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
    844         return (false);
     754      psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     755      return (false);
    845756    }
    846757
    … …  
    857768    y = Y1/Sum;
    858769    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    859         psTrace (".psModules.pmSourceMoments", 5,
    860                  "large centroid swing; invalid peak %d, %d\n",
    861                  source->peak->x, source->peak->y);
    862         return (false);
    863     }
    864 
    865     source->moments->x = x + ((psF32) source->peak->x);
    866     source->moments->y = y + ((psF32) source->peak->y);
    867 
     770      psTrace (".psModules.pmSourceMoments", 5,
     771               "large centroid swing; invalid peak %d, %d\n",
     772               source->peak->x, source->peak->y);
     773      return (false);
     774    }
     775
     776    source->moments->x = x + xPeak;
     777    source->moments->y = y + yPeak;
     778
     779    // XXX EAM : Sxy needs to have x*y subtracted
    868780    source->moments->Sxy = XY/Sum - x*y;
    869781    source->moments->Sum = Sum;
    … …  
    919831*****************************************************************************/
    920832
     833// XXX EAM include a S/N cutoff in selecting the sources?
    921834pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
    922835{
    … …  
    11231036        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
    11241037        if ((Nsatpix > 1) && big) {
    1125             tmpSrc->type |= PS_SOURCE_SATSTAR;
     1038            tmpSrc->type |= PM_SOURCE_SATSTAR;
    11261039            Nsatstar ++;
    11271040            continue;
    … …  
    11301043        // saturated object (not a star, eg bleed trails, hot pixels)
    11311044        if (Nsatpix > 1) {
    1132             tmpSrc->type |= PS_SOURCE_SATURATED;
     1045            tmpSrc->type |= PM_SOURCE_SATURATED;
    11331046            Nsat ++;
    11341047            continue;
    … …  
    11391052        // only set candidate defects if
    11401053        if ((sigX < 0.05) || (sigY < 0.05)) {
    1141             tmpSrc->type |= PS_SOURCE_DEFECT;
     1054            tmpSrc->type |= PM_SOURCE_DEFECT;
    11421055            Ncr ++;
    11431056            continue;
    … …  
    11461059        // likely unsaturated galaxy (too large to be stellar)
    11471060        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    1148             tmpSrc->type |= PS_SOURCE_GALAXY;
     1061            tmpSrc->type |= PM_SOURCE_GALAXY;
    11491062            Ngal ++;
    11501063            continue;
    … …  
    11591072        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    11601073        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    1161             tmpSrc->type |= PS_SOURCE_PSFSTAR;
     1074            tmpSrc->type |= PM_SOURCE_PSFSTAR;
    11621075            Npsf ++;
    11631076            continue;
    … …  
    11651078
    11661079        // random type of star
    1167         tmpSrc->type |= PS_SOURCE_OTHER;
     1080        tmpSrc->type |= PM_SOURCE_OTHER;
    11681081    }
    11691082
    … …  
    11971110XXX: The circle will have a diameter of (1+radius).  This is different from
    11981111     the pmSourceSetLocal() function.
     1112
     1113XXX: this should probably use the psImageMaskCircle Function for consistency
    11991114*****************************************************************************/
    12001115bool pmSourceSetPixelsCircle(pmSource *source,
    … …  
    15371452XXX EAM : fit the specified model (not necessarily the one in source)
    15381453*****************************************************************************/
    1539 bool pmSourceFitModel(pmSource *source,
    1540                       pmModel *model,
    1541                       const bool PSF)
     1454bool pmSourceFitModel_v5(pmSource *source,
     1455                        pmModel *model,
     1456                        const bool PSF)
    15421457{
    15431458    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    15541469    //           tests below could be conditions (!NULL)
    15551470
    1556     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1557 
    15581471    psVector *params = model->params;
    15591472    psVector *dparams = model->dparams;
    15601473    psVector *paramMask = NULL;
     1474
     1475    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    15611476
    15621477    int nParams = PSF ? params->n - 4 : params->n;
    … …  
    16151530    }       
    16161531
    1617     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     1532    // XXX EAM : covar must be F64?
     1533    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    16181534
    16191535    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    … …  
    16211537    for (int i = 0; i < dparams->n; i++) {
    16221538        if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
    1623         dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     1539        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    16241540    }
    16251541 
    … …  
    16401556    // XXX EAM get the Gauss-Newton distance for fixed model parameters
    16411557    if (paramMask != NULL) {
    1642         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
     1558        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    16431559        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
    16441560        for (int i = 0; i < dparams->n; i++) {
    16451561            if (!paramMask->data.U8[i]) continue;
    1646             dparams->data.F32[i] = delta->data.F32[i];
     1562            dparams->data.F32[i] = delta->data.F64[i];
    16471563        }
    16481564        psFree (delta);
    … …  
    16551571    psFree(covar);
    16561572    psFree(paramMask);
     1573
     1574    rc = (onPic && fitStatus);
     1575    return(rc);
     1576}
     1577
     1578// XXX EAM : new version with parameter range limits and noise enhancement
     1579bool pmSourceFitModel (pmSource *source,
     1580                       pmModel *model,
     1581                       const bool PSF)
     1582{
     1583    PS_ASSERT_PTR_NON_NULL(source, false);
     1584    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     1585    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1586    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     1587    PS_ASSERT_PTR_NON_NULL(source->mask, false);
     1588    PS_ASSERT_PTR_NON_NULL(source->noise, false);
     1589
     1590    // XXX EAM : is it necessary for the mask & noise to exist?  the
     1591    //           tests below could be conditions (!NULL)
     1592
     1593    psBool fitStatus = true;
     1594    psBool onPic     = true;
     1595    psBool rc        = true;
     1596    psF32  Ro, ymodel;
     1597
     1598    psVector *params = model->params;
     1599    psVector *dparams = model->dparams;
     1600    psVector *paramMask = NULL;
     1601
     1602    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     1603
     1604    // XXX EAM : I need to use the sky value to constrain the noise model
     1605    int nParams = PSF ? params->n - 4 : params->n;
     1606    psF32 So = params->data.F32[0];
     1607
     1608    // find the number of valid pixels
     1609    // XXX EAM : this loop and the loop below could just be one pass
     1610    //           using the psArrayAdd and psVectorExtend functions
     1611    psS32 count = 0;
     1612    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1613        for (psS32 j = 0; j < source->pixels->numCols; j++) {
     1614            if (source->mask->data.U8[i][j] == 0) {
     1615                count++;
     1616            }
     1617        }
     1618    }
     1619    if (count <  nParams + 1) {
     1620        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
     1621        return(false);
     1622    }
     1623
     1624    // construct the coordinate and value entries
     1625    psArray *x = psArrayAlloc(count);
     1626    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
     1627    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
     1628    psS32 tmpCnt = 0;
     1629    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1630        for (psS32 j = 0; j < source->pixels->numCols; j++) {
     1631            if (source->mask->data.U8[i][j] == 0) {
     1632                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1633                // XXX: Convert i/j to image space:
     1634                // XXX EAM: coord order is (x,y) == (col,row)
     1635                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1636                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     1637                x->data[tmpCnt] = (psPtr *) coord;
     1638                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
     1639
     1640                // compare observed flux to model flux to adjust weight
     1641                ymodel = modelFunc (NULL, model->params, coord);
     1642               
     1643                // this test enhances the noise based on deviation from the model flux
     1644                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
     1645
     1646                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
     1647                //           the minimization function calculates sq().
     1648                //           should psMinimizeLMChi2 take dY^2?
     1649                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro);
     1650                tmpCnt++;
     1651            }
     1652        }
     1653    }
     1654
     1655    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
     1656                                                PM_SOURCE_FIT_MODEL_TOLERANCE);
     1657
     1658    // PSF model only fits first 4 parameters, FLT model fits all
     1659    if (PSF) {
     1660        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
     1661        for (int i = 0; i < 4; i++) {
     1662            paramMask->data.U8[i] = 0;
     1663        }
     1664        for (int i = 4; i < paramMask->n; i++) {
     1665            paramMask->data.U8[i] = 1;
     1666        }
     1667    } 
     1668
     1669    // XXX EAM : I've added three types of parameter range checks
     1670    // XXX EAM : this requires my new psMinimization functions
     1671    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
     1672    psVector *beta_lim = NULL;
     1673    psVector *params_min = NULL;
     1674    psVector *params_max = NULL;
     1675
     1676    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
     1677    //           in the SDRS, I define a new psMinimization which will take these in
     1678    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
     1679    modelLimits (&beta_lim, &params_min, &params_max);
     1680    for (int i = 0; i < params->n; i++) {
     1681        covar->data.F64[0][i] = beta_lim->data.F32[i];
     1682        covar->data.F64[1][i] = params_min->data.F32[i];
     1683        covar->data.F64[2][i] = params_max->data.F32[i];
     1684    }
     1685
     1686    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
     1687    fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
     1688    for (int i = 0; i < dparams->n; i++) {
     1689        if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
     1690        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     1691    }
     1692 
     1693    // XXX EAM: we need to do something (give an error?) if rc is false
     1694    // XXX EAM: psMinimizeLMChi2 does not check convergence
     1695
     1696    // XXX models can go insane: reject these
     1697    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1698    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1699    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1700    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1701
     1702    // XXX EAM: save the resulting chisq, nDOF, nIter
     1703    model->chisq = myMin->value;
     1704    model->nIter = myMin->iter;
     1705    model->nDOF  = y->n - nParams;
     1706
     1707    // XXX EAM get the Gauss-Newton distance for fixed model parameters
     1708    if (paramMask != NULL) {
     1709        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
     1710        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
     1711        for (int i = 0; i < dparams->n; i++) {
     1712            if (!paramMask->data.U8[i]) continue;
     1713            dparams->data.F32[i] = delta->data.F64[i];
     1714        }
     1715    }
     1716
     1717    psFree(paramMask);
     1718    psFree(x);
     1719    psFree(y);
     1720    psFree(myMin);
    16571721
    16581722    rc = (onPic && fitStatus);

  • trunk/psphot/src/pmObjects_EAM.h

    r4983 r5048  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-09-09 18:31:29 $
     7 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-09-14 01:35:20 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    210210
    211211*****************************************************************************/
    212 pmSource *pmSourceLocalSky(
    213     const psImage *image,               ///< The input image (float)
    214     const pmPeak *peak,                 ///< The peak for which the pmSource struct is created.
     212bool pmSourceLocalSky(
     213    pmSource *source,                   ///< The input image (float)
    215214    psStatsOptions statsOptions,        ///< The statistic used in calculating the background sky
    216     float innerRadius,                  ///< The inner radius of the suqare annulus for calculating sky
    217     float outerRadius                   ///< The outer radius of the suqare annulus for calculating sky
    218 );
     215    float Radius                        ///< The inner radius of the square annulus to exclude
     216    );
    219217
    220218/******************************************************************************
    … …  
    298296);
    299297
     298bool pmSourceFitModel_v5(
     299    pmSource *source,                   ///< The input pmSource
     300    pmModel *model,                     ///< model to be fitted
     301    const bool PSF                      ///< Treat model as PSF or FLT?
     302);
     303
     304bool pmSourceFitModel_v7(
     305    pmSource *source,                   ///< The input pmSource
     306    pmModel *model,                     ///< model to be fitted
     307    const bool PSF                      ///< Treat model as PSF or FLT?
     308);
     309
    300310#endif

  • trunk/psphot/src/psLibUtils.c

    r4977 r5048  
    147147}
    148148
     149# if (0)
    149150// alternate implementation of this function from pmObjects.c
    150 // XXX EAM : move this to pmObjects_EAM.c
     151// now defined in psLib SDRS as psImageRow
    151152psVector *psGetRowVectorFromImage(psImage *image,
    152153                                  psU32 row)
    … …  
    159160    return(tmpVector);
    160161}
     162# endif
    161163
    162164// XXX EAM : this is now in psLib

  • trunk/psphot/src/psLibUtils.h

    r4983 r5048  
    3939int          psArgumentRemove (int N, int *argc, char **argv); // added to SDRS
    4040psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type); // added to SDRS
    41 psVector    *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow)
     41// psVector    *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow)
    4242int          psImageCountPixelMask (psImage *mask, psU8 value); // added to SDRS
    4343
    … …  
    5656
    5757// polynomial functions -- all in SDRS, not done for v.15
     58// XXX the pslib 0.7.0 polynomials are *still* nTerm, not nOrder!!!
    5859psPolynomial1D *Polynomial1DAlloc_EAM(psPolynomialType type, psS32 nOrder);
    5960void            Polynomial1DDump_EAM(psPolynomial1D *poly);

  • trunk/psphot/src/psMinimize.c

    r4977 r5048  
    9595
    9696    // iterate until the tolerance is reached, or give up
    97     while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
     97    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
    9898
    9999        // set a new guess for Alpha, Beta, Params

  • trunk/psphot/src/psModulesUtils.c

    r4977 r5048  
    22
    33// extract config informatin from config data or from image header, if specified
     4// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme
    45psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
    56
    … …  
    4849}
    4950
    50 // XXX EAM : these are my alternate implementations of psModule functions
    51 
    52 // this sets and clears bit 0x80
    53 bool pmSourceLocalSky_EAM (pmSource *source,
    54                            psStatsOptions statsOptions,
    55                            psF32 Radius)
    56 {
    57 
    58     psImage *image = source->pixels;
    59     psImage *mask  = source->mask;
    60     pmPeak *peak  = source->peak;
    61     psRegion srcRegion;
    62 
    63     // XXX EAM : psRegionXXX funcs need value not ptr
    64 
    65     srcRegion = psRegionForSquare (peak->x, peak->y, Radius);
    66     srcRegion = psRegionForImage (mask, &srcRegion);
    67 
    68     psImageMaskRegion (mask, &srcRegion, "OR", PSPHOT_MASK_MARKED);
    69     psStats *myStats = psStatsAlloc(statsOptions);
    70     myStats = psImageStats(myStats, image, mask, 0xff);
    71     psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED);
    72 
    73     psF64 tmpF64;
    74     p_psGetStatValue(myStats, &tmpF64);
    75     psFree (myStats);
    76 
    77     if (isnan(tmpF64)) return (false);
    78     source->moments = pmMomentsAlloc();
    79     source->moments->Sky = (psF32) tmpF64;
    80     return (true);
    81 }
    82 
    83 # define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
    84 # define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    85 bool pmSourceFitModel_EAM (pmSource *source,
    86                              pmModel *model,
    87                              const bool PSF)
    88 {
    89     // PS_PTR_CHECK_NULL(source, false);
    90     // PS_PTR_CHECK_NULL(source->moments, false);
    91     // PS_PTR_CHECK_NULL(source->peak, false);
    92     // PS_PTR_CHECK_NULL(source->pixels, false);
    93     // PS_PTR_CHECK_NULL(source->mask, false);
    94     // PS_PTR_CHECK_NULL(source->noise, false);
    95     psBool fitStatus = true;
    96     psBool onPic     = true;
    97     psBool rc        = true;
    98     psF32  Ro, ymodel;
    99 
    100 
    101     // XXX EAM : is it necessary for the mask & noise to exist?  the
    102     //           tests below could be conditions (!NULL)
    103 
    104     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    105     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    106 
    107     psVector *params = model->params;
    108     psVector *dparams = model->dparams;
    109     psVector *paramMask = NULL;
    110 
    111     psVector *beta_lim = NULL;
    112     psVector *params_min = NULL;
    113     psVector *params_max = NULL;
    114 
    115     int nParams = PSF ? params->n - 4 : params->n;
    116     psF32 So = params->data.F32[0];
    117 
    118     // find the number of valid pixels
    119     // XXX EAM : this loop and the loop below could just be one pass
    120     //           using the psArrayAdd and psVectorExtend functions
    121     psS32 count = 0;
    122     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    123         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    124             if (source->mask->data.U8[i][j] == 0) {
    125                 count++;
    126             }
    127         }
    128     }
    129     if (count <  nParams + 1) {
    130       psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    131       return(false);
    132     }
    133 
    134     // construct the coordinate and value entries
    135     psArray *x = psArrayAlloc(count);
    136     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    137     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    138     psS32 tmpCnt = 0;
    139     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    140         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    141             if (source->mask->data.U8[i][j] == 0) {
    142                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    143                 // XXX: Convert i/j to image space:
    144                 // XXX EAM: coord order is (x,y) == (col,row)
    145                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    146                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    147                 x->data[tmpCnt] = (psPtr *) coord;
    148                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    149 
    150                 ymodel = modelFunc (NULL, model->params, coord);
    151                
    152                 // this test enhances the noise based on deviation from the model flux
    153                 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
    154                 yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro);
    155                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
    156                 //           the minimization function calculates sq()
    157                 tmpCnt++;
    158             }
    159         }
    160     }
    161 
    162     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    163                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    164 
    165     // PSF model only fits first 4 parameters, FLT model fits all
    166     if (PSF) {
    167       paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    168       for (int i = 0; i < 4; i++) {
    169         paramMask->data.U8[i] = 0;
    170       }
    171       for (int i = 4; i < paramMask->n; i++) {
    172         paramMask->data.U8[i] = 1;
    173       }
    174     } 
    175 
    176     psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
    177     modelLimits (&beta_lim, &params_min, &params_max);
    178     for (int i = 0; i < params->n; i++) {
    179         covar->data.F64[0][i] = beta_lim->data.F32[i];
    180         covar->data.F64[1][i] = params_min->data.F32[i];
    181         covar->data.F64[2][i] = params_max->data.F32[i];
    182     }
    183 
    184     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    185     fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
    186     for (int i = 0; i < dparams->n; i++) {
    187         if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
    188         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    189     }
    190  
    191     // XXX EAM: we need to do something (give an error?) if rc is false
    192     // XXX EAM: psMinimizeLMChi2 does not check convergence
    193 
    194     // XXX models can go insane: reject these
    195     onPic &= (params->data.F32[2] >= source->pixels->col0);
    196     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    197     onPic &= (params->data.F32[3] >= source->pixels->row0);
    198     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    199 
    200     // XXX EAM: save the resulting chisq, nDOF, nIter
    201     model->chisq = myMin->value;
    202     model->nIter = myMin->iter;
    203     model->nDOF  = y->n - nParams;
    204 
    205     // XXX EAM get the Gauss-Newton distance for fixed model parameters
    206     if (paramMask != NULL) {
    207         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    208         psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
    209         for (int i = 0; i < dparams->n; i++) {
    210             if (!paramMask->data.U8[i]) continue;
    211             dparams->data.F32[i] = delta->data.F64[i];
    212         }
    213     }
    214 
    215     psFree(paramMask);
    216     psFree(x);
    217     psFree(y);
    218     psFree(myMin);
    219 
    220     rc = (onPic && fitStatus);
    221     return(rc);
    222 }
    223 
    224 /******************************************************************************
    225 pmSourceMoments(source, radius): this function takes a subImage defined in the
    226 pmSource data structure, along with the peak location, and determines the
    227 various moments associated with that peak.
    228  
    229 Requires the following to have been created:
    230     pmSource
    231     pmSource->peak
    232     pmSource->pixels
    233  
    234 XXX: The peak calculations are done in image coords, not subImage coords.
    235  
    236 XXX: mask values?
    237 
    238 XXX EAM : this version clips input pixels on S/N
    239 *****************************************************************************/
    240 # define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    241 
    242 bool pmSourceMoments_EAM(pmSource *source,
    243                          psF32 radius)
    244 {
    245     // PS_PTR_CHECK_NULL(source, NULL);
    246     // PS_PTR_CHECK_NULL(source->peak, NULL);
    247     // PS_PTR_CHECK_NULL(source->pixels, NULL);
    248     // PS_PTR_CHECK_NULL(source->noise, NULL);
    249     PS_FLOAT_COMPARE(0.0, radius, NULL);
    250 
    251     //
    252     // XXX: Verify the setting for sky if source->moments == NULL.
    253     //
    254     psF32 sky = 0.0;
    255     if (source->moments == NULL) {
    256         source->moments = pmMomentsAlloc();
    257     } else {
    258         sky = source->moments->Sky;
    259     }
    260 
    261     //
    262     // Sum = SUM (z - sky)
    263     // X1  = SUM (x - xc)*(z - sky)
    264     // X2  = SUM (x - xc)^2 * (z - sky)
    265     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    266     //
    267     psF32 peakPixel = -PS_MAX_F32;
    268     psS32 numPixels = 0;
    269     psF32 Sum = 0.0;
    270     psF32 X1 = 0.0;
    271     psF32 Y1 = 0.0;
    272     psF32 X2 = 0.0;
    273     psF32 Y2 = 0.0;
    274     psF32 XY = 0.0;
    275     psF32 x  = 0;
    276     psF32 y  = 0;
    277     psF32 R2 = PS_SQR(radius);
    278 
    279     psF32 xPeak = source->peak->x;
    280     psF32 yPeak = source->peak->y;
    281 
    282     // XXX why do I get different results for these two methods of finding Sx?
    283     // XXX Sx, Sy would be better measured if we clip pixels close to sky
    284     // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
    285     // We loop through all pixels in this subimage (source->pixels), and for each
    286     // pixel that is not masked, AND within the radius of the peak pixel, we
    287     // proceed with the moments calculation.  need to do two loops for a
    288     // numerically stable result.  first loop: get the sums.
    289     // XXX EAM : mask == 0 is valid
    290 
    291     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    292         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    293             if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue;
    294 
    295             psF32 xDiff = col + source->pixels->col0 - xPeak;
    296             psF32 yDiff = row + source->pixels->row0 - yPeak;
    297 
    298             if (!VALID_RADIUS(xDiff, yDiff)) continue;
    299 
    300             psF32 pDiff = source->pixels->data.F32[row][col] - sky;
    301 
    302             // XXX EAM : check for valid S/N in pixel
    303             if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue;
    304            
    305             Sum += pDiff;
    306             X1  += xDiff * pDiff;
    307             Y1  += yDiff * pDiff;
    308             XY  += xDiff * yDiff * pDiff;
    309            
    310             X2  += PS_SQR(xDiff) * pDiff;
    311             Y2  += PS_SQR(yDiff) * pDiff;
    312            
    313             peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
    314             numPixels++;
    315         }
    316     }
    317     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    318     if ((numPixels < 3) || (Sum <= 0)) {
    319       psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
    320       return (false);
    321     }
    322 
    323     psTrace (".psModules.pmSourceMoments", 5,
    324              "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
    325              sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    326 
    327     //
    328     // first moment X  = X1/Sum + xc
    329     // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    330     // Sxy             = XY / Sum
    331     //
    332     x = X1/Sum;
    333     y = Y1/Sum;
    334     if ((fabs(x) > radius) || (fabs(y) > radius)) {
    335       psTrace (".psModules.pmSourceMoments", 5,
    336                "large centroid swing; invalid peak %d, %d\n",
    337                source->peak->x, source->peak->y);
    338       return (false);
    339     }
    340 
    341     source->moments->x = x + xPeak;
    342     source->moments->y = y + yPeak;
    343 
    344     source->moments->Sxy = XY/Sum - x*y;
    345     source->moments->Sum = Sum;
    346     source->moments->Peak = peakPixel;
    347     source->moments->nPixels = numPixels;
    348 
    349     // XXX EAM : these values can be negative, so we need to limit the range
    350     source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    351     source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    352 
    353     psTrace (".psModules.pmSourceMoments", 4,
    354              "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
    355              sky, Sum, source->moments->x, source->moments->y,
    356              source->moments->Sx, source->moments->Sy, source->moments->Sxy);
    357 
    358     return(true);
    359 }
    360 
    36151bool pmModelFitStatus (pmModel *model) {
    36252

  • trunk/psphot/src/psModulesUtils.h

    r4983 r5048  
    1010
    1111// psModule extra utilities
    12 bool         pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius);
    13 bool         pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF);
    14 bool         pmSourceMoments_EAM(pmSource *source, psF32 radius);
     12// bool      pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF);
    1513bool         pmModelFitStatus (pmModel *model);
    1614int          pmSourceDophotType (pmSource *source);

  • trunk/psphot/src/psphotApplyPSF.c

    r4977 r5048  
    3939        // skip non-astronomical objects (very likely defects)
    4040        // XXX EAM : should we try these anyway?
    41         if (source->type == PS_SOURCE_DEFECT) continue;
    42         if (source->type == PS_SOURCE_SATURATED) continue;
     41        if (source->type == PM_SOURCE_DEFECT) continue;
     42        if (source->type == PM_SOURCE_SATURATED) continue;
    4343
    4444        // use the source moments, etc to guess basic model parameters
    … …  
    6565        }
    6666
    67         if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5);
     67        // if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5);
    6868           
    6969        // fit PSF model (set/unset the pixel mask)
    … …  
    7373        if (!status || (model->params->data.F32[1] < 0)) {
    7474          psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
    75           source->type = PS_SOURCE_FAIL_FIT_PSF;  // better choice?
     75          source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
    7676          psFree (model);
    7777          continue;

  • trunk/psphot/src/psphotChoosePSF.c

    r4977 r5048  
    1818    for (int i = 0; i < sources->n; i++) {
    1919        pmSource *source = sources->data[i];
    20         if (source->type != PS_SOURCE_PSFSTAR) continue;
     20        if (source->type != PM_SOURCE_PSFSTAR) continue;
    2121        psArrayAdd (stars, 200, source);
    2222    }

  • trunk/psphot/src/psphotFitGalaxies.c

    r4977 r5048  
    3131        // sources which should not be fitted
    3232        // skip all valid stars
    33         if (source->type == PS_SOURCE_PSFSTAR) continue;
    34         if (source->type == PS_SOURCE_SATSTAR) continue;
    35         if (source->type == PS_SOURCE_GOODSTAR) continue;
     33        if (source->type == PM_SOURCE_PSFSTAR) continue;
     34        if (source->type == PM_SOURCE_SATSTAR) continue;
     35        if (source->type == PM_SOURCE_GOODSTAR) continue;
    3636
    3737        // skip all likely defects
    38         if (source->type == PS_SOURCE_DEFECT) continue;
    39         if (source->type == PS_SOURCE_SATURATED) continue;
     38        if (source->type == PM_SOURCE_DEFECT) continue;
     39        if (source->type == PM_SOURCE_SATURATED) continue;
    4040
    4141        // skip poorly fitted stars
    42         if (source->type == PS_SOURCE_FAINTSTAR) continue;
    43         if (source->type == PS_SOURCE_POOR_FIT_PSF) continue;
     42        if (source->type == PM_SOURCE_FAINTSTAR) continue;
     43        if (source->type == PM_SOURCE_POOR_FIT_PSF) continue;
    4444
    4545        // XXX EAM when do we pick these up again?
    4646        if (source->moments->SN < GAL_MIN_SN) {
    47           source->type = PS_SOURCE_FAINT_GALAXY;  // better choice?
     47          source->type = PM_SOURCE_FAINT_GALAXY;  // better choice?
    4848          continue;
    4949        }
    5050
    5151        // recalculate the source moments using the larger galaxy moments radius
    52         // XXX EAM : 0.5.0 code used _EAM version
    53         // status = pmSourceMoments (source, GAL_MOMENTS_RAD);
    54         status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD);
     52        status = pmSourceMoments (source, GAL_MOMENTS_RAD);
    5553        if (!status) {
    56           source->type = PS_SOURCE_DROP_GALAXY;  // better choice?
     54          source->type = PM_SOURCE_DROP_GALAXY;  // better choice?
    5755          continue;
    5856        }
    … …  
    8482          // if the fit fails, we need to change the classification
    8583          psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
    86           source->type = PS_SOURCE_FAIL_FIT_GAL;  // better choice?
     84          source->type = PM_SOURCE_FAIL_FIT_GAL;  // better choice?
    8785          source->modelFLT = model;
    8886          Nfail ++;
    … …  
    9593          // if the fit fails, we need to change the classification
    9694          psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
    97           source->type = PS_SOURCE_POOR_FIT_GAL;  // better choice?
     95          source->type = PM_SOURCE_POOR_FIT_GAL;  // better choice?
    9896          source->modelFLT = model;
    9997          Nfail ++;
    … …  
    102100
    103101        // accept the model
    104         source->type = PS_SOURCE_GALAXY;
     102        source->type = PM_SOURCE_GALAXY;
    105103        source->modelFLT = model;
    106104        Niter += model[0].nIter;

  • trunk/psphot/src/psphotMarkPSF.c

    r4949 r5048  
    1818// we also reject objects with S/N too low or ChiSquare to high
    1919
    20 // any object which meets the criterion is marked as PS_SOURCE_BRIGHTSTAR
     20// any object which meets the criterion is marked as PM_SOURCE_BRIGHTSTAR
    2121
    2222// floor for DS value
    … …  
    3535    // remember: fit does not use saturated pixels (masked)
    3636    if (source->modelPSF->params->data.F32[1] >= SATURATE) {
    37         if (source->type == PS_SOURCE_PSFSTAR) {
     37        if (source->type == PM_SOURCE_PSFSTAR) {
    3838            psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n");
    3939        }
    40         source->type = PS_SOURCE_SATSTAR;
     40        source->type = PM_SOURCE_SATSTAR;
    4141        return (true);
    4242    }
    43     if (source->type == PS_SOURCE_SATSTAR) {
     43    if (source->type == PM_SOURCE_SATSTAR) {
    4444        psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
    45         source->type = PS_SOURCE_GOODSTAR;
     45        source->type = PM_SOURCE_GOODSTAR;
    4646    }
    4747
    … …  
    5757    nSy = dSY / hypot (MIN_DS, 1 / (SY * SN));
    5858
    59     // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[]
     59    // assign PM_SOURCE_GOODSTAR to bright objects within PSF region of dparams[]
    6060    keep = TRUE;
    6161    keep &= (fabs(nSx) < shapeNsigma);
    … …  
    6464    keep &= (Chi < maxChi);
    6565    if (keep) {
    66         if (source->type == PS_SOURCE_PSFSTAR) return (true);
    67         source->type = PS_SOURCE_GOODSTAR;
     66        if (source->type == PM_SOURCE_PSFSTAR) return (true);
     67        source->type = PM_SOURCE_GOODSTAR;
    6868        return (true);
    6969    }
    7070   
    71     if (source->type == PS_SOURCE_PSFSTAR) {
     71    if (source->type == PM_SOURCE_PSFSTAR) {
    7272        psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n");
    7373    }
    … …  
    7575    // object appears to be small, suspected defect
    7676    if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
    77         source->type = PS_SOURCE_DEFECT;
     77        source->type = PM_SOURCE_DEFECT;
    7878        return (false);
    7979    }
    … …  
    8181    // object appears to be large, suspected galaxy
    8282    if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
    83         source->type = PS_SOURCE_GALAXY;
     83        source->type = PM_SOURCE_GALAXY;
    8484        return (false);
    8585    }
    … …  
    8787    // object appears to be extremely faint: what is this?
    8888    if (SN < minSN) {
    89         source->type = PS_SOURCE_FAINTSTAR;
     89        source->type = PM_SOURCE_FAINTSTAR;
    9090        return (false);
    9191    }
    9292
    9393    // these are pooly fitted, probable stars near other stars?
    94     source->type = PS_SOURCE_POOR_FIT_PSF;
     94    source->type = PM_SOURCE_POOR_FIT_PSF;
    9595    return (false);
    9696}       

  • trunk/psphot/src/psphotOutput.c

    r4954 r5048  
    9191
    9292            // use PSF model with stars
    93           case PS_SOURCE_PSFSTAR:
    94           case PS_SOURCE_SATSTAR:
    95           case PS_SOURCE_GOODSTAR:
     93          case PM_SOURCE_PSFSTAR:
     94          case PM_SOURCE_SATSTAR:
     95          case PM_SOURCE_GOODSTAR:
    9696            model = source->modelPSF;
    9797            break;
    9898
    9999            // use FLT model with galaxies
    100           case PS_SOURCE_GALAXY:
     100          case PM_SOURCE_GALAXY:
    101101            model = source->modelFLT;
    102102            break;
    103103           
    104104            // skip defects and poorly fitted stars or galaxies
    105           case PS_SOURCE_DEFECT:
    106           case PS_SOURCE_SATURATED:
    107           case PS_SOURCE_FAINTSTAR:
    108           case PS_SOURCE_POOR_FIT_PSF:
    109           case PS_SOURCE_POOR_FIT_GAL:
    110           case PS_SOURCE_FAIL_FIT_GAL:
     105          case PM_SOURCE_DEFECT:
     106          case PM_SOURCE_SATURATED:
     107          case PM_SOURCE_FAINTSTAR:
     108          case PM_SOURCE_POOR_FIT_PSF:
     109          case PM_SOURCE_POOR_FIT_GAL:
     110          case PM_SOURCE_FAIL_FIT_GAL:
    111111          default:
    112112            model = NULL;
    … …  
    183183
    184184            // use PSF model with stars
    185           case PS_SOURCE_SATSTAR:
     185          case PM_SOURCE_SATSTAR:
    186186            flags = 4;
    187187            model = source->modelPSF;
    188188            break;
    189189
    190           case PS_SOURCE_PSFSTAR:
    191           case PS_SOURCE_GOODSTAR:
     190          case PM_SOURCE_PSFSTAR:
     191          case PM_SOURCE_GOODSTAR:
    192192            model = source->modelPSF;
    193193            break;
    194194
    195195            // use FLT model with galaxies
    196           case PS_SOURCE_GALAXY:
     196          case PM_SOURCE_GALAXY:
    197197            model = source->modelFLT;
    198198            break;
    199199           
    200200            // skip defects and poorly fitted stars or galaxies
    201           case PS_SOURCE_DEFECT:
    202           case PS_SOURCE_SATURATED:
    203           case PS_SOURCE_FAINTSTAR:
    204           case PS_SOURCE_POOR_FIT_PSF:
    205           case PS_SOURCE_POOR_FIT_GAL:
    206           case PS_SOURCE_FAIL_FIT_GAL:
     201          case PM_SOURCE_DEFECT:
     202          case PM_SOURCE_SATURATED:
     203          case PM_SOURCE_FAINTSTAR:
     204          case PM_SOURCE_POOR_FIT_PSF:
     205          case PM_SOURCE_POOR_FIT_GAL:
     206          case PM_SOURCE_FAIL_FIT_GAL:
    207207          default:
    208208            model = NULL;
    … …  
    288288
    289289            // use PSF model with stars
    290           case PS_SOURCE_PSFSTAR:
    291           case PS_SOURCE_SATSTAR:
    292           case PS_SOURCE_GOODSTAR:
     290          case PM_SOURCE_PSFSTAR:
     291          case PM_SOURCE_SATSTAR:
     292          case PM_SOURCE_GOODSTAR:
    293293            model = source->modelPSF;
    294294            break;
    295295
    296296            // use FLT model with galaxies
    297           case PS_SOURCE_GALAXY:
     297          case PM_SOURCE_GALAXY:
    298298            model = source->modelFLT;
    299299            break;
    300300           
    301301            // skip defects and poorly fitted stars or galaxies
    302           case PS_SOURCE_DEFECT:
    303           case PS_SOURCE_SATURATED:
    304           case PS_SOURCE_FAINTSTAR:
    305           case PS_SOURCE_POOR_FIT_PSF:
    306           case PS_SOURCE_POOR_FIT_GAL:
    307           case PS_SOURCE_FAIL_FIT_GAL:
     302          case PM_SOURCE_DEFECT:
     303          case PM_SOURCE_SATURATED:
     304          case PM_SOURCE_FAINTSTAR:
     305          case PM_SOURCE_POOR_FIT_PSF:
     306          case PM_SOURCE_POOR_FIT_GAL:
     307          case PM_SOURCE_FAIL_FIT_GAL:
    308308          default:
    309309            model = NULL;
    … …  
    445445
    446446        // valid source types for this function
    447         if (source->type == PS_SOURCE_SATSTAR) goto valid;
    448         if (source->type == PS_SOURCE_PSFSTAR) goto valid;
    449         if (source->type == PS_SOURCE_GOODSTAR) goto valid;
     447        if (source->type == PM_SOURCE_SATSTAR) goto valid;
     448        if (source->type == PM_SOURCE_PSFSTAR) goto valid;
     449        if (source->type == PM_SOURCE_GOODSTAR) goto valid;
    450450        continue;
    451451
    … …  
    509509        model = (pmModel  *) source->modelFLT;
    510510        if (model == NULL) continue;
    511         if (source->type == PS_SOURCE_GALAXY) goto valid;
     511        if (source->type == PM_SOURCE_GALAXY) goto valid;
    512512        continue;
    513513
    … …  
    570570
    571571        // skip these sources (in PSF or FLT)
    572         if (source->type == PS_SOURCE_GALAXY) continue;
    573         if (source->type == PS_SOURCE_PSFSTAR) continue;
    574         if (source->type == PS_SOURCE_SATSTAR) continue;
    575         if (source->type == PS_SOURCE_GOODSTAR) continue;
     572        if (source->type == PM_SOURCE_GALAXY) continue;
     573        if (source->type == PM_SOURCE_PSFSTAR) continue;
     574        if (source->type == PM_SOURCE_SATSTAR) continue;
     575        if (source->type == PM_SOURCE_GOODSTAR) continue;
    576576
    577577        if (source->moments == NULL) {
    … …  
    602602    switch (source->type) {
    603603
    604       case PS_SOURCE_DEFECT:
    605       case PS_SOURCE_SATURATED:
     604      case PM_SOURCE_DEFECT:
     605      case PM_SOURCE_SATURATED:
    606606        return (8);
    607607
    608       case PS_SOURCE_SATSTAR:
     608      case PM_SOURCE_SATSTAR:
    609609        return (10);
    610610
    611       case PS_SOURCE_PSFSTAR:
    612       case PS_SOURCE_GOODSTAR:
     611      case PM_SOURCE_PSFSTAR:
     612      case PM_SOURCE_GOODSTAR:
    613613        return (1);
    614614
    615       case PS_SOURCE_POOR_FIT_PSF:
     615      case PM_SOURCE_POOR_FIT_PSF:
    616616        return (7);
    617617
    618       case PS_SOURCE_FAIL_FIT_PSF:
    619       case PS_SOURCE_FAINTSTAR:
     618      case PM_SOURCE_FAIL_FIT_PSF:
     619      case PM_SOURCE_FAINTSTAR:
    620620        return (4);
    621621
    622       case PS_SOURCE_GALAXY:
    623       case PS_SOURCE_FAINT_GALAXY:
    624       case PS_SOURCE_DROP_GALAXY:
    625       case PS_SOURCE_FAIL_FIT_GAL:
    626       case PS_SOURCE_POOR_FIT_GAL:
     622      case PM_SOURCE_GALAXY:
     623      case PM_SOURCE_FAINT_GALAXY:
     624      case PM_SOURCE_DROP_GALAXY:
     625      case PM_SOURCE_FAIL_FIT_GAL:
     626      case PM_SOURCE_POOR_FIT_GAL:
    627627        return (2);
    628628
    629       case PS_SOURCE_OTHER:
     629      case PM_SOURCE_OTHER:
    630630      default:
    631631        return (0);

  • trunk/psphot/src/psphotSourceStats.c

    r4977 r5048  
    2828        // measure a local sky value
    2929        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
    30         status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     30        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
    3131        if (!status) {
    3232          psFree (source);
    … …  
    3636        // measure basic source moments
    3737        // XXX EAM : choose between these two versions
    38         status = pmSourceMoments_EAM (source, RADIUS);
     38        status = pmSourceMoments (source, RADIUS);
    3939        if (!status) {
    4040          psFree (source);

  • trunk/psphot/src/psphotSubtractPSF.c

    r4949 r5048  
    1111
    1212  // only subtract successful fits
    13   if (source->type == PS_SOURCE_SATSTAR) goto valid;
    14   if (source->type == PS_SOURCE_PSFSTAR) goto valid;
    15   if (source->type == PS_SOURCE_GOODSTAR) goto valid;
     13  if (source->type == PM_SOURCE_SATSTAR) goto valid;
     14  if (source->type == PM_SOURCE_PSFSTAR) goto valid;
     15  if (source->type == PM_SOURCE_GOODSTAR) goto valid;
    1616  return (false);
    1717
Note: See TracChangeset for help on using the changeset viewer.