agg/inc/agg_math.h

   1 //----------------------------------------------------------------------------
   2 // Anti-Grain Geometry - Version 2.3
   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
  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     AGG_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     AGG_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     AGG_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     AGG_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     AGG_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     AGG_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     AGG_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     AGG_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