songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Packages/Distance_2/web/sqdistance_2.fw

2225 lines
64 KiB
Plaintext
Raw Blame History

@! $RCSfile$
@! $Revision$
@! $Date$
@! Author: Geert-Jan Giezeman (geert@cs.ruu.nl)
@i cgal_util.fwi
@A@<Squared distance computations@>
The routines in this chapter compute the square of the euclidean distance
between objects.
By computing the square of the distance between two objects instead of the
distance itself, often the computation of a square root is avoided.
This is not only important for efficiency reasons.
There may be number types (especially when an exact representation is used)
where taking the square root is not possible.
@B@<Utility routines@>
While computing square distances I need a number of routines that are not
supplied as basic facilities in CGAL, but that are pretty basic and not tied
to distance computations. Perhaps they should be added to more basic
components.
@C@<Low level routines that avoid divisions@>
In order to avoid division (which is important for homogeneous types), we
sometimes need to compute with different routines than the natural ones.
We prefix the natural name of such a routine with a @{w@} to indicate this.
For example, the dot product of two vectors @{u@} and @{v@} is computed as
@{dot(u,v) == (ux.vx+uy.vy)/(uw.vw)@}.
The routine @{wdot@} is computed as follows: @{wdot(u,v) == ux.vx+uy.vy@}.
We use a short notation here where @{ux@} stands for @{u.hx()@}, @{uw@} stands
for @{u.hw()@} and so on.
Often we only need the sign of a routine. Then multiplication by homogeneous
factors has no influence, as they are positive by definition.
Also, when we need to compare the dot product of two different vectors, we can
use multiplications instead of divisions. Instead of comparing @{dot(p,q)@}
with @{dot(u,v)@}, we can compare @{uw.vw.wdot(p,q)@} with
@{pw.qw.wdot(u,v)@}. The homogenising factors are brought to the other side of
the comparison operator which can be done, again, because they are positive.
The homogeneous factors that are multiplied differ from case to case. This
makes those routines pretty low level and care should be taken when they are
used.
In the case of cartesian coordinates, we know that the homogenising factors
are 1 by definition. This leads to ample possibilities for optimisation.
Therefore, most routines have a specialisation for the cartesian case that is
more efficient.
@D@<Multiplication with homogeneous factors (wmult)@>
As we have seen in the example where dot products are compared, we sometimes
need to multiply by homogenising factors. Of course, the multiplication
operator can be used for this. Still, we provide a special function, @{wmult@}
for this. The reason to do that is that now we can define a specialisation for
cartesian types. In that case, the multiplication is a noop and can be
optimised away.
There are versions of @{wmult@} taking a different number of parameters. This
depends on the number of factors that are to be multiplied.
Of the multiplying factors the first may an arbitrary number.
The other (1 or 2) factors must be homogenising factors.
This allows the optimisation in the cartesian case.
The routines take a, quite unnatural, first parameter.
This is the mechanism that allows to decide between the cartesian or general
case.
@O@<../include/CGAL/wmult.h@>==@{@-
@<cgal_heading@>@(@-
include/CGAL/wmult.h@,@-
sqdistance_2.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_WMULT_H
#define CGAL_WMULT_H
CGAL_BEGIN_NAMESPACE
#if defined CGAL_HOMOGENEOUS_H
template <class RT>
inline RT wmult(Homogeneous<RT>*, const RT &a, const RT & w)
{
return a*w;
}
template <class RT>
inline RT wmult(Homogeneous<RT>*,
const RT &a, const RT & w1, const RT & w2)
{
return a*w1*w2;
}
template <class RT>
inline RT wmult(Homogeneous<RT>*, const RT &a,
const RT & w1, const RT & w2, const RT & w3)
{
return a*w1*w2*w3;
}
template <class RT>
inline RT wmult(Homogeneous<RT>*, const RT &a,
const RT & w1, const RT & w2, const RT & w3, const RT & w4)
{
return a*w1*w2*w3*w4;
}
#endif // CGAL_HOMOGENEOUS_H
#if defined CGAL_SIMPLE_HOMOGENEOUS_H
template <class RT>
inline RT wmult(Simple_homogeneous<RT>*, const RT &a, const RT & w)
{
return a*w;
}
template <class RT>
inline RT wmult(Simple_homogeneous<RT>*,
const RT &a, const RT & w1, const RT & w2)
{
return a*w1*w2;
}
template <class RT>
inline RT wmult(Simple_homogeneous<RT>*, const RT &a,
const RT & w1, const RT & w2, const RT & w3)
{
return a*w1*w2*w3;
}
template <class RT>
inline RT wmult(Simple_homogeneous<RT>*, const RT &a,
const RT & w1, const RT & w2, const RT & w3, const RT & w4)
{
return a*w1*w2*w3*w4;
}
#endif // CGAL_SIMPLE_HOMOGENEOUS_H
#if defined CGAL_CARTESIAN_H
template <class RT>
inline RT wmult(Cartesian<RT> *, const RT &a, const RT & )
{
return a;
}
template <class RT>
inline RT wmult(Cartesian<RT> *,
const RT &a, const RT & , const RT & )
{
return a;
}
template <class RT>
inline RT wmult(Cartesian<RT> *, const RT &a,
const RT & , const RT & , const RT & )
{
return a;
}
template <class RT>
inline RT wmult(Cartesian<RT> *, const RT &a,
const RT & , const RT & , const RT & , const RT & )
{
return a;
}
#endif // CGAL_CARTESIAN_H
#if defined CGAL_SIMPLE_CARTESIAN_H
template <class RT>
inline RT wmult(Simple_cartesian<RT> *, const RT &a, const RT & )
{
return a;
}
template <class RT>
inline RT wmult(Simple_cartesian<RT> *,
const RT &a, const RT & , const RT & )
{
return a;
}
template <class RT>
inline RT wmult(Simple_cartesian<RT> *, const RT &a,
const RT & , const RT & , const RT & )
{
return a;
}
template <class RT>
inline RT wmult(Simple_cartesian<RT> *, const RT &a,
const RT & , const RT & , const RT & , const RT & )
{
return a;
}
#endif // CGAL_SIMPLE_CARTESIAN_H
CGAL_END_NAMESPACE
#endif // CGAL_WMULT_H
@}
@D@<Testing for null vector@>
A routine to test whether a vector is a null vector.
This code does not belong here.
@$@<2D squared distance general utilities@>+=@{
template <class R>
bool is_null(const Vector_2<R> &v)
{
typedef typename R::RT RT;
return v.hx()==RT(0) && v.hy()==RT(0);
}
@}
@D@<The dot product@>
The dot product is useful for projecting a point on a line, and for deciding
whether an angle is acute, straight or obtuse.
Apart from the @{wdot@} routine there are also routines that use @{wdot@} to
decide on the acuteness of an angle.
Of all the routines there are different flavours, taking points and
vectors as arguments respectively. The points @{p,q@} and @{r@} stand
for the vectors @{p-q@} and @{r-q@}. The homogenising factors that are used
for multiplication are different though.
For vectors @{u@} and @{v@} we have @{wdot(u,v) == wu.wv.dot(u,v)@}
For points @{p,q@} and @{r@} we have @{wdot(p,q,r) == wp.wq.wq.wr.dot(p-q,r-q)@}
Note the double occurrence of the factor@{wq@} in the latter definition.
@$@<2D squared distance general utilities@>+=@{
template <class R>
typename R::RT
wdot(const Vector_2<R> &u, const Vector_2<R> &v)
{
return (u.hx()*v.hx() + u.hy()*v.hy());
}
template <class R>
typename R::RT wdot(const Point_2< R > &p,
const Point_2< R > &q,
const Point_2< R > &r)
{
R* pR = 0;
return (wmult(pR, p.hx(),q.hw()) - wmult(pR, q.hx(),p.hw()))
* (wmult(pR, r.hx(),q.hw()) - wmult(pR, q.hx(),r.hw()))
+ (wmult(pR, p.hy(),q.hw()) - wmult(pR, q.hy(),p.hw()))
* (wmult(pR, r.hy(),q.hw()) - wmult(pR, q.hy(),r.hw()));
}
@}
@D@<The cross product@>
The cross product of two vectors can be used for various purposes.
The sign tells whether the two vectors make a left or right turn or are
collinear (have the same or opposite direction). The absolute value gives the
area of the parallelogram spanned by the two vectors. This can also be used to
compute the (signed) distance from a point to a line.
Like with the dot product, we have versions for two vectors and for three
points.
We define @{wcross(u,v) == wu.wv.cross(u,v)@} and
@{wcross(p,q,r) == @}@{wp.wq.wr.cross(q-p,r-q)@}.
Different from the dot case, the factor @{wq@} occurs only once in the
definition of the wcross product for points.
@$@<2D squared distance general utilities@>+=@{
template <class R>
typename R::RT
wcross(const Vector_2<R> &u,
const Vector_2<R> &v)
{
return (typename R::RT)(u.hx()*v.hy() - u.hy()*v.hx());
}
#if defined CGAL_HOMOGENEOUS_H
template <class RT>
inline
RT wcross_impl(const Homogeneous<RT>*,const Point_2< Homogeneous<RT> > &p,
const Point_2< Homogeneous<RT> > &q,
const Point_2< Homogeneous<RT> > &r)
{
return p.hx() * (q.hy()*r.hw() - q.hw()*r.hy() )
+ p.hy() * (q.hw()*r.hx() - q.hx()*r.hw() )
+ p.hw() * (q.hx()*r.hy() - q.hy()*r.hx() );
}
#endif // CGAL_HOMOGENEOUS_H
#if defined CGAL_SIMPLE_HOMOGENEOUS_H
template <class RT>
inline
RT wcross_impl(const Simple_homogeneous<RT>*,
const Point_2< Simple_homogeneous<RT> > &p,
const Point_2< Simple_homogeneous<RT> > &q,
const Point_2< Simple_homogeneous<RT> > &r)
{
return det3x3_by_formula(
p.hx(), q.hx(), r.hx(),
p.hy(), q.hy(), r.hy(),
p.hw(), q.hw(), r.hw());
}
#endif // CGAL_SIMPLE_HOMOGENEOUS_H
#if defined CGAL_CARTESIAN_H
template <class FT>
inline
FT wcross_impl(const Cartesian<FT> *, const Point_2< Cartesian<FT> > &p,
const Point_2< Cartesian<FT> > &q,
const Point_2< Cartesian<FT> > &r)
{
return (q.x()-p.x())*(r.y()-q.y()) - (q.y()-p.y())*(r.x()-q.x());
}
#endif // CGAL_CARTESIAN_H
#if defined CGAL_SIMPLE_CARTESIAN_H
template <class FT>
inline
FT wcross_impl(const Simple_cartesian<FT> *,
const Point_2< Simple_cartesian<FT> > &p,
const Point_2< Simple_cartesian<FT> > &q,
const Point_2< Simple_cartesian<FT> > &r)
{
return (q.x()-p.x())*(r.y()-q.y()) - (q.y()-p.y())*(r.x()-q.x());
}
#endif // CGAL_SIMPLE_CARTESIAN_H
template <class R>
typename R::RT wcross(const Point_2< R > &p,
const Point_2< R > &q,
const Point_2< R > &r)
{
return wcross_impl(static_cast<R*>(0), p, q, r);
}
@}
@C@<Routines that characterise angles@>
Often one needs to know whether an angle is sharp or obtuse or whether it is a
left turn or a right turn. An angle can be characterised by three points
@{p,q@} and @{r@} (the angle formed at @{q@} when going from @{p@} to @{r@}
via @{q@}), or by two vectors @{u@} and @{v@} (the angle formed when their
tails are placed in the origin).
Then we can use the routines @{wdot@} and @{wcross@} for this purpose.
@D@<Sharp or obtuse@>
Here are routines to decide whether an angle is sharp, straight or obtuse.
There should also be an enum with three values (@{ACUTE@},
@{STRAIGHT@},@{OBTUSE@}) and a routine that returns such a value.
@$@<2D squared distance general utilities@>+=@{
template <class R>
inline bool is_acute_angle(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wdot(u, v)) > RT(0) ;
}
template <class R>
inline bool is_straight_angle(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wdot(u, v)) == RT(0) ;
}
template <class R>
inline bool is_obtuse_angle(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wdot(u, v)) < RT(0) ;
}
template <class R>
inline bool is_acute_angle(const Point_2<R> &p,
const Point_2<R> &q, const Point_2<R> &r)
{
typedef typename R::RT RT;
return RT(wdot(p, q, r)) > RT(0) ;
}
template <class R>
inline bool is_straight_angle(const Point_2<R> &p,
const Point_2<R> &q, const Point_2<R> &r)
{
typedef typename R::RT RT;
return RT(wdot(p, q, r)) == RT(0) ;
}
template <class R>
inline bool is_obtuse_angle(const Point_2<R> &p,
const Point_2<R> &q, const Point_2<R> &r)
{
typedef typename R::RT RT;
return RT(wdot(p, q, r)) < RT(0) ;
}
@}
@D@<Clockwise or counterclockwise@>
The following routines decide whether the smallest angle between two vectors
(when their tails are in the same point) is clockwise or counterclockwise.
If both angles are equal (zero or pi), the vectors are collinear.
If the angle is characterised by three points @{p@}, @{q@} and @{r@}, it is
more natural to think of making a left or right turn at @{q@} when going from
@{p@} to @{r@}. Those routines are implemented elsewhere, though.
@$@<2D squared distance general utilities@>+=@{
template <class R>
Orientation orientation(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
RT wcr = wcross(u,v);
return (wcr > RT(0)) ? COUNTERCLOCKWISE :
(wcr < RT(0)) ? CLOCKWISE
: COLLINEAR;
}
template <class R>
inline bool counterclockwise(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wcross(u,v)) > RT(0);
}
template <class R>
inline bool leftturn(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wcross(u,v)) > RT(0);
}
template <class R>
inline bool clockwise(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wcross(u,v)) < RT(0);
}
template <class R>
inline bool rightturn(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wcross(u,v)) < RT(0);
}
template <class R>
inline bool collinear(const Vector_2<R> &u,
const Vector_2<R> &v)
{
typedef typename R::RT RT;
return RT(wcross(u,v)) == RT(0);
}
/*
the ordertype, rightturn, leftturn and collinear routines for points are
defined elsewhere.
*/
@}
@B@<Point Point distance@>
The squared distance between two points is computed by taking the dot product
of the difference vector.
@$@<2D squared distance header declarations 1@>+=@{@-
template <class R>
inline typename R::FT
squared_distance(
const Point_2<R> & pt1,
const Point_2<R> & pt2)
{
Vector_2<R> vec(pt1-pt2);
return (typename R::FT)(vec*vec);
}
@}
@B@<Point Line distance@>
@C@<Declarations@>
@$@<2D squared distance header declarations 1@>+=@{@-
template <class R>
class Squared_distance_to_line {
typename R::RT a, b, c, sqnorm;
public:
Squared_distance_to_line(typename R::Line_2 const &line)
: a(line.a()), b(line.b()), c(line.c())
{
sqnorm = a*a+b*b;
}
typename R::FT operator()(typename R::Point_2 const &pt) const
{
typedef typename R::RT RT;
RT w = pt.hw();
RT n = a*pt.hx() + b*pt.hy() + wmult((R*)0, c, w);
RT d = wmult((R*)0, sqnorm, w, w);
return R::make_FT(n*n, d);
}
};
@<Point Line distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Line_2<R> & line,
const Point_2<R> & pt)
{
return squared_distance(pt, line);
}
@}
@C@<Implementation@>
We compute a vector perpendicular to the line(@{normal@}).
The dot product of a vector from the point to an arbitrary point on the line
with the @{normal@} is the (signed) distance from the point to the line,
multiplied with the length of @{normal@}. Taking the square of this value and
dividing by the square of the length of the @{normal@} gives the square
distance between the line and the point.
@$@<Point Line distance implementation@>==@{
template <class R>
typename R::FT
squared_distance(
const Point_2<R> &pt,
const Line_2<R> &line)
{
typedef typename R::RT RT;
RT a = line.a();
RT b = line.b();
RT w = pt.hw();
RT n = a*pt.hx() + b*pt.hy() + wmult((R*)0, line.c(), w);
RT d = wmult((R*)0, a*a+b*b, w, w);
return R::make_FT(n*n, d);
}
@}
@B@<Point Ray distance@>
@C@<Declarations@>
@$@<2D squared distance header declarations 1@>+=@{@-
template <class R>
class Squared_distance_to_ray {
typename R::Vector_2 ray_dir;
typename R::Point_2 ray_source;
Squared_distance_to_line<R> supline_dist;
public:
Squared_distance_to_ray(typename R::Ray_2 const &ray)
: ray_dir(ray.direction().vector()),
ray_source(ray.source()),
supline_dist(ray.supporting_line())
{ }
typename R::FT operator()(typename R::Point_2 const &pt) const
{
Vector_2<R> diff = pt-ray_source;
if (!is_acute_angle(ray_dir,diff) )
return (typename R::FT)(diff*diff);
return supline_dist(pt);
}
};
@<Point Ray distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Ray_2<R> & ray,
const Point_2<R> & pt)
{
return squared_distance(pt, ray);
}
@}
@C@<Implementation@>
We test the angle between the direction vector of the ray and the vector from
the start point of the ray to the point @{pt@}. If this angle is acute,
a line perpendicular to the ray through @{pt@} passes through the ray, so we
the distance from the point to the ray is the same as the distance from the
point to the supporting line of the ray. Otherwise we compute the squared
distance from @{pt@} to the start point of the ray, as this is the closest
point.
@$@<Point Ray distance implementation@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Point_2<R> &pt,
const Ray_2<R> &ray)
{
Vector_2<R> diff = pt-ray.source();
const Vector_2<R> &dir = ray.direction().vector();
if (!is_acute_angle(dir,diff) )
return (typename R::FT)(diff*diff);
return squared_distance(pt, ray.supporting_line());
}
@}
@B@<Point Ray distance discerners@>
@$@<2D squared distance header declarations 1@>+=@{
@<Point Ray distance discerners implementation@>
@}
@C@<Implementation@>
@$@<Point Ray distance discerners implementation@>==@{
template <class R>
extern void
distance_index(
int &ind,
const Point_2<R> &pt,
const Ray_2<R> &ray)
{
if (!is_acute_angle(ray.direction().vector(),pt-ray.source())) {
ind = 0;
return;
}
ind = -1;
}
template <class R>
typename R::FT
squared_distance_indexed(const Point_2<R> &pt,
const Ray_2<R> &ray, int ind)
{
if (ind == 0)
return squared_distance(pt, ray.source());
return squared_distance(pt, ray.supporting_line());
}
@}
@B@<Point Segment distance@>
@$@<2D squared distance header declarations 1@>+=@{@-
template <class R>
class Squared_distance_to_segment {
typename R::Point_2 seg_source, seg_target;
Squared_distance_to_line<R> supline_dist;
typename R::Vector_2 segvec;
typename R::RT e;
public:
Squared_distance_to_segment(typename R::Segment_2 const &seg)
: seg_source(seg.source()), seg_target(seg.target()),
supline_dist(seg.supporting_line())
{
segvec = seg_target-seg_source;
e = wdot(segvec,segvec);
}
typename R::FT operator()(typename R::Point_2 const &pt) const
{
typedef typename R::RT RT;
// assert that the segment is valid (non zero length).
Vector_2<R> diff = pt-seg_source;
RT d = wdot(diff,segvec);
if (d <= (RT)0)
return (typename R::FT)(diff*diff);
if (wmult((R*)0 ,d, segvec.hw()) > wmult((R*)0, e, diff.hw()))
return squared_distance(pt, seg_target);
return supline_dist(pt);
}
};
@<Point Segment distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Segment_2<R> & seg,
const Point_2<R> & pt)
{
return squared_distance(pt, seg);
}
@}
@C@<Implementation@>
To compute the distance of a point to a segment, first we divide the plane in
three parts. We draw two lines perpendicular to the segment, one through the
start point, the other through the end point of the segment. If the query
point lies between those two lines, we need to compute the distance of the
query point to the supporting line of the segment. Otherwise we need to
compute the distance from the point to either the endpoint or the startpoint.
To decide which case applies, we compute the dot poduct of the vector from
segment start to segment end and the vector from segment start to the point.
If this product is negative, the distance from the point to the start of the
segment is taken.
@$@<Point Segment distance implementation@>==@{
template <class R>
typename R::FT
squared_distance(
const Point_2<R> &pt,
const Segment_2<R> &seg)
{
typedef typename R::RT RT;
// assert that the segment is valid (non zero length).
Vector_2<R> diff = pt-seg.source();
Vector_2<R> segvec = seg.target()-seg.source();
RT d = wdot(diff,segvec);
if (d <= (RT)0)
return (typename R::FT)(diff*diff);
RT e = wdot(segvec,segvec);
if (wmult((R*)0 ,d, segvec.hw()) > wmult((R*)0, e, diff.hw()))
return squared_distance(pt, seg.target());
return squared_distance(pt, seg.supporting_line());
}
@}
@B@<Point Segment distance discerners@>
@$@<2D squared distance header declarations 1@>+=@{
@<Point Segment distance discerners implementation@>
@}
@C@<Implementation@>
@$@<Point Segment distance discerners implementation@>==@{
template <class R>
extern void
distance_index(
int &ind,
const Point_2<R> &pt,
const Segment_2<R> &seg)
{
if (!is_acute_angle(seg.target(),seg.source(),pt)) {
ind = 0;
return;
}
if (!is_acute_angle(seg.source(),seg.target(),pt)) {
ind = 1;
return;
}
ind = -1;
}
template <class R>
typename R::FT
squared_distance_indexed(const Point_2<R> &pt,
const Segment_2<R> &seg, int ind)
{
if (ind == 0)
return squared_distance(pt, seg.source());
if (ind == 1)
return squared_distance(pt, seg.target());
return squared_distance(pt, seg.supporting_line());
}
@}
@B@<Segment Segment distance@>
@C@<Declarations@>
@$@<2D squared distance header declarations 1@>+=@{@-
@<Segment Segment distance implementation@>
@}
@C@<Implementation@>
If two segments do not cross, the distance of the two segments can be found by
comparing the distance of the endpoints to the other segment. The minimum
of those four distances is the distance between the two segments. This is the
fact that is exploited in the implementation.
@D@<Distance computation@>
To compute the squared distance of two segments, the following steps are
taken.
First we check if one of the segments is degenerate (coinciding start and
end point). If so, we solve the problem for a point and a segment.
Then we decide if the one segment crosses the supporting line of the other
segment. The boolean @{crossing1@} tells whether @{seg1@} crosses the
supporting line of @{seg2@}, and vice versa for @{crossing2@}.
To decide whether they cross, we look at the sign of the cross product. If an
endpoint lies on the supporting line, we consider this as crossing. The value
of the cross product is stored in a variable because it is needed later on.
This is also the reason why we do not use the ordertype routine to decide on
which side the points lie. That would mean double work.
A special case arises if both segments lie on the same supporting line.
This is checked separately.
@$@<2D segment segment squared distance main routine@>+=@{
template <class R>
typename R::FT
squared_distance(
const Segment_2<R> &seg1,
const Segment_2<R> &seg2)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
bool crossing1, crossing2;
RT c1s, c1e, c2s, c2e;
if (seg1.source() == seg1.target())
return squared_distance(seg1.source(), seg2);
if (seg2.source() == seg2.target())
return squared_distance(seg2.source(), seg1);
c1s = wcross(seg2.source(), seg2.target(), seg1.source());
c1e = wcross(seg2.source(), seg2.target(), seg1.target());
c2s = wcross(seg1.source(), seg1.target(), seg2.source());
c2e = wcross(seg1.source(), seg1.target(), seg2.target());
if (c1s < RT(0)) {
crossing1 = (c1e >= RT(0));
} else {
if (c1e <= RT(0)) {
if (c1s == RT(0) && c1e == RT(0))
return squared_distance_parallel(seg1, seg2);
crossing1 = true;
} else {
crossing1 = (c1s == RT(0));
}
}
if (c2s < RT(0)) {
crossing2 = (c2e >= RT(0));
} else {
if (c2e <= RT(0)) {
if (c2s == RT(0) && c2e == RT(0))
return squared_distance_parallel(seg1, seg2);
crossing2 = true;
} else {
crossing2 = (c2s == RT(0));
}
}
@}
Depending on the values of the crossing variables we now proceed. If both
segments cross the supporting line of the other segment, the segments
themselves must cross, so they have zero distance.
If exactly one of the segments does not cross, the distance is the
distance of one of the endpoints of this (non crossing) segment to the
other (crossing) segment. Which endpoint? The one that is closest to the
supporting lines of the other segment. This latter property can easily be
computed with the value of the cross product that was computed earlier.
The function @{_distance_measure_sub@} makes this comparison (it
subtracts the two absolute values of the cross product, first corrected by
homogenising factors).
Note that the computation of the endpoint closest to a line may be
susceptible to roundoff error. The result can be incorrect in that
case. This must still be solved, for instance by computing the squared
distance from both endpoints to the other segment and taking the minimum.
@$@<unused implementation@>@Z==@{
template <class R, class RT>
typename R::FT
_squared_distance_of_endpoint_exact(RT start_dm, RT end_dm,
const Point_2<R> &start, const Point_2<R> &end,
const Segment<R> &seg)
{
typedef typename R::RT RT;
RT dm;
dm = _distance_measure_sub(start_dm, end_dm, start, end);
if (dm < RT(0)) {
return squared_distance(start, seg);
} else {
return squared_distance(end, seg);
}
}
template <class R, class RT>
typename R::FT
_squared_distance_of_endpoint_inexact<R, RT>(
const Point_2<R> &start, const Point_2<R> &end,
const Segment<R> &seg)
{
typedef typename R::RT RT;
RT d1, d2;
d1 = squared_distance(start, seg);
d2 = squared_distance(end, seg);
return (d1 < d2) ? (typename R::FT)(d1) : (typename_R::FT)(d2);
}
template <class R, class RT>
inline typename R::FT
_squared_distance_of_endpoint<false, R, RT>(RT start_dm, RT end_dm,
const Point_2<R> &start, const Point_2<R> &end,
const Segment<R> &seg)
{
if (numeric_limits<R>::is_exact)
return _squared_distance_of_endpoint_exact(
start_dm, end_dm, start, end, seg);
else
return _squared_distance_of_endpoint_inexact(start, end, seg);
}
@}
@$@<2D segment segment squared distance main routine@>+=@{
if (crossing1) {
if (crossing2)
return (FT)0;
RT dm;
dm = _distance_measure_sub(c2s,c2e, seg2.source(), seg2.target());
if (dm < RT(0)) {
return squared_distance(seg2.source(), seg1);
} else {
if (dm > RT(0)) {
return squared_distance(seg2.target(), seg1);
} else {
// parallel, should not happen (no crossing)
return squared_distance_parallel(seg1, seg2);
}
}
} else {
if (crossing2) {
RT dm;
dm =
_distance_measure_sub(c1s, c1e,seg1.source(),seg1.target());
if (dm < RT(0)) {
return squared_distance(seg1.source(), seg2);
} else {
if (dm > RT(0)) {
return squared_distance(seg1.target(), seg2);
} else {
// parallel, should not happen (no crossing)
return squared_distance_parallel(seg1, seg2);
}
}
} else {
@}
Finally, neither of the segments may cross. This may be because the two
segments are parallel, which we deal with as a separate case.
Otherwise we have to compare the distance of two endpoints to the other
segment and take the minimum.
@$@<2D segment segment squared distance main routine@>+=@{
FT min1, min2;
RT dm = _distance_measure_sub(
c1s, c1e, seg1.source(), seg1.target());
if (dm == RT(0))
return squared_distance_parallel(seg1, seg2);
min1 = (dm < RT(0)) ?
squared_distance(seg1.source(), seg2):
squared_distance(seg1.target(), seg2);
dm = _distance_measure_sub(
c2s, c2e, seg2.source(), seg2.target());
if (dm == RT(0)) // should not happen.
return squared_distance_parallel(seg1, seg2);
min2 = (dm < RT(0)) ?
squared_distance(seg2.source(), seg1):
squared_distance(seg2.target(), seg1);
return (min1 < min2) ? min1 : min2;
}
}
}
@}
We need a small helper function to compare the distance from a point to the
intersection point of the supporting lines of the two segments. With the help
of the cross product (that is known), this can be computed cheaply.
@$@<2D segment segment squared distance helper routine@>==@{
template <class RT, class R>
RT _distance_measure_sub(RT startwcross, RT endwcross,
const Point_2<R> &start, const Point_2<R> &end
)
{
return CGAL_NTS abs(wmult((R*)0, startwcross, end.hw())) -
CGAL_NTS abs(wmult((R*)0, endwcross, start.hw()));
}
@}
We still need to deal with the case where both segments are parallel. The
segments may not be degenerate (this was tested before), but they may have
opposite directions. This is checked and stored in the boolean variable
@{same_direction@}.
The distance between two parallel segments is either the
distance between two endpoints or the distance between an arbitrary point on
one segment to the supporting line of the other segment. In order to decide
which case applies, we can project the second segment perpendicularly on the
supporting line of the first segment. If the resulting segments overlap, the
second case apllies, otherwise the former case. The projections are just a way
of explaining, they are not actually done. The sign of the dot product is used
to make the necessary decisions.
@$@<2D segment segment squared distance parallel case@>+=@{
template <class R>
typename R::FT
squared_distance_parallel(
const Segment_2<R> &seg1,
const Segment_2<R> &seg2)
{
bool same_direction;
const Vector_2<R> &dir1 = seg1.direction().vector();
const Vector_2<R> &dir2 = seg2.direction().vector();
if (CGAL_NTS abs(dir1.hx()) > CGAL_NTS abs(dir1.hy())) {
same_direction = (CGAL_NTS sign(dir1.hx()) == CGAL_NTS sign(dir2.hx()));
} else {
same_direction = (CGAL_NTS sign(dir1.hy()) == CGAL_NTS sign(dir2.hy()));
}
if (same_direction) {
if (!is_acute_angle(seg1.source(), seg1.target(), seg2.source()))
return squared_distance(seg1.target(), seg2.source());
if (!is_acute_angle(seg1.target(), seg1.source(), seg2.target()))
return squared_distance(seg1.source(), seg2.target());
} else {
if (!is_acute_angle(seg1.source(), seg1.target(), seg2.target()))
return squared_distance(seg1.target(), seg2.target());
if (!is_acute_angle(seg1.target(), seg1.source(), seg2.source()))
return squared_distance(seg1.source(), seg2.source());
}
return squared_distance(seg2.source(), seg1.supporting_line());
}
@}
@$@<Segment Segment distance implementation@>==@{
@<2D segment segment squared distance parallel case@>
@<2D segment segment squared distance helper routine@>
@<2D segment segment squared distance main routine@>
@}
@B@<Segment Ray distance@>
@C@<Declarations@>
@$@<2D squared distance header declarations 1@>+=@{@-
@<Segment Ray distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Ray_2<R> & ray,
const Segment_2<R> & seg)
{
return squared_distance(seg, ray);
}
@}
@C@<Implementation@>
The computation of the squared distance between a segment and a ray goes along
the same lines as the computation of the square distance of two segments.
We give only short comments here.
First we check if the segment is degenerate (coinciding start and
end point). If so, we solve the problem for a point and a ray.
Then we decide if the one segment crosses the supporting line of the other
segment. The boolean @{crossing1@} tells whether the segment crosses the
supporting line of the ray, @{crossing2@} tells whether the ray crosses the
supporting line of the segment.
For the ray we look at the direction and on which side its starting point
lies.
@$@<2D ray segment squared distance computation@>+=@{
template <class R>
typename R::FT
squared_distance(
const Segment_2<R> &seg,
const Ray_2<R> &ray)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const Vector_2<R> &raydir = ray.direction().vector();
Vector_2<R> startvec(seg.source()-ray.source());
Vector_2<R> endvec(seg.target()-ray.source());
bool crossing1, crossing2;
RT c1s, c1e;
Orientation ray_s_side;
if (seg.source() == seg.target())
return squared_distance(seg.source(), ray);
c1s = wcross(raydir, startvec);
c1e = wcross(raydir, endvec);
if (c1s < RT(0)) {
crossing1 = (c1e >= RT(0));
} else {
if (c1e <= RT(0)) {
if (c1s == RT(0) && c1e == RT(0))
return squared_distance_parallel(seg, ray);
crossing1 = true;
} else {
crossing1 = (c1s == RT(0));
}
}
ray_s_side = orientation(seg.source(), seg.target(), ray.source());
switch (ray_s_side) {
case LEFTTURN:
crossing2 = rightturn(seg.target()-seg.source(), raydir);
break;
case RIGHTTURN:
crossing2 = leftturn(seg.target()-seg.source(), raydir);
break;
case COLLINEAR:
crossing2 = true;
break;
}
@}
Depending on the values of the crossing variables we now proceed. If both
segments cross the supporting line of the other segment, the segments
themselves must cross, so they have zero distance.
If exactly one of the segments does not cross, the distance is the
distance of one of the endpoints of this (non crossing) segment to the
other (crossing) segment. Which endpoint? The one that is closest to the
intersection point of the supporting lines of the two segments. This
latter property can easily be computed with the value of the cross product
that was computed earlier. The function @{_distance_measure_sub@}
makes this comparison (it subtracts the two absolute values of the cross
product, first corrected by homogenising factors).
@$@<2D ray segment squared distance computation@>+=@{
if (crossing1) {
if (crossing2)
return FT(0);
return squared_distance(ray.source(), seg);
} else {
if (crossing2) {
RT dm;
dm = _distance_measure_sub(c1s, c1e, startvec, endvec);
if (dm < RT(0)) {
return squared_distance(seg.source(), ray);
} else {
if (dm > RT(0)) {
return squared_distance(seg.target(), ray);
} else {
// parallel, should not happen (no crossing)
return squared_distance_parallel(seg, ray);
}
}
} else {
@}
Finally, neither of the segments may cross. This may be because the two
segments are parallel, which we deal with as a separate case.
Otherwise we have to compare the distance of two endpoints to the other
segment and take the minimum.
@$@<2D ray segment squared distance computation@>+=@{
FT min1, min2;
RT dm;
dm = _distance_measure_sub(c1s, c1e, startvec, endvec);
if (dm == RT(0))
return squared_distance_parallel(seg, ray);
min1 = (dm < RT(0))
? squared_distance(seg.source(), ray)
: squared_distance(seg.target(), ray);
min2 = squared_distance(ray.source(), seg);
return (min1 < min2) ? min1 : min2;
}
}
}
@}
@$@<2D ray segment squared distance parallel case@>==@{
template <class R>
typename R::FT
squared_distance_parallel(
const Segment_2<R> &seg,
const Ray_2<R> &ray)
{
bool same_direction;
const Vector_2<R> &dir1 = seg.direction().vector();
const Vector_2<R> &dir2 = ray.direction().vector();
if (CGAL_NTS abs(dir1.hx()) > CGAL_NTS abs(dir1.hy())) {
same_direction = (CGAL_NTS sign(dir1.hx()) == CGAL_NTS sign(dir2.hx()));
} else {
same_direction = (CGAL_NTS sign(dir1.hy()) == CGAL_NTS sign(dir2.hy()));
}
if (same_direction) {
if (!is_acute_angle(seg.source(), seg.target(), ray.source()))
return squared_distance(seg.target(), ray.source());
} else {
if (!is_acute_angle(seg.target(), seg.source(), ray.source()))
return squared_distance(seg.source(), ray.source());
}
return squared_distance(ray.source(), seg.supporting_line());
}
@}
@$@<2D ray segment squared distance helper routine@>==@{
template <class RT, class R>
RT _distance_measure_sub(RT startwcross, RT endwcross,
const Vector_2<R> &start, const Vector_2<R> &end
)
{
return CGAL_NTS abs(wmult((R*)0, startwcross, end.hw())) -
CGAL_NTS abs(wmult((R*)0, endwcross, start.hw()));
}
@}
@$@<Segment Ray distance implementation@>==@{
@<2D ray segment squared distance helper routine@>
@<2D ray segment squared distance parallel case@>
@<2D ray segment squared distance computation@>
@}
@B@<Segment Line distance@>
@C@<Declarations@>
@$@<2D squared distance header declarations 1@>+=@{@-
@<Segment Line distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Line_2<R> & line,
const Segment_2<R> & seg)
{
return squared_distance(seg, line);
}
@}
@C@<Implementation@>
The computation of the squared distance between a segment and a line goes along
the same lines as the computation of the square distance of two segments.
First we check if the segment is degenerate (coinciding start and
end point). If so, we solve the problem for a point and a line.
Then we decide if the segment crosses the line.
The boolean @{crossing1@} tells whether the segment crosses the line.
@$@<2D line segment squared distance computation@>+=@{
template <class R>
typename R::FT
squared_distance(
const Segment_2<R> &seg,
const Line_2<R> &line)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const Vector_2<R> &linedir = line.direction().vector();
const Point_2<R> &linepoint = line.point();
Vector_2<R> startvec(seg.source()-linepoint);
Vector_2<R> endvec(seg.target()-linepoint);
bool crossing1;
RT c1s, c1e;
if (seg.source() == seg.target())
return squared_distance(seg.source(), line);
c1s = wcross(linedir, startvec);
c1e = wcross(linedir, endvec);
if (c1s < RT(0)) {
crossing1 = (c1e >= RT(0));
} else {
if (c1e <= RT(0)) {
crossing1 = true;
} else {
crossing1 = (c1s == RT(0));
}
}
@}
If the segment crosses the line, they have zero distance. Otherwise, the
distance is the distance of one of the endpoints of the segment to the
line. Which endpoint? The one that is closest to the supporting line of
the other segment. This latter property can easily be computed with the
value of the cross product that was computed earlier. The function
@{_distance_measure_sub@} makes this comparison (it subtracts the two
absolute values of the cross product, first corrected by homogenising
factors).
@$@<2D line segment squared distance computation@>+=@{
if (crossing1) {
return (FT)0;
} else {
RT dm;
dm = _distance_measure_sub(c1s, c1e, startvec, endvec);
if (dm <= RT(0)) {
return _sqd_to_line(startvec, c1s, linedir);
} else {
return _sqd_to_line(endvec, c1e, linedir);
}
}
}
@}
The following routine returns the square distance from a point to a line.
It is an optimised routine, that makes use of the intermediate results that
are available. It is based on the observation that given @{diff@}, the
difference vector between the point and a point on the line, and @{dir@}, the
(not necessarily normalized) direction vector of the line, the squared
distance can be computed as: @{cross(dir,diff)^2 / dot(dir,dir)@}.
We can write this in terms of wcross and wdot operations and multiplications
with hw factors. The hw factors of @{dir@} cancel.
@$@<2D line segment squared distance helper routine@>==@{
template <class RT, class R>
typename R::FT
_sqd_to_line(const Vector_2<R> &diff,
const RT & wcross, const Vector_2<R> &dir )
{
typedef typename R::FT FT;
RT numerator = wcross*wcross;
RT denominator = wmult((R*)0, RT(wdot(dir,dir)),
diff.hw(), diff.hw());
FT result = R::make_FT(numerator, denominator);
return result;
}
@}
@$@<Segment Line distance implementation@>==@{
@<2D line segment squared distance helper routine@>
@<2D line segment squared distance computation@>
@}
@B@<Ray Ray distance@>
@$@<2D squared distance header declarations 1@>+=@{@-
@<Ray Ray distance@>
@}
@C@<Implementation@>
The computation of the squared distance two rays goes along the same lines as
the computation of the square distance between a segment and a ray.
First we decide if the one ray crosses the supporting line of the other ray.
The boolean @{crossing1@} tells whether @{ray1@} crosses the supporting line
of the @{ray2@}, @{crossing2@} tells whether @{ray2@} crosses the
supporting line of @{ray1@}.
We look at the direction of the ray and on which side its starting point lies.
@$@<2D ray ray squared distance computation@>+=@{
template <class R>
typename R::FT
squared_distance(
const Ray_2<R> &ray1,
const Ray_2<R> &ray2)
{
typedef typename R::FT FT;
const Vector_2<R> &ray1dir = ray1.direction().vector();
const Vector_2<R> &ray2dir = ray2.direction().vector();
Vector_2<R> diffvec(ray2.source()-ray1.source());
bool crossing1, crossing2;
Orientation dirorder;
dirorder = orientation(ray1dir, ray2dir);
switch (dirorder) {
case COUNTERCLOCKWISE:
crossing1 = !clockwise(diffvec, ray2dir);
crossing2 = !counterclockwise(ray1dir, diffvec);
break;
case CLOCKWISE:
crossing1 = !counterclockwise(diffvec, ray2dir);
crossing2 = !clockwise(ray1dir, diffvec);
break;
case COLLINEAR:
return ray_ray_squared_distance_parallel(ray1dir,ray2dir,diffvec);
}
@}
Depending on the values of the crossing variables we now proceed. If both
rays cross the supporting line of the other ray, the rays
themselves must cross, so they have zero distance.
If exactly one of the rays does not cross, the distance is the distance of
the startpoint of this (non crossing) ray to the other (crossing)
ray.
@$@<2D ray ray squared distance computation@>+=@{
if (crossing1) {
if (crossing2)
return (FT)0;
return squared_distance(ray2.source(), ray1);
} else {
if (crossing2) {
return squared_distance(ray1.source(), ray2);
} else {
@}
Finally, neither of the rays may cross. This may be because the two
rays are parallel, which we deal with as a separate case.
Otherwise we have to compare the distance of two endpoints to the other
ray and take the minimum.
@$@<2D ray ray squared distance computation@>+=@{
FT min1, min2;
min1 = squared_distance(ray1.source(), ray2);
min2 = squared_distance(ray2.source(), ray1);
return (min1 < min2) ? min1 : min2;
}
}
}
@}
The following routine computes the squared distance between two rays if they
are parallel. The parameters are the directions of the rays and the vector
from the start point of the first ray to the start point of the second.
The distance is the distance between the two start points if there is no
line perpendicular to the rays that passes through both rays. Otherwise we can
take the distance between the starting point of the one ray to the supporting
line of the other.
@$@<2D ray ray squared distance parallel case@>==@{
template <class R>
typename R::FT
ray_ray_squared_distance_parallel(
const Vector_2<R> &ray1dir,
const Vector_2<R> &ray2dir,
const Vector_2<R> &from1to2)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
if (!is_acute_angle(ray1dir, from1to2)) {
bool same_direction;
if (CGAL_NTS abs(ray1dir.hx()) > CGAL_NTS abs(ray1dir.hy())) {
same_direction =
(CGAL_NTS sign(ray1dir.hx()) == CGAL_NTS sign(ray2dir.hx()));
} else {
same_direction =
(CGAL_NTS sign(ray1dir.hy()) == CGAL_NTS sign(ray2dir.hy()));
}
if (!same_direction)
return (typename R::FT)(from1to2*from1to2);
}
RT wcr, w;
wcr = wcross(ray1dir, from1to2);
w = from1to2.hw();
return (typename R::FT)(FT(wcr*wcr)
/ FT(wmult((R*)0, RT(wdot(ray1dir, ray1dir)), w, w)));
}
@}
@$@<Ray Ray distance@>==@{
@<2D ray ray squared distance parallel case@>
@<2D ray ray squared distance computation@>
@}
@B@<Line Ray distance@>
@$@<2D squared distance header declarations 1@>+=@{@-
@<Line Ray distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Ray_2<R> & ray,
const Line_2<R> & line)
{
return squared_distance(line, ray);
}
@}
@C@<Implementation@>
First we decide if the ray and the line intersect.
The normal to the line is taken.
If the start point of the ray lies on the normal side of the line and the
direction of the ray has a component in the opposite direction of the normal
(or the other way around), the line and ray intersect and 0 is returned.
Otherwise the squared distance from the start point to the line is reported.
Because we have intermediate results for this computation, we do not call
that routine but compute it using primitive operations.
It would be nice to have a routine to get the normal vector of a line (or of a
direction).
@$@<Line Ray distance implementation@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Line_2<R> &line,
const Ray_2<R> &ray)
{
typedef typename R::FT FT;
Vector_2<R> normalvec(line.a(), line.b());
Vector_2<R> diff = ray.source()-line.point();
FT sign_dist = diff*normalvec;
if (sign_dist < FT(0)) {
if (is_acute_angle(normalvec, ray.direction().vector()) )
return (FT)0;
} else {
if (is_obtuse_angle(normalvec, ray.direction().vector()) )
return (FT)0;
}
return (typename R::FT)((sign_dist*sign_dist)/(normalvec*normalvec));
}
@}
@B@<Line Line distance@>
@$@<2D squared distance header declarations 1@>+=@{@-
@<Line Line distance@>
@}
@C@<Implementation@>
The distance between two lines is zero unless the two lines are paralel. In
that case we can take the distance of an arbitrary point on the one line to
the other line.
@$@<Line Line distance@>==@{
template <class R>
bool
_are_parallel(
const Line_2<R> &line1,
const Line_2<R> &line2)
{
return line1.a()*line2.b() == line2.a()*line1.b();
}
template <class R>
typename R::FT
squared_distance(
const Line_2<R> &line1,
const Line_2<R> &line2)
{
typedef typename R::FT FT;
if (_are_parallel(line1,line2))
return squared_distance(line1.point(), line2);
else
return (FT)0;
}
@}
@B@<Point Triangle distance helper function@>
The @{distance_index@} function computes to which feature of the
triangle the point is closest.
If the point lies inside the triangle, both indices are set to -1.
If the point is closest to the interior of a triangle edge, the indices
are set to the indices of the vertices of that edge, such that the point
lies to the right of the line through vertex(ind1) and vertex(ind2).
Otherwise, the point is closest to a vertex of the triangle. Then ind2 is
set to -1 and vertex(ind1) is the closest vertex.
@$@<2D squared distance header declarations 2@>+=@{
@<Point Triangle distance implementation 1@>
@}
@$@<2D squared distance header declarations 2@>+=@{
@<Point Triangle distance implementation 2@>
@}
@C@<Implementation@>
The implementation discerns between a whole lot of cases.
The outermost if then else clause discerns between a positive oriented
triangle (vertices in counterclockwise order) and a negative oriented.
The code in the else part is identical, except that the role of vertex 1 and 2
of the triangle are swapped.
A point can either lie inside the triangle, be closest to the interior of an
edge or closest to a vertex.
It lies inside the triangle iff it lies not to the right of any supporting
line of the triangle edges (for a positive oriented triangle).
It is closest to an edge if it lies to the right of the edge and if both the
angle edge.source, edge.target, point and edge.target, edge.source, point are
sharp.
Otherwise the point lies closest to a vertex.
@$@<Point Triangle distance implementation 1@>==@{
template <class R>
extern void
distance_index(
int &ind1,
int &ind2,
const Point_2<R> &pt,
const Triangle_2<R> &triangle)
{
const Point_2<R> &vt0 = triangle.vertex(0);
const Point_2<R> &vt1 = triangle.vertex(1);
const Point_2<R> &vt2 = triangle.vertex(2);
if (leftturn(vt0, vt1, vt2)) {
if (rightturn(vt0, vt1, pt)) {
if (!is_acute_angle(vt0, vt1, pt)) {
if (rightturn(vt1, vt2, pt)) {
if (!is_acute_angle(vt1, vt2, pt)) {
ind1 = 2; ind2 = -1;
return;
}
if (!is_acute_angle(vt2, vt1, pt)) {
ind1 = 1; ind2 = -1;
return;
}
ind1 = 1; ind2 = 2;
return;
}
ind1 = 1; ind2 = -1;
return;
}
if (!is_acute_angle(vt1, vt0, pt)) {
if (rightturn(vt2, vt0, pt)) {
if (!is_acute_angle(vt0, vt2, pt)) {
ind1 = 2; ind2 = -1;
return;
}
if (!is_acute_angle(vt2, vt0, pt)) {
ind1 = 0; ind2 = -1;
return;
}
ind1 = 2; ind2 = 0;
return;
}
ind1 = 0; ind2 = -1;
return;
}
ind1 = 0; ind2 = 1;
return;
} else {
if (rightturn(vt1, vt2, pt)) {
if (!is_acute_angle(vt1, vt2, pt)) {
if (rightturn(vt2, vt0, pt)) {
if (!is_acute_angle(vt0, vt2, pt)) {
ind1 = 2; ind2 = -1;
return;
}
if (!is_acute_angle(vt2, vt0, pt)) {
ind1 = 0; ind2 = -1;
return;
}
ind1 = 2; ind2 = 0;
return;
}
ind1 = 0; ind2 = -1;
return;
}
if (!is_acute_angle(vt2, vt1, pt)) {
ind1 = 1; ind2 = -1;
return;
}
ind1 = 1; ind2 = 2;
return;
} else {
if (rightturn(vt2, vt0, pt)) {
if (!is_acute_angle(vt2, vt0, pt)) {
ind1 = 0; ind2 = -1;
return;
}
if (!is_acute_angle(vt0, vt2, pt)) {
ind1 = 2; ind2 = -1;
return;
}
ind1 = 2; ind2 = 0;
return;
} else {
ind1 = -1; ind2 = -1; // point inside or on boundary.
return;
}
}
}
} else {
if (rightturn(vt0, vt2, pt)) {
if (!is_acute_angle(vt0, vt2, pt)) {
if (rightturn(vt2, vt1, pt)) {
if (!is_acute_angle(vt2, vt1, pt)) {
ind1 = 1; ind2 = -1;
return;
}
if (!is_acute_angle(vt1, vt2, pt)) {
ind1 = 2; ind2 = -1;
return;
}
ind1 = 2; ind2 = 1;
return;
}
ind1 = 2; ind2 = -1;
return;
}
if (!is_acute_angle(vt2, vt0, pt)) {
if (rightturn(vt1, vt0, pt)) {
if (!is_acute_angle(vt0, vt1, pt)) {
ind1 = 1; ind2 = -1;
return;
}
if (!is_acute_angle(vt1, vt0, pt)) {
ind1 = 0; ind2 = -1;
return;
}
ind1 = 1; ind2 = 0;
return;
}
ind1 = 0; ind2 = -1;
return;
}
ind1 = 0; ind2 = 2;
return;
} else {
if (rightturn(vt2, vt1, pt)) {
if (!is_acute_angle(vt2, vt1, pt)) {
if (rightturn(vt1, vt0, pt)) {
if (!is_acute_angle(vt0, vt1, pt)) {
ind1 = 1; ind2 = -1;
return;
}
if (!is_acute_angle(vt1, vt0, pt)) {
ind1 = 0; ind2 = -1;
return;
}
ind1 = 1; ind2 = 0;
return;
}
ind1 = 0; ind2 = -1;
return;
}
if (!is_acute_angle(vt1, vt2, pt)) {
ind1 = 2; ind2 = -1;
return;
}
ind1 = 2; ind2 = 1;
return;
} else {
if (rightturn(vt1, vt0, pt)) {
if (!is_acute_angle(vt1, vt0, pt)) {
ind1 = 0; ind2 = -1;
return;
}
if (!is_acute_angle(vt0, vt1, pt)) {
ind1 = 1; ind2 = -1;
return;
}
ind1 = 1; ind2 = 0;
return;
} else {
ind1 = -1; ind2 = -1; // point inside or on boundary.
return;
}
}
}
}
}
@}
@$@<Point Triangle distance implementation 2@>==@{
template <class R>
extern typename R::FT
squared_distance_indexed(const Point_2<R> &pt,
const Triangle_2<R> &triangle, int ind1, int ind2)
{
typedef typename R::FT FT;
if (ind1 == -1)
return FT(0);
if (ind2 == -1)
return squared_distance(pt, triangle.vertex(ind1));
return squared_distance(pt,
Line_2<R>(triangle.vertex(ind1), triangle.vertex(ind2)));
}
@}
@B@<Point Triangle distance@>
@$@<2D squared distance header declarations 2@>+=@{@-
@<Point Triangle distance implementation 3@>
template <class R>
inline typename R::FT
squared_distance(
const Triangle_2<R> & triangle,
const Point_2<R> & pt)
{
return squared_distance(pt, triangle);
}
@}
@C@<Implementation@>
@$@<Point Triangle distance implementation 3@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Point_2<R> &pt,
const Triangle_2<R> &triangle)
{
int ind1,ind2;
distance_index(ind1, ind2, pt, triangle);
return squared_distance_indexed(pt, triangle, ind1, ind2);
}
@}
@B@<Line Triangle distance@>
@$@<2D squared distance header declarations 2@>+=@{@-
@<Line Triangle distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Triangle_2<R> & triangle,
const Line_2<R> & line)
{
return squared_distance(line, triangle);
}
@}
@C@<Implementation@>
@$@<Line Triangle distance implementation@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Line_2<R> &line,
const Triangle_2<R> &triangle)
{
typedef typename R::FT FT;
Oriented_side side0;
side0 = line.oriented_side(triangle.vertex(0));
if (line.oriented_side(triangle.vertex(1)) != side0)
return FT(0);
if (line.oriented_side(triangle.vertex(2)) != side0)
return FT(0);
FT mindist, dist;
int i;
mindist = squared_distance(triangle.vertex(0),line);
for (i=1; i<3; i++) {
dist = squared_distance(triangle.vertex(i),line);
if (dist < mindist)
mindist = dist;
}
return mindist;
}
@}
@B@<Ray Triangle distance@>
@$@<2D squared distance header declarations 2@>+=@{@-
@<Ray Triangle distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Triangle_2<R> & triangle,
const Ray_2<R> & ray)
{
return squared_distance(ray, triangle);
}
@}
@C@<Implementation@>
@$@<Ray Triangle distance implementation@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Ray_2<R> &ray,
const Triangle_2<R> &triangle)
{
typedef typename R::FT FT;
int i, ind_tr1, ind_tr2, ind_ray = 0, ind1;
FT mindist, dist;
distance_index(ind_tr1, ind_tr2, ray.source(), triangle);
mindist =
squared_distance_indexed(ray.source(), triangle, ind_tr1,ind_tr2);
for (i=0; i<3; i++) {
const Point_2<R>& pt = triangle.vertex(i);
distance_index(ind1, pt, ray);
dist = squared_distance_indexed(pt, ray, ind1);
if (dist < mindist) {
ind_ray = ind1;
ind_tr1 = i; ind_tr2 = -1;
mindist = dist;
}
}
// now check if all vertices are on the right side of the separating line.
// In case of vertex-vertex smallest distance this is the case.
if (ind_tr2 == -1 && ind_ray != -1)
return mindist;
if (ind_tr2 != -1) {
// Check if all the segment vertices lie at the same side of
// the triangle segment.
const Point_2<R> &vt1 = triangle.vertex(ind_tr1);
const Point_2<R> &vt2 = triangle.vertex(ind_tr2);
if (clockwise(ray.direction().vector(), vt2-vt1)) {
mindist = FT(0);
}
} else {
// Check if all the triangle vertices lie at the same side of the segment.
const Line_2<R> &sl = ray.supporting_line();
Oriented_side or_s = sl.oriented_side(triangle.vertex(0));
for (i=1; i<3; i++) {
if (sl.oriented_side(triangle.vertex(i)) != or_s) {
mindist = FT(0);
break;
}
}
}
return mindist;
}
@}
@B@<Segment Triangle distance@>
@$@<2D squared distance header declarations 2@>+=@{@-
@<Segment Triangle distance implementation@>
template <class R>
inline typename R::FT
squared_distance(
const Triangle_2<R> & triangle,
const Segment_2<R> & seg)
{
return squared_distance(seg, triangle);
}
@}
@C@<Implementation@>
@$@<Segment Triangle distance implementation@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Segment_2<R> &seg,
const Triangle_2<R> &triangle)
{
typedef typename R::FT FT;
int i, ind_tr1 = 0, ind_tr2 = -1, ind_seg = 0, ind1, ind2;
FT mindist, dist;
mindist = squared_distance(seg.source(), triangle.vertex(0));
for (i=0; i<2; i++) {
const Point_2<R> &pt = seg.vertex(i);
distance_index(ind1, ind2, pt, triangle);
dist = squared_distance_indexed(pt, triangle, ind1, ind2);
if (dist < mindist) {
ind_seg = i;
ind_tr1 = ind1; ind_tr2 = ind2;
mindist = dist;
}
}
for (i=0; i<3; i++) {
const Point_2<R>& pt = triangle.vertex(i);
distance_index(ind1, pt, seg);
dist = squared_distance_indexed(pt, seg, ind1);
if (dist < mindist) {
ind_seg = ind1;
ind_tr1 = i; ind_tr2 = -1;
mindist = dist;
}
}
// now check if all vertices are on the right side of the separating line.
// In case of vertex-vertex smallest distance this is the case.
if (ind_tr2 == -1 && ind_seg != -1)
return mindist;
if (ind_tr2 != -1) {
// Check if all the segment vertices lie at the same side of
// the triangle segment.
const Point_2<R> &vt1 = triangle.vertex(ind_tr1);
const Point_2<R> &vt2 = triangle.vertex(ind_tr2);
Orientation or_s = orientation(vt1, vt2, seg.source());
if (orientation(vt1, vt2, seg.target()) != or_s) {
mindist = FT(0);
}
} else {
// Check if all the triangle vertices lie at the same side of the segment.
const Point_2<R> &vt1 = seg.source();
const Point_2<R> &vt2 = seg.target();
Orientation or_s = orientation(vt1, vt2, triangle.vertex(0));
for (i=1; i<3; i++) {
if (orientation(vt1, vt2, triangle.vertex(i)) != or_s) {
mindist = FT(0);
break;
}
}
}
return mindist;
}
@}
@B@<Triangle Triangle distance@>
@$@<2D squared distance header declarations 2@>+=@{@-
@<Triangle Triangle distance implementation@>
@}
@C@<Implementation@>
@$@<Triangle Triangle distance implementation@>==@{
template <class R>
extern typename R::FT
squared_distance(
const Triangle_2<R> &triangle1,
const Triangle_2<R> &triangle2)
{
typedef typename R::FT FT;
int i, ind1_1 = 0,ind1_2 = -1, ind2_1 = 0, ind2_2 = -1, ind1, ind2;
FT mindist, dist;
mindist =
squared_distance(triangle1.vertex(0), triangle2.vertex(0));
for (i=0; i<3; i++) {
const Point_2<R>& pt = triangle1.vertex(i);
distance_index(ind1, ind2, pt, triangle2);
dist = squared_distance_indexed(pt, triangle2, ind1, ind2);
if (dist < mindist) {
ind1_1 = i; ind1_2 = -1;
ind2_1 = ind1; ind2_2 = ind2;
mindist = dist;
}
}
for (i=0; i<3; i++) {
const Point_2<R>& pt = triangle2.vertex(i);
distance_index(ind1, ind2, pt, triangle1);
dist = squared_distance_indexed(pt, triangle1, ind1, ind2);
if (dist < mindist) {
ind1_1 = ind1; ind1_2 = ind2;
ind2_1 = i; ind2_2 = -1;
mindist = dist;
}
}
// now check if all vertices are on the right side of the separating line.
if (ind1_2 == -1 && ind2_2 == -1)
return mindist;
// In case of point-segment closest distance, there is still the possibility
// of overlapping triangles.
// Check if all the vertices lie at the same side of the segment.
if (ind1_2 != -1) {
const Point_2<R> &vt1 = triangle1.vertex(ind1_1);
const Point_2<R> &vt2 = triangle1.vertex(ind1_2);
Orientation or_s = orientation(vt1, vt2, triangle2.vertex(0));
for (i=1; i<3; i++) {
if (orientation(vt1, vt2, triangle2.vertex(i)) != or_s) {
mindist = FT(0);
break;
}
}
} else {
const Point_2<R> &vt1 = triangle2.vertex(ind2_1);
const Point_2<R> &vt2 = triangle2.vertex(ind2_2);
Orientation or_s = orientation(vt1, vt2, triangle1.vertex(0));
for (i=1; i<3; i++) {
if (orientation(vt1, vt2, triangle1.vertex(i)) != or_s) {
mindist = FT(0);
break;
}
}
}
return mindist;
}
@}
@B@<The header files@>
@O@<../include/CGAL/squared_distance_utils.h@>==@{@-
@<cgal_heading@>@(@-
include/CGAL/squared_distance_utils.h@,@-
sqdistance_2.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_SQUARED_DISTANCE_UTILS_H
#define CGAL_SQUARED_DISTANCE_UTILS_H
#include <CGAL/determinant.h>
#include <CGAL/wmult.h>
CGAL_BEGIN_NAMESPACE
@<2D squared distance general utilities@>
CGAL_END_NAMESPACE
#endif // CGAL_SQUARED_DISTANCE_UTILS_H
@}
@O@<../include/CGAL/squared_distance_2_1.h@>==@{@-
@<cgal_heading@>@(@-
include/CGAL/squared_distance_2_1.h@,@-
sqdistance_2.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_SQUARED_DISTANCE_2_1_H
#define CGAL_SQUARED_DISTANCE_2_1_H
#include <CGAL/user_classes.h>
#include <CGAL/utils.h>
#include <CGAL/Point_2.h>
#include <CGAL/Segment_2.h>
#include <CGAL/Line_2.h>
#include <CGAL/Ray_2.h>
#include <CGAL/enum.h>
#include <CGAL/wmult.h>
#include <CGAL/squared_distance_utils.h>
CGAL_BEGIN_NAMESPACE
@<2D squared distance header declarations 1@>
CGAL_END_NAMESPACE
#endif
@}
@O@<../include/CGAL/squared_distance_2_2.h@>==@{@-
@<cgal_heading@>@(@-
include/CGAL/squared_distance_2_2.h@,@-
sqdistance_2.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_SQUARED_DISTANCE_2_2_H
#define CGAL_SQUARED_DISTANCE_2_2_H
#include <CGAL/user_classes.h>
#include <CGAL/utils.h>
#include <CGAL/Point_2.h>
#include <CGAL/Segment_2.h>
#include <CGAL/Line_2.h>
#include <CGAL/Ray_2.h>
#include <CGAL/Triangle_2.h>
#include <CGAL/enum.h>
#include <CGAL/wmult.h>
#include <CGAL/squared_distance_utils.h>
#include <CGAL/squared_distance_2_1.h>
CGAL_BEGIN_NAMESPACE
@<2D squared distance header declarations 2@>
CGAL_END_NAMESPACE
#endif
@}
@O@<../include/CGAL/squared_distance_2.h@>==@{@-
@<cgal_heading@>@(@-
include/CGAL/squared_distance_2.h@,@-
sqdistance_2.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_SQUARED_DISTANCE_2_H
#define CGAL_SQUARED_DISTANCE_2_H
#include <CGAL/squared_distance_2_1.h>
#include <CGAL/squared_distance_2_2.h>
#endif
@}
Reference in New Issue View Git Blame Copy Permalink