1 //---------------------------------------------------------------------------- 2 // Anti-Grain Geometry - Version 2.2 3 // Copyright (C) 2002-2004 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 16 #ifndef AGG_MATH_INCLUDED 17 #define AGG_MATH_INCLUDED 18 19 #include <math.h> 20 #include "agg_basics.h" 21 22 namespace agg 23 { 24 25 const double intersection_epsilon = 1.0e-8; 26 27 //------------------------------------------------------calc_point_location 28 inline double calc_point_location(double x1, double y1, 29 double x2, double y2, 30 double x, double y) 31 { 32 return (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1); 33 } 34 35 36 //--------------------------------------------------------point_in_triangle 37 inline bool point_in_triangle(double x1, double y1, 38 double x2, double y2, 39 double x3, double y3, 40 double x, double y) 41 { 42 bool cp1 = calc_point_location(x1, y1, x2, y2, x, y) < 0.0; 43 bool cp2 = calc_point_location(x2, y2, x3, y3, x, y) < 0.0; 44 bool cp3 = calc_point_location(x3, y3, x1, y1, x, y) < 0.0; 45 return cp1 == cp2 && cp2 == cp3 && cp3 == cp1; 46 } 47 48 49 //-----------------------------------------------------------calc_distance 50 inline double calc_distance(double x1, double y1, double x2, double y2) 51 { 52 double dx = x2-x1; 53 double dy = y2-y1; 54 return sqrt(dx * dx + dy * dy); 55 } 56 57 58 //------------------------------------------------calc_point_line_distance 59 inline double calc_point_line_distance(double x1, double y1, 60 double x2, double y2, 61 double x, double y) 62 { 63 double dx = x2-x1; 64 double dy = y2-y1; 65 return ((x - x2) * dy - (y - y2) * dx) / sqrt(dx * dx + dy * dy); 66 } 67 68 69 //-------------------------------------------------------calc_intersection 70 inline bool calc_intersection(double ax, double ay, double bx, double by, 71 double cx, double cy, double dx, double dy, 72 double* x, double* y) 73 { 74 double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy); 75 double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx); 76 if(fabs(den) < intersection_epsilon) return false; 77 double r = num / den; 78 *x = ax + r * (bx-ax); 79 *y = ay + r * (by-ay); 80 return true; 81 } 82 83 84 //--------------------------------------------------------calc_orthogonal 85 inline void calc_orthogonal(double thickness, 86 double x1, double y1, 87 double x2, double y2, 88 double* x, double* y) 89 { 90 double dx = x2 - x1; 91 double dy = y2 - y1; 92 double d = sqrt(dx*dx + dy*dy); 93 *x = thickness * dy / d; 94 *y = thickness * dx / d; 95 } 96 97 98 //--------------------------------------------------------dilate_triangle 99 inline void dilate_triangle(double x1, double y1, 100 double x2, double y2, 101 double x3, double y3, 102 double *x, double* y, 103 double d) 104 { 105 double dx1=0.0; 106 double dy1=0.0; 107 double dx2=0.0; 108 double dy2=0.0; 109 double dx3=0.0; 110 double dy3=0.0; 111 double loc = calc_point_location(x1, y1, x2, y2, x3, y3); 112 if(fabs(loc) > intersection_epsilon) 113 { 114 if(calc_point_location(x1, y1, x2, y2, x3, y3) > 0.0) 115 { 116 d = -d; 117 } 118 calc_orthogonal(d, x1, y1, x2, y2, &dx1, &dy1); 119 calc_orthogonal(d, x2, y2, x3, y3, &dx2, &dy2); 120 calc_orthogonal(d, x3, y3, x1, y1, &dx3, &dy3); 121 } 122 *x++ = x1 + dx1; *y++ = y1 - dy1; 123 *x++ = x2 + dx1; *y++ = y2 - dy1; 124 *x++ = x2 + dx2; *y++ = y2 - dy2; 125 *x++ = x3 + dx2; *y++ = y3 - dy2; 126 *x++ = x3 + dx3; *y++ = y3 - dy3; 127 *x++ = x1 + dx3; *y++ = y1 - dy3; 128 } 129 130 //-------------------------------------------------------calc_polygon_area 131 template<class Storage> double calc_polygon_area(const Storage& st) 132 { 133 unsigned i; 134 double sum = 0.0; 135 double x = st[0].x; 136 double y = st[0].y; 137 double xs = x; 138 double ys = y; 139 140 for(i = 1; i < st.size(); i++) 141 { 142 const typename Storage::value_type& v = st[i]; 143 sum += x * v.y - y * v.x; 144 x = v.x; 145 y = v.y; 146 } 147 return (sum + x * ys - y * xs) * 0.5; 148 } 149 150 //------------------------------------------------------------------------ 151 // Tables for fast sqrt 152 extern int16u g_sqrt_table[1024]; 153 extern int8 g_elder_bit_table[256]; 154 155 156 //---------------------------------------------------------------fast_sqrt 157 //Fast integer Sqrt - really fast: no cycles, divisions or multiplications 158 #if defined(_MSC_VER) 159 #pragma warning(push) 160 #pragma warning(disable : 4035) //Disable warning "no return value" 161 #endif 162 inline unsigned fast_sqrt(unsigned val) 163 { 164 #if defined(_M_IX86) && defined(_MSC_VER) && !defined(AGG_NO_ASM) 165 //For Ix86 family processors this assembler code is used. 166 //The key command here is bsr - determination the number of the most 167 //significant bit of the value. For other processors 168 //(and maybe compilers) the pure C "#else" section is used. 169 __asm 170 { 171 mov ebx, val 172 mov edx, 11 173 bsr ecx, ebx 174 sub ecx, 9 175 jle less_than_9_bits 176 shr ecx, 1 177 adc ecx, 0 178 sub edx, ecx 179 shl ecx, 1 180 shr ebx, cl 181 less_than_9_bits: 182 xor eax, eax 183 mov ax, g_sqrt_table[ebx*2] 184 mov ecx, edx 185 shr eax, cl 186 } 187 #else 188 189 //This code is actually pure C and portable to most 190 //arcitectures including 64bit ones. 191 unsigned t = val; 192 int bit=0; 193 unsigned shift = 11; 194 195 //The following piece of code is just an emulation of the 196 //Ix86 assembler command "bsr" (see above). However on old 197 //Intels (like Intel MMX 233MHz) this code is about twice 198 //faster (sic!) then just one "bsr". On PIII and PIV the 199 //bsr is optimized quite well. 200 bit = t >> 24; 201 if(bit) 202 { 203 bit = g_elder_bit_table[bit] + 24; 204 } 205 else 206 { 207 bit = (t >> 16) & 0xFF; 208 if(bit) 209 { 210 bit = g_elder_bit_table[bit] + 16; 211 } 212 else 213 { 214 bit = (t >> 8) & 0xFF; 215 if(bit) 216 { 217 bit = g_elder_bit_table[bit] + 8; 218 } 219 else 220 { 221 bit = g_elder_bit_table[t]; 222 } 223 } 224 } 225 226 //This is calculation sqrt itself. 227 bit -= 9; 228 if(bit > 0) 229 { 230 bit = (bit >> 1) + (bit & 1); 231 shift -= bit; 232 val >>= (bit << 1); 233 } 234 return g_sqrt_table[val] >> shift; 235 #endif 236 } 237 #if defined(_MSC_VER) 238 #pragma warning(pop) 239 #endif 240 241 242 243 244 } 245 246 247 #endif 248