cgal/Packages/Distance_3/web/sqdistance_3.fw

1756 lines
53 KiB
Plaintext

@! $RCSfile$
@! $Revision$
@! $Date$
@! Author: Geert-Jan Giezeman
@i cgal_util.fwi
@A@<Three dimensional squared distances@>
The distance between zero and one dimensional objects (points, segments, rays
and lines) in three dimensions is much the same as the computation in two
dimensions.
If one of the objects is a point, this is not surprising. Then there is a plane
that contains both the point and the other object.
If both objects are one dimensional, we can use the following trick.
There is a normal vector perpendicular to the direction vectors of both objects.
We can project the two objects on a plane perpendicular to this normal and
compute the distance of the projected objects (now a two dimensional problem).
This projection doesn't actually take place in the code. Instead, the one
dimensional object is extended to a two dimensional object in the direction of
the normal. If the one dimensional object is a segment, the result is a strip
formed by the points @{p+s*normal@}, where @{p@} is an arbitrary point on the
segment and @{s@} is an arbitrary real value (positive, negative or zero). In
the same way a plane can be formed by extending a line in the direction of the
normal. If we extend a ray, we also speak of the resulting strip.
In the sequel, when we talk about the extension strip of a segment or ray, or
the extension plane of a line, we mean the strip or plane as defined above.
The extension vector will not be mentioned explicitly, but is always understood
to be the normal vector. The phrase 'extension plane of a segment (or ray)'
will be used as an abbreviation for the 'extension plane of the supporting line
of a segment (or ray)'.
@B@<Utility routines@>
While computing square distances we 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@<Testing for null vector@>
A routine to test whether a vector is a null vector.
This code does not belong in this chapter. It should be moved to the vector
class.
@$@<3D squared distance general utilities@>+=@{
template <class R>
bool CGAL_is_null(const CGAL_Vector_3<R> &v)
{
typedef typename R::RT RT;
return v.hx()==RT(0) && v.hy()==RT(0) && v.hz()==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 or
vectors as arguments.
The homogenising factors that are used for multiplication are different.
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.
@$@<3D squared distance general utilities@>+=@{
template <class R>
R_RT_return(R)
CGAL_wdot(const CGAL_Vector_3<R> &u,
const CGAL_Vector_3<R> &v)
{
return (R_RT_return(R))(u.hx()*v.hx() + u.hy()*v.hy() + u.hz()*v.hz());
}
#ifdef CGAL_HOMOGENEOUS_H
template <class RT>
RT CGAL_wdot(const CGAL_Point_3< CGAL_Homogeneous<RT> > &p,
const CGAL_Point_3<CGAL_Homogeneous<RT> > &q,
const CGAL_Point_3<CGAL_Homogeneous<RT> > &r)
{
return (q.hw()*p.hx() - p.hw()*q.hx())*(q.hw()*r.hx() - r.hw()*q.hx()) +
(q.hw()*p.hy() - p.hw()*q.hy())*(q.hw()*r.hy() - r.hw()*q.hy()) +
(q.hw()*p.hz() - p.hw()*q.hz())*(q.hw()*r.hz() - r.hw()*q.hz());
}
#endif // CGAL_HOMOGENEOUS_H
#ifdef CGAL_CARTESIAN_H
template <class FT>
FT CGAL_wdot(const CGAL_Point_3< CGAL_Cartesian<FT> > &p,
const CGAL_Point_3< CGAL_Cartesian<FT> > &q,
const CGAL_Point_3< CGAL_Cartesian<FT> > &r)
{
return (p.hx() - q.hx())*(r.hx() - q.hx()) +
(p.hy() - q.hy())*(r.hy() - q.hy()) +
(p.hz() - q.hz())*(r.hz() - q.hz());
}
#endif // CGAL_CARTESIAN_H
@}
@D@<The cross product@>
The cross product of two vectors is a vector perpendicular to both.
The absolute value gives the area of the parallelogram spanned by the two
vectors.
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.
@$@<3D squared distance general utilities@>+=@{
template <class R>
CGAL_Vector_3<R> CGAL_wcross(const CGAL_Vector_3<R> &u,
const CGAL_Vector_3<R> &v)
{
return CGAL_Vector_3<R>(
u.hy()*v.hz() - u.hz()*v.hy(),
u.hz()*v.hx() - u.hx()*v.hz(),
u.hx()*v.hy() - u.hy()*v.hx());
}
#ifdef CGAL_HOMOGENEOUS_H
template <class RT>
CGAL_Vector_3< CGAL_Homogeneous<RT> >
CGAL_wcross(const CGAL_Point_3< CGAL_Homogeneous<RT> > &p,
const CGAL_Point_3< CGAL_Homogeneous<RT> > &q,
const CGAL_Point_3< CGAL_Homogeneous<RT> > &r)
{
RT x,y,z;
x = p.hy() * (q.hz()*r.hw() - q.hw()*r.hz() )
+ p.hz() * (q.hw()*r.hy() - q.hy()*r.hw() )
+ p.hw() * (q.hy()*r.hz() - q.hz()*r.hy() );
y = p.hz() * (q.hx()*r.hw() - q.hw()*r.hx() )
+ p.hx() * (q.hw()*r.hz() - q.hz()*r.hw() )
+ p.hw() * (q.hz()*r.hx() - q.hx()*r.hz() );
z = 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() );
return CGAL_Vector_3< CGAL_Homogeneous<RT> >(x, y, z);
}
#endif // CGAL_HOMOGENEOUS_H
#ifdef CGAL_CARTESIAN_H
template <class FT>
CGAL_Vector_3< CGAL_Cartesian<FT> >
CGAL_wcross(const CGAL_Point_3< CGAL_Cartesian<FT> > &p,
const CGAL_Point_3< CGAL_Cartesian<FT> > &q,
const CGAL_Point_3< CGAL_Cartesian<FT> > &r)
{
FT x,y,z;
x = (q.y()-p.y())*(r.z()-q.z()) - (q.z()-p.z())*(r.y()-q.y());
y = (q.z()-p.z())*(r.x()-q.x()) - (q.x()-p.x())*(r.z()-q.z());
z = (q.x()-p.x())*(r.y()-q.y()) - (q.y()-p.y())*(r.x()-q.x());
return CGAL_Vector_3< CGAL_Cartesian<FT> >(x, y, z);
}
#endif // CGAL_CARTESIAN_H
@}
@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 (@{CGAL_ACUTE@},
@{CGAL_STRAIGHT@},@{CGAL_OBTUSE@}) and a routine that returns such a value.
@$@<3D squared distance general utilities@>+=@{
template <class R>
inline bool CGAL_is_acute_angle(const CGAL_Vector_3<R> &u,
const CGAL_Vector_3<R> &v)
{
typedef typename R::RT RT;
return RT(CGAL_wdot(u, v)) > RT(0) ;
}
template <class R>
inline bool CGAL_is_straight_angle(const CGAL_Vector_3<R> &u,
const CGAL_Vector_3<R> &v)
{
typedef typename R::RT RT;
return RT(CGAL_wdot(u, v)) == RT(0) ;
}
template <class R>
inline bool CGAL_is_obtuse_angle(const CGAL_Vector_3<R> &u,
const CGAL_Vector_3<R> &v)
{
typedef typename R::RT RT;
return RT(CGAL_wdot(u, v)) < RT(0) ;
}
template <class R>
inline bool CGAL_is_acute_angle(const CGAL_Point_3<R> &p,
const CGAL_Point_3<R> &q, const CGAL_Point_3<R> &r)
{
typedef typename R::RT RT;
return RT(CGAL_wdot(p, q, r)) > RT(0) ;
}
template <class R>
inline bool CGAL_is_straight_angle(const CGAL_Point_3<R> &p,
const CGAL_Point_3<R> &q, const CGAL_Point_3<R> &r)
{
typedef typename R::RT RT;
return RT(CGAL_wdot(p, q, r)) == RT(0) ;
}
template <class R>
inline bool CGAL_is_obtuse_angle(const CGAL_Point_3<R> &p,
const CGAL_Point_3<R> &q, const CGAL_Point_3<R> &r)
{
typedef typename R::RT RT;
return RT(CGAL_wdot(p, q, r)) < RT(0) ;
}
@}
@B@<Point Point distance@>
The squared distance between two points is computed by taking the dot product
of the difference vector.
@$@<3D squared distance 0 header declarations@>+=@{@-
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> & pt1,
const CGAL_Point_3<R> & pt2)
{
CGAL_Vector_3<R> vec(pt1-pt2);
return vec*vec;
}
@}
@B@<Point Plane distance@>
@C@<Declarations@>
The regular squared distance routines for the point plane case are given
further on. They are defined in another header file.
Here we define a specialised routine that computes this squared distance with
the following parameters: the normal vector to the plane and the difference
vector between the point and an arbitrary point on the plane.
This specialised routine is used extensively in other routines.
@$@<3D squared distance 0 header declarations@>+=@{@-
template <class R>
R_FT_return(R)
CGAL_squared_distance_to_plane(
const CGAL_Vector_3<R> & normal,
const CGAL_Vector_3<R> & diff);
@}
@C@<Implementation@>
@$@<3D squared distance 0 outline definitions@>+=@{@-
template <class R>
R_FT_return(R)
CGAL_squared_distance_to_plane(
const CGAL_Vector_3<R> & normal,
const CGAL_Vector_3<R> & diff)
{
typedef typename R::RT RT;
RT dot, squared_length;
dot = CGAL_wdot(normal, diff);
squared_length = CGAL_wdot(normal, normal);
return (R_FT_return(R))
((dot*dot) / CGAL_wmult((R*)0, squared_length, diff.hw(), diff.hw()));
}
@}
@B@<Point Line distance@>
@C@<Declarations@>
The third function is meant for internal use. It computes the squared
distance from a point to a line, given a direction vector (not necessarily
normalised) of the line and the vector between the point and an arbitrary
point on the line.
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Line_3<R> &line);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> & line,
const CGAL_Point_3<R> & pt)
{
return CGAL_squared_distance(pt, line);
}
@}
@$@<3D squared distance 0 header declarations@>+=@{@-
template <class R>
R_FT_return(R)
CGAL_squared_distance_to_line(
const CGAL_Vector_3<R> & dir,
const CGAL_Vector_3<R> & diff);
@}
@C@<Implementation@>
The cross product of two vectors is a vector whose length is equal to the area
of the parallellogram spanned by the two vectors. This area is length of base
times height.
We can use this to compute the distance from a point to a line. As base vector
we take the direction vector of the line. The other vector is the vector
connecting a point on the line to the point. Now, the length of the cross
product vector divided by the length of the base vector is the distance from
the point to the line. Squaring everything, we get the square distance.
@$@<3D squared distance 0 outline definitions@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance_to_line(
const CGAL_Vector_3<R> & dir,
const CGAL_Vector_3<R> & diff)
{
typedef typename R::RT RT;
CGAL_Vector_3<R> wcr = CGAL_wcross(dir, diff);
return (R_FT_return(R))((wcr*wcr)/
CGAL_wmult((R*)0, RT(CGAL_wdot(dir, dir)), diff.hw(), diff.hw()));
}
@}
@$@<3D squared distance outline definitions@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Line_3<R> &line)
{
CGAL_Vector_3<R> dir(line.direction().vector());
CGAL_Vector_3<R> diff = pt - line.point();
return CGAL_squared_distance_to_line(dir, diff);
}
@}
@B@<Point Ray distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Ray_3<R> &ray);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> & ray,
const CGAL_Point_3<R> & pt)
{
return CGAL_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,
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.
@$@<3D squared distance outline definitions@>+=@{
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Ray_3<R> &ray)
{
CGAL_Vector_3<R> diff = pt-ray.start();
const CGAL_Vector_3<R> &dir = ray.direction().vector();
if (!CGAL_is_acute_angle(dir,diff) )
return (R_FT_return(R))(diff*diff);
return CGAL_squared_distance_to_line(dir, diff);
}
@}
@B@<Point Segment distance@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Segment_3<R> &seg);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> & seg,
const CGAL_Point_3<R> & pt)
{
return CGAL_squared_distance(pt, seg);
}
@}
@C@<Implementation@>
To compute the distance of a point to a segment, first we divide the space in
three parts. We draw two planes 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 planes, 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.
@$@<3D squared distance outline definitions@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Segment_3<R> &seg)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
// assert that the segment is valid (non zero length).
CGAL_Vector_3<R> diff = pt-seg.start();
CGAL_Vector_3<R> segvec = seg.end()-seg.start();
RT d = CGAL_wdot(diff,segvec);
if (d <= (RT)0)
return (R_FT_return(R))(FT(diff*diff));
RT e = CGAL_wdot(segvec,segvec);
if (CGAL_wmult((R*)0 ,d, segvec.hw()) > CGAL_wmult((R*)0, e, diff.hw()))
return CGAL_squared_distance(pt, seg.end());
return CGAL_squared_distance_to_line(segvec, diff);
}
@}
@B@<Segment Segment distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &,
const CGAL_Segment_3<R> &);
@}
@C@<Implementation@>
The computation of the distance between two segments in 3D goes along the same
lines as the computation in 2D. An important role is played by the @{normal@}
vector, which is a vector perpendicular to both segments.
@D@<Distance computation@>
To compute the squared distance of two segments, the following steps are
taken.
@$@<3D seg seg squared distance main routine start@>==@{@-
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg1,
const CGAL_Segment_3<R> &seg2)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const CGAL_Point_3<R> &start1 = seg1.start();
const CGAL_Point_3<R> &start2 = seg2.start();
const CGAL_Point_3<R> &end1 = seg1.end();
const CGAL_Point_3<R> &end2 = seg2.end();
@}
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 compute a vector that is normal to the direction vectors of both
segments. This is done by taking the cross product of the two vectors.
If the normal is the null vector, the segments must be parallel to each other.
This is dealt with as a special case.
@$@<3D seg seg squared distance main routine degeneracy checking@>==@{@-
if (start1 == end1)
return CGAL_squared_distance(start1, seg2);
if (start2 == end2)
return CGAL_squared_distance(start2, seg1);
CGAL_Vector_3<R> dir1, dir2, normal;
dir1 = seg1.direction().vector();
dir2 = seg2.direction().vector();
normal = CGAL_wcross(dir1, dir2);
if (CGAL_is_null(normal))
return CGAL_squared_distance_parallel(seg1, seg2);
@}
We decide if the segments intersects the extension plane of the other segment.
See the top of this chapter on three dimensional squared distances for the
definition of the extension plane of a segment. The boolean @{crossing1@}
tells whether @{seg1@} intersects the extension plane of @{seg2@}, and vice
versa for @{crossing2@}. If an endpoint lies in the plane, we consider this as
an intersection.
To decide whether a segment and a plane intersect, we compute the signed
distance of the endpoints of the segment to the plane. We first compute the
vectors perpendicular to the planes (@{perpend1@} and @{perpend2@}). The signed
distance is then computed by taking the dot product between this vector and the
vectors connecting endpoints of the segment to a point in the plane.
As the @{perpend@} vectors are not normalised to length 1, this is not the
exact distance. What we compute is a measure related to the signed distance.
For the moment we are only interested in the sign of the distances (positive
means on one side of the plane, negative on the other side and 0 means in the
plane), and this is retained in this measure.
We store the signed distance measures in variables which we will need later
on. The names are a bit cryptic: @{sdm_s1to2@} stands for the signed distance
measure from the start point of @{seg1@} to the extension plane of @{seg2@}.
@$@<3D seg seg squared distance main routine crossing computation@>==@{@-
bool crossing1, crossing2;
RT sdm_s1to2, sdm_e1to2, sdm_s2to1, sdm_e2to1;
CGAL_Vector_3<R> perpend1, perpend2, s2mins1, e2mins1, e1mins2;
perpend1 = CGAL_wcross(dir1, normal);
perpend2 = CGAL_wcross(dir2, normal);
s2mins1 = start2-start1;
e2mins1 = end2-start1;
e1mins2 = end1-start2;
sdm_s1to2 = -RT(CGAL_wdot(perpend2, s2mins1));
sdm_e1to2 = CGAL_wdot(perpend2, e1mins2);
sdm_s2to1 = CGAL_wdot(perpend1, s2mins1);
sdm_e2to1 = CGAL_wdot(perpend1, e2mins1);
if (sdm_s1to2 < RT(0)) {
crossing1 = (sdm_e1to2 >= RT(0));
} else {
if (sdm_e1to2 <= RT(0)) {
crossing1 = true;
} else {
crossing1 = (sdm_s1to2 == RT(0));
}
}
if (sdm_s2to1 < RT(0)) {
crossing2 = (sdm_e2to1 >= RT(0));
} else {
if (sdm_e2to1 <= RT(0)) {
crossing2 = true;
} else {
crossing2 = (sdm_s2to1 == RT(0));
}
}
@}
Depending on the values of the crossing variables we now proceed. If both
segments intersect the extension plane of the other segment, the extension
strips of the segments themselves must cross. In this case we compute the
distance from an arbitrary point of segment 1 to the plane perpendicular to
@{normal@} that goes through segment 2.
We use the special routine that computes this squared distance from a point to
a plane, given the normal to the plane and the vector from the point to an
arbitrary point in the plane.
@$@<3D seg seg squared distance main routine two crossings case@>==@{@-
if (crossing1) {
if (crossing2) {
return CGAL_squared_distance_to_plane(normal, s2mins1);
}
@}
If exactly one of the segments does not intersect, the distance is the
distance of one of the endpoints of this (non intersecting) segment to the
other (intersecting) segment. Which endpoint? The one that is closest to
the extension plane of the other segment. This latter property can easily
be computed with the value of the signed distance measure that was
computed earlier. The function @{CGAL__distance_measure_sub@} makes this
comparison (it subtracts the two absolute values of the signed distance
measure, first corrected by homogenising factors).
@$@<3D seg seg squared distance main routine one crossing case@>==@{@-
RT dm;
dm = CGAL__distance_measure_sub(
sdm_s2to1, sdm_e2to1, s2mins1, e2mins1);
if (dm < RT(0)) {
return CGAL_squared_distance(start2, seg1);
} else {
if (dm > RT(0)) {
return CGAL_squared_distance(end2, seg1);
} else {
// should not happen with exact arithmetic.
return CGAL_squared_distance_parallel(seg1, seg2);
}
}
} else {
if (crossing2) {
RT dm;
dm =CGAL__distance_measure_sub(
sdm_s1to2, sdm_e1to2, s2mins1, e1mins2);
if (dm < RT(0)) {
return CGAL_squared_distance(start1, seg2);
} else {
if (dm > RT(0)) {
return CGAL_squared_distance(end1, seg2);
} else {
// should not happen with exact arithmetic.
return CGAL_squared_distance_parallel(seg1, seg2);
}
}
} else {
@}
Finally, neither of the segments may cross. This is the case when the two
segments are parallel, which we deal with as a separate case. This case should
have been detected before, when the normal was checked for being null, but we
check anyway to cope with approximation errors due to possible inexact
arithmetic.
Otherwise we have to compare the squared distance of two endpoints to the
other segment and take the minimum.
@$@<3D seg seg squared distance main routine no crossing case@>==@{@-
FT min1, min2;
RT dm;
dm = CGAL__distance_measure_sub(
sdm_s1to2, sdm_e1to2, s2mins1, e1mins2);
if (dm == RT(0)) // should not happen with exact arithmetic.
return CGAL_squared_distance_parallel(seg1, seg2);
min1 = (dm < RT(0)) ?
CGAL_squared_distance(seg1.start(), seg2):
CGAL_squared_distance(end1, seg2);
dm = CGAL__distance_measure_sub(
sdm_s2to1, sdm_e2to1, s2mins1, e2mins1);
if (dm == RT(0)) // should not happen with exact arithmetic.
return CGAL_squared_distance_parallel(seg1, seg2);
min2 = (dm < RT(0)) ?
CGAL_squared_distance(start2, seg1):
CGAL_squared_distance(end2, seg1);
return (min1 < min2)
? (R_FT_return(R))(min1)
: (R_FT_return(R))(min2);
}
}
@}
We need a small helper function to compare the distance from a point to the
the extension planes of the two segments. With the help
of the signed distance measure (that is known), this can be computed cheaply.
The previously computed signed distance measure differ from the actual signed
distance by several factors. First, they they are computed by taking the dot
product with a vector that was not normalised to 1. This is not important,
because this multiplies both signed distance measures by the same factor.
The other thing is that we took the wdot product.
For this we have to compensate, which accounts for the last two parameters.
@$@<3D segment segment squared distance helper routine@>==@{
template <class RT, class R>
RT CGAL__distance_measure_sub(RT startwdist, RT endwdist,
const CGAL_Vector_3<R> &start, const CGAL_Vector_3<R> &end
)
{
return CGAL_abs(CGAL_wmult((R*)0, startwdist, end.hw())) -
CGAL_abs(CGAL_wmult((R*)0, endwdist, 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 first case. The projections are just a way
of explaining, they are not actually done. The acuteness of angles is used
to make the necessary decisions.
@$@<3D segment segment squared distance parallel case@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance_parallel(
const CGAL_Segment_3<R> &seg1,
const CGAL_Segment_3<R> &seg2)
{
bool same_direction;
const CGAL_Vector_3<R> &dir1 = seg1.direction().vector();
const CGAL_Vector_3<R> &dir2 = seg2.direction().vector();
if (CGAL_abs(dir1.hx()) > CGAL_abs(dir1.hy())) {
if (CGAL_abs(dir1.hx()) > CGAL_abs(dir1.hz())) {
same_direction = (CGAL_sign(dir1.hx()) == CGAL_sign(dir2.hx()));
} else {
same_direction = (CGAL_sign(dir1.hz()) == CGAL_sign(dir2.hz()));
}
} else {
if (CGAL_abs(dir1.hy()) > CGAL_abs(dir1.hz())) {
same_direction = (CGAL_sign(dir1.hy()) == CGAL_sign(dir2.hy()));
} else {
same_direction = (CGAL_sign(dir1.hz()) == CGAL_sign(dir2.hz()));
}
}
if (same_direction) {
if (!CGAL_is_acute_angle(seg1.start(), seg1.end(), seg2.start()))
return CGAL_squared_distance(seg1.end(), seg2.start());
if (!CGAL_is_acute_angle(seg1.end(), seg1.start(), seg2.end()))
return CGAL_squared_distance(seg1.start(), seg2.end());
} else {
if (!CGAL_is_acute_angle(seg1.start(), seg1.end(), seg2.end()))
return CGAL_squared_distance(seg1.end(), seg2.end());
if (!CGAL_is_acute_angle(seg1.end(), seg1.start(), seg2.start()))
return CGAL_squared_distance(seg1.start(), seg2.start());
}
return CGAL_squared_distance(seg2.start(), seg1.supporting_line());
}
@}
@$@<3D segment segment squared distance main routine@>==@{@-
@<3D seg seg squared distance main routine start@>
@<3D seg seg squared distance main routine degeneracy checking@>
@<3D seg seg squared distance main routine crossing computation@>
@<3D seg seg squared distance main routine two crossings case@>
@<3D seg seg squared distance main routine one crossing case@>
@<3D seg seg squared distance main routine no crossing case@>
}
@}
@$@<3D squared distance outline definitions@>+=@{
@<3D segment segment squared distance parallel case@>
@<3D segment segment squared distance helper routine@>
@<3D segment segment squared distance main routine@>
@}
@B@<Segment Ray distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg,
const CGAL_Ray_3<R> &ray);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> & ray,
const CGAL_Segment_3<R> & seg)
{
return CGAL_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 compute a vector that is normal to the direction vectors of both the
segment and the ray. This is done by taking the cross product of the two
vectors. If the normal is the null vector, the objects must be parallel to
each other. This is dealt with as a special case.
@$@<3D ray segment squared distance computation@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg,
const CGAL_Ray_3<R> &ray)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const CGAL_Point_3<R> & ss (seg.start());
const CGAL_Point_3<R> & se (seg.end());
if (ss == se)
return CGAL_squared_distance(ss, ray);
CGAL_Vector_3<R> raydir, segdir, normal;
raydir = ray.direction().vector();
segdir = seg.direction().vector();
normal = CGAL_wcross(segdir, raydir);
if (CGAL_is_null(normal))
return CGAL_squared_distance_parallel(seg, ray);
@}
We decide if the objects intersect the extension plane of the other object.
The boolean @{crossing1@} tells whether @{seg@} intersects the extension plane
of @{ray@}, and vice versa for @{crossing2@}. If an endpoint lies in the
plane, we consider this as an intersection.
This intersection computation goes in the same way as the two segments case.
The only difference is that the ray only has a start point, not an end point,
so we can not compute the signed distance measure of the end point to the
extension plane of the segment. Instead, we can take the signed distance
measure of a point that is computed by adding the direction vector of the ray
to an arbitrary point of the extension plane of the segment.
As for naming: @{ss_min_rs@} stands for the segment start point minus the ray
start point. @{sdm_se2r@} stands for the signed distance measure(@{sdm@}) from
the segment end point(@{se@}) to the extension plane of the ray(@{r@}).
@$@<3D ray segment squared distance computation@>+=@{
bool crossing1, crossing2;
RT sdm_ss2r, sdm_se2r, sdm_rs2s, sdm_re2s;
CGAL_Vector_3<R> perpend2seg, perpend2ray, ss_min_rs, se_min_rs;
perpend2seg = CGAL_wcross(segdir, normal);
perpend2ray = CGAL_wcross(raydir, normal);
ss_min_rs = ss-ray.start();
se_min_rs = se-ray.start();
sdm_ss2r = CGAL_wdot(perpend2ray, ss_min_rs);
sdm_se2r = CGAL_wdot(perpend2ray, se_min_rs);
if (sdm_ss2r < RT(0)) {
crossing1 = (sdm_se2r >= RT(0));
} else {
if (sdm_se2r <= RT(0)) {
crossing1 = true;
} else {
crossing1 = (sdm_ss2r == RT(0));
}
}
sdm_rs2s = -RT(CGAL_wdot(perpend2seg, ss_min_rs));
sdm_re2s = CGAL_wdot(perpend2seg, raydir);
if (sdm_rs2s < RT(0)) {
crossing2 = (sdm_re2s >= RT(0));
} else {
if (sdm_re2s <= RT(0)) {
crossing2 = true;
} else {
crossing2 = (sdm_rs2s == RT(0));
}
}
@}
Depending on the values of the crossing variables we now proceed. If both
the segment and the ray intersect the extension plane of the other object,
the extension strips of the objects intersect.
In this case we compute the
distance from an arbitrary point of the segment to the plane perpendicular to
@{normal@} that goes through the ray.
Also the other cases go very similar to the segment segment case, so we omit
further commentary.
@$@<3D ray segment squared distance computation@>+=@{
if (crossing1) {
if (crossing2) {
return CGAL_squared_distance_to_plane(normal, ss_min_rs);
}
return CGAL_squared_distance(ray.start(), seg);
} else {
if (crossing2) {
RT dm;
dm = CGAL__distance_measure_sub(
sdm_ss2r, sdm_se2r, ss_min_rs, se_min_rs);
if (dm < RT(0)) {
return CGAL_squared_distance(ss, ray);
} else {
if (dm > RT(0)) {
return CGAL_squared_distance(se, ray);
} else {
// parallel, should not happen (no crossing)
return CGAL_squared_distance_parallel(seg, ray);
}
}
} else {
FT min1, min2;
RT dm;
dm = CGAL__distance_measure_sub(
sdm_ss2r, sdm_se2r, ss_min_rs, se_min_rs);
if (dm == RT(0))
return CGAL_squared_distance_parallel(seg, ray);
min1 = (dm < RT(0))
? CGAL_squared_distance(ss, ray)
: CGAL_squared_distance(se, ray);
min2 = CGAL_squared_distance(ray.start(), seg);
return (min1 < min2)
? (R_FT_return(R))(min1)
: (R_FT_return(R))(min2);
}
}
}
@}
@$@<3D ray segment squared distance parallel case@>==@{
template <class R>
R_FT_return(R)
CGAL_squared_distance_parallel(
const CGAL_Segment_3<R> &seg,
const CGAL_Ray_3<R> &ray)
{
bool same_direction;
const CGAL_Vector_3<R> &dir1 = seg.direction().vector();
const CGAL_Vector_3<R> &dir2 = ray.direction().vector();
if (CGAL_abs(dir1.hx()) > CGAL_abs(dir1.hy())) {
same_direction = (CGAL_sign(dir1.hx()) == CGAL_sign(dir2.hx()));
} else {
same_direction = (CGAL_sign(dir1.hy()) == CGAL_sign(dir2.hy()));
}
if (same_direction) {
if (!CGAL_is_acute_angle(seg.start(), seg.end(), ray.start()))
return CGAL_squared_distance(seg.end(), ray.start());
} else {
if (!CGAL_is_acute_angle(seg.end(), seg.start(), ray.start()))
return CGAL_squared_distance(seg.start(), ray.start());
}
return CGAL_squared_distance(ray.start(), seg.supporting_line());
}
@}
@$@<3D squared distance outline definitions@>+=@{
@! @<3D ray segment squared distance helper routine@>
@<3D ray segment squared distance parallel case@>
@<3D ray segment squared distance computation@>
@}
@B@<Segment Line distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg,
const CGAL_Line_3<R> &line);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> & line,
const CGAL_Segment_3<R> & seg)
{
return CGAL_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 squared distance of two segments.
We do not repeat the comments of the code that is the same.
@$@<3D line segment squared distance computation@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg,
const CGAL_Line_3<R> &line)
{
typedef typename R::RT RT;
const CGAL_Point_3<R> &linepoint = line.point();
const CGAL_Point_3<R> &start = seg.start();
const CGAL_Point_3<R> &end = seg.end();
@}
We start with checking for degenerate cases.
@$@<3D line segment squared distance computation@>+=@{
if (start == end)
return CGAL_squared_distance(start, line);
CGAL_Vector_3<R> linedir = line.direction().vector();
CGAL_Vector_3<R> segdir = seg.direction().vector();
CGAL_Vector_3<R> normal = CGAL_wcross(segdir, linedir);
if (CGAL_is_null(normal))
return CGAL_squared_distance_to_line(linedir, start-linepoint);
@}
Then we test for intersection with the extension planes. Things are a bit
easier in this case. We do not have to check if the line intersects the
extension plane of the segment. It always does, unless the segment and the line
are parallel to each other. This is the reason why there is only one
@{crossing@} variable.
@$@<3D line segment squared distance computation@>+=@{
bool crossing;
RT sdm_ss2l, sdm_se2l;
CGAL_Vector_3<R> perpend2line, start_min_lp, end_min_lp;
perpend2line = CGAL_wcross(linedir, normal);
start_min_lp = start-linepoint;
end_min_lp = end-linepoint;
sdm_ss2l = CGAL_wdot(perpend2line, start_min_lp);
sdm_se2l = CGAL_wdot(perpend2line, end_min_lp);
if (sdm_ss2l < RT(0)) {
crossing = (sdm_se2l >= RT(0));
} else {
if (sdm_se2l <= RT(0)) {
crossing = true;
} else {
crossing = (sdm_ss2l == RT(0));
}
}
@}
Finally, we compute the squared distance, depending on which case applies.
@$@<3D line segment squared distance computation@>+=@{
if (crossing) {
return CGAL_squared_distance_to_plane(normal, start_min_lp);
} else {
RT dm;
dm = CGAL__distance_measure_sub(
sdm_ss2l, sdm_se2l, start_min_lp, end_min_lp);
if (dm <= RT(0)) {
return CGAL_squared_distance_to_line(linedir, start_min_lp);
} else {
return CGAL_squared_distance_to_line(linedir, end_min_lp);
}
}
}
@}
@$@<3D squared distance outline definitions@>+=@{
@<3D line segment squared distance computation@>
@}
@B@<Ray Ray distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> &ray1,
const CGAL_Ray_3<R> &ray2);
@}
@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.
Look at that routine for comments.
@$@<3D ray ray squared distance computation@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> &ray1,
const CGAL_Ray_3<R> &ray2)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const CGAL_Point_3<R> & s1 (ray1.start());
const CGAL_Point_3<R> & s2 (ray2.start());
CGAL_Vector_3<R> dir1, dir2, normal;
dir1 = ray1.direction().vector();
dir2 = ray2.direction().vector();
normal = CGAL_wcross(dir1, dir2);
CGAL_Vector_3<R> s1_min_s2 = s1-s2;
if (CGAL_is_null(normal))
return CGAL_ray_ray_squared_distance_parallel(dir1, dir2, s1_min_s2);
bool crossing1, crossing2;
RT sdm_s1_2, sdm_s2_1;
CGAL_Vector_3<R> perpend1, perpend2;
perpend1 = CGAL_wcross(dir1, normal);
perpend2 = CGAL_wcross(dir2, normal);
sdm_s1_2 = CGAL_wdot(perpend2, s1_min_s2);
if (sdm_s1_2 < RT(0)) {
crossing1 = (RT(CGAL_wdot(perpend2, dir1)) >= RT(0));
} else {
if (RT(CGAL_wdot(perpend2, dir1)) <= RT(0)) {
crossing1 = true;
} else {
crossing1 = (sdm_s1_2 == RT(0));
}
}
sdm_s2_1 = -RT(CGAL_wdot(perpend1, s1_min_s2));
if (sdm_s2_1 < RT(0)) {
crossing2 = (RT(CGAL_wdot(perpend1, dir2)) >= RT(0));
} else {
if (RT(CGAL_wdot(perpend1, dir2)) <= RT(0)) {
crossing2 = true;
} else {
crossing2 = (sdm_s2_1 == RT(0));
}
}
if (crossing1) {
if (crossing2)
return CGAL_squared_distance_to_plane(normal, s1_min_s2);
return CGAL_squared_distance(ray2.start(), ray1);
} else {
if (crossing2) {
return CGAL_squared_distance(ray1.start(), ray2);
} else {
FT min1, min2;
min1 = CGAL_squared_distance(ray1.start(), ray2);
min2 = CGAL_squared_distance(ray2.start(), ray1);
return (min1 < min2)
? (R_FT_return(R))(min1)
: (R_FT_return(R))(min2);
}
}
}
@}
The following routine computes the squared distance between two rays if they
are parallel. The parameters are the two 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.
@$@<3D ray ray squared distance parallel case@>==@{
template <class R>
R_FT_return(R)
CGAL_ray_ray_squared_distance_parallel(
const CGAL_Vector_3<R> &ray1dir,
const CGAL_Vector_3<R> &ray2dir,
const CGAL_Vector_3<R> &s1_min_s2)
{
if (!CGAL_is_acute_angle(ray2dir, s1_min_s2)) {
bool same_direction;
if (CGAL_abs(ray1dir.hx()) > CGAL_abs(ray1dir.hy())) {
if (CGAL_abs(ray1dir.hx()) > CGAL_abs(ray1dir.hz()))
same_direction =
(CGAL_sign(ray1dir.hx()) == CGAL_sign(ray2dir.hx()));
else
same_direction =
(CGAL_sign(ray1dir.hz()) == CGAL_sign(ray2dir.hz()));
} else {
if (CGAL_abs(ray1dir.hy()) > CGAL_abs(ray1dir.hz()))
same_direction =
(CGAL_sign(ray1dir.hy()) == CGAL_sign(ray2dir.hy()));
else
same_direction =
(CGAL_sign(ray1dir.hz()) == CGAL_sign(ray2dir.hz()));
}
if (!same_direction)
return (R_FT_return(R))(s1_min_s2*s1_min_s2);
}
return CGAL_squared_distance_to_line(ray1dir, s1_min_s2);
}
@}
@$@<3D squared distance outline definitions@>+=@{
@<3D ray ray squared distance parallel case@>
@<3D ray ray squared distance computation@>
@}
@B@<Line Ray distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> &line,
const CGAL_Ray_3<R> &ray);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> & ray,
const CGAL_Line_3<R> & line)
{
return CGAL_squared_distance(line, ray);
}
@}
@C@<Implementation@>
The computation of the squared distance a line and a ray goes along the same
lines as the computation of the square distance between two segments. With
fewer endpoints, things get ever simpler to compute.
@$@<3D squared distance outline definitions@>+=@{
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> &line,
const CGAL_Ray_3<R> &ray)
{
typedef typename R::RT RT;
const CGAL_Point_3<R> & rs (ray.start());
CGAL_Vector_3<R> raydir, linedir, normal;
linedir = line.direction().vector();
raydir = ray.direction().vector();
normal = CGAL_wcross(raydir, linedir);
CGAL_Vector_3<R> rs_min_lp = rs-line.point();
if (CGAL_is_null(normal))
return CGAL_squared_distance_to_line(linedir, rs_min_lp);
bool crossing;
RT sdm_sr_l;
CGAL_Vector_3<R> perpend2l;
perpend2l = CGAL_wcross(linedir, normal);
sdm_sr_l = CGAL_wdot(perpend2l, rs_min_lp);
if (sdm_sr_l < RT(0)) {
crossing = (RT(CGAL_wdot(perpend2l, raydir)) >= RT(0));
} else {
if (RT(CGAL_wdot(perpend2l, raydir)) <= RT(0)) {
crossing = true;
} else {
crossing = (sdm_sr_l == RT(0));
}
}
if (crossing)
return CGAL_squared_distance_to_plane(normal, rs_min_lp);
else
return CGAL_squared_distance_to_line(linedir, rs_min_lp);
}
@}
@B@<Line Line distance@>
@C@<Declarations@>
@$@<3D squared distance header declarations@>+=@{
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> &line1,
const CGAL_Line_3<R> &line2);
@}
@C@<Implementation@>
We discern two cases: parallel and non-parallel lines.
For parallel lines, we can take the squared distance from an arbitrary point
on one line to the other line.
Otherwise, the problem can be reduced to the squared distance from a point to a
plane.
@$@<3D squared distance outline definitions@>+=@{
template <class R>
R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> &line1,
const CGAL_Line_3<R> &line2)
{
CGAL_Vector_3<R> dir1, dir2, normal, diff;
dir1 = line1.direction().vector();
dir2 = line2.direction().vector();
normal = CGAL_wcross(dir1, dir2);
diff = line2.point() - line1.point();
if (CGAL_is_null(normal))
return CGAL_squared_distance_to_line(dir2, diff);
return CGAL_squared_distance_to_plane(normal, diff);
}
@}
@B@<Point Plane distance@>
@C@<Declarations@>
@$@<3D squared distance 2 header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> &pt,
const CGAL_Plane_3<R> &plane);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Plane_3<R> & plane,
const CGAL_Point_3<R> & pt)
{
return CGAL_squared_distance(pt, plane);
}
@}
@C@<Implementation@>
@$@<3D squared distance 2 outline definitions@>+=@{@-
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Point_3<R> & pt,
const CGAL_Plane_3<R> & plane)
{
CGAL_Vector_3<R> diff(pt-plane.point());
return CGAL_squared_distance_to_plane(plane.orthogonal_vector(), diff);
}
@}
@B@<3D Plane Line distance@>
The following routine decides whether a vector lies in a plane. This is useful
for detecting whether a line is parallel to a plane.
@$@<3D squared distance general utilities 2@>==@{
template <class R>
bool
contains_vector(const CGAL_Plane_3<R> &pl, const CGAL_Vector_3<R> &vec)
{
typedef typename R::RT RT;
return pl.a()*vec.hx() + pl.b()*vec.hy() + pl.c() * vec.hz() == RT(0);
}
@}
@C@<Declarations@>
@$@<3D squared distance 2 header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> &line,
const CGAL_Plane_3<R> &plane);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Plane_3<R> & p,
const CGAL_Line_3<R> & line)
{
return CGAL_squared_distance(line, p);
}
@}
@C@<Implementation@>
If the line is parallel to the plane, we can take the squared distance of an
arbitrary point of the line to the plane. Otherwise the distance is zero.
@$@<3D squared distance 2 outline definitions@>+=@{
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Line_3<R> &line,
const CGAL_Plane_3<R> &plane)
{
typedef typename R::FT FT;
if (contains_vector(plane, line.direction().vector() ))
return CGAL_squared_distance(plane, line.point());
return (R_FT_return(R))(FT(0));
}
@}
@B@<3D Plane Ray distance@>
@C@<Declarations@>
@$@<3D squared distance 2 header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> &ray,
const CGAL_Plane_3<R> &plane);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Plane_3<R> & plane,
const CGAL_Ray_3<R> & ray)
{
return CGAL_squared_distance(ray, plane);
}
@}
@C@<Implementation@>
@$@<3D squared distance 2 outline definitions@>+=@{
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Ray_3<R> &ray,
const CGAL_Plane_3<R> &plane)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const CGAL_Point_3<R> &start = ray.start();
// const CGAL_Vector_3<R> &end = ray.direction().vector();
const CGAL_Point_3<R> &planepoint = plane.point();
CGAL_Vector_3<R> start_min_pp = start - planepoint;
CGAL_Vector_3<R> end_min_pp = ray.direction().vector();
const CGAL_Vector_3<R> &normal = plane.orthogonal_vector();
RT sdm_rs2pp = CGAL_wdot(normal, start_min_pp);
RT sdm_re2pp = CGAL_wdot(normal, end_min_pp);
switch (CGAL_sign(sdm_rs2pp)) {
case -1:
if (sdm_re2pp > RT(0))
return (R_FT_return(R))(FT(0));
return CGAL_squared_distance_to_plane(normal, start_min_pp);
case 0:
default:
return (R_FT_return(R))(FT(0));
case 1:
if (sdm_re2pp < RT(0))
return (R_FT_return(R))(FT(0));
return CGAL_squared_distance_to_plane(normal, start_min_pp);
}
}
@}
@B@<3D Plane Segment distance@>
@C@<Declarations@>
@$@<3D squared distance 2 header declarations@>+=@{@-
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg,
const CGAL_Plane_3<R> &plane);
template <class R>
inline R_FT_return(R)
CGAL_squared_distance(
const CGAL_Plane_3<R> & plane,
const CGAL_Segment_3<R> & seg)
{
return CGAL_squared_distance(seg, plane);
}
@}
@C@<Implementation@>
@$@<3D squared distance 2 outline definitions@>+=@{
template <class R>
extern R_FT_return(R)
CGAL_squared_distance(
const CGAL_Segment_3<R> &seg,
const CGAL_Plane_3<R> &plane)
{
typedef typename R::RT RT;
typedef typename R::FT FT;
const CGAL_Point_3<R> &start = seg.start();
const CGAL_Point_3<R> &end = seg.end();
if (start == end)
return CGAL_squared_distance(start, plane);
const CGAL_Point_3<R> &planepoint = plane.point();
CGAL_Vector_3<R> start_min_pp = start - planepoint;
CGAL_Vector_3<R> end_min_pp = end - planepoint;
const CGAL_Vector_3<R> &normal = plane.orthogonal_vector();
RT sdm_ss2pp = CGAL_wdot(normal, start_min_pp);
RT sdm_se2pp = CGAL_wdot(normal, end_min_pp);
switch (CGAL_sign(sdm_ss2pp)) {
case -1:
if (sdm_se2pp >= RT(0))
return (R_FT_return(R))(FT(0));
if (sdm_ss2pp >= sdm_se2pp)
return CGAL_squared_distance_to_plane(normal, start_min_pp);
else
return CGAL_squared_distance_to_plane(normal, end_min_pp);
case 0:
default:
return (R_FT_return(R))(FT(0));
case 1:
if (sdm_se2pp <= RT(0))
return (R_FT_return(R))(FT(0));
if (sdm_ss2pp <= sdm_se2pp)
return CGAL_squared_distance_to_plane(normal, start_min_pp);
else
return CGAL_squared_distance_to_plane(normal, end_min_pp);
}
}
@}
@B@<The header files@>
Here we collect all code in the appropriate header files.
@O@<../include/CGAL/squared_distance_3_0.h@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3_0.h@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_DISTANCE_3_0_H
#define CGAL_DISTANCE_3_0_H
#ifndef CGAL_POINT_3_H
#include <CGAL/Point_3.h>
#endif // CGAL_POINT_3_H
@<3D squared distance 0 header declarations@>
#ifdef CGAL_CFG_NO_AUTOMATIC_TEMPLATE_INCLUSION
#include <CGAL/squared_distance_3_0.C>
#endif
#endif
@}
@O@<../include/CGAL/squared_distance_3_0.C@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3_0.C@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_UTILS_H
#include <CGAL/utils.h>
#endif // CGAL_UTILS_H
#ifndef CGAL_ENUM_H
#include <CGAL/enum.h>
#endif // CGAL_ENUM_H
#ifndef CGAL_WMULT_H
#include <CGAL/wmult.h>
#endif // CGAL_WMULT_H
@<3D squared distance general utilities@>
@<3D squared distance 0 outline definitions@>
@}
@O@<../include/CGAL/squared_distance_3_1.h@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3_1.h@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_DISTANCE_3_1_H
#define CGAL_DISTANCE_3_1_H
#ifndef CGAL_SEGMENT_3_H
#include <CGAL/Segment_3.h>
#endif // CGAL_SEGMENT_3_H
#ifndef CGAL_LINE_3_H
#include <CGAL/Line_3.h>
#endif // CGAL_LINE_3_H
#ifndef CGAL_RAY_3_H
#include <CGAL/Ray_3.h>
#endif // CGAL_RAY_3_H
@<3D squared distance header declarations@>
#ifdef CGAL_CFG_NO_AUTOMATIC_TEMPLATE_INCLUSION
#include <CGAL/squared_distance_3_1.C>
#endif
#endif
@}
@O@<../include/CGAL/squared_distance_3_1.C@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3_1.C@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_UTILS_H
#include <CGAL/utils.h>
#endif // CGAL_UTILS_H
#ifndef CGAL_POINT_3_H
#include <CGAL/Point_3.h>
#endif // CGAL_POINT_3_H
#ifndef CGAL_ENUM_H
#include <CGAL/enum.h>
#endif // CGAL_ENUM_H
#ifndef CGAL_WMULT_H
#include <CGAL/wmult.h>
#endif // CGAL_WMULT_H
#ifndef CGAL_SQUARED_DISTANCE_3_0_H
#include <CGAL/squared_distance_3_0.h>
#endif // CGAL_SQUARED_DISTANCE_3_0_H
@<3D squared distance outline definitions@>
@}
@O@<../include/CGAL/squared_distance_3_2.h@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3_2.h@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_DISTANCE_3_2_H
#define CGAL_DISTANCE_3_2_H
#ifndef CGAL_SEGMENT_3_H
#include <CGAL/Segment_3.h>
#endif // CGAL_SEGMENT_3_H
#ifndef CGAL_LINE_3_H
#include <CGAL/Line_3.h>
#endif // CGAL_LINE_3_H
#ifndef CGAL_RAY_3_H
#include <CGAL/Ray_3.h>
#endif // CGAL_RAY_3_H
#ifndef CGAL_PLANE_3_H
#include <CGAL/Plane_3.h>
#endif // CGAL_PLANE_3_H
@<3D squared distance 2 header declarations@>
#ifdef CGAL_CFG_NO_AUTOMATIC_TEMPLATE_INCLUSION
#include <CGAL/squared_distance_3_2.C>
#endif
#endif
@}
@O@<../include/CGAL/squared_distance_3_2.C@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3_2.C@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_UTILS_H
#include <CGAL/utils.h>
#endif // CGAL_UTILS_H
#ifndef CGAL_POINT_3_H
#include <CGAL/Point_3.h>
#endif // CGAL_POINT_3_H
#ifndef CGAL_PLANE_3_H
#include <CGAL/Plane_3.h>
#endif // CGAL_PLANE_3_H
#ifndef CGAL_ENUM_H
#include <CGAL/enum.h>
#endif // CGAL_ENUM_H
#ifndef CGAL_WMULT_H
#include <CGAL/wmult.h>
#endif // CGAL_WMULT_H
#ifndef CGAL_SQUARED_DISTANCE_3_0_H
#include <CGAL/squared_distance_3_0.h>
#endif // CGAL_SQUARED_DISTANCE_3_0_H
@<3D squared distance general utilities 2@>
@<3D squared distance 2 outline definitions@>
@}
@O@<../include/CGAL/squared_distance_3.h@>==@{@-
@<CGAL_heading@>@(@-
include/CGAL/squared_distance_3.h@,@-
sqdistance_3.fw@,@-
Geert-Jan Giezeman@,@-
Saarbruecken@)
#ifndef CGAL_DISTANCE_3_H
#define CGAL_DISTANCE_3_H
#ifndef CGAL_SQUARED_DISTANCE_3_0_H
#include <CGAL/squared_distance_3_0.h>
#endif // CGAL_SQUARED_DISTANCE_3_0_H
#ifndef CGAL_SQUARED_DISTANCE_3_1_H
#include <CGAL/squared_distance_3_1.h>
#endif // CGAL_SQUARED_DISTANCE_3_1_H
#ifndef CGAL_SQUARED_DISTANCE_3_2_H
#include <CGAL/squared_distance_3_2.h>
#endif // CGAL_SQUARED_DISTANCE_3_2_H
#endif
@}