IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1823


Ignore:
Timestamp:
Sep 16, 2004, 4:14:52 PM (22 years ago)
Author:
gusciora
Message:

Added vector bin disect table lookup type methods.

Location:
trunk/psLib/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psFunctions.c

    r1784 r1823  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-10 23:20:29 $
     9 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-17 02:14:52 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424#include "psMemory.h"
    2525#include "psVector.h"
     26#include "psScalar.h"
    2627#include "psTrace.h"
    2728#include "psError.h"
    … …  
    16821683
    16831684/*****************************************************************************
    1684 VectorBinDisect(): This is a private function which takes as input a vector
    1685 of floating point data as well as a single floating point values.  The input
    1686 vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i).
    1687 This routine does a binary disection of the vector and returns "i" such
    1688 that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
     1685VectorBinDisectF32(): This is a private function which takes as input a
     1686vector of floating point data as well as a single floating point values.
     1687The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
     1688all j>=i).  This routine does a binary disection of the vector and returns
     1689"i" such that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
    16891690then this routine prints a warning message and returns -1.
     1691 
     1692XXX: Macro this for a few different types.
    16901693 *****************************************************************************/
    1691 int VectorBinDisect(float *bins,
    1692                     int numBins,
    1693                     float x)
     1694int VectorBinDisectF32(float *bins,
     1695                       int numBins,
     1696                       float x)
    16941697{
    16951698    int min;
    … …  
    16971700    int mid;
    16981701
    1699     psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1700             "---- Calling VectorBinDisect(%f)\n", x);
     1702    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1703            "---- Calling VectorBinDisectF32(%f)\n", x);
    17011704
    17021705    if ((x < bins[0]) ||
    17031706            (x > bins[numBins-1])) {
    17041707        psLogMsg(__func__, PS_LOG_WARN,
    1705                  "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).",
     1708                 "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
    17061709                 x, bins[0], bins[numBins-1]);
    17071710        return(-1);
    … …  
    17131716
    17141717    while (min != max) {
    1715         psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1718        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
    17161719                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
    17171720                min, mid, max, x, bins[mid]);
    17181721
    17191722        if (x == bins[mid]) {
    1720             psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1721                     "---- Exiting VectorBinDisect(): bin %d\n", min);
     1723            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1724                    "---- Exiting VectorBinDisectF32(): bin %d\n", mid);
     1725            return(mid);
     1726        } else if (x < bins[mid]) {
     1727            max = mid-1;
     1728        } else {
     1729            min = mid;
     1730        }
     1731        mid = ((max+1)+min)/2;
     1732    }
     1733
     1734    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1735            "---- Exiting VectorBinDisectF32(): bin %d\n", min);
     1736    return(min);
     1737}
     1738
     1739/*****************************************************************************
     1740VectorBinDisectS32(): integer version of above.
     1741 *****************************************************************************/
     1742int VectorBinDisectS32(int *bins,
     1743                       int numBins,
     1744                       int x)
     1745{
     1746    int min;
     1747    int max;
     1748    int mid;
     1749
     1750    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1751            "---- Calling VectorBinDisectS32(%f)\n", x);
     1752
     1753    if ((x < bins[0]) ||
     1754            (x > bins[numBins-1])) {
     1755        psLogMsg(__func__, PS_LOG_WARN,
     1756                 "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
     1757                 x, bins[0], bins[numBins-1]);
     1758        return(-1);
     1759    }
     1760
     1761    min = 0;
     1762    max = numBins-2;
     1763    mid = ((max+1)-min)/2;
     1764
     1765    while (min != max) {
     1766        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1767                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
     1768                min, mid, max, x, bins[mid]);
     1769
     1770        if (x == bins[mid]) {
     1771            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1772                    "---- Exiting VectorBinDisectS32(): bin %d\n", min);
    17221773            return(min);
    17231774        } else if (x < bins[mid]) {
    … …  
    17291780    }
    17301781
    1731     psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1732             "---- Exiting VectorBinDisect(): bin %d\n", min);
     1782    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1783            "---- Exiting VectorBinDisectS32(): bin %d\n", min);
    17331784    return(min);
    17341785}
     1786
     1787/*****************************************************************************
     1788p_psVectorBinDisect(): A wrapper to the above VectorBinDisect().
     1789 *****************************************************************************/
     1790int p_psVectorBinDisect(psVector *bins,
     1791                        psScalar *x)
     1792{
     1793    if (x->type.type != bins->type.type) {
     1794        // XXX: Generate error message.
     1795        return(-1);
     1796    }
     1797
     1798    if (x->type.type == PS_TYPE_S32) {
     1799        return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     1800    } else if (x->type.type == PS_TYPE_F32) {
     1801        return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     1802    } else {
     1803        // XXX: Generate error message.
     1804    }
     1805    return(-1);
     1806}
     1807
     1808/*****************************************************************************
     1809p_psInterpolate1D(): This routine will take as input n-element floating
     1810point arrays domain and range, and the x value, assumed to lie with the
     1811domain vector.  It produces as output the n-order LaGrange interpolated
     1812value of x.
     1813 
     1814XXX: do we error check for non-distinct domain values?
     1815 *****************************************************************************/
     1816float p_psFullInterpolate1DF32(float *domain,
     1817                               float *range,
     1818                               int n,
     1819                               float x)
     1820{
     1821    int i;
     1822    int m;
     1823    static psVector *p = NULL;
     1824    psVectorRecycle(p, n, PS_TYPE_F32);
     1825    p_psMemSetPersistent(p, true);
     1826
     1827    // From NR, during each iteration of the m loop, we are computing the
     1828    // p_{i ... i+m} terms.
     1829    for (m=0;m<n;m++) {
     1830        for (i=0;i<n-m;i++) {
     1831            if (m == 0) {
     1832                p->data.F32[i] = range[i];
     1833            } else {
     1834                // From NR: we are computing P_{i ... i+m}
     1835                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
     1836                                  ((range[i]-x) * p->data.F32[i+1])) /
     1837                                 (domain[i] - domain[i+m]);
     1838            }
     1839        }
     1840    }
     1841    return(p->data.F32[0]);
     1842}
     1843
     1844
     1845// This is a
     1846float p_psInterpolate1DF32(float *domain,
     1847                           float *range,
     1848                           int n,
     1849                           int order,
     1850                           float x)
     1851{
     1852    int binNum;
     1853    int numIntPoints = order+1;
     1854    int origin;
     1855
     1856    binNum = VectorBinDisectF32(domain, n, x);
     1857    if (0 == numIntPoints%2) {
     1858        origin = binNum - ((numIntPoints/2) - 1);
     1859    } else {
     1860        origin = binNum - (numIntPoints/2);
     1861        if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
     1862            // x is closer to binNum+1.
     1863            origin = 1 + (binNum - (numIntPoints/2));
     1864        }
     1865    }
     1866    if (origin < 0) {
     1867        origin = 0;
     1868    }
     1869    if ((origin + numIntPoints) > n) {
     1870        origin = n - numIntPoints;
     1871    }
     1872
     1873    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
     1874}
     1875
     1876/*****************************************************************************
     1877p_psInterpolate1D(): This routine will take as input psVectors domain and
     1878range, and the x value, assumed to lie with the domain vector.  It produces
     1879as output the LaGrange interpolated value of a polynomial of the specified
     1880order around the point x.
     1881 
     1882XXX: This stuff does not work with a mask.
     1883 *****************************************************************************/
     1884float p_psInterpolate1D(psVector *domain,
     1885                        psVector *range,
     1886                        int order,
     1887                        psScalar *x)
     1888{
     1889    if (domain->type.type != range->type.type != x->type.type) {
     1890        // XXX psError
     1891    }
     1892    if (domain->n != range->n) {
     1893        // XXX psError
     1894    }
     1895    if (order > (domain->n - 1)) {
     1896        // XXX psError: not enough data points for order-order interpolation.
     1897    }
     1898
     1899    if (x->type.type == PS_TYPE_F32) {
     1900        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
     1901                                    domain->n, order, x->data.F32));
     1902    } else {
     1903        // XXX psError: type not supported
     1904    }
     1905
     1906    return(-1.0);
     1907}
     1908
    17351909
    17361910float psSpline1DEval(const psSpline1D *spline,
    … …  
    17411915
    17421916    n = spline->n;
    1743     binNum = VectorBinDisect(spline->domains, (spline->n)+1, x);
     1917    binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x);
    17441918    if (binNum == -1) {
    17451919        psLogMsg(__func__, PS_LOG_WARN,

  • trunk/psLib/src/dataManip/psFunctions.h

    r1775 r1823  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-10 02:52:02 $
     14*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-17 02:14:52 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424
    2525#include "psVector.h"
     26#include "psScalar.h"
    2627
    2728/** \addtogroup Stats
    … …  
    405406                          psVector *data);
    406407
     408int p_psVectorBinDisect(psVector *bins,
     409                        psScalar *x);
    407410
    408411/* \} */// End of MathGroup Functions
    409412
    410413#endif
     414

  • trunk/psLib/src/dataManip/psStats.c

    r1722 r1823  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-08 06:04:30 $
     11 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-17 02:14:52 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    14491449void p_ps_FitTheGaussian()
    14501450{
    1451  
    14521451    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    14531452    // NOTE: This code uses the psMinimize.c functions to perform
    … …  
    14551454    // I am commenting this code out, and replacing it with code which
    14561455    // calculates the mean directly on the robustHistogram.
     1456    //
     1457    // XXX: Replace this with the LaGrange interpolation stuff.
    14571458 
    14581459    domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
    … …  
    16711672psHistogram* psVectorHistogram(psHistogram* out,
    16721673                               const psVector* restrict in,
    1673                                const psVector* restrict mask, unsigned int maskVal)
     1674                               const psVector* restrict mask,
     1675                               unsigned int maskVal)
    16741676{
    16751677    int i = 0;                  // Loop index variable
    1676     int j = 0;                  // Loop index variable
    16771678    float binSize = 0.0;        // Histogram bin size
    16781679    int binNum = 0;             // A temporary bin number
    16791680    int numBins = 0;            // The total number of bins
     1681    int tmp = 0;
     1682    psScalar tmpScalar;
     1683    tmpScalar.type.type = PS_TYPE_F32;
    16801684
    16811685    // NOTE: Verify that this is the correct action.
    … …  
    17081712        }
    17091713    }
     1714
    17101715    // NOTE: determine the correct action for a variety of other cases:
    17111716    // in vector has 0 elements, and histogram structure has zero bins.
    … …  
    17421747                    // bin number requires a bit more work.
    17431748                } else {
    1744                     // NOTE: This is slow.  Put a smarter algorithm here to
    1745                     // find the correct bin number (bin search, probably)
    1746                     for (j = 0; j < (out->bounds->n) - 1; j++) {
    1747                         if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
    1748                                 (in->data.F32[i] <= out->bounds->data.S32[j + 1])) {
    1749                             (out->nums->data.S32[j])++;
    1750                         }
     1749                    tmpScalar.data.F32 = in->data.F32[i];
     1750                    tmp = p_psVectorBinDisect(out->bounds,
     1751                                              &tmpScalar);
     1752                    if (tmp == -1) {
     1753                        //                      XXX: Generate warning message
     1754                    } else {
     1755                        (out->nums->data.S32[tmp])++;
    17511756                    }
    17521757                }

  • trunk/psLib/src/math/psPolynomial.c

    r1784 r1823  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-10 23:20:29 $
     9 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-17 02:14:52 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424#include "psMemory.h"
    2525#include "psVector.h"
     26#include "psScalar.h"
    2627#include "psTrace.h"
    2728#include "psError.h"
    … …  
    16821683
    16831684/*****************************************************************************
    1684 VectorBinDisect(): This is a private function which takes as input a vector
    1685 of floating point data as well as a single floating point values.  The input
    1686 vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i).
    1687 This routine does a binary disection of the vector and returns "i" such
    1688 that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
     1685VectorBinDisectF32(): This is a private function which takes as input a
     1686vector of floating point data as well as a single floating point values.
     1687The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
     1688all j>=i).  This routine does a binary disection of the vector and returns
     1689"i" such that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
    16891690then this routine prints a warning message and returns -1.
     1691 
     1692XXX: Macro this for a few different types.
    16901693 *****************************************************************************/
    1691 int VectorBinDisect(float *bins,
    1692                     int numBins,
    1693                     float x)
     1694int VectorBinDisectF32(float *bins,
     1695                       int numBins,
     1696                       float x)
    16941697{
    16951698    int min;
    … …  
    16971700    int mid;
    16981701
    1699     psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1700             "---- Calling VectorBinDisect(%f)\n", x);
     1702    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1703            "---- Calling VectorBinDisectF32(%f)\n", x);
    17011704
    17021705    if ((x < bins[0]) ||
    17031706            (x > bins[numBins-1])) {
    17041707        psLogMsg(__func__, PS_LOG_WARN,
    1705                  "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).",
     1708                 "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
    17061709                 x, bins[0], bins[numBins-1]);
    17071710        return(-1);
    … …  
    17131716
    17141717    while (min != max) {
    1715         psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1718        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
    17161719                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
    17171720                min, mid, max, x, bins[mid]);
    17181721
    17191722        if (x == bins[mid]) {
    1720             psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1721                     "---- Exiting VectorBinDisect(): bin %d\n", min);
     1723            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1724                    "---- Exiting VectorBinDisectF32(): bin %d\n", mid);
     1725            return(mid);
     1726        } else if (x < bins[mid]) {
     1727            max = mid-1;
     1728        } else {
     1729            min = mid;
     1730        }
     1731        mid = ((max+1)+min)/2;
     1732    }
     1733
     1734    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1735            "---- Exiting VectorBinDisectF32(): bin %d\n", min);
     1736    return(min);
     1737}
     1738
     1739/*****************************************************************************
     1740VectorBinDisectS32(): integer version of above.
     1741 *****************************************************************************/
     1742int VectorBinDisectS32(int *bins,
     1743                       int numBins,
     1744                       int x)
     1745{
     1746    int min;
     1747    int max;
     1748    int mid;
     1749
     1750    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1751            "---- Calling VectorBinDisectS32(%f)\n", x);
     1752
     1753    if ((x < bins[0]) ||
     1754            (x > bins[numBins-1])) {
     1755        psLogMsg(__func__, PS_LOG_WARN,
     1756                 "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
     1757                 x, bins[0], bins[numBins-1]);
     1758        return(-1);
     1759    }
     1760
     1761    min = 0;
     1762    max = numBins-2;
     1763    mid = ((max+1)-min)/2;
     1764
     1765    while (min != max) {
     1766        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1767                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
     1768                min, mid, max, x, bins[mid]);
     1769
     1770        if (x == bins[mid]) {
     1771            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1772                    "---- Exiting VectorBinDisectS32(): bin %d\n", min);
    17221773            return(min);
    17231774        } else if (x < bins[mid]) {
    … …  
    17291780    }
    17301781
    1731     psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1732             "---- Exiting VectorBinDisect(): bin %d\n", min);
     1782    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1783            "---- Exiting VectorBinDisectS32(): bin %d\n", min);
    17331784    return(min);
    17341785}
     1786
     1787/*****************************************************************************
     1788p_psVectorBinDisect(): A wrapper to the above VectorBinDisect().
     1789 *****************************************************************************/
     1790int p_psVectorBinDisect(psVector *bins,
     1791                        psScalar *x)
     1792{
     1793    if (x->type.type != bins->type.type) {
     1794        // XXX: Generate error message.
     1795        return(-1);
     1796    }
     1797
     1798    if (x->type.type == PS_TYPE_S32) {
     1799        return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     1800    } else if (x->type.type == PS_TYPE_F32) {
     1801        return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     1802    } else {
     1803        // XXX: Generate error message.
     1804    }
     1805    return(-1);
     1806}
     1807
     1808/*****************************************************************************
     1809p_psInterpolate1D(): This routine will take as input n-element floating
     1810point arrays domain and range, and the x value, assumed to lie with the
     1811domain vector.  It produces as output the n-order LaGrange interpolated
     1812value of x.
     1813 
     1814XXX: do we error check for non-distinct domain values?
     1815 *****************************************************************************/
     1816float p_psFullInterpolate1DF32(float *domain,
     1817                               float *range,
     1818                               int n,
     1819                               float x)
     1820{
     1821    int i;
     1822    int m;
     1823    static psVector *p = NULL;
     1824    psVectorRecycle(p, n, PS_TYPE_F32);
     1825    p_psMemSetPersistent(p, true);
     1826
     1827    // From NR, during each iteration of the m loop, we are computing the
     1828    // p_{i ... i+m} terms.
     1829    for (m=0;m<n;m++) {
     1830        for (i=0;i<n-m;i++) {
     1831            if (m == 0) {
     1832                p->data.F32[i] = range[i];
     1833            } else {
     1834                // From NR: we are computing P_{i ... i+m}
     1835                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
     1836                                  ((range[i]-x) * p->data.F32[i+1])) /
     1837                                 (domain[i] - domain[i+m]);
     1838            }
     1839        }
     1840    }
     1841    return(p->data.F32[0]);
     1842}
     1843
     1844
     1845// This is a
     1846float p_psInterpolate1DF32(float *domain,
     1847                           float *range,
     1848                           int n,
     1849                           int order,
     1850                           float x)
     1851{
     1852    int binNum;
     1853    int numIntPoints = order+1;
     1854    int origin;
     1855
     1856    binNum = VectorBinDisectF32(domain, n, x);
     1857    if (0 == numIntPoints%2) {
     1858        origin = binNum - ((numIntPoints/2) - 1);
     1859    } else {
     1860        origin = binNum - (numIntPoints/2);
     1861        if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
     1862            // x is closer to binNum+1.
     1863            origin = 1 + (binNum - (numIntPoints/2));
     1864        }
     1865    }
     1866    if (origin < 0) {
     1867        origin = 0;
     1868    }
     1869    if ((origin + numIntPoints) > n) {
     1870        origin = n - numIntPoints;
     1871    }
     1872
     1873    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
     1874}
     1875
     1876/*****************************************************************************
     1877p_psInterpolate1D(): This routine will take as input psVectors domain and
     1878range, and the x value, assumed to lie with the domain vector.  It produces
     1879as output the LaGrange interpolated value of a polynomial of the specified
     1880order around the point x.
     1881 
     1882XXX: This stuff does not work with a mask.
     1883 *****************************************************************************/
     1884float p_psInterpolate1D(psVector *domain,
     1885                        psVector *range,
     1886                        int order,
     1887                        psScalar *x)
     1888{
     1889    if (domain->type.type != range->type.type != x->type.type) {
     1890        // XXX psError
     1891    }
     1892    if (domain->n != range->n) {
     1893        // XXX psError
     1894    }
     1895    if (order > (domain->n - 1)) {
     1896        // XXX psError: not enough data points for order-order interpolation.
     1897    }
     1898
     1899    if (x->type.type == PS_TYPE_F32) {
     1900        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
     1901                                    domain->n, order, x->data.F32));
     1902    } else {
     1903        // XXX psError: type not supported
     1904    }
     1905
     1906    return(-1.0);
     1907}
     1908
    17351909
    17361910float psSpline1DEval(const psSpline1D *spline,
    … …  
    17411915
    17421916    n = spline->n;
    1743     binNum = VectorBinDisect(spline->domains, (spline->n)+1, x);
     1917    binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x);
    17441918    if (binNum == -1) {
    17451919        psLogMsg(__func__, PS_LOG_WARN,

  • trunk/psLib/src/math/psPolynomial.h

    r1775 r1823  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-10 02:52:02 $
     14*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-17 02:14:52 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424
    2525#include "psVector.h"
     26#include "psScalar.h"
    2627
    2728/** \addtogroup Stats
    … …  
    405406                          psVector *data);
    406407
     408int p_psVectorBinDisect(psVector *bins,
     409                        psScalar *x);
    407410
    408411/* \} */// End of MathGroup Functions
    409412
    410413#endif
     414

  • trunk/psLib/src/math/psSpline.c

    r1784 r1823  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-10 23:20:29 $
     9 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-17 02:14:52 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424#include "psMemory.h"
    2525#include "psVector.h"
     26#include "psScalar.h"
    2627#include "psTrace.h"
    2728#include "psError.h"
    … …  
    16821683
    16831684/*****************************************************************************
    1684 VectorBinDisect(): This is a private function which takes as input a vector
    1685 of floating point data as well as a single floating point values.  The input
    1686 vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i).
    1687 This routine does a binary disection of the vector and returns "i" such
    1688 that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
     1685VectorBinDisectF32(): This is a private function which takes as input a
     1686vector of floating point data as well as a single floating point values.
     1687The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
     1688all j>=i).  This routine does a binary disection of the vector and returns
     1689"i" such that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
    16891690then this routine prints a warning message and returns -1.
     1691 
     1692XXX: Macro this for a few different types.
    16901693 *****************************************************************************/
    1691 int VectorBinDisect(float *bins,
    1692                     int numBins,
    1693                     float x)
     1694int VectorBinDisectF32(float *bins,
     1695                       int numBins,
     1696                       float x)
    16941697{
    16951698    int min;
    … …  
    16971700    int mid;
    16981701
    1699     psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1700             "---- Calling VectorBinDisect(%f)\n", x);
     1702    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1703            "---- Calling VectorBinDisectF32(%f)\n", x);
    17011704
    17021705    if ((x < bins[0]) ||
    17031706            (x > bins[numBins-1])) {
    17041707        psLogMsg(__func__, PS_LOG_WARN,
    1705                  "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).",
     1708                 "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
    17061709                 x, bins[0], bins[numBins-1]);
    17071710        return(-1);
    … …  
    17131716
    17141717    while (min != max) {
    1715         psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1718        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
    17161719                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
    17171720                min, mid, max, x, bins[mid]);
    17181721
    17191722        if (x == bins[mid]) {
    1720             psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1721                     "---- Exiting VectorBinDisect(): bin %d\n", min);
     1723            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1724                    "---- Exiting VectorBinDisectF32(): bin %d\n", mid);
     1725            return(mid);
     1726        } else if (x < bins[mid]) {
     1727            max = mid-1;
     1728        } else {
     1729            min = mid;
     1730        }
     1731        mid = ((max+1)+min)/2;
     1732    }
     1733
     1734    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
     1735            "---- Exiting VectorBinDisectF32(): bin %d\n", min);
     1736    return(min);
     1737}
     1738
     1739/*****************************************************************************
     1740VectorBinDisectS32(): integer version of above.
     1741 *****************************************************************************/
     1742int VectorBinDisectS32(int *bins,
     1743                       int numBins,
     1744                       int x)
     1745{
     1746    int min;
     1747    int max;
     1748    int mid;
     1749
     1750    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1751            "---- Calling VectorBinDisectS32(%f)\n", x);
     1752
     1753    if ((x < bins[0]) ||
     1754            (x > bins[numBins-1])) {
     1755        psLogMsg(__func__, PS_LOG_WARN,
     1756                 "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
     1757                 x, bins[0], bins[numBins-1]);
     1758        return(-1);
     1759    }
     1760
     1761    min = 0;
     1762    max = numBins-2;
     1763    mid = ((max+1)-min)/2;
     1764
     1765    while (min != max) {
     1766        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1767                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
     1768                min, mid, max, x, bins[mid]);
     1769
     1770        if (x == bins[mid]) {
     1771            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1772                    "---- Exiting VectorBinDisectS32(): bin %d\n", min);
    17221773            return(min);
    17231774        } else if (x < bins[mid]) {
    … …  
    17291780    }
    17301781
    1731     psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
    1732             "---- Exiting VectorBinDisect(): bin %d\n", min);
     1782    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
     1783            "---- Exiting VectorBinDisectS32(): bin %d\n", min);
    17331784    return(min);
    17341785}
     1786
     1787/*****************************************************************************
     1788p_psVectorBinDisect(): A wrapper to the above VectorBinDisect().
     1789 *****************************************************************************/
     1790int p_psVectorBinDisect(psVector *bins,
     1791                        psScalar *x)
     1792{
     1793    if (x->type.type != bins->type.type) {
     1794        // XXX: Generate error message.
     1795        return(-1);
     1796    }
     1797
     1798    if (x->type.type == PS_TYPE_S32) {
     1799        return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     1800    } else if (x->type.type == PS_TYPE_F32) {
     1801        return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     1802    } else {
     1803        // XXX: Generate error message.
     1804    }
     1805    return(-1);
     1806}
     1807
     1808/*****************************************************************************
     1809p_psInterpolate1D(): This routine will take as input n-element floating
     1810point arrays domain and range, and the x value, assumed to lie with the
     1811domain vector.  It produces as output the n-order LaGrange interpolated
     1812value of x.
     1813 
     1814XXX: do we error check for non-distinct domain values?
     1815 *****************************************************************************/
     1816float p_psFullInterpolate1DF32(float *domain,
     1817                               float *range,
     1818                               int n,
     1819                               float x)
     1820{
     1821    int i;
     1822    int m;
     1823    static psVector *p = NULL;
     1824    psVectorRecycle(p, n, PS_TYPE_F32);
     1825    p_psMemSetPersistent(p, true);
     1826
     1827    // From NR, during each iteration of the m loop, we are computing the
     1828    // p_{i ... i+m} terms.
     1829    for (m=0;m<n;m++) {
     1830        for (i=0;i<n-m;i++) {
     1831            if (m == 0) {
     1832                p->data.F32[i] = range[i];
     1833            } else {
     1834                // From NR: we are computing P_{i ... i+m}
     1835                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
     1836                                  ((range[i]-x) * p->data.F32[i+1])) /
     1837                                 (domain[i] - domain[i+m]);
     1838            }
     1839        }
     1840    }
     1841    return(p->data.F32[0]);
     1842}
     1843
     1844
     1845// This is a
     1846float p_psInterpolate1DF32(float *domain,
     1847                           float *range,
     1848                           int n,
     1849                           int order,
     1850                           float x)
     1851{
     1852    int binNum;
     1853    int numIntPoints = order+1;
     1854    int origin;
     1855
     1856    binNum = VectorBinDisectF32(domain, n, x);
     1857    if (0 == numIntPoints%2) {
     1858        origin = binNum - ((numIntPoints/2) - 1);
     1859    } else {
     1860        origin = binNum - (numIntPoints/2);
     1861        if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
     1862            // x is closer to binNum+1.
     1863            origin = 1 + (binNum - (numIntPoints/2));
     1864        }
     1865    }
     1866    if (origin < 0) {
     1867        origin = 0;
     1868    }
     1869    if ((origin + numIntPoints) > n) {
     1870        origin = n - numIntPoints;
     1871    }
     1872
     1873    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
     1874}
     1875
     1876/*****************************************************************************
     1877p_psInterpolate1D(): This routine will take as input psVectors domain and
     1878range, and the x value, assumed to lie with the domain vector.  It produces
     1879as output the LaGrange interpolated value of a polynomial of the specified
     1880order around the point x.
     1881 
     1882XXX: This stuff does not work with a mask.
     1883 *****************************************************************************/
     1884float p_psInterpolate1D(psVector *domain,
     1885                        psVector *range,
     1886                        int order,
     1887                        psScalar *x)
     1888{
     1889    if (domain->type.type != range->type.type != x->type.type) {
     1890        // XXX psError
     1891    }
     1892    if (domain->n != range->n) {
     1893        // XXX psError
     1894    }
     1895    if (order > (domain->n - 1)) {
     1896        // XXX psError: not enough data points for order-order interpolation.
     1897    }
     1898
     1899    if (x->type.type == PS_TYPE_F32) {
     1900        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
     1901                                    domain->n, order, x->data.F32));
     1902    } else {
     1903        // XXX psError: type not supported
     1904    }
     1905
     1906    return(-1.0);
     1907}
     1908
    17351909
    17361910float psSpline1DEval(const psSpline1D *spline,
    … …  
    17411915
    17421916    n = spline->n;
    1743     binNum = VectorBinDisect(spline->domains, (spline->n)+1, x);
     1917    binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x);
    17441918    if (binNum == -1) {
    17451919        psLogMsg(__func__, PS_LOG_WARN,

  • trunk/psLib/src/math/psSpline.h

    r1775 r1823  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-10 02:52:02 $
     14*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-17 02:14:52 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424
    2525#include "psVector.h"
     26#include "psScalar.h"
    2627
    2728/** \addtogroup Stats
    … …  
    405406                          psVector *data);
    406407
     408int p_psVectorBinDisect(psVector *bins,
     409                        psScalar *x);
    407410
    408411/* \} */// End of MathGroup Functions
    409412
    410413#endif
     414

  • trunk/psLib/src/math/psStats.c

    r1722 r1823  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-08 06:04:30 $
     11 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-17 02:14:52 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    14491449void p_ps_FitTheGaussian()
    14501450{
    1451  
    14521451    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    14531452    // NOTE: This code uses the psMinimize.c functions to perform
    … …  
    14551454    // I am commenting this code out, and replacing it with code which
    14561455    // calculates the mean directly on the robustHistogram.
     1456    //
     1457    // XXX: Replace this with the LaGrange interpolation stuff.
    14571458 
    14581459    domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
    … …  
    16711672psHistogram* psVectorHistogram(psHistogram* out,
    16721673                               const psVector* restrict in,
    1673                                const psVector* restrict mask, unsigned int maskVal)
     1674                               const psVector* restrict mask,
     1675                               unsigned int maskVal)
    16741676{
    16751677    int i = 0;                  // Loop index variable
    1676     int j = 0;                  // Loop index variable
    16771678    float binSize = 0.0;        // Histogram bin size
    16781679    int binNum = 0;             // A temporary bin number
    16791680    int numBins = 0;            // The total number of bins
     1681    int tmp = 0;
     1682    psScalar tmpScalar;
     1683    tmpScalar.type.type = PS_TYPE_F32;
    16801684
    16811685    // NOTE: Verify that this is the correct action.
    … …  
    17081712        }
    17091713    }
     1714
    17101715    // NOTE: determine the correct action for a variety of other cases:
    17111716    // in vector has 0 elements, and histogram structure has zero bins.
    … …  
    17421747                    // bin number requires a bit more work.
    17431748                } else {
    1744                     // NOTE: This is slow.  Put a smarter algorithm here to
    1745                     // find the correct bin number (bin search, probably)
    1746                     for (j = 0; j < (out->bounds->n) - 1; j++) {
    1747                         if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
    1748                                 (in->data.F32[i] <= out->bounds->data.S32[j + 1])) {
    1749                             (out->nums->data.S32[j])++;
    1750                         }
     1749                    tmpScalar.data.F32 = in->data.F32[i];
     1750                    tmp = p_psVectorBinDisect(out->bounds,
     1751                                              &tmpScalar);
     1752                    if (tmp == -1) {
     1753                        //                      XXX: Generate warning message
     1754                    } else {
     1755                        (out->nums->data.S32[tmp])++;
    17511756                    }
    17521757                }
Note: See TracChangeset for help on using the changeset viewer.