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 25namespace 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 438