cartesian_predicates_3.h

很多二维 三维几何计算算法 C++ 类库

// Copyright (c) 2005  Stanford University (USA).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Kinetic_data_structures/include/CGAL/Kinetic/internal/Kernel/cartesian_predicates_3.h $
// $Id: cartesian_predicates_3.h 37180 2007-03-17 09:00:08Z afabri $
// 
//
// Author(s)     : Daniel Russel <[e-mail]>

#ifndef CGAL_KINETIC_CARTESIAN_PREDICATES_3_H
#define CGAL_KINETIC_CARTESIAN_PREDICATES_3_H
#include <CGAL/Kinetic/basic.h>
#include <CGAL/determinant.h>

//#include <CGAL/Kinetic/kernel/Cartesian_moving_sphere_3.h>

CGAL_KINETIC_BEGIN_INTERNAL_NAMESPACE



template <class KK>
struct Cartesian_side_of_oriented_sphere_3
{
  Cartesian_side_of_oriented_sphere_3(){}

  typedef typename KK::Certificate_function result_type;
  typedef typename KK::Point_3 first_argument_type;
  typedef typename KK::Point_3 second_argument_type;
  typedef typename KK::Point_3 third_argument_type;
  typedef typename KK::Point_3 fourth_argument_type;
  typedef typename KK::Point_3 fifth_argument_type;
  result_type operator()(const first_argument_type &a,
			 const second_argument_type &b,
			 const third_argument_type &c,
			 const fourth_argument_type &d,
			 const fourth_argument_type &e)const
  {
    typedef typename KK::Motion_function FT;
    // We translate the points so that T becomes the origin.
    FT dpx = (a.x()) - (e.x());
    FT dpy = (a.y()) - (e.y());
    FT dpz = (a.z()) - (e.z());
    FT dpt = dpx*dpx + dpy*dpy + dpz*dpz;
    FT dqx = (b.x()) - (e.x());
    FT dqy = (b.y()) - (e.y());
    FT dqz = (b.z()) - (e.z());
    FT dqt = dqx*dqx + dqy*dqy + dqz*dqz;
    FT drx = (c.x()) - (e.x());
    FT dry = (c.y()) - (e.y());
    FT drz = (c.z()) - (e.z());
    FT drt = drx*drx + dry*dry + drz*drz;
    FT dsx = (d.x()) - (e.x());
    FT dsy = (d.y()) - (e.y());
    FT dsz = (d.z()) - (e.z());
    FT dst = dsx*dsx + dsy*dsy + dsz*dsz;
    FT ret= CGAL::det4x4_by_formula(dpx, dpy, dpz, dpt,
				    dqx, dqy, dqz, dqt,
				    drx, dry, drz, drt,
				    dsx, dsy, dsz, dst);
    /*  CGAL_KINETIC_MAPLE_LOG( << std::endl << std::endl;
	CGAL_KINETIC_MAPLE_LOG( << "pt3\n";
	CGAL_KINETIC_MAPLE_LOG( << "m:=matrix(5,5,[[";
	CGAL_KINETIC_MAPLE_LOG(
	<< a.x() << ", " << a.y() << ", " << a.z() << ", " << a.x()*a.x()+a.y()*a.y()+a.z()*a.z() << ", 1], ["
	<< a.x() << ", " << b.y() << ", " << b.z() << ", " << b.x()*b.x()+b.y()*b.y()+b.z()*b.z() << ", 1], ["
	<< c.x() << ", " << c.y() << ", " << c.z() << ", " << c.x()*c.x()+c.y()*c.y()+c.z()*c.z() << ", 1], ["
	<< d.x() << ", " << d.y() << ", " << d.z() << ", " << d.x()*d.x()+d.y()*d.y()+d.z()*a.z() << ", 1], ["
	<< e.x() << ", " << e.y() << ", " << e.z() << ", " << e.x()*e.x()+e.y()*e.y()+e.z()*e.z() << ", 1]]);\n";
	CGAL_KINETIC_MAPLE_LOG( << "det(m)-( " << ret << ");\n";
	::CGAL::Kinetic::log()->stream(::CGAL::KDS::Log::MAPLE) << "m2:= matrix(4,4,[[";
	::CGAL::Kinetic::log()->stream(::CGAL::KDS::Log::MAPLE) << dpx << "," << dpy << "," << dpz << "," << dpt << "], [";
	::CGAL::Kinetic::log()->stream(::CGAL::KDS::Log::MAPLE) << dqx << "," << dqy << "," << dqz << "," << dqt << "], [";
	::CGAL::Kinetic::log()->stream(::CGAL::KDS::Log::MAPLE) << drx << "," << dry << "," << drz << "," << drt << "], [";
	::CGAL::Kinetic::log()->stream(::CGAL::KDS::Log::MAPLE) << dsx << "," << dsy << "," << dsz << "," << dst << "]]);\n";
    */
    return ret;
  }
};

template <class KK>
struct Cartesian_power_test_3
{
  typedef typename KK::Certificate_function result_type;
  typedef typename KK::Weighted_point_3 first_argument_type;
  typedef typename KK::Weighted_point_3 second_argument_type;
  typedef typename KK::Weighted_point_3 third_argument_type;
  typedef typename KK::Weighted_point_3 fourth_argument_type;
  typedef typename KK::Weighted_point_3 fifth_argument_type;
  Cartesian_power_test_3(){}
  result_type operator()(const first_argument_type &a,
			 const second_argument_type &b,
			 const third_argument_type &c,
			 const fourth_argument_type &d,
			 const fourth_argument_type &e)const
  {
    typedef typename KK::Motion_function FT;
    // We translate the points so that T becomes the origin.
    FT dpx = (a.point().x()) - (e.point().x());
    FT dpy = (a.point().y()) - (e.point().y());
    FT dpz = (a.point().z()) - (e.point().z());
    FT dpt = dpx*dpx + dpy*dpy
      + dpz*dpz - (a.weight()) + (e.weight());
    FT dqx = (b.point().x()) - (e.point().x());
    FT dqy = (b.point().y()) - (e.point().y());
    FT dqz = (b.point().z()) - (e.point().z());
    FT dqt = dqx*dqx + dqy*dqy
      + dqz*dqz - (b.weight()) + (e.weight());
    FT drx = (c.point().x()) - (e.point().x());
    FT dry = (c.point().y()) - (e.point().y());
    FT drz = (c.point().z()) - (e.point().z());
    FT drt = drx*drx + dry*dry
      + drz*drz - (c.weight()) + (e.weight());
    FT dsx = (d.point().x()) - (e.point().x());
    FT dsy = (d.point().y()) - (e.point().y());
    FT dsz = (d.point().z()) - (e.point().z());
    FT dst = dsx*dsx + dsy*dsy
      + dsz*dsz - (d.weight()) + (e.weight());

    return ::CGAL::det4x4_by_formula(dpx, dpy, dpz, dpt,
				     dqx, dqy, dqz, dqt,
				     drx, dry, drz, drt,
				     dsx, dsy, dsz, dst);
  }

 result_type operator()(const first_argument_type &,
			const second_argument_type &,
			const third_argument_type &,
			const fourth_argument_type &)const
  {
    CGAL_assertion(0);
    return result_type();
  }
  result_type operator()(const first_argument_type &,
			 const second_argument_type &,
			 const third_argument_type &)const
  {
    CGAL_assertion(0);
    return result_type();
  }
  result_type operator()(const first_argument_type &,
			 const second_argument_type &)const
  {
    CGAL_assertion(0);
    return result_type();
  }
};

// Works for unweighted points.
template <class KK>
struct Cartesian_lifted_power_test_3
{
  Cartesian_lifted_power_test_3(){}
  typedef typename KK::Certificate_function result_type;
  typedef typename KK::Weighted_point_3 first_argument_type;
  typedef typename KK::Weighted_point_3 second_argument_type;
  typedef typename KK::Weighted_point_3 third_argument_type;
  typedef typename KK::Weighted_point_3 fourth_argument_type;
  typedef typename KK::Weighted_point_3 fifth_argument_type;
  result_type operator()(const first_argument_type &a,
			 const second_argument_type &b,
			 const third_argument_type &c,
			 const fourth_argument_type &d,
			 const fourth_argument_type &e)const
  {
    typedef typename KK::Motion_function FT;

    FT dpx = (a.point().x()) - (e.point().x());
    FT dpy = (a.point().y()) - (e.point().y());
    FT dpz = (a.point().z()) - (e.point().z());
    FT dpt = (a.lifted()) - (e.lifted());
    FT dqx = (b.point().x()) - (e.point().x());
    FT dqy = (b.point().y()) - (e.point().y());
    FT dqz = (b.point().z()) - (e.point().z());
    FT dqt = (b.lifted()) - (e.lifted());
    FT drx = (c.point().x()) - (e.point().x());
    FT dry = (c.point().y()) - (e.point().y());
    FT drz = (c.point().z()) - (e.point().z());
    FT drt = (c.lifted()) - (e.lifted());
    FT dsx = (d.point().x()) - (e.point().x());
    FT dsy = (d.point().y()) - (e.point().y());
    FT dsz = (d.point().z()) - (e.point().z());
    FT dst = (d.lifted()) - (e.lifted());

    return CGAL::det4x4_by_formula(dpx, dpy, dpz, dpt,
				   dqx, dqy, dqz, dqt,
				   drx, dry, drz, drt,
				   dsx, dsy, dsz, dst);
  }
};

template <class KK>
struct Cartesian_linear_lifted_power_test_3
{
  Cartesian_linear_lifted_power_test_3(){}
  typedef typename KK::Certificate_function result_type;
  typedef typename KK::Weighted_point_3 first_argument_type;
  typedef typename KK::Weighted_point_3 second_argument_type;
  typedef typename KK::Weighted_point_3 third_argument_type;
  typedef typename KK::Weighted_point_3 fourth_argument_type;
  typedef typename KK::Weighted_point_3 fifth_argument_type;
  result_type operator()(first_argument_type a,
			 second_argument_type b,
			 third_argument_type c,
			 fourth_argument_type d,
			 fourth_argument_type e)const
  {
    typedef typename KK::Motion_function FT;
    typedef typename FT::NT NT;
    /*if (a.is_constant()){
      std::swap(a,e);
      std::swap(b,c);
      }*/

    //typedef typename RT::NT FT;
    // We translate the points so that T becomes the origin.
    FT safe_ret;

    {
      FT dpx = (a.point().x()) - (e.point().x());
      FT dpy = (a.point().y()) - (e.point().y());
      FT dpz = (a.point().z()) - (e.point().z());
      FT dpt = (a.lifted()) - (e.lifted());
      FT dqx = (b.point().x()) - (e.point().x());
      FT dqy = (b.point().y()) - (e.point().y());
      FT dqz = (b.point().z()) - (e.point().z());
      FT dqt = (b.lifted()) - (e.lifted());
      FT drx = (c.point().x()) - (e.point().x());
      FT dry = (c.point().y()) - (e.point().y());
      FT drz = (c.point().z()) - (e.point().z());
      FT drt = (c.lifted()) - (e.lifted());
      FT dsx = (d.point().x()) - (e.point().x());
      FT dsy = (d.point().y()) - (e.point().y());
      FT dsz = (d.point().z()) - (e.point().z());
      FT dst = (d.lifted()) - (e.lifted());

      safe_ret = CGAL::det4x4_by_formula(dpx, dpy, dpz, dpt,
					 dqx, dqy, dqz, dqt,
					 drx, dry, drz, drt,
					 dsx, dsy, dsz, dst);
    }
    FT ret;
    {
      if (!b.is_constant()) {
	std::swap(a,b);
	std::swap(c,d);
      }
      else if (! c.is_constant()) {
	std::swap(a,c);
	std::swap(b,d);
      }
      else if (! d.is_constant()) {
	std::swap(a,d);
	std::swap(b,c);
      }
      else if (! e.is_constant()) {
	std::swap(a,e);
	std::swap(b,c);
      }
      CGAL_precondition(b.is_constant());
      CGAL_precondition(c.is_constant());
      CGAL_precondition(d.is_constant());
      CGAL_precondition(e.is_constant());

      FT a00 = a.point().x() - e.point().x();
      FT a01 = a.point().y() - e.point().y();
      FT a02 = a.point().z() - e.point().z();
      FT a03 = a.lifted() - e.lifted();
      NT a10 = b.point().x()[0] - e.point().x()[0];
      NT a11 = b.point().y()[0] - e.point().y()[0];
      NT a12 = b.point().z()[0] - e.point().z()[0];
      NT a13 = b.lifted()[0] - e.lifted()[0];
      NT a20 = c.point().x()[0] - e.point().x()[0];
      NT a21 = c.point().y()[0] - e.point().y()[0];
      NT a22 = c.point().z()[0] - e.point().z()[0];
      NT a23 = c.lifted()[0] - e.lifted()[0];
      NT a30 = d.point().x()[0] - e.point().x()[0];
      NT a31 = d.point().y()[0] - e.point().y()[0];
      NT a32 = d.point().z()[0] - e.point().z()[0];
      NT a33 = d.lifted()[0] - e.lifted()[0];
      // First ompute the det2x2
      const FT m01 = a10*a01 - a00*a11;
      const FT m02 = a20*a01 - a00*a21;
      const FT m03 = a30*a01 - a00*a31;
      const NT m12 = a20*a11 - a10*a21;
      const NT m13 = a30*a11 - a10*a31;
      const NT m23 = a30*a21 - a20*a31;
      // Now compute the minors of rank 3
      const FT m012 = m12*a02 - m02*a12 + m01*a22;
      const FT m013 = m13*a02 - m03*a12 + m01*a32;
      const FT m023 = m23*a02 - m03*a22 + m02*a32;
      const NT m123 = m23*a12 - m13*a22 + m12*a32;
      // Now compute the minors of rank 4
      const FT m0123 = m123*a03 - m023*a13 + m013*a23 - m012*a33;
      ret= m0123;
    }
    if (ret != safe_ret) {
      ret.print();
      safe_ret.print();
      std::cout << a << std::endl;
      std::cout << b << std::endl;
      std::cout << c << std::endl;
      std::cout << d << std::endl;
      std::cout << e << std::endl;
    }
    CGAL_exactness_postcondition(ret==safe_ret);
    return ret;
  }
};

template <class Pt>
typename Pt::Coordinate co3(const Pt &a, const Pt &b, const Pt &c, const Pt &d)
{
  //std::cout << "Computing orientation of matrix(4,4, [[" << a << ",1], [" << b << ",1], [" << c << ",1], [" << d << ",1]]);\n";
  typedef typename Pt::Coordinate RT;

  RT px= (a.x());
  RT py= (a.y());
  RT pz= (a.z());
  RT qx= (b.x());
  RT qy= (b.y());
  RT qz= (b.z());
  RT rx= (c.x());
  RT ry= (c.y());
  RT rz= (c.z());
  RT sx= (d.x());
  RT sy= (d.y());
  RT sz= (d.z());
  if (d.is_constant()) {
    std::swap(px, sx);
    std::swap(py, sy);
    std::swap(pz, sz);
    std::swap(rx, qx);
    std::swap(ry, qy);
    std::swap(rz, qz);
  }

  RT a00= qx-px;
  RT a01= rx-px;
  RT a02= sx-px;
  RT a10= qy-py;
  RT a11= ry-py;
  RT a12= sy-py;
  RT a20= qz-pz;
  RT a21= rz-pz;
  RT a22= sz-pz;
  RT ret= CGAL::det3x3_by_formula(a00, a01, a02,
				  a10, a11, a12,
				  a20, a21, a22);
  /*CGAL_KINETIC_LOG_MAPLE(std::endl << std::endl);
    CGAL_KINETIC_LOG_MAPLE( "co3\n");
    CGAL_KINETIC_LOG_MAPLE( "m:= matrix(4,4, [[1," << a.x()<<","<<a.y() <<","<<a.z() << "], [1," << b.x()<<","<<b.y()<<","<<b.z());
    CGAL_KINETIC_LOG_MAPLE( "], [1," << c.x() <<","<<c.y()<<","<<c.z() << "], [1,");
    CGAL_KINETIC_LOG_MAPLE( d.x() << ","<<d.y()<<","<<d.z() << "]]);"<< std::endl);
    CGAL_KINETIC_LOG_MAPLE( "det(m)- (" << ret << ");" << std::endl);
    CGAL_KINETIC_LOG_MAPLE( "m2:= matrix(3,3,[[,");
    CGAL_KINETIC_LOG_MAPLE( a00 << "," << a01 << "," << a02 << "], [");
    CGAL_KINETIC_LOG_MAPLE( a10 << "," << a11 << "," << a12 << "], [");
    CGAL_KINETIC_LOG_MAPLE( a20 << "," << a21 << "," << a22 << "]]);");*/

  //std::cout << "returning " << ret << std::endl;
  return ret;
}




template <class KK>
struct Cartesian_orientation_3
{
  Cartesian_orientation_3(){}

  typedef typename KK::Certificate_function result_type;

  typedef typename KK::Point_3 first_argument_type;
  typedef typename KK::Point_3 second_argument_type;
  typedef typename KK::Point_3 third_argument_type;
  typedef typename KK::Point_3 fourth_argument_type;
  result_type operator()(const first_argument_type &a,
			 const second_argument_type &b,
			 const third_argument_type &c,
			 const fourth_argument_type &d) const
  {
    return co3(a, b, c, d);
  }
};

template <class KK>