intersection_3_1_impl.h

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

// Copyright (c) 1997-2004  Utrecht University (The Netherlands),
// ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany),
// INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg
// (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria),
// and Tel-Aviv University (Israel).  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/Intersections_3/include/CGAL/Intersections_3/intersection_3_1_impl.h $
// $Id: intersection_3_1_impl.h 36691 2007-02-28 16:52:55Z spion $
// 
//
// Author(s)     : Geert-Jan Giezeman <[e-mail]>


#include <CGAL/wmult.h>

CGAL_BEGIN_NAMESPACE

namespace CGALi {

template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Plane_3  &plane, 
	     const typename CGAL_WRAP(K)::Line_3 &line, 
	     const K&)
{
    typedef typename K::Point_3 Point_3;
    typedef typename K::Direction_3 Direction_3;
    typedef typename K::RT RT;

    const Point_3 &line_pt = line.point();
    const Direction_3 &line_dir = line.direction();

    RT num = plane.a()*line_pt.hx() + plane.b()*line_pt.hy()
           + plane.c()*line_pt.hz() + wmult_hw((K*)0, plane.d(), line_pt);
    RT den = plane.a()*line_dir.dx() + plane.b()*line_dir.dy()
           + plane.c()*line_dir.dz();
    if (den == 0) {
        if (num == 0) {
            // all line
            return make_object(line);
        } else {
            // no intersection
            return Object();
        }
    }
    return make_object(Point_3(
        den*line_pt.hx()-num*line_dir.dx(),
        den*line_pt.hy()-num*line_dir.dy(),
        den*line_pt.hz()-num*line_dir.dz(),
        wmult_hw((K*)0, den, line_pt)));
}

template <class K>
inline
Object
intersection(const typename CGAL_WRAP(K)::Line_3 &line, 
	     const typename CGAL_WRAP(K)::Plane_3  &plane, 
	     const K& k)
{
  return intersection(plane, line, k);
}

template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Plane_3 &plane1, 
	     const typename CGAL_WRAP(K)::Plane_3 &plane2, 
	     const K&)
{
  typedef typename K::Point_3 Point_3;
  typedef typename K::Direction_3 Direction_3;
  typedef typename K::Line_3 Line_3;

    typedef typename K::RT RT;
    const RT &a = plane1.a();
    const RT &b = plane1.b();
    const RT &c = plane1.c();
    const RT &d = plane1.d();
    const RT &p = plane2.a();
    const RT &q = plane2.b();
    const RT &r = plane2.c();
    const RT &s = plane2.d();

    RT det = a*q-p*b;
    if (det != 0) {
        Point_3 is_pt = Point_3(b*s-d*q, p*d-a*s, 0, det);
        Direction_3 is_dir = Direction_3(b*r-c*q, p*c-a*r, det);
        return make_object(Line_3(is_pt, is_dir));
    }
    det = a*r-p*c;
    if (det != 0) {
        Point_3 is_pt = Point_3(c*s-d*r, 0, p*d-a*s, det);
        Direction_3 is_dir = Direction_3(c*q-b*r, det, p*b-a*q);
        return make_object(Line_3(is_pt, is_dir));
    }
    det = b*r-c*q;
    if (det != 0) {
        Point_3 is_pt = Point_3(0, c*s-d*r, d*q-b*s, det);
        Direction_3 is_dir = Direction_3(det, c*p-a*r, a*q-b*p);
        return make_object(Line_3(is_pt, is_dir));
    }
// degenerate case
    if (a!=0 || p!=0) {
        if (a*s == p*d)
            return make_object(plane1);
        else
            return Object();
    }
    if (b!=0 || q!=0) {
        if (b*s == q*d)
            return make_object(plane1);
        else
            return Object();
    }
    if (c!=0 || r!=0) {
        if (c*s == r*d)
            return make_object(plane1);
        else
            return Object();
    }
    return make_object(plane1);
}

template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Plane_3 &plane1,
	     const typename CGAL_WRAP(K)::Plane_3 &plane2,
	     const typename CGAL_WRAP(K)::Plane_3 &plane3,
	     const K& k)
{
    typedef typename K::Line_3       Line_3;
    typedef typename K::Plane_3      Plane_3;

    // Intersection between plane1 and plane2 can either be
    // a line, a plane, or empty.
    Object o12 = CGALi::intersection(plane1, plane2, k);

    if (const Line_3 *l = object_cast<Line_3>(&o12))
        return CGALi::intersection(plane3, *l, k);

    if (const Plane_3 *pl = object_cast<Plane_3>(&o12))
        return CGALi::intersection(plane3, *pl, k);

    return Object();
}


template <class K>
bool
do_intersect(const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const typename CGAL_WRAP(K)::Line_3 &line,
	     const K&)
{
    typedef typename K::Point_3 Point_3;
    typedef typename K::Direction_3 Direction_3;
    typedef typename K::RT RT;
    const Point_3 &line_pt = line.point();
    const Direction_3 &line_dir = line.direction();

    RT den = plane.a()*line_dir.dx() + plane.b()*line_dir.dy()
           + plane.c()*line_dir.dz();
    if (den != 0)
        return true;
    RT num = plane.a()*line_pt.hx() + plane.b()*line_pt.hy()
           + plane.c()*line_pt.hz() + wmult_hw((K*)0, plane.d(), line_pt);
    if (num == 0) {
        // all line
        return true;
    } else {
        // no intersection
        return false;
    }
}

template <class K>
inline
bool
do_intersect(const typename CGAL_WRAP(K)::Line_3 &line, 
	     const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const K& k)
{
  return do_intersect(plane, line, k);
}


template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const typename CGAL_WRAP(K)::Ray_3 &ray, 
	     const K& k)
{
    typedef typename K::Point_3 Point_3;
    const Object line_intersection =
            intersection(plane, ray.supporting_line(), k);
    if (const Point_3 *isp = object_cast<Point_3>(&line_intersection)) {
        if (ray.collinear_has_on(*isp))
            return line_intersection;
        else
            return Object();
    }
    if (line_intersection.is_empty())
        return line_intersection;
    return make_object(ray);
}


template <class K>
inline
Object
intersection(const typename CGAL_WRAP(K)::Ray_3 &ray, 
	     const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const K& k)
{
  return intersection(plane, ray, k);
}



template <class K>
bool
do_intersect(const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const typename CGAL_WRAP(K)::Ray_3 &ray, 
	     const K& k)
{
    typedef typename K::Point_3 Point_3;
    const Object line_intersection =
            intersection(plane, ray.supporting_line(), k);
    if (line_intersection.is_empty())
        return false;
    if (const Point_3 *isp = object_cast<Point_3>(&line_intersection))
        return ray.collinear_has_on(*isp);
    return true;
}


template <class K>
inline
bool
do_intersect(const typename CGAL_WRAP(K)::Ray_3 &ray, 
	     const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const K& k)
{
  return do_intersect(plane, ray, k);
}


template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const typename CGAL_WRAP(K)::Segment_3 &seg, 
	     const K& k)
{
    typedef typename K::Point_3 Point_3;
    const Point_3 &source = seg.source();
    const Point_3 &target = seg.target();

    Oriented_side source_side = plane.oriented_side(source);
    Oriented_side target_side = plane.oriented_side(target);

    switch (source_side) {
    case ON_ORIENTED_BOUNDARY:
        if (target_side == ON_ORIENTED_BOUNDARY)
            return make_object(seg);
        else
            return make_object(source);
    case ON_POSITIVE_SIDE:
        switch (target_side) {
        case ON_ORIENTED_BOUNDARY:
            return make_object(target);
        case ON_POSITIVE_SIDE:
            return Object();
        case ON_NEGATIVE_SIDE:
            return intersection(plane, seg.supporting_line(), k);
        }
    case ON_NEGATIVE_SIDE:
        switch (target_side) {
        case ON_ORIENTED_BOUNDARY:
            return make_object(target);
        case ON_POSITIVE_SIDE:
            return intersection(plane, seg.supporting_line(), k);
        case ON_NEGATIVE_SIDE:
            return Object();
        }
    }
    CGAL_kernel_assertion_msg(false, "Supposedly unreachable code.");
    return Object();
}


template <class K>
inline
Object
intersection(const typename CGAL_WRAP(K)::Segment_3 &seg, 
	     const typename CGAL_WRAP(K)::Plane_3 &plane, 
	     const K& k)
{
  return intersection(plane, seg, k);
}


template <class K>
bool
do_intersect(const typename CGAL_WRAP(K)::Plane_3  &plane, 
	     const typename CGAL_WRAP(K)::Segment_3 &seg, 
	     const K&)
{
    typedef typename K::Point_3 Point_3;
    const Point_3 &source = seg.source();
    const Point_3 &target = seg.target();

    Oriented_side source_side = plane.oriented_side(source);
    Oriented_side target_side = plane.oriented_side(target);

    if ( source_side == target_side
       && target_side != ON_ORIENTED_BOUNDARY) {
        return false;
    }
    return true;
}


template <class K>
inline
bool
do_intersect(const typename CGAL_WRAP(K)::Segment_3 &seg, 
	     const typename CGAL_WRAP(K)::Plane_3  &plane, 
	     const K& k)
{
  return do_intersect(plane, seg, k);
}


template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Line_3 &line,
	     const Bbox_3 &box, 
	     const K&)
{
    typedef typename K::Point_3 Point_3;
    typedef typename K::Direction_3 Direction_3;
    const Point_3 &linepoint = line.point();
    const Direction_3 &linedir = line.direction();
    return intersection_bl(box,
        CGAL::to_double(linepoint.x()),
        CGAL::to_double(linepoint.y()),
        CGAL::to_double(linepoint.z()),
        CGAL::to_double(linedir.dx()),
        CGAL::to_double(linedir.dy()),
        CGAL::to_double(linedir.dz()),
        true, true
        );
}


template <class K>
inline
Object
intersection(const Bbox_3 &box, 
	     const typename CGAL_WRAP(K)::Line_3 &line, 
	     const K& k)
{
  return intersection(line, box, k);
}


template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Ray_3 &ray,
	     const Bbox_3 &box, 
	     const K&)
{
    typedef typename K::Point_3 Point_3;
    typedef typename K::Direction_3 Direction_3;
    const Point_3 &linepoint = ray.source();
    const Direction_3 &linedir = ray.direction();
    return intersection_bl(box,
        CGAL::to_double(linepoint.x()),
        CGAL::to_double(linepoint.y()),
        CGAL::to_double(linepoint.z()),
        CGAL::to_double(linedir.dx()),
        CGAL::to_double(linedir.dy()),
        CGAL::to_double(linedir.dz()),
        false, true
        );
}


template <class K>
inline
Object
intersection(const Bbox_3 &box, 
	     const typename CGAL_WRAP(K)::Ray_3 &ray, 
	     const K& k)
{
  return intersection(ray, box, k);
}



template <class K>
Object
intersection(const typename CGAL_WRAP(K)::Segment_3 &seg, 
	     const Bbox_3 &box, 
	     const K&)
{
    typedef typename K::Point_3 Point_3;
    typedef typename K::Vector_3 Vector_3;
    const Point_3 &linepoint = seg.source();
    const Vector_3 &diffvec = seg.target()-linepoint;
    return intersection_bl(box,
        CGAL::to_double(linepoint.x()),
        CGAL::to_double(linepoint.y()),
        CGAL::to_double(linepoint.z()),
        CGAL::to_double(diffvec.x()),
        CGAL::to_double(diffvec.y()),
        CGAL::to_double(diffvec.z()),
        false, false
        );
}


template <class K>
inline
Object
intersection(const Bbox_3 &box, 
	     const typename CGAL_WRAP(K)::Segment_3 &seg,