IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5727


Ignore:
Timestamp:
Dec 7, 2005, 10:01:05 AM (21 years ago)
Author:
magnier
Message:
  • added a 'status' element to the pmModel structure, which carries the resulting status of the fit. this gets set in pmSourceModelFit.
  • added config information MOMENTS_SX_MAX and MOMENTS_SY_MAX to carry in the size of the sigma x, sigma y search plane in pmSourcePSFClump.
  • fixed pmSourceRoughClass to use the default full array for psImageCountPixelMask
  • changed the |= to = in pmSourceRoughClass source assignments
  • added the sky offset option to p_pmSourceAddOrSubModel, inverted logic on 'flag' and renamed it 'add'.
  • modified pmSourceAddModel, pmSourceSubModel to reflect the above
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5435 r5727  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-20 23:06:24 $
     8 *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-12-07 20:01:05 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    8181    tmp->chisq = 0.0;
    8282    tmp->nIter = 0;
     83    tmp->radius = 0;
     84    tmp->status = PM_MODEL_UNTRIED;
     85
    8386    psS32 Nparams = pmModelParameterCount(type);
    8487    if (Nparams == 0) {
    … …  
    769772    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    770773        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    771             if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
     774            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
    772775                continue;
     776            }
    773777
    774778            psF32 xDiff = col + source->pixels->col0 - xPeak;
    … …  
    777781            // XXX EAM : calculate xDiff, yDiff up front;
    778782            //           radius is just a function of (xDiff, yDiff)
    779             if (!VALID_RADIUS(xDiff, yDiff, R2))
     783            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
    780784                continue;
     785            }
    781786
    782787            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
    … …  
    784789            // XXX EAM : check for valid S/N in pixel
    785790            // XXX EAM : should this limit be user-defined?
    786             if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1)
     791            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
    787792                continue;
     793            }
    788794
    789795            Sum += pDiff;
    … …  
    913919        psImage *splane = NULL;
    914920        int binX, binY;
     921        bool status;
     922
     923        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
     924        if (!status)
     925            SX_MAX = 10.0;
     926        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
     927        if (!status)
     928            SY_MAX = 10.0;
    915929
    916930        // construct a sigma-plane image
    917931        // psImageAlloc does zero the data
    918         splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
     932        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
    919933        for (int i = 0; i < splane->numRows; i++)
    920934        {
    … …  
    10651079    int Ncr      = 0;
    10661080    int Nsatstar = 0;
     1081    psRegion allArray = psRegionSet (0, 0, 0, 0);
    10671082
    10681083    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    … …  
    10961111
    10971112        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1098         //
    1099         // XXX: Must verify this region (the region argument was added to psImageCountPixelMask()
    1100         // after EAM wrote this code.
    1101         //
    1102         psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows);
    1103         int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED);
     1113        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
    11041114
    11051115        // saturated star (size consistent with PSF or larger)
    … …  
    11071117        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
    11081118        if ((Nsatpix > 1) && big) {
    1109             tmpSrc->type |= PM_SOURCE_SATSTAR;
     1119            tmpSrc->type = PM_SOURCE_SATSTAR;
    11101120            Nsatstar ++;
    11111121            continue;
    … …  
    11141124        // saturated object (not a star, eg bleed trails, hot pixels)
    11151125        if (Nsatpix > 1) {
    1116             tmpSrc->type |= PM_SOURCE_SATURATED;
     1126            tmpSrc->type = PM_SOURCE_SATURATED;
    11171127            Nsat ++;
    11181128            continue;
    … …  
    11231133        // only set candidate defects if
    11241134        if ((sigX < 0.05) || (sigY < 0.05)) {
    1125             tmpSrc->type |= PM_SOURCE_DEFECT;
     1135            tmpSrc->type = PM_SOURCE_DEFECT;
    11261136            Ncr ++;
    11271137            continue;
    … …  
    11301140        // likely unsaturated galaxy (too large to be stellar)
    11311141        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    1132             tmpSrc->type |= PM_SOURCE_GALAXY;
     1142            tmpSrc->type = PM_SOURCE_GALAXY;
    11331143            Ngal ++;
    11341144            continue;
    … …  
    11431153        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    11441154        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    1145             tmpSrc->type |= PM_SOURCE_PSFSTAR;
     1155            tmpSrc->type = PM_SOURCE_PSFSTAR;
    11461156            Npsf ++;
    11471157            continue;
    … …  
    11491159
    11501160        // random type of star
    1151         tmpSrc->type |= PM_SOURCE_OTHER;
     1161        tmpSrc->type = PM_SOURCE_OTHER;
    11521162    }
    11531163
    … …  
    16031613        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    16041614        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     1615        model->status = PM_MODEL_BADARGS;
    16051616        return(false);
    16061617    }
    … …  
    16541665    }
    16551666
    1656     // XXX EAM: we need to do something (give an error?) if rc is false
    1657     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1658 
    1659     // XXX models can go insane: reject these
    1660     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1661     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1662     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1663     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1664 
    1665     // XXX EAM: save the resulting chisq, nDOF, nIter
     1667    // save the resulting chisq, nDOF, nIter
    16661668    model->chisq = myMin->value;
    16671669    model->nIter = myMin->iter;
    16681670    model->nDOF  = y->n - nParams;
    16691671
    1670     // XXX EAM get the Gauss-Newton distance for fixed model parameters
     1672    // get the Gauss-Newton distance for fixed model parameters
    16711673    if (paramMask != NULL) {
    16721674        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    … …  
    16781680        }
    16791681        psFree (delta);
     1682    }
     1683
     1684    // set the model success or failure status
     1685    if (!fitStatus) {
     1686        model->status = PM_MODEL_NONCONVERGE;
     1687    } else {
     1688        model->status = PM_MODEL_SUCCESS;
     1689    }
     1690
     1691    // models can go insane: reject these
     1692    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1693    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1694    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1695    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1696    if (!onPic) {
     1697        model->status = PM_MODEL_OFFIMAGE;
    16801698    }
    16811699
    … …  
    17371755        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    17381756        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     1757        model->status = PM_MODEL_BADARGS;
    17391758        return(false);
    17401759    }
    … …  
    18151834    // XXX EAM: psMinimizeLMChi2 does not check convergence
    18161835
    1817     // XXX models can go insane: reject these
    1818     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1819     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1820     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1821     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1822 
    18231836    // XXX EAM: save the resulting chisq, nDOF, nIter
    18241837    model->chisq = myMin->value;
    … …  
    18371850    }
    18381851
    1839     psFree(paramMask);
     1852    // set the model success or failure status
     1853    if (!fitStatus) {
     1854        model->status = PM_MODEL_NONCONVERGE;
     1855    } else {
     1856        model->status = PM_MODEL_SUCCESS;
     1857    }
     1858
     1859    // XXX models can go insane: reject these
     1860    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1861    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1862    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1863    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1864    if (!onPic) {
     1865        model->status = PM_MODEL_OFFIMAGE;
     1866    }
     1867
    18401868    psFree(x);
    18411869    psFree(y);
     1870    psFree(yErr);
    18421871    psFree(myMin);
     1872    psFree(covar);
     1873    psFree(paramMask);
    18431874
    18441875    rc = (onPic && fitStatus);
    … …  
    18511882                             pmModel *model,
    18521883                             bool center,
    1853                              psS32 flag)
     1884                             bool sky,
     1885                             bool add
     1886                                )
    18541887{
    18551888    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    … …  
    18641897    psS32 imageCol;
    18651898    psS32 imageRow;
     1899    psF32 skyValue = params->data.F32[0];
     1900    psF32 pixelValue;
    18661901
    18671902    for (psS32 i = 0; i < image->numRows; i++) {
    … …  
    18691904            if ((mask != NULL) && mask->data.U8[i][j])
    18701905                continue;
    1871             psF32 pixelValue;
     1906
    18721907            // XXX: Should you be adding the pixels for the entire subImage,
    18731908            // or a radius of pixels around it?
    … …  
    18871922            x->data.F32[0] = (float) imageCol;
    18881923            x->data.F32[1] = (float) imageRow;
    1889             pixelValue = modelFunc (NULL, params, x);
    1890             // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
    1891 
    1892             if (flag == 1) {
    1893                 pixelValue = -pixelValue;
    1894             }
    1895 
    1896             // XXX: Must figure out how to calculate the image coordinates and
    1897             // how to use the boolean "center" flag.
    1898 
    1899             image->data.F32[i][j]+= pixelValue;
     1924
     1925            // set the appropriate pixel value for this coordinate
     1926            if (sky) {
     1927                pixelValue = modelFunc (NULL, params, x);
     1928            } else {
     1929                pixelValue = modelFunc (NULL, params, x) - skyValue;
     1930            }
     1931
     1932
     1933            // add or subtract the value
     1934            if (add
     1935               ) {
     1936                image->data.F32[i][j] += pixelValue;
     1937            }
     1938            else {
     1939                image->data.F32[i][j] -= pixelValue;
     1940            }
    19001941        }
    19011942    }
    … …  
    19121953                      psImage *mask,
    19131954                      pmModel *model,
    1914                       bool center)
    1915 {
    1916     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1917     psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
     1955                      bool center,
     1956                      bool sky)
     1957{
     1958    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1959    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
    19181960    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    19191961    return(rc);
    … …  
    19251967                      psImage *mask,
    19261968                      pmModel *model,
    1927                       bool center)
    1928 {
    1929     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1930     psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
     1969                      bool center,
     1970                      bool sky)
     1971{
     1972    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1973    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
    19311974    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    19321975    return(rc);
Note: See TracChangeset for help on using the changeset viewer.