headers/libs/agg/agg_math.h

   1 //----------------------------------------------------------------------------
   2 // Anti-Grain Geometry - Version 2.4
   3 // Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com)
   4 //
   5 // Permission to copy, use, modify, sell and distribute this software
   6 // is granted provided this copyright notice appears in all copies.
   7 // This software is provided "as is" without express or implied
   8 // warranty, and with no claim as to its suitability for any purpose.
   9 //
  10 //----------------------------------------------------------------------------
  11 // Contact: mcseem@antigrain.com
  12 //          mcseemagg@yahoo.com
  13 //          http://www.antigrain.com
  14 //----------------------------------------------------------------------------
  15 // Bessel function (besj) was adapted for use in AGG library by Andy Wilk
  16 // Contact: castor.vulgaris@gmail.com
  17 //----------------------------------------------------------------------------
  18
  19 #ifndef AGG_MATH_INCLUDED
  20 #define AGG_MATH_INCLUDED
  21
  22 #include <math.h>
  23 #include "agg_basics.h"
  24
  25 namespace agg
  26 {
  27
  28     //------------------------------------------------------vertex_dist_epsilon
  29     // Coinciding points maximal distance (Epsilon)
  30     const double vertex_dist_epsilon = 1e-14;
  31
  32     //-----------------------------------------------------intersection_epsilon
  33     // See calc_intersection
  34     const double intersection_epsilon = 1.0e-30;
  35
  36     //------------------------------------------------------------cross_product
  37     AGG_INLINE double cross_product(double x1, double y1,
  38                                     double x2, double y2,
  39                                     double x,  double y)
  40     {
  41         return (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1);
  42     }
  43
  44     //--------------------------------------------------------point_in_triangle
  45     AGG_INLINE bool point_in_triangle(double x1, double y1,
  46                                       double x2, double y2,
  47                                       double x3, double y3,
  48                                       double x,  double y)
  49     {
  50         bool cp1 = cross_product(x1, y1, x2, y2, x, y) < 0.0;
  51         bool cp2 = cross_product(x2, y2, x3, y3, x, y) < 0.0;
  52         bool cp3 = cross_product(x3, y3, x1, y1, x, y) < 0.0;
  53         return cp1 == cp2 && cp2 == cp3 && cp3 == cp1;
  54     }
  55
  56     //-----------------------------------------------------------calc_distance
  57     AGG_INLINE double calc_distance(double x1, double y1, double x2, double y2)
  58     {
  59         double dx = x2-x1;
  60         double dy = y2-y1;
  61         return sqrt(dx * dx + dy * dy);
  62     }
  63
  64     //--------------------------------------------------------calc_sq_distance
  65     AGG_INLINE double calc_sq_distance(double x1, double y1, double x2, double y2)
  66     {
  67         double dx = x2-x1;
  68         double dy = y2-y1;
  69         return dx * dx + dy * dy;
  70     }
  71
  72     //------------------------------------------------calc_line_point_distance
  73     AGG_INLINE double calc_line_point_distance(double x1, double y1,
  74                                                double x2, double y2,
  75                                                double x,  double y)
  76     {
  77         double dx = x2-x1;
  78         double dy = y2-y1;
  79         double d = sqrt(dx * dx + dy * dy);
  80         if(d < vertex_dist_epsilon)
  81         {
  82             return calc_distance(x1, y1, x, y);
  83         }
  84         return ((x - x2) * dy - (y - y2) * dx) / d;
  85     }
  86
  87     //-------------------------------------------------------calc_line_point_u
  88     AGG_INLINE double calc_segment_point_u(double x1, double y1,
  89                                            double x2, double y2,
  90                                            double x,  double y)
  91     {
  92         double dx = x2 - x1;
  93         double dy = y2 - y1;
  94
  95         if(dx == 0 && dy == 0)
  96         {
  97                 return 0;
  98         }
  99
 100         double pdx = x - x1;
 101         double pdy = y - y1;
 102
 103         return (pdx * dx + pdy * dy) / (dx * dx + dy * dy);
 104     }
 105
 106     //---------------------------------------------calc_line_point_sq_distance
 107     AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1,
 108                                                      double x2, double y2,
 109                                                      double x,  double y,
 110                                                      double u)
 111     {
 112         if(u <= 0)
 113         {
 114                 return calc_sq_distance(x, y, x1, y1);
 115         }
 116         else
 117         if(u >= 1)
 118         {
 119                 return calc_sq_distance(x, y, x2, y2);
 120         }
 121         return calc_sq_distance(x, y, x1 + u * (x2 - x1), y1 + u * (y2 - y1));
 122     }
 123
 124     //---------------------------------------------calc_line_point_sq_distance
 125     AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1,
 126                                                      double x2, double y2,
 127                                                      double x,  double y)
 128     {
 129         return
 130             calc_segment_point_sq_distance(
 131                 x1, y1, x2, y2, x, y,
 132                 calc_segment_point_u(x1, y1, x2, y2, x, y));
 133     }
 134
 135     //-------------------------------------------------------calc_intersection
 136     AGG_INLINE bool calc_intersection(double ax, double ay, double bx, double by,
 137                                       double cx, double cy, double dx, double dy,
 138                                       double* x, double* y)
 139     {
 140         double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy);
 141         double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx);
 142         if(fabs(den) < intersection_epsilon) return false;
 143         double r = num / den;
 144         *x = ax + r * (bx-ax);
 145         *y = ay + r * (by-ay);
 146         return true;
 147     }
 148
 149     //-----------------------------------------------------intersection_exists
 150     AGG_INLINE bool intersection_exists(double x1, double y1, double x2, double y2,
 151                                         double x3, double y3, double x4, double y4)
 152     {
 153         // It's less expensive but you can't control the
 154         // boundary conditions: Less or LessEqual
 155         double dx1 = x2 - x1;
 156         double dy1 = y2 - y1;
 157         double dx2 = x4 - x3;
 158         double dy2 = y4 - y3;
 159         return ((x3 - x2) * dy1 - (y3 - y2) * dx1 < 0.0) !=
 160                ((x4 - x2) * dy1 - (y4 - y2) * dx1 < 0.0) &&
 161                ((x1 - x4) * dy2 - (y1 - y4) * dx2 < 0.0) !=
 162                ((x2 - x4) * dy2 - (y2 - y4) * dx2 < 0.0);
 163
 164         // It's is more expensive but more flexible
 165         // in terms of boundary conditions.
 166         //--------------------
 167         //double den  = (x2-x1) * (y4-y3) - (y2-y1) * (x4-x3);
 168         //if(fabs(den) < intersection_epsilon) return false;
 169         //double nom1 = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
 170         //double nom2 = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
 171         //double ua = nom1 / den;
 172         //double ub = nom2 / den;
 173         //return ua >= 0.0 && ua <= 1.0 && ub >= 0.0 && ub <= 1.0;
 174     }
 175
 176     //--------------------------------------------------------calc_orthogonal
 177     AGG_INLINE void calc_orthogonal(double thickness,
 178                                     double x1, double y1,
 179                                     double x2, double y2,
 180                                     double* x, double* y)
 181     {
 182         double dx = x2 - x1;
 183         double dy = y2 - y1;
 184         double d = sqrt(dx*dx + dy*dy);
 185         *x =  thickness * dy / d;
 186         *y = -thickness * dx / d;
 187     }
 188
 189     //--------------------------------------------------------dilate_triangle
 190     AGG_INLINE void dilate_triangle(double x1, double y1,
 191                                     double x2, double y2,
 192                                     double x3, double y3,
 193                                     double *x, double* y,
 194                                     double d)
 195     {
 196         double dx1=0.0;
 197         double dy1=0.0;
 198         double dx2=0.0;
 199         double dy2=0.0;
 200         double dx3=0.0;
 201         double dy3=0.0;
 202         double loc = cross_product(x1, y1, x2, y2, x3, y3);
 203         if(fabs(loc) > intersection_epsilon)
 204         {
 205             if(cross_product(x1, y1, x2, y2, x3, y3) > 0.0)
 206             {
 207                 d = -d;
 208             }
 209             calc_orthogonal(d, x1, y1, x2, y2, &dx1, &dy1);
 210             calc_orthogonal(d, x2, y2, x3, y3, &dx2, &dy2);
 211             calc_orthogonal(d, x3, y3, x1, y1, &dx3, &dy3);
 212         }
 213         *x++ = x1 + dx1;  *y++ = y1 + dy1;
 214         *x++ = x2 + dx1;  *y++ = y2 + dy1;
 215         *x++ = x2 + dx2;  *y++ = y2 + dy2;
 216         *x++ = x3 + dx2;  *y++ = y3 + dy2;
 217         *x++ = x3 + dx3;  *y++ = y3 + dy3;
 218         *x++ = x1 + dx3;  *y++ = y1 + dy3;
 219     }
 220
 221     //------------------------------------------------------calc_triangle_area
 222     AGG_INLINE double calc_triangle_area(double x1, double y1,
 223                                          double x2, double y2,
 224                                          double x3, double y3)
 225     {
 226         return (x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3) * 0.5;
 227     }
 228
 229     //-------------------------------------------------------calc_polygon_area
 230     template<class Storage> double calc_polygon_area(const Storage& st)
 231     {
 232         unsigned i;
 233         double sum = 0.0;
 234         double x  = st[0].x;
 235         double y  = st[0].y;
 236         double xs = x;
 237         double ys = y;
 238
 239         for(i = 1; i < st.size(); i++)
 240         {
 241             const typename Storage::value_type& v = st[i];
 242             sum += x * v.y - y * v.x;
 243             x = v.x;
 244             y = v.y;
 245         }
 246         return (sum + x * ys - y * xs) * 0.5;
 247     }
 248
 249     //------------------------------------------------------------------------
 250     // Tables for fast sqrt
 251     extern int16u g_sqrt_table[1024];
 252     extern int8   g_elder_bit_table[256];
 253
 254
 255     //---------------------------------------------------------------fast_sqrt
 256     //Fast integer Sqrt - really fast: no cycles, divisions or multiplications
 257     #if defined(_MSC_VER)
 258     #pragma warning(push)
 259     #pragma warning(disable : 4035) //Disable warning "no return value"
 260     #endif
 261     AGG_INLINE unsigned fast_sqrt(unsigned val)
 262     {
 263     #if defined(_M_IX86) && defined(_MSC_VER) && !defined(AGG_NO_ASM)
 264         //For Ix86 family processors this assembler code is used.
 265         //The key command here is bsr - determination the number of the most
 266         //significant bit of the value. For other processors
 267         //(and maybe compilers) the pure C "#else" section is used.
 268         __asm
 269         {
 270             mov ebx, val
 271             mov edx, 11
 272             bsr ecx, ebx
 273             sub ecx, 9
 274             jle less_than_9_bits
 275             shr ecx, 1
 276             adc ecx, 0
 277             sub edx, ecx
 278             shl ecx, 1
 279             shr ebx, cl
 280     less_than_9_bits:
 281             xor eax, eax
 282             mov  ax, g_sqrt_table[ebx*2]
 283             mov ecx, edx
 284             shr eax, cl
 285         }
 286     #else
 287
 288         //This code is actually pure C and portable to most
 289         //arcitectures including 64bit ones.
 290         unsigned t = val;
 291         int bit=0;
 292         unsigned shift = 11;
 293
 294         //The following piece of code is just an emulation of the
 295         //Ix86 assembler command "bsr" (see above). However on old
 296         //Intels (like Intel MMX 233MHz) this code is about twice
 297         //faster (sic!) then just one "bsr". On PIII and PIV the
 298         //bsr is optimized quite well.
 299         bit = t >> 24;
 300         if(bit)
 301         {
 302             bit = g_elder_bit_table[bit] + 24;
 303         }
 304         else
 305         {
 306             bit = (t >> 16) & 0xFF;
 307             if(bit)
 308             {
 309                 bit = g_elder_bit_table[bit] + 16;
 310             }
 311             else
 312             {
 313                 bit = (t >> 8) & 0xFF;
 314                 if(bit)
 315                 {
 316                     bit = g_elder_bit_table[bit] + 8;
 317                 }
 318                 else
 319                 {
 320                     bit = g_elder_bit_table[t];
 321                 }
 322             }
 323         }
 324
 325         //This is calculation sqrt itself.
 326         bit -= 9;
 327         if(bit > 0)
 328         {
 329             bit = (bit >> 1) + (bit & 1);
 330             shift -= bit;
 331             val >>= (bit << 1);
 332         }
 333         return g_sqrt_table[val] >> shift;
 334     #endif
 335     }
 336     #if defined(_MSC_VER)
 337     #pragma warning(pop)
 338     #endif
 339
 340
 341
 342
 343     //--------------------------------------------------------------------besj
 344     // Function BESJ calculates Bessel function of first kind of order n
 345     // Arguments:
 346     //     n - an integer (>=0), the order
 347     //     x - value at which the Bessel function is required
 348     //--------------------
 349     // C++ Mathematical Library
 350     // Convereted from equivalent FORTRAN library
 351     // Converetd by Gareth Walker for use by course 392 computational project
 352     // All functions tested and yield the same results as the corresponding
 353     // FORTRAN versions.
 354     //
 355     // If you have any problems using these functions please report them to
 356     // M.Muldoon@UMIST.ac.uk
 357     //
 358     // Documentation available on the web
 359     // http://www.ma.umist.ac.uk/mrm/Teaching/392/libs/392.html
 360     // Version 1.0   8/98
 361     // 29 October, 1999
 362     //--------------------
 363     // Adapted for use in AGG library by Andy Wilk (castor.vulgaris@gmail.com)
 364     //------------------------------------------------------------------------
 365     inline double besj(double x, int n)
 366     {
 367         if(n < 0)
 368         {
 369             return 0;
 370         }
 371         double d = 1E-6;
 372         double b = 0;
 373         if(fabs(x) <= d)
 374         {
 375             if(n != 0) return 0;
 376             return 1;
 377         }
 378         double b1 = 0; // b1 is the value from the previous iteration
 379         // Set up a starting order for recurrence
 380         int m1 = (int)fabs(x) + 6;
 381         if(fabs(x) > 5)
 382         {
 383             m1 = (int)(fabs(1.4 * x + 60 / x));
 384         }
 385         int m2 = (int)(n + 2 + fabs(x) / 4);
 386         if (m1 > m2)
 387         {
 388             m2 = m1;
 389         }
 390
 391         // Apply recurrence down from curent max order
 392         for(;;)
 393         {
 394             double c3 = 0;
 395             double c2 = 1E-30;
 396             double c4 = 0;
 397             int m8 = 1;
 398             if (m2 / 2 * 2 == m2)
 399             {
 400                 m8 = -1;
 401             }
 402             int imax = m2 - 2;
 403             for (int i = 1; i <= imax; i++)
 404             {
 405                 double c6 = 2 * (m2 - i) * c2 / x - c3;
 406                 c3 = c2;
 407                 c2 = c6;
 408                 if(m2 - i - 1 == n)
 409                 {
 410                     b = c6;
 411                 }
 412                 m8 = -1 * m8;
 413                 if (m8 > 0)
 414                 {
 415                     c4 = c4 + 2 * c6;
 416                 }
 417             }
 418             double c6 = 2 * c2 / x - c3;
 419             if(n == 0)
 420             {
 421                 b = c6;
 422             }
 423             c4 += c6;
 424             b /= c4;
 425             if(fabs(b - b1) < d)
 426             {
 427                 return b;
 428             }
 429             b1 = b;
 430             m2 += 3;
 431         }
 432     }
 433
 434 }
 435
 436
 437 #endif