vigra - vigra/bessel.hxx Source File

[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/bessel.hxx |

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 2010-2011 by Ullrich Koethe                  */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    The VIGRA Website is                                              */
    00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00008 /*    Please direct questions, bug reports, and contributions to        */
    00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00010 /*        vigra@informatik.uni-hamburg.de                               */
    00011 /*                                                                      */
    00012 /*    Permission is hereby granted, free of charge, to any person       */
    00013 /*    obtaining a copy of this software and associated documentation    */
    00014 /*    files (the "Software"), to deal in the Software without           */
    00015 /*    restriction, including without limitation the rights to use,      */
    00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00017 /*    sell copies of the Software, and to permit persons to whom the    */
    00018 /*    Software is furnished to do so, subject to the following          */
    00019 /*    conditions:                                                       */
    00020 /*                                                                      */
    00021 /*    The above copyright notice and this permission notice shall be    */
    00022 /*    included in all copies or substantial portions of the             */
    00023 /*    Software.                                                         */
    00024 /*                                                                      */
    00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00033 /*                                                                      */
    00034 /************************************************************************/
    00035
    00036 #ifndef VIGRA_BESSEL_HXX
    00037 #define VIGRA_BESSEL_HXX
    00038
    00039 #include "mathutil.hxx"
    00040
    00041 #ifdef HasBoostMath
    00042 #include <boost/math/special_functions/bessel.hpp>
    00043 #endif
    00044
    00045 namespace vigra {
    00046
    00047 /** \addtogroup MathFunctions
    00048 */
    00049 //@{
    00050
    00051 namespace detail {
    00052
    00053 template <class REAL>
    00054 int msta1(REAL x, int mp)
    00055 {
    00056     double a0,f0,f1,f;
    00057     int i,n0,n1,nn;
    00058
    00059     a0 = abs(x);
    00060     n0 = (int)(1.1*a0)+1;
    00061     f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-mp;
    00062     n1 = n0+5;
    00063     f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-mp;
    00064     for(i=0;i<20;i++)
    00065     {
    00066         nn = int(n1-(n1-n0)/(1.0-f0/f1));
    00067         f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-mp;
    00068         if(abs(nn-n1) < 1)
    00069             break;
    00070         n0 = n1;
    00071         f0 = f1;
    00072         n1 = nn;
    00073         f1 = f;
    00074     }
    00075     return nn;
    00076 }
    00077
    00078 template <class REAL>
    00079 int msta2(REAL x, int n, int mp)
    00080 {
    00081     double a0,ejn,hmp,f0,f1,f,obj;
    00082     int i,n0,n1,nn;
    00083
    00084     a0 = abs(x);
    00085     hmp = 0.5*mp;
    00086     ejn = 0.5*std::log10(6.28*n) - n*std::log10(1.36*a0/n);
    00087     if (ejn <= hmp)
    00088     {
    00089         obj = mp;
    00090         n0 = (int)(1.1*a0);
    00091         if (n0 < 1)
    00092             n0 = 1;
    00093     }
    00094     else
    00095     {
    00096         obj = hmp+ejn;
    00097         n0 = n;
    00098     }
    00099     f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-obj;
    00100     n1 = n0+5;
    00101     f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-obj;
    00102     for (i=0;i<20;i++)
    00103     {
    00104         nn = int(n1-(n1-n0)/(1.0-f0/f1));
    00105         f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-obj;
    00106         if (abs(nn-n1) < 1)
    00107             break;
    00108         n0 = n1;
    00109         f0 = f1;
    00110         n1 = nn;
    00111         f1 = f;
    00112     }
    00113     return nn+10;
    00114 }
    00115
    00116 //
    00117 //  INPUT:
    00118 //  double x    -- argument of Bessel function of 1st and 2nd kind.
    00119 //  int n       -- order
    00120 //
    00121 //  OUPUT:
    00122 //
    00123 //  int nm      -- highest order actually computed (nm <= n)
    00124 //  double jn[] -- Bessel function of 1st kind, orders from 0 to nm
    00125 //  double yn[] -- Bessel function of 2nd kind, orders from 0 to nm
    00126 //
    00127 //  Computes Bessel functions of all order up to 'n' using recurrence
    00128 //  relations. If 'nm' < 'n' only 'nm' orders are returned.
    00129 //
    00130 // code has been adapted from C.R. Bond's implementation
    00131 // see http://www.crbond.com/math.htm
    00132 //
    00133 template <class REAL>
    00134 void bessjyn(int n, REAL x,int &nm, double *jn, double *yn)
    00135 {
    00136     double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
    00137     double ec,bs,byk,p0,p1,q0,q1;
    00138     static double a[] = {
    00139         -0.7031250000000000e-1,
    00140          0.1121520996093750,
    00141         -0.5725014209747314,
    00142          6.074042001273483};
    00143     static double b[] = {
    00144          0.7324218750000000e-1,
    00145         -0.2271080017089844,
    00146          1.727727502584457,
    00147         -2.438052969955606e1};
    00148     static double a1[] = {
    00149          0.1171875,
    00150         -0.1441955566406250,
    00151          0.6765925884246826,
    00152         -6.883914268109947};
    00153     static double b1[] = {
    00154        -0.1025390625,
    00155         0.2775764465332031,
    00156        -1.993531733751297,
    00157         2.724882731126854e1};
    00158
    00159     int i,k,m;
    00160     nm = n;
    00161     if (x < 1e-15)
    00162     {
    00163         for (i=0;i<=n;i++)
    00164         {
    00165             jn[i] = 0.0;
    00166             yn[i] = -1e308;
    00167         }
    00168         jn[0] = 1.0;
    00169         return;
    00170     }
    00171     if (x <= 300.0 || n > (int)(0.9*x))
    00172     {
    00173         if (n == 0)
    00174             nm = 1;
    00175         m = msta1(x,200);
    00176         if (m < nm)
    00177             nm = m;
    00178         else
    00179             m = msta2(x,nm,15);
    00180         bs = 0.0;
    00181         su = 0.0;
    00182         sv = 0.0;
    00183         f2 = 0.0;
    00184         f1 = 1.0e-100;
    00185         for (k = m;k>=0;k--)
    00186         {
    00187             f = 2.0*(k+1.0)/x*f1 - f2;
    00188             if (k <= nm)
    00189                 jn[k] = f;
    00190             if ((k == 2*(int)(k/2)) && (k != 0))
    00191             {
    00192                 bs += 2.0*f;
    00193                 su += (-1)*((k & 2)-1)*f/(double)k;
    00194             }
    00195             else if (k > 1)
    00196             {
    00197                 sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0);
    00198             }
    00199             f2 = f1;
    00200             f1 = f;
    00201         }
    00202         s0 = bs+f;
    00203         for (k=0;k<=nm;k++)
    00204         {
    00205             jn[k] /= s0;
    00206         }
    00207         ec = std::log(0.5*x) + M_EULER_GAMMA;
    00208         by0 = M_2_PI*(ec*jn[0]-4.0*su/s0);
    00209         yn[0] = by0;
    00210         by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0);
    00211         yn[1] = by1;
    00212     }
    00213     else
    00214     {
    00215         t1 = x-M_PI_4;
    00216         p0 = 1.0;
    00217         q0 = -0.125/x;
    00218         for (k=0;k<4;k++)
    00219         {
    00220             p0 += a[k]*std::pow(x,-2*k-2);
    00221             q0 += b[k]*std::pow(x,-2*k-3);
    00222         }
    00223         cu = std::sqrt(M_2_PI/x);
    00224         bj0 = cu*(p0*std::cos(t1)-q0*std::sin(t1));
    00225         by0 = cu*(p0*std::sin(t1)+q0*std::cos(t1));
    00226         jn[0] = bj0;
    00227         yn[0] = by0;
    00228         t2 = x-0.75*M_PI;
    00229         p1 = 1.0;
    00230         q1 = 0.375/x;
    00231         for (k=0;k<4;k++)
    00232         {
    00233             p1 += a1[k]*std::pow(x,-2*k-2);
    00234             q1 += b1[k]*std::pow(x,-2*k-3);
    00235         }
    00236         bj1 = cu*(p1*std::cos(t2)-q1*std::sin(t2));
    00237         by1 = cu*(p1*std::sin(t2)+q1*std::cos(t2));
    00238         jn[1] = bj1;
    00239         yn[1] = by1;
    00240         for (k=2;k<=nm;k++)
    00241         {
    00242             bjk = 2.0*(k-1.0)*bj1/x-bj0;
    00243             jn[k] = bjk;
    00244             bj0 = bj1;
    00245             bj1 = bjk;
    00246         }
    00247     }
    00248     for (k=2;k<=nm;k++)
    00249     {
    00250         byk = 2.0*(k-1.0)*by1/x-by0;
    00251         yn[k] = byk;
    00252         by0 = by1;
    00253         by1 = byk;
    00254     }
    00255 }
    00256
    00257
    00258
    00259 } // namespace detail
    00260
    00261     /*! Bessel function of the first kind.
    00262
    00263         Computes the value of BesselJ of integer order <tt>n</tt> and argument <tt>x</tt>.
    00264         Negative <tt>x</tt> are unsupported and will result in a <tt>std::domain_error</tt>.
    00265
    00266         This function wraps a number of existing implementations and falls back to
    00267         a rather slow algorithm if none of them is available. In particular,
    00268         it uses boost::math when <tt>HasBoostMath</tt> is \#defined, or native
    00269         implementations on gcc and MSVC otherwise.
    00270
    00271         <b>\#include</b> <vigra/bessel.hxx><br>
    00272         Namespace: vigra
    00273     */
    00274 inline double besselJ(int n, double x)
    00275 {
    00276     if(x < 0.0)
    00277         throw std::domain_error("besselJ(n, x): x cannot be negative");
    00278     if(x < 1e-15)
    00279         return n == 0 ? 1.0 : 0.0;
    00280 #if defined(HasBoostMath)
    00281     return boost::math::cyl_bessel_j((double)n, x);
    00282 #elif defined(__GNUC__)
    00283     return ::jn(n, x);
    00284 #elif defined(_MSC_VER)
    00285     return _jn(n, x);
    00286 #else
    00287     int an = abs(n), nr = n, s = an+2;
    00288     ArrayVector<double> t(2*s);
    00289     detail::bessjyn(an, x, nr, &t[0], &t[s]);
    00290     if(n < 0 && odd(an))
    00291         return -t[an];
    00292     else
    00293         return  t[an];
    00294 #endif
    00295 }
    00296
    00297     /*! Bessel function of the second kind.
    00298
    00299         Computes the value of BesselY of integer order <tt>n</tt> and argument <tt>x</tt>.
    00300         Negative <tt>x</tt> are unsupported and will result in a <tt>std::domain_error</tt>.
    00301
    00302         This function wraps a number of existing implementations and falls back to
    00303         a rather slow algorithm if none of them is available. In particular,
    00304         it uses boost::math when <tt>HasBoostMath</tt> is \#defined, or native
    00305         implementations on gcc and MSVC otherwise.
    00306
    00307         <b>\#include</b> <vigra/bessel.hxx><br>
    00308         Namespace: vigra
    00309     */
    00310 inline double besselY(int n, double x)
    00311 {
    00312     if(x < 0.0)
    00313         throw std::domain_error("besselY(n, x): x cannot be negative");
    00314     if(x == 0.0 )
    00315         return -std::numeric_limits<double>::infinity();
    00316 #if defined(HasBoostMath)
    00317     return boost::math::cyl_neumann((double)n, x);
    00318 #elif defined(__GNUC__)
    00319     return ::yn(n, x);
    00320 #elif defined(_MSC_VER)
    00321     return _yn(n, x);
    00322 #else
    00323     int an = abs(n), nr = n, s = an+2;
    00324     ArrayVector<double> t(2*s);
    00325     detail::bessjyn(an, x, nr, &t[0], &t[s]);
    00326     if(n < 0.0 && odd(n))
    00327         return -t[an+s];
    00328     else
    00329         return  t[an+s];
    00330 #endif
    00331 }
    00332
    00333 //@}
    00334
    00335 } // namespace vigra
    00336
    00337 #endif // VIGRA_BESSEL_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)