IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 9, 2006, 6:46:03 PM (20 years ago)
Author:
magnier
Message:

updated objects code with ApTrend concept (was in eam_rel9_b0)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_rel9_p0/psModules/src/objects/pmObjects.c

    r5765 r5958  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-12 21:14:38 $
     8 *  @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-01-10 04:46:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    271271    tmp->chisq = 0.0;
    272272    tmp->nIter = 0;
     273    tmp->radius = 0;
     274    tmp->status = PM_MODEL_UNTRIED;
     275
    273276    psS32 Nparams = pmModelParameterCount(type);
    274277    if (Nparams == 0) {
    … …  
    307310    tmp->mask = NULL;
    308311    tmp->moments = NULL;
     312    tmp->blends = NULL;
    309313    tmp->modelPSF = NULL;
    310314    tmp->modelFLT = NULL;
    311     tmp->type = 0;
     315    tmp->type = PM_SOURCE_UNKNOWN;
     316    tmp->mode = PM_SOURCE_DEFAULT;
    312317    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
    313318
    … …  
    827832    psS32 numPixels = 0;
    828833    psF32 Sum = 0.0;
     834    psF32 Var = 0.0;
    829835    psF32 X1 = 0.0;
    830836    psF32 Y1 = 0.0;
    … …  
    850856    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    851857        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    852             if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
     858            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
    853859                continue;
     860            }
    854861
    855862            psF32 xDiff = col + source->pixels->col0 - xPeak;
    … …  
    858865            // XXX EAM : calculate xDiff, yDiff up front;
    859866            //           radius is just a function of (xDiff, yDiff)
    860             if (!VALID_RADIUS(xDiff, yDiff, R2))
     867            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
    861868                continue;
     869            }
    862870
    863871            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     872            psF32 wDiff = source->weight->data.F32[row][col];
    864873
    865874            // XXX EAM : check for valid S/N in pixel
    866875            // XXX EAM : should this limit be user-defined?
    867             if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1)
     876            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
    868877                continue;
    869 
     878            }
     879
     880            Var += wDiff;
    870881            Sum += pDiff;
    871882            X1  += xDiff * pDiff;
    … …  
    880891        }
    881892    }
     893
     894    // if we have less than (1/4) of the possible pixels, force a retry
    882895    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    883     if ((numPixels < 3) || (Sum <= 0)) {
    884         psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     896    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
     897        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
    885898        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    886899        return (false);
    … …  
    899912    y = Y1/Sum;
    900913    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    901         psTrace (".psModules.pmSourceMoments", 5,
     914        psTrace (".psModules.pmSourceMoments", 3,
    902915                 "large centroid swing; invalid peak %d, %d\n",
    903916                 source->peak->x, source->peak->y);
    … …  
    912925    source->moments->Sxy = XY/Sum - x*y;
    913926    source->moments->Sum = Sum;
     927    source->moments->SN  = Sum / sqrt(Var);
    914928    source->moments->Peak = peakPixel;
    915929    source->moments->nPixels = numPixels;
    … …  
    9941008        psImage *splane = NULL;
    9951009        int binX, binY;
     1010        bool status;
     1011
     1012        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
     1013        if (!status)
     1014            SX_MAX = 10.0;
     1015        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
     1016        if (!status)
     1017            SY_MAX = 10.0;
    9961018
    9971019        // construct a sigma-plane image
    9981020        // psImageAlloc does zero the data
    999         splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
     1021        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
    10001022        for (int i = 0; i < splane->numRows; i++)
    10011023        {
    … …  
    11311153XXX: How can this function ever return FALSE?
    11321154 
    1133 XXX EAM : add the saturated mask value to metadata
     1155EAM: I moved S/N calculation to pmSourceMoments, using weight image
    11341156*****************************************************************************/
    11351157
    … …  
    11461168    int Ncr      = 0;
    11471169    int Nsatstar = 0;
    1148 
     1170    // psRegion allArray = psRegionSet (0, 0, 0, 0);
     1171    psRegion inner;
     1172
     1173    // report stats on S/N values for star-like objects
    11491174    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    11501175    starsn->n = 0;
    … …  
    11521177    // check return status value (do these exist?)
    11531178    bool status;
    1154     psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
    1155     psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
    11561179    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    1157     // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
    11581180
    11591181    // XXX allow clump size to be scaled relative to sigmas?
    … …  
    11681190        psF32 sigY = tmpSrc->moments->Sy;
    11691191
    1170         // calculate and save signal-to-noise estimates
    1171         psF32 S  = tmpSrc->moments->Sum;
    1172         psF32 A  = 4 * M_PI * sigX * sigY;
    1173         psF32 B  = tmpSrc->moments->Sky;
    1174         psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
    1175         psF32 SN = (S * sqrt(GAIN) / RT);
    1176         tmpSrc->moments->SN = SN;
    1177 
    11781192        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1179         //
    1180         // XXX: Must verify this region (the region argument was added to psImageCountPixelMask()
    1181         // after EAM wrote this code.
    1182         //
    1183         psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows);
    1184         int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED);
     1193        inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
     1194        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
    11851195
    11861196        // saturated star (size consistent with PSF or larger)
    11871197        // Nsigma should be user-configured parameter
    11881198        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
     1199        big = true;
    11891200        if ((Nsatpix > 1) && big) {
    1190             tmpSrc->type = PM_SOURCE_SATSTAR;
     1201            tmpSrc->type = PM_SOURCE_STAR;
     1202            tmpSrc->mode = PM_SOURCE_SATSTAR;
    11911203            Nsatstar ++;
    11921204            continue;
    … …  
    11961208        if (Nsatpix > 1) {
    11971209            tmpSrc->type = PM_SOURCE_SATURATED;
     1210            tmpSrc->mode = PM_SOURCE_DEFAULT;
    11981211            Nsat ++;
    11991212            continue;
    … …  
    12051218        if ((sigX < 0.05) || (sigY < 0.05)) {
    12061219            tmpSrc->type = PM_SOURCE_DEFECT;
     1220            tmpSrc->mode = PM_SOURCE_DEFAULT;
    12071221            Ncr ++;
    12081222            continue;
    … …  
    12121226        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    12131227            tmpSrc->type = PM_SOURCE_GALAXY;
     1228            tmpSrc->mode = PM_SOURCE_DEFAULT;
    12141229            Ngal ++;
    12151230            continue;
    … …  
    12171232
    12181233        // the rest are probable stellar objects
    1219         starsn->data.F32[starsn->n] = SN;
     1234        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
    12201235        starsn->n ++;
    12211236        Nstar ++;
    12221237
    1223         // PSF star (within 1.5 sigma of clump center
     1238        // PSF star (within 1.5 sigma of clump center, S/N > limit)
    12241239        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    1225         if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    1226             tmpSrc->type = PM_SOURCE_PSFSTAR;
     1240        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
     1241            tmpSrc->type = PM_SOURCE_STAR;
     1242            tmpSrc->mode = PM_SOURCE_PSFSTAR;
    12271243            Npsf ++;
    12281244            continue;
    … …  
    12301246
    12311247        // random type of star
    1232         tmpSrc->type = PM_SOURCE_OTHER;
     1248        tmpSrc->type = PM_SOURCE_STAR;
     1249        tmpSrc->mode = PM_SOURCE_DEFAULT;
    12331250    }
    12341251
    … …  
    15511568#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
    15521569#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    1553 /******************************************************************************
    1554 pmSourceFitModel(source, model): must create the appropiate arguments to the
    1555 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
    1556  
    1557 XXX: should there be a mask value?
    1558 XXX EAM : fit the specified model (not necessarily the one in source)
    1559 *****************************************************************************/
    1560 bool pmSourceFitModel_v5(pmSource *source,
    1561                          pmModel *model,
    1562                          const bool PSF)
     1570bool pmSourceFitModel (pmSource *source,
     1571                       pmModel *model,
     1572                       const bool PSF)
    15631573{
    15641574    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    … …  
    15691579    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    15701580    PS_ASSERT_PTR_NON_NULL(source->weight, false);
     1581
     1582    // XXX EAM : is it necessary for the mask & weight to exist?  the
     1583    //           tests below could be conditions (!NULL)
     1584
    15711585    psBool fitStatus = true;
    15721586    psBool onPic     = true;
    15731587    psBool rc        = true;
    15741588
    1575     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1576     //           tests below could be conditions (!NULL)
    1577 
    15781589    psVector *params = model->params;
    15791590    psVector *dparams = model->dparams;
    … …  
    15821593    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    15831594
    1584     int nParams = PSF ? params->n - 4 : params->n;
    1585 
    1586     // find the number of valid pixels
    1587     // XXX EAM : this loop and the loop below could just be one pass
    1588     //           using the psArrayAdd and psVectorExtend functions
    1589     psS32 count = 0;
     1595    int nParams = PSF ? 4 : params->n;
     1596
     1597    // maximum number of valid pixels
     1598    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
     1599
     1600    // construct the coordinate and value entries
     1601    psArray *x = psArrayAlloc(nPix);
     1602    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
     1603    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
     1604
     1605    nPix = 0;
    15901606    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    15911607        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1592             if (source->mask->data.U8[i][j] == 0) {
    1593                 count++;
    1594             }
    1595         }
    1596     }
    1597     if (count <  nParams + 1) {
     1608            if (source->mask->data.U8[i][j]) {
     1609                continue;
     1610            }
     1611            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1612
     1613            // Convert i/j to image space:
     1614            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1615            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     1616            x->data[nPix] = (psPtr *) coord;
     1617            y->data.F32[nPix] = source->pixels->data.F32[i][j];
     1618
     1619            // psMinimizeLMChi2 takes wt = 1/dY^2
     1620            if (source->weight->data.F32[i][j] == 0) {
     1621                continue;
     1622            }
     1623            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     1624            nPix++;
     1625        }
     1626    }
     1627    if (nPix <  nParams + 1) {
    15981628        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    15991629        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     1630        model->status = PM_MODEL_BADARGS;
     1631        psFree (x);
     1632        psFree (y);
     1633        psFree (yErr);
    16001634        return(false);
    16011635    }
    1602 
    1603     // construct the coordinate and value entries
    1604     psArray *x = psArrayAlloc(count);
    1605     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1606     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1607     psS32 tmpCnt = 0;
    1608     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1609         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1610             if (source->mask->data.U8[i][j] == 0) {
    1611                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1612                 // XXX: Convert i/j to image space:
    1613                 // XXX EAM: coord order is (x,y) == (col,row)
    1614                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1615                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1616                 x->data[tmpCnt] = (psPtr *) coord;
    1617                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1618                 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
    1619                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
    1620                 //           the minimization function calculates sq()
    1621                 tmpCnt++;
    1622             }
    1623         }
    1624     }
     1636    x->n = nPix;
     1637    y->n = nPix;
     1638    yErr->n = nPix;
    16251639
    16261640    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    … …  
    16381652    }
    16391653
    1640     // XXX EAM : covar must be F64?
    1641     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    1642 
    1643     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    1644     fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
    1645     for (int i = 0; i < dparams->n; i++) {
    1646         if ((paramMask != NULL) && paramMask->data.U8[i])
    1647             continue;
    1648         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    1649     }
    1650 
    1651     // XXX EAM: we need to do something (give an error?) if rc is false
    1652     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1653 
    1654     // XXX models can go insane: reject these
    1655     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1656     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1657     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1658     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1659 
    1660     // XXX EAM: save the resulting chisq, nDOF, nIter
    1661     model->chisq = myMin->value;
    1662     model->nIter = myMin->iter;
    1663     model->nDOF  = y->n - nParams;
    1664 
    1665     // XXX EAM get the Gauss-Newton distance for fixed model parameters
    1666     if (paramMask != NULL) {
    1667         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    1668         psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
    1669         for (int i = 0; i < dparams->n; i++) {
    1670             if (!paramMask->data.U8[i])
    1671                 continue;
    1672             dparams->data.F32[i] = delta->data.F64[i];
    1673         }
    1674         psFree (delta);
    1675     }
    1676 
    1677     psFree(x);
    1678     psFree(y);
    1679     psFree(yErr);
    1680     psFree(myMin);
    1681     psFree(covar);
    1682     psFree(paramMask);
    1683 
    1684     rc = (onPic && fitStatus);
    1685     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1686     return(rc);
    1687 }
    1688 
    1689 // XXX EAM : new version with parameter range limits and weight enhancement
    1690 bool pmSourceFitModel (pmSource *source,
    1691                        pmModel *model,
    1692                        const bool PSF)
    1693 {
    1694     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1695     PS_ASSERT_PTR_NON_NULL(source, false);
    1696     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1697     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1698     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1699     PS_ASSERT_PTR_NON_NULL(source->mask, false);
    1700     PS_ASSERT_PTR_NON_NULL(source->weight, false);
    1701 
    1702     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1703     //           tests below could be conditions (!NULL)
    1704 
    1705     psBool fitStatus = true;
    1706     psBool onPic     = true;
    1707     psBool rc        = true;
    1708     psF32  Ro, ymodel;
    1709 
    1710     psVector *params = model->params;
    1711     psVector *dparams = model->dparams;
    1712     psVector *paramMask = NULL;
    1713 
    1714     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1715 
    1716     // XXX EAM : I need to use the sky value to constrain the weight model
    1717     int nParams = PSF ? params->n - 4 : params->n;
    1718     psF32 So = params->data.F32[0];
    1719 
    1720     // find the number of valid pixels
    1721     // XXX EAM : this loop and the loop below could just be one pass
    1722     //           using the psArrayAdd and psVectorExtend functions
    1723     psS32 count = 0;
    1724     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1725         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1726             if (source->mask->data.U8[i][j] == 0) {
    1727                 count++;
    1728             }
    1729         }
    1730     }
    1731     if (count <  nParams + 1) {
    1732         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    1733         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    1734         return(false);
    1735     }
    1736 
    1737     // construct the coordinate and value entries
    1738     psArray *x = psArrayAlloc(count);
    1739     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1740     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1741     psS32 tmpCnt = 0;
    1742     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1743         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1744             if (source->mask->data.U8[i][j] == 0) {
    1745                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1746                 // XXX: Convert i/j to image space:
    1747                 // XXX EAM: coord order is (x,y) == (col,row)
    1748                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1749                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1750                 x->data[tmpCnt] = (psPtr *) coord;
    1751                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1752 
    1753                 // compare observed flux to model flux to adjust weight
    1754                 ymodel = modelFunc (NULL, model->params, coord);
    1755 
    1756                 // this test enhances the weight based on deviation from the model flux
    1757                 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
    1758 
    1759                 // psMinimizeLMChi2_EAM takes wt = 1/dY^2
    1760                 if (source->weight->data.F32[i][j] == 0) {
    1761                     yErr->data.F32[tmpCnt] = 0.0;
    1762                 } else {
    1763                     yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
    1764                 }
    1765                 tmpCnt++;
    1766             }
    1767         }
    1768     }
    1769 
    1770     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    1771                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    1772 
    1773     // PSF model only fits first 4 parameters, FLT model fits all
    1774     if (PSF) {
    1775         paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    1776         for (int i = 0; i < 4; i++) {
    1777             paramMask->data.U8[i] = 0;
    1778         }
    1779         for (int i = 4; i < paramMask->n; i++) {
    1780             paramMask->data.U8[i] = 1;
    1781         }
    1782     }
    1783 
    1784     // XXX EAM : I've added three types of parameter range checks
    1785     // XXX EAM : this requires my new psMinimization functions
     1654    // Set the parameter range checks
    17861655    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    17871656    psVector *beta_lim = NULL;
    … …  
    18071676    }
    18081677
    1809     // XXX EAM: we need to do something (give an error?) if rc is false
    1810     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1811 
    1812     // XXX models can go insane: reject these
    1813     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1814     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1815     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1816     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1817 
    1818     // XXX EAM: save the resulting chisq, nDOF, nIter
     1678    // save the resulting chisq, nDOF, nIter
    18191679    model->chisq = myMin->value;
    18201680    model->nIter = myMin->iter;
    18211681    model->nDOF  = y->n - nParams;
    18221682
    1823     // XXX EAM get the Gauss-Newton distance for fixed model parameters
     1683    // get the Gauss-Newton distance for fixed model parameters
    18241684    if (paramMask != NULL) {
    18251685        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    … …  
    18321692    }
    18331693
    1834     psFree(paramMask);
     1694    // set the model success or failure status
     1695    if (!fitStatus) {
     1696        model->status = PM_MODEL_NONCONVERGE;
     1697    } else {
     1698        model->status = PM_MODEL_SUCCESS;
     1699    }
     1700
     1701    // models can go insane: reject these
     1702    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1703    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1704    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1705    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1706    if (!onPic) {
     1707        model->status = PM_MODEL_OFFIMAGE;
     1708    }
     1709
     1710    source->mode |= PM_SOURCE_FITTED;
     1711
    18351712    psFree(x);
    18361713    psFree(y);
     1714    psFree(yErr);
    18371715    psFree(myMin);
     1716    psFree(covar);
     1717    psFree(paramMask);
    18381718
    18391719    rc = (onPic && fitStatus);
    … …  
    18461726                             pmModel *model,
    18471727                             bool center,
    1848                              psS32 flag)
     1728                             bool sky,
     1729                             bool add
     1730                                )
    18491731{
    18501732    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    … …  
    18591741    psS32 imageCol;
    18601742    psS32 imageRow;
     1743    psF32 skyValue = params->data.F32[0];
     1744    psF32 pixelValue;
    18611745
    18621746    for (psS32 i = 0; i < image->numRows; i++) {
    … …  
    18641748            if ((mask != NULL) && mask->data.U8[i][j])
    18651749                continue;
    1866             psF32 pixelValue;
     1750
    18671751            // XXX: Should you be adding the pixels for the entire subImage,
    18681752            // or a radius of pixels around it?
    … …  
    18821766            x->data.F32[0] = (float) imageCol;
    18831767            x->data.F32[1] = (float) imageRow;
    1884             pixelValue = modelFunc (NULL, params, x);
    1885             // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
    1886 
    1887             if (flag == 1) {
    1888                 pixelValue = -pixelValue;
    1889             }
    1890 
    1891             // XXX: Must figure out how to calculate the image coordinates and
    1892             // how to use the boolean "center" flag.
    1893 
    1894             image->data.F32[i][j]+= pixelValue;
     1768
     1769            // set the appropriate pixel value for this coordinate
     1770            if (sky) {
     1771                pixelValue = modelFunc (NULL, params, x);
     1772            } else {
     1773                pixelValue = modelFunc (NULL, params, x) - skyValue;
     1774            }
     1775
     1776
     1777            // add or subtract the value
     1778            if (add
     1779               ) {
     1780                image->data.F32[i][j] += pixelValue;
     1781            }
     1782            else {
     1783                image->data.F32[i][j] -= pixelValue;
     1784            }
    18951785        }
    18961786    }
    … …  
    19071797                      psImage *mask,
    19081798                      pmModel *model,
    1909                       bool center)
    1910 {
    1911     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1912     psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
     1799                      bool center,
     1800                      bool sky)
     1801{
     1802    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1803    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
    19131804    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    19141805    return(rc);
    … …  
    19201811                      psImage *mask,
    19211812                      pmModel *model,
    1922                       bool center)
    1923 {
    1924     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1925     psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
     1813                      bool center,
     1814                      bool sky)
     1815{
     1816    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1817    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
    19261818    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    19271819    return(rc);
Note: See TracChangeset for help on using the changeset viewer.