Anti-Grain Geometry - AGG (libagg)
2.5
|
00001 //---------------------------------------------------------------------------- 00002 // Anti-Grain Geometry (AGG) - Version 2.5 00003 // A high quality rendering engine for C++ 00004 // Copyright (C) 2002-2006 Maxim Shemanarev 00005 // Contact: mcseem@antigrain.com 00006 // mcseemagg@yahoo.com 00007 // http://antigrain.com 00008 // 00009 // AGG is free software; you can redistribute it and/or 00010 // modify it under the terms of the GNU General Public License 00011 // as published by the Free Software Foundation; either version 2 00012 // of the License, or (at your option) any later version. 00013 // 00014 // AGG is distributed in the hope that it will be useful, 00015 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00017 // GNU General Public License for more details. 00018 // 00019 // You should have received a copy of the GNU General Public License 00020 // along with AGG; if not, write to the Free Software 00021 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 00022 // MA 02110-1301, USA. 00023 //---------------------------------------------------------------------------- 00024 // Bessel function (besj) was adapted for use in AGG library by Andy Wilk 00025 // Contact: castor.vulgaris@gmail.com 00026 //---------------------------------------------------------------------------- 00027 00028 #ifndef AGG_MATH_INCLUDED 00029 #define AGG_MATH_INCLUDED 00030 00031 #include <math.h> 00032 #include "agg_basics.h" 00033 00034 namespace agg 00035 { 00036 00037 //------------------------------------------------------vertex_dist_epsilon 00038 // Coinciding points maximal distance (Epsilon) 00039 const double vertex_dist_epsilon = 1e-14; 00040 00041 //-----------------------------------------------------intersection_epsilon 00042 // See calc_intersection 00043 const double intersection_epsilon = 1.0e-30; 00044 00045 //------------------------------------------------------------cross_product 00046 AGG_INLINE double cross_product(double x1, double y1, 00047 double x2, double y2, 00048 double x, double y) 00049 { 00050 return (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1); 00051 } 00052 00053 //--------------------------------------------------------point_in_triangle 00054 AGG_INLINE bool point_in_triangle(double x1, double y1, 00055 double x2, double y2, 00056 double x3, double y3, 00057 double x, double y) 00058 { 00059 bool cp1 = cross_product(x1, y1, x2, y2, x, y) < 0.0; 00060 bool cp2 = cross_product(x2, y2, x3, y3, x, y) < 0.0; 00061 bool cp3 = cross_product(x3, y3, x1, y1, x, y) < 0.0; 00062 return cp1 == cp2 && cp2 == cp3 && cp3 == cp1; 00063 } 00064 00065 //-----------------------------------------------------------calc_distance 00066 AGG_INLINE double calc_distance(double x1, double y1, double x2, double y2) 00067 { 00068 double dx = x2-x1; 00069 double dy = y2-y1; 00070 return sqrt(dx * dx + dy * dy); 00071 } 00072 00073 //--------------------------------------------------------calc_sq_distance 00074 AGG_INLINE double calc_sq_distance(double x1, double y1, double x2, double y2) 00075 { 00076 double dx = x2-x1; 00077 double dy = y2-y1; 00078 return dx * dx + dy * dy; 00079 } 00080 00081 //------------------------------------------------calc_line_point_distance 00082 AGG_INLINE double calc_line_point_distance(double x1, double y1, 00083 double x2, double y2, 00084 double x, double y) 00085 { 00086 double dx = x2-x1; 00087 double dy = y2-y1; 00088 double d = sqrt(dx * dx + dy * dy); 00089 if(d < vertex_dist_epsilon) 00090 { 00091 return calc_distance(x1, y1, x, y); 00092 } 00093 return ((x - x2) * dy - (y - y2) * dx) / d; 00094 } 00095 00096 //-------------------------------------------------------calc_line_point_u 00097 AGG_INLINE double calc_segment_point_u(double x1, double y1, 00098 double x2, double y2, 00099 double x, double y) 00100 { 00101 double dx = x2 - x1; 00102 double dy = y2 - y1; 00103 00104 if(dx == 0 && dy == 0) 00105 { 00106 return 0; 00107 } 00108 00109 double pdx = x - x1; 00110 double pdy = y - y1; 00111 00112 return (pdx * dx + pdy * dy) / (dx * dx + dy * dy); 00113 } 00114 00115 //---------------------------------------------calc_line_point_sq_distance 00116 AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1, 00117 double x2, double y2, 00118 double x, double y, 00119 double u) 00120 { 00121 if(u <= 0) 00122 { 00123 return calc_sq_distance(x, y, x1, y1); 00124 } 00125 else 00126 if(u >= 1) 00127 { 00128 return calc_sq_distance(x, y, x2, y2); 00129 } 00130 return calc_sq_distance(x, y, x1 + u * (x2 - x1), y1 + u * (y2 - y1)); 00131 } 00132 00133 //---------------------------------------------calc_line_point_sq_distance 00134 AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1, 00135 double x2, double y2, 00136 double x, double y) 00137 { 00138 return 00139 calc_segment_point_sq_distance( 00140 x1, y1, x2, y2, x, y, 00141 calc_segment_point_u(x1, y1, x2, y2, x, y)); 00142 } 00143 00144 //-------------------------------------------------------calc_intersection 00145 AGG_INLINE bool calc_intersection(double ax, double ay, double bx, double by, 00146 double cx, double cy, double dx, double dy, 00147 double* x, double* y) 00148 { 00149 double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy); 00150 double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx); 00151 if(fabs(den) < intersection_epsilon) return false; 00152 double r = num / den; 00153 *x = ax + r * (bx-ax); 00154 *y = ay + r * (by-ay); 00155 return true; 00156 } 00157 00158 //-----------------------------------------------------intersection_exists 00159 AGG_INLINE bool intersection_exists(double x1, double y1, double x2, double y2, 00160 double x3, double y3, double x4, double y4) 00161 { 00162 // It's less expensive but you can't control the 00163 // boundary conditions: Less or LessEqual 00164 double dx1 = x2 - x1; 00165 double dy1 = y2 - y1; 00166 double dx2 = x4 - x3; 00167 double dy2 = y4 - y3; 00168 return ((x3 - x2) * dy1 - (y3 - y2) * dx1 < 0.0) != 00169 ((x4 - x2) * dy1 - (y4 - y2) * dx1 < 0.0) && 00170 ((x1 - x4) * dy2 - (y1 - y4) * dx2 < 0.0) != 00171 ((x2 - x4) * dy2 - (y2 - y4) * dx2 < 0.0); 00172 00173 // It's is more expensive but more flexible 00174 // in terms of boundary conditions. 00175 //-------------------- 00176 //double den = (x2-x1) * (y4-y3) - (y2-y1) * (x4-x3); 00177 //if(fabs(den) < intersection_epsilon) return false; 00178 //double nom1 = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3); 00179 //double nom2 = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3); 00180 //double ua = nom1 / den; 00181 //double ub = nom2 / den; 00182 //return ua >= 0.0 && ua <= 1.0 && ub >= 0.0 && ub <= 1.0; 00183 } 00184 00185 //--------------------------------------------------------calc_orthogonal 00186 AGG_INLINE void calc_orthogonal(double thickness, 00187 double x1, double y1, 00188 double x2, double y2, 00189 double* x, double* y) 00190 { 00191 double dx = x2 - x1; 00192 double dy = y2 - y1; 00193 double d = sqrt(dx*dx + dy*dy); 00194 *x = thickness * dy / d; 00195 *y = -thickness * dx / d; 00196 } 00197 00198 //--------------------------------------------------------dilate_triangle 00199 AGG_INLINE void dilate_triangle(double x1, double y1, 00200 double x2, double y2, 00201 double x3, double y3, 00202 double *x, double* y, 00203 double d) 00204 { 00205 double dx1=0.0; 00206 double dy1=0.0; 00207 double dx2=0.0; 00208 double dy2=0.0; 00209 double dx3=0.0; 00210 double dy3=0.0; 00211 double loc = cross_product(x1, y1, x2, y2, x3, y3); 00212 if(fabs(loc) > intersection_epsilon) 00213 { 00214 if(cross_product(x1, y1, x2, y2, x3, y3) > 0.0) 00215 { 00216 d = -d; 00217 } 00218 calc_orthogonal(d, x1, y1, x2, y2, &dx1, &dy1); 00219 calc_orthogonal(d, x2, y2, x3, y3, &dx2, &dy2); 00220 calc_orthogonal(d, x3, y3, x1, y1, &dx3, &dy3); 00221 } 00222 *x++ = x1 + dx1; *y++ = y1 + dy1; 00223 *x++ = x2 + dx1; *y++ = y2 + dy1; 00224 *x++ = x2 + dx2; *y++ = y2 + dy2; 00225 *x++ = x3 + dx2; *y++ = y3 + dy2; 00226 *x++ = x3 + dx3; *y++ = y3 + dy3; 00227 *x++ = x1 + dx3; *y++ = y1 + dy3; 00228 } 00229 00230 //------------------------------------------------------calc_triangle_area 00231 AGG_INLINE double calc_triangle_area(double x1, double y1, 00232 double x2, double y2, 00233 double x3, double y3) 00234 { 00235 return (x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3) * 0.5; 00236 } 00237 00238 //-------------------------------------------------------calc_polygon_area 00239 template<class Storage> double calc_polygon_area(const Storage& st) 00240 { 00241 unsigned i; 00242 double sum = 0.0; 00243 double x = st[0].x; 00244 double y = st[0].y; 00245 double xs = x; 00246 double ys = y; 00247 00248 for(i = 1; i < st.size(); i++) 00249 { 00250 const typename Storage::value_type& v = st[i]; 00251 sum += x * v.y - y * v.x; 00252 x = v.x; 00253 y = v.y; 00254 } 00255 return (sum + x * ys - y * xs) * 0.5; 00256 } 00257 00258 //------------------------------------------------------------------------ 00259 // Tables for fast sqrt 00260 extern int16u g_sqrt_table[1024]; 00261 extern int8 g_elder_bit_table[256]; 00262 00263 00264 //---------------------------------------------------------------fast_sqrt 00265 //Fast integer Sqrt - really fast: no cycles, divisions or multiplications 00266 #if defined(_MSC_VER) 00267 #pragma warning(push) 00268 #pragma warning(disable : 4035) //Disable warning "no return value" 00269 #endif 00270 AGG_INLINE unsigned fast_sqrt(unsigned val) 00271 { 00272 #if defined(_M_IX86) && defined(_MSC_VER) && !defined(AGG_NO_ASM) 00273 //For Ix86 family processors this assembler code is used. 00274 //The key command here is bsr - determination the number of the most 00275 //significant bit of the value. For other processors 00276 //(and maybe compilers) the pure C "#else" section is used. 00277 __asm 00278 { 00279 mov ebx, val 00280 mov edx, 11 00281 bsr ecx, ebx 00282 sub ecx, 9 00283 jle less_than_9_bits 00284 shr ecx, 1 00285 adc ecx, 0 00286 sub edx, ecx 00287 shl ecx, 1 00288 shr ebx, cl 00289 less_than_9_bits: 00290 xor eax, eax 00291 mov ax, g_sqrt_table[ebx*2] 00292 mov ecx, edx 00293 shr eax, cl 00294 } 00295 #else 00296 00297 //This code is actually pure C and portable to most 00298 //arcitectures including 64bit ones. 00299 unsigned t = val; 00300 int bit=0; 00301 unsigned shift = 11; 00302 00303 //The following piece of code is just an emulation of the 00304 //Ix86 assembler command "bsr" (see above). However on old 00305 //Intels (like Intel MMX 233MHz) this code is about twice 00306 //faster (sic!) then just one "bsr". On PIII and PIV the 00307 //bsr is optimized quite well. 00308 bit = t >> 24; 00309 if(bit) 00310 { 00311 bit = g_elder_bit_table[bit] + 24; 00312 } 00313 else 00314 { 00315 bit = (t >> 16) & 0xFF; 00316 if(bit) 00317 { 00318 bit = g_elder_bit_table[bit] + 16; 00319 } 00320 else 00321 { 00322 bit = (t >> 8) & 0xFF; 00323 if(bit) 00324 { 00325 bit = g_elder_bit_table[bit] + 8; 00326 } 00327 else 00328 { 00329 bit = g_elder_bit_table[t]; 00330 } 00331 } 00332 } 00333 00334 //This code calculates the sqrt. 00335 bit -= 9; 00336 if(bit > 0) 00337 { 00338 bit = (bit >> 1) + (bit & 1); 00339 shift -= bit; 00340 val >>= (bit << 1); 00341 } 00342 return g_sqrt_table[val] >> shift; 00343 #endif 00344 } 00345 #if defined(_MSC_VER) 00346 #pragma warning(pop) 00347 #endif 00348 00349 00350 00351 00352 //--------------------------------------------------------------------besj 00353 // Function BESJ calculates Bessel function of first kind of order n 00354 // Arguments: 00355 // n - an integer (>=0), the order 00356 // x - value at which the Bessel function is required 00357 //-------------------- 00358 // C++ Mathematical Library 00359 // Convereted from equivalent FORTRAN library 00360 // Converetd by Gareth Walker for use by course 392 computational project 00361 // All functions tested and yield the same results as the corresponding 00362 // FORTRAN versions. 00363 // 00364 // If you have any problems using these functions please report them to 00365 // M.Muldoon@UMIST.ac.uk 00366 // 00367 // Documentation available on the web 00368 // http://www.ma.umist.ac.uk/mrm/Teaching/392/libs/392.html 00369 // Version 1.0 8/98 00370 // 29 October, 1999 00371 //-------------------- 00372 // Adapted for use in AGG library by Andy Wilk (castor.vulgaris@gmail.com) 00373 //------------------------------------------------------------------------ 00374 inline double besj(double x, int n) 00375 { 00376 if(n < 0) 00377 { 00378 return 0; 00379 } 00380 double d = 1E-6; 00381 double b = 0; 00382 if(fabs(x) <= d) 00383 { 00384 if(n != 0) return 0; 00385 return 1; 00386 } 00387 double b1 = 0; // b1 is the value from the previous iteration 00388 // Set up a starting order for recurrence 00389 int m1 = (int)fabs(x) + 6; 00390 if(fabs(x) > 5) 00391 { 00392 m1 = (int)(fabs(1.4 * x + 60 / x)); 00393 } 00394 int m2 = (int)(n + 2 + fabs(x) / 4); 00395 if (m1 > m2) 00396 { 00397 m2 = m1; 00398 } 00399 00400 // Apply recurrence down from curent max order 00401 for(;;) 00402 { 00403 double c3 = 0; 00404 double c2 = 1E-30; 00405 double c4 = 0; 00406 int m8 = 1; 00407 if (m2 / 2 * 2 == m2) 00408 { 00409 m8 = -1; 00410 } 00411 int imax = m2 - 2; 00412 for (int i = 1; i <= imax; i++) 00413 { 00414 double c6 = 2 * (m2 - i) * c2 / x - c3; 00415 c3 = c2; 00416 c2 = c6; 00417 if(m2 - i - 1 == n) 00418 { 00419 b = c6; 00420 } 00421 m8 = -1 * m8; 00422 if (m8 > 0) 00423 { 00424 c4 = c4 + 2 * c6; 00425 } 00426 } 00427 double c6 = 2 * c2 / x - c3; 00428 if(n == 0) 00429 { 00430 b = c6; 00431 } 00432 c4 += c6; 00433 b /= c4; 00434 if(fabs(b - b1) < d) 00435 { 00436 return b; 00437 } 00438 b1 = b; 00439 m2 += 3; 00440 } 00441 } 00442 00443 } 00444 00445 00446 #endif