agg_math.h

windows ce 下的画各种b样条曲线

//----------------------------------------------------------------------------
// Anti-Grain Geometry (AGG) - Version 2.5
// A high quality rendering engine for C++
// Copyright (C) 2002-2006 Maxim Shemanarev
// Contact: [e-mail]
//          [e-mail]
//          http://antigrain.com
// 
// AGG is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
// 
// AGG is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with AGG; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 
// MA 02110-1301, USA.
//----------------------------------------------------------------------------
// Bessel function (besj) was adapted for use in AGG library by Andy Wilk 
// Contact: [e-mail]
//----------------------------------------------------------------------------

#ifndef AGG_MATH_INCLUDED
#define AGG_MATH_INCLUDED

#include <math.h>
#include "agg_basics.h"

namespace agg
{

    //------------------------------------------------------vertex_dist_epsilon
    // Coinciding points maximal distance (Epsilon)
    const double vertex_dist_epsilon = 1e-14;

    //-----------------------------------------------------intersection_epsilon
    // See calc_intersection
    const double intersection_epsilon = 1.0e-30;

    //------------------------------------------------------------cross_product
    AGG_INLINE double cross_product(double x1, double y1, 
                                    double x2, double y2, 
                                    double x,  double y)
    {
        return (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1);
    }

    //--------------------------------------------------------point_in_triangle
    AGG_INLINE bool point_in_triangle(double x1, double y1, 
                                      double x2, double y2, 
                                      double x3, double y3, 
                                      double x,  double y)
    {
        bool cp1 = cross_product(x1, y1, x2, y2, x, y) < 0.0;
        bool cp2 = cross_product(x2, y2, x3, y3, x, y) < 0.0;
        bool cp3 = cross_product(x3, y3, x1, y1, x, y) < 0.0;
        return cp1 == cp2 && cp2 == cp3 && cp3 == cp1;
    }

    //-----------------------------------------------------------calc_distance
    AGG_INLINE double calc_distance(double x1, double y1, double x2, double y2)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        return sqrt(dx * dx + dy * dy);
    }

    //--------------------------------------------------------calc_sq_distance
    AGG_INLINE double calc_sq_distance(double x1, double y1, double x2, double y2)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        return dx * dx + dy * dy;
    }

    //------------------------------------------------calc_line_point_distance
    AGG_INLINE double calc_line_point_distance(double x1, double y1, 
                                               double x2, double y2, 
                                               double x,  double y)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        double d = sqrt(dx * dx + dy * dy);
        if(d < vertex_dist_epsilon)
        {
            return calc_distance(x1, y1, x, y);
        }
        return ((x - x2) * dy - (y - y2) * dx) / d;
    }

    //-------------------------------------------------------calc_line_point_u
    AGG_INLINE double calc_segment_point_u(double x1, double y1, 
                                           double x2, double y2, 
                                           double x,  double y)
    {
        double dx = x2 - x1;
        double dy = y2 - y1;

        if(dx == 0 && dy == 0)
        {
	        return 0;
        }

        double pdx = x - x1;
        double pdy = y - y1;

        return (pdx * dx + pdy * dy) / (dx * dx + dy * dy);
    }

    //---------------------------------------------calc_line_point_sq_distance
    AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1, 
                                                     double x2, double y2, 
                                                     double x,  double y,
                                                     double u)
    {
        if(u <= 0)
        {
	        return calc_sq_distance(x, y, x1, y1);
        }
        else 
        if(u >= 1)
        {
	        return calc_sq_distance(x, y, x2, y2);
        }
        return calc_sq_distance(x, y, x1 + u * (x2 - x1), y1 + u * (y2 - y1));
    }

    //---------------------------------------------calc_line_point_sq_distance
    AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1, 
                                                     double x2, double y2, 
                                                     double x,  double y)
    {
        return 
            calc_segment_point_sq_distance(
                x1, y1, x2, y2, x, y,
                calc_segment_point_u(x1, y1, x2, y2, x, y));
    }

    //-------------------------------------------------------calc_intersection
    AGG_INLINE bool calc_intersection(double ax, double ay, double bx, double by,
                                      double cx, double cy, double dx, double dy,
                                      double* x, double* y)
    {
        double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy);
        double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx);
        if(fabs(den) < intersection_epsilon) return false;
        double r = num / den;
        *x = ax + r * (bx-ax);
        *y = ay + r * (by-ay);
        return true;
    }

    //-----------------------------------------------------intersection_exists
    AGG_INLINE bool intersection_exists(double x1, double y1, double x2, double y2,
                                        double x3, double y3, double x4, double y4)
    {
        // It's less expensive but you can't control the 
        // boundary conditions: Less or LessEqual
        double dx1 = x2 - x1;
        double dy1 = y2 - y1;
        double dx2 = x4 - x3;
        double dy2 = y4 - y3;
        return ((x3 - x2) * dy1 - (y3 - y2) * dx1 < 0.0) != 
               ((x4 - x2) * dy1 - (y4 - y2) * dx1 < 0.0) &&
               ((x1 - x4) * dy2 - (y1 - y4) * dx2 < 0.0) !=
               ((x2 - x4) * dy2 - (y2 - y4) * dx2 < 0.0);

        // It's is more expensive but more flexible 
        // in terms of boundary conditions.
        //--------------------
        //double den  = (x2-x1) * (y4-y3) - (y2-y1) * (x4-x3);
        //if(fabs(den) < intersection_epsilon) return false;
        //double nom1 = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
        //double nom2 = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
        //double ua = nom1 / den;
        //double ub = nom2 / den;
        //return ua >= 0.0 && ua <= 1.0 && ub >= 0.0 && ub <= 1.0;
    }

    //--------------------------------------------------------calc_orthogonal
    AGG_INLINE void calc_orthogonal(double thickness,
                                    double x1, double y1,
                                    double x2, double y2,
                                    double* x, double* y)
    {
        double dx = x2 - x1;
        double dy = y2 - y1;
        double d = sqrt(dx*dx + dy*dy); 
        *x =  thickness * dy / d;
        *y = -thickness * dx / d;
    }

    //--------------------------------------------------------dilate_triangle
    AGG_INLINE void dilate_triangle(double x1, double y1,
                                    double x2, double y2,
                                    double x3, double y3,
                                    double *x, double* y,
                                    double d)
    {
        double dx1=0.0;
        double dy1=0.0; 
        double dx2=0.0;
        double dy2=0.0; 
        double dx3=0.0;
        double dy3=0.0; 
        double loc = cross_product(x1, y1, x2, y2, x3, y3);
        if(fabs(loc) > intersection_epsilon)
        {
            if(cross_product(x1, y1, x2, y2, x3, y3) > 0.0) 
            {
                d = -d;
            }
            calc_orthogonal(d, x1, y1, x2, y2, &dx1, &dy1);
            calc_orthogonal(d, x2, y2, x3, y3, &dx2, &dy2);
            calc_orthogonal(d, x3, y3, x1, y1, &dx3, &dy3);
        }
        *x++ = x1 + dx1;  *y++ = y1 + dy1;
        *x++ = x2 + dx1;  *y++ = y2 + dy1;
        *x++ = x2 + dx2;  *y++ = y2 + dy2;