IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5829


Ignore:
Timestamp:
Dec 22, 2005, 4:23:00 PM (21 years ago)
Author:
magnier
Message:

cleanups, dropped pmSourceFitModel_v5, dropped double loop in pmSourceFitModel, skip masked points

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5806 r5829  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.3.4.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-17 20:41:29 $
     8 *  @version $Revision: 1.3.4.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-12-23 02:23:00 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    132132    tmp->mask = NULL;
    133133    tmp->moments = NULL;
     134    tmp->blends = NULL;
    134135    tmp->modelPSF = NULL;
    135136    tmp->modelFLT = NULL;
    … …  
    11151116        // Nsigma should be user-configured parameter
    11161117        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
     1118        big = true;
    11171119        if ((Nsatpix > 1) && big) {
    11181120            tmpSrc->type = PM_SOURCE_STAR;
    … …  
    15711573#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
    15721574#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    1573 /******************************************************************************
    1574 pmSourceFitModel(source, model): must create the appropiate arguments to the
    1575 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
    1576  
    1577 XXX: should there be a mask value?
    1578 XXX EAM : fit the specified model (not necessarily the one in source)
    1579 *****************************************************************************/
    1580 bool pmSourceFitModel_v5(pmSource *source,
    1581                          pmModel *model,
    1582                          const bool PSF)
     1575bool pmSourceFitModel (pmSource *source,
     1576                       pmModel *model,
     1577                       const bool PSF)
    15831578{
    15841579    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    … …  
    15891584    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    15901585    PS_ASSERT_PTR_NON_NULL(source->weight, false);
     1586
     1587    // XXX EAM : is it necessary for the mask & weight to exist?  the
     1588    //           tests below could be conditions (!NULL)
     1589
    15911590    psBool fitStatus = true;
    15921591    psBool onPic     = true;
    15931592    psBool rc        = true;
    15941593
    1595     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1596     //           tests below could be conditions (!NULL)
    1597 
    15981594    psVector *params = model->params;
    15991595    psVector *dparams = model->dparams;
    … …  
    16021598    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    16031599
    1604     int nParams = PSF ? params->n - 4 : params->n;
    1605 
    1606     // find the number of valid pixels
    1607     // XXX EAM : this loop and the loop below could just be one pass
    1608     //           using the psArrayAdd and psVectorExtend functions
    1609     psS32 count = 0;
     1600    int nParams = PSF ? 4 : params->n;
     1601
     1602    // maximum number of valid pixels
     1603    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
     1604
     1605    // construct the coordinate and value entries
     1606    psArray *x = psArrayAlloc(nPix);
     1607    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
     1608    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
     1609
     1610    nPix = 0;
    16101611    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    16111612        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1612             if (source->mask->data.U8[i][j] == 0) {
    1613                 count++;
    1614             }
    1615         }
    1616     }
    1617     if (count <  nParams + 1) {
     1613            if (source->mask->data.U8[i][j]) {
     1614                continue;
     1615            }
     1616            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1617
     1618            // Convert i/j to image space:
     1619            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1620            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     1621            x->data[nPix] = (psPtr *) coord;
     1622            y->data.F32[nPix] = source->pixels->data.F32[i][j];
     1623
     1624            // psMinimizeLMChi2 takes wt = 1/dY^2
     1625            if (source->weight->data.F32[i][j] == 0) {
     1626                continue;
     1627            }
     1628            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     1629            nPix++;
     1630        }
     1631    }
     1632    if (nPix <  nParams + 1) {
    16181633        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    16191634        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    16201635        model->status = PM_MODEL_BADARGS;
     1636        psFree (x);
     1637        psFree (y);
     1638        psFree (yErr);
    16211639        return(false);
    16221640    }
    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                 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
    1640                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
    1641                 //           the minimization function calculates sq()
    1642                 tmpCnt++;
    1643             }
    1644         }
    1645     }
     1641    x->n = nPix;
     1642    y->n = nPix;
     1643    yErr->n = nPix;
    16461644
    16471645    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    … …  
    16591657    }
    16601658
    1661     // XXX EAM : covar must be F64?
    1662     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    1663 
    1664     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    1665     fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
    1666     for (int i = 0; i < dparams->n; i++) {
    1667         if ((paramMask != NULL) && paramMask->data.U8[i])
    1668             continue;
    1669         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    1670     }
    1671 
    1672     // save the resulting chisq, nDOF, nIter
    1673     model->chisq = myMin->value;
    1674     model->nIter = myMin->iter;
    1675     model->nDOF  = y->n - nParams;
    1676 
    1677     // get the Gauss-Newton distance for fixed model parameters
    1678     if (paramMask != NULL) {
    1679         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    1680         psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
    1681         for (int i = 0; i < dparams->n; i++) {
    1682             if (!paramMask->data.U8[i])
    1683                 continue;
    1684             dparams->data.F32[i] = delta->data.F64[i];
    1685         }
    1686         psFree (delta);
    1687     }
    1688 
    1689     // set the model success or failure status
    1690     if (!fitStatus) {
    1691         model->status = PM_MODEL_NONCONVERGE;
    1692     } else {
    1693         model->status = PM_MODEL_SUCCESS;
    1694     }
    1695 
    1696     // 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     if (!onPic) {
    1702         model->status = PM_MODEL_OFFIMAGE;
    1703     }
    1704 
    1705     source->mode |= PM_SOURCE_FITTED;
    1706 
    1707     psFree(x);
    1708     psFree(y);
    1709     psFree(yErr);
    1710     psFree(myMin);
    1711     psFree(covar);
    1712     psFree(paramMask);
    1713 
    1714     rc = (onPic && fitStatus);
    1715     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1716     return(rc);
    1717 }
    1718 
    1719 // XXX EAM : new version with parameter range limits and weight enhancement
    1720 bool pmSourceFitModel (pmSource *source,
    1721                        pmModel *model,
    1722                        const bool PSF)
    1723 {
    1724     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1725     PS_ASSERT_PTR_NON_NULL(source, false);
    1726     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1727     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1728     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1729     PS_ASSERT_PTR_NON_NULL(source->mask, false);
    1730     PS_ASSERT_PTR_NON_NULL(source->weight, false);
    1731 
    1732     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1733     //           tests below could be conditions (!NULL)
    1734 
    1735     psBool fitStatus = true;
    1736     psBool onPic     = true;
    1737     psBool rc        = true;
    1738     // psF32  Ro, ymodel;
    1739 
    1740     psVector *params = model->params;
    1741     psVector *dparams = model->dparams;
    1742     psVector *paramMask = NULL;
    1743 
    1744     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1745 
    1746     // XXX EAM : I need to use the sky value to constrain the weight model
    1747     int nParams = PSF ? params->n - 4 : params->n;
    1748     // psF32 So = params->data.F32[0];
    1749 
    1750     // find the number of valid pixels
    1751     // XXX EAM : this loop and the loop below could just be one pass
    1752     //           using the psArrayAdd and psVectorExtend functions
    1753     psS32 count = 0;
    1754     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1755         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1756             if (source->mask->data.U8[i][j] == 0) {
    1757                 count++;
    1758             }
    1759         }
    1760     }
    1761     if (count <  nParams + 1) {
    1762         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    1763         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    1764         model->status = PM_MODEL_BADARGS;
    1765         return(false);
    1766     }
    1767 
    1768     // construct the coordinate and value entries
    1769     psArray *x = psArrayAlloc(count);
    1770     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1771     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1772     psS32 tmpCnt = 0;
    1773     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1774         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1775             if (source->mask->data.U8[i][j] == 0) {
    1776                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1777                 // XXX: Convert i/j to image space:
    1778                 // XXX EAM: coord order is (x,y) == (col,row)
    1779                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1780                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1781                 x->data[tmpCnt] = (psPtr *) coord;
    1782                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1783 
    1784                 // compare observed flux to model flux to adjust weight
    1785                 // ymodel = modelFunc (NULL, model->params, coord);
    1786 
    1787                 // this test enhances the weight based on deviation from the model flux
    1788                 // Ro = 1.0 + sqrt (fabs (y->data.F32[tmpCnt] - ymodel) / (PS_SQR(ymodel - So) + 1.0));
    1789 
    1790                 // psMinimizeLMChi2_EAM takes wt = 1/dY^2
    1791                 if (source->weight->data.F32[i][j] == 0) {
    1792                     yErr->data.F32[tmpCnt] = 0.0;
    1793                 } else {
    1794                     // yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
    1795                     yErr->data.F32[tmpCnt] = 1.0 / source->weight->data.F32[i][j];
    1796                 }
    1797                 tmpCnt++;
    1798             }
    1799         }
    1800     }
    1801 
    1802     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    1803                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    1804 
    1805     // PSF model only fits first 4 parameters, FLT model fits all
    1806     if (PSF) {
    1807         paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    1808         for (int i = 0; i < 4; i++) {
    1809             paramMask->data.U8[i] = 0;
    1810         }
    1811         for (int i = 4; i < paramMask->n; i++) {
    1812             paramMask->data.U8[i] = 1;
    1813         }
    1814     }
    1815 
    1816     // XXX EAM : I've added three types of parameter range checks
    1817     // XXX EAM : this requires my new psMinimization functions
     1659    // Set the parameter range checks
    18181660    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    18191661    psVector *beta_lim = NULL;
    … …  
    18391681    }
    18401682
    1841     // XXX EAM: we need to do something (give an error?) if rc is false
    1842     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1843 
    1844     // XXX EAM: save the resulting chisq, nDOF, nIter
     1683    // save the resulting chisq, nDOF, nIter
    18451684    model->chisq = myMin->value;
    18461685    model->nIter = myMin->iter;
    18471686    model->nDOF  = y->n - nParams;
    18481687
    1849     // XXX EAM get the Gauss-Newton distance for fixed model parameters
     1688    // get the Gauss-Newton distance for fixed model parameters
    18501689    if (paramMask != NULL) {
    18511690        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    … …  
    18651704    }
    18661705
    1867     // XXX models can go insane: reject these
     1706    // models can go insane: reject these
    18681707    onPic &= (params->data.F32[2] >= source->pixels->col0);
    18691708    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
Note: See TracChangeset for help on using the changeset viewer.