Geometry.h 17.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425
//============================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//
//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
//  Copyright 2014 UT-Battelle, LLC.
//  Copyright 2014 Los Alamos National Security.
//
//  Under the terms of Contract DE-NA0003525 with NTESS,
//  the U.S. Government retains certain rights in this software.
//
//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
//  Laboratory (LANL), the U.S. Government retains certain rights in
//  this software.
//============================================================================

#ifndef vtk_m_Geometry_h
#define vtk_m_Geometry_h

#include <vtkm/VectorAnalysis.h>

namespace vtkm
{

// Forward declarations of geometric types:
template <typename CoordType, int Dim, bool IsTwoSided>
struct Ray;
template <typename CoordType, int Dim>
struct LineSegment;
template <typename CoordType>
struct Plane;
template <typename CoordType, int Dim>
struct Sphere;

/// Represent an infinite or semi-infinite line segment with a point and a direction
///
/// The \a IsTwoSided template parameter indicates whether the class represents
/// an infinite line extending in both directions from the base point or a
/// semi-infinite ray extending only in the positive direction from the base points.
template <typename CoordType = vtkm::FloatDefault, int Dim = 3, bool IsTwoSided = false>
struct Ray
{
  static constexpr int Dimension = Dim;
  using Vector = vtkm::Vec<CoordType, Dim>;
  static constexpr bool TwoSided = IsTwoSided;

  Vector Origin;
  Vector Direction; // always stored as a unit length.

  /// Construct a default 2-D ray, from (0,0) pointing along the +x axis.
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
  VTKM_EXEC_CONT Ray();

  /// Construct a default 3-D ray from (0,0,0) pointing along the +x axis.
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 3, int>::type = 0>
  VTKM_EXEC_CONT Ray();

  /// Construct a ray from a point and direction.
  VTKM_EXEC_CONT
  Ray(const Vector& point, const Vector& direction);

  /// Construct a ray from a line segment.
  VTKM_EXEC_CONT
  Ray(const LineSegment<CoordType, Dim>& segment);

  /// Return whether the ray is valid or not.
  ///
  /// It is possible for an invalid direction (zero length)
  /// to be passed to the constructor. When this happens,
  /// the constructor divides by zero, leaving Inf in all
  /// components.
  VTKM_EXEC_CONT
  bool IsValid() const;

  /// Compute a point along the line. \a param values > 0 lie on the ray.
  VTKM_EXEC_CONT
  Vector Evaluate(CoordType param) const;

  /// Return the minmum distance from \a point to this line/ray.
  ///
  /// Note that when the direction has zero length, this simplifies
  /// the distance between \a point and the ray's origin.
  /// Otherwise, the distance returned is either the perpendicular
  /// distance from \a point to the line or the distance
  /// to the ray's origin.
  VTKM_EXEC_CONT
  CoordType DistanceTo(const Vector& point) const;

  /// Return the minimum distance between the ray/line and \a point.
  VTKM_EXEC_CONT
  CoordType DistanceTo(const Vector& point, CoordType& param, Vector& projectedPoint) const;

  /// Compute the non-degenerate point where two 2-D rays intersect, or return false.
  ///
  /// If true is returned, then the rays intersect in a unique point and
  /// \a point is set to that location.
  ///
  /// If false is returned, then either
  /// (1) the rays are parallel and may be
  /// coincident (intersect everywhere) or offset (intersect nowhere); or
  /// (2) the lines intersect but not the rays (because the the intersection
  /// occurs in the negative parameter space of one or both rays).
  /// In the latter case (2), the \a point is still set to the location of the
  /// intersection.
  ///
  /// The tolerance \a tol is the minimum acceptable denominator used to
  /// compute the intersection point coordinates and thus dictates the
  /// maximum distance from the segments at which intersections will be
  /// reported as valid.
  template <bool OtherTwoSided, int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
  VTKM_EXEC_CONT bool Intersect(const Ray<CoordType, Dim, OtherTwoSided>& other,
                                Vector& point,
                                CoordType tol = 0.f);
};

/// Represent a finite line segment with a pair of points
template <typename CoordType = vtkm::FloatDefault, int Dim = 3>
struct LineSegment
{
  static constexpr int Dimension = Dim;
  using Vector = vtkm::Vec<CoordType, Dim>;
  Vector Endpoints[2];

  /// Construct a default segment from (0,0) to (1,0).
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
  VTKM_EXEC_CONT LineSegment();

  /// Construct a default segment from (0,0,0) to (1,0,0).
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 3, int>::type = 0>
  VTKM_EXEC_CONT LineSegment();

  /// Construct a segment spanning points \a p0 and \a p1.
  VTKM_EXEC_CONT
  LineSegment(const Vector& p0, const Vector& p1);

  /// Return whether this line segment has an infinitesimal extent (i.e., whether the endpoints are coincident)
  VTKM_EXEC_CONT
  bool IsSingular(CoordType tol2 = static_cast<CoordType>(1.0e-6f)) const;

  /// Construct a plane bisecting this line segment (only when Dimension is 3).
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 3, int>::type = 0>
  VTKM_EXEC_CONT Plane<CoordType> PerpendicularBisector() const;

  /// Construct a perpendicular bisector to this line segment (only when Dimension is 2).
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
  VTKM_EXEC_CONT Ray<CoordType, Dim, true> PerpendicularBisector() const;

  /// Return the midpoint of the line segment
  VTKM_EXEC_CONT
  Vector Center() const { return this->Evaluate(0.5f); }

  /// Return the vector pointing to endpoint 1 from endpoint 0.
  /// This vector is not of unit length and in the case of
  /// degenerate lines, may have zero length.
  /// Call vtkm::Normal() on the return value if you want a normalized result.
  VTKM_EXEC_CONT
  Vector Direction() const { return this->Endpoints[1] - this->Endpoints[0]; }

  /// Compute a point along the line. \a param values in [0,1] lie on the line segment.
  VTKM_EXEC_CONT
  Vector Evaluate(CoordType param) const;

  /// Return the minmum distance from \a point to this line segment.
  ///
  /// Note that when the endpoints are coincident, this simplifies
  /// the distance between \a point and either endpoint.
  /// Otherwise, the distance returned is either the perpendicular
  /// distance from \a point to the line segment or the distance
  /// to the nearest endpoint (whichever is smaller).
  VTKM_EXEC_CONT
  CoordType DistanceTo(const Vector& point) const;

  /// Return the minimum distance between the line segment and \a point.
  VTKM_EXEC_CONT
  CoordType DistanceTo(const Vector& point, CoordType& param, Vector& projectedPoint) const;

  /// Compute the non-degenerate point where two (infinite) 2-D line segments intersect, or return false.
  ///
  /// If true is returned, then the lines intersect in a unique point and
  /// \a point is set to that location.
  ///
  /// If false is returned, then the lines are parallel and either they are
  /// coincident (intersect everywhere) or offset (intersect nowhere).
  ///
  /// The tolerance \a tol is the minimum acceptable denominator used to
  /// compute the intersection point coordinates and thus dictates the
  /// maximum distance from the segments at which intersections will be
  /// reported as valid.
  template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
  VTKM_EXEC_CONT bool IntersectInfinite(const LineSegment<CoordType, Dim>& other,
                                        Vector& point,
                                        CoordType tol = 0.f);
};

/// Represent a plane with a base point (origin) and normal vector.
template <typename CoordType = vtkm::FloatDefault>
struct Plane
{
  using Vector = vtkm::Vec<CoordType, 3>;
  Vector Origin;
  Vector Normal;

  /// Construct a default plane whose base point is the origin and whose normal is (0,0,1)
  VTKM_EXEC_CONT
  Plane();

  /// Construct a plane with the given \a origin and \a normal.
  VTKM_EXEC_CONT
  Plane(const Vector& origin, const Vector& normal, CoordType tol2 = static_cast<CoordType>(1e-8f));

  /// Return true if the plane's normal is well-defined to within the given tolerance.
  VTKM_EXEC_CONT
  bool IsValid() const { return !vtkm::IsInf(this->Normal[0]); }

  /// Return the **signed** distance from the plane to the point.
  VTKM_EXEC_CONT
  CoordType DistanceTo(const Vector& point) const;

  /// Return the closest point in the plane to the given point.
  VTKM_EXEC_CONT
  Vector ClosestPoint(const Vector& point) const;

  /// Intersect this plane with the ray (or line if the ray is two-sided).
  ///
  /// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection).
  /// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the
  /// entire ray/line lies in the plane).
  /// In the latter case, \a lineInPlane will be true upon exit.
  ///
  /// If this method returns true, then \a parameter will be set to a number indicating
  /// where along the ray/line  the plane hits and \a point will be set to that location.
  /// If the input is a ray, the \a parameter will be non-negative.
  template <bool IsTwoSided>
  VTKM_EXEC_CONT bool Intersect(const Ray<CoordType, 3, IsTwoSided>& ray,
                                CoordType& parameter,
                                Vector& point,
                                bool& lineInPlane,
                                CoordType tol = CoordType(1e-6f)) const;

  /// Intersect this plane with the line \a segment.
  ///
  /// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection).
  /// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the
  /// entire line segment lies in the plane).
  /// In the latter case, \a lineInPlane will be true upon exit.
  ///
  /// If this method returns true, then \a parameter will be set to a number in [0,1] indicating
  /// where along the line segment the plane hits.
  VTKM_EXEC_CONT
  bool Intersect(const LineSegment<CoordType>& segment,
                 CoordType& parameter,
                 bool& lineInPlane) const;

  /// Intersect this plane with the line \a segment.
  ///
  /// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection).
  /// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the
  /// entire line segment lines in the plane).
  /// In the latter case, \a lineInPlane will be true upon exit.
  ///
  /// If this method returns true, then \a parameter will be set to a number in [0,1] indicating
  /// where along the line segment the plane hits and \a point will be set to that location.
  VTKM_EXEC_CONT
  bool Intersect(const LineSegment<CoordType>& segment,
                 CoordType& parameter,
                 Vector& point,
                 bool& lineInPlane) const;

  /// Intersect this plane with another plane.
  ///
  /// Returns true if there is a non-degenrate intersection (i.e., a line of intersection).
  /// Returns false if there is no intersection *or* if the intersection is degenerate
  /// (i.e., the planes are coincident).
  /// In the latter case, \a coincident will be true upon exit and \a segment will
  /// unmodified.
  ///
  /// If this method returns true, then the resulting \a segment will have its
  /// base point on the line of intersection and its second point will be a unit
  /// length away in the direction of the cross produce of the input plane normals
  /// (this plane crossed with the \a other).
  ///
  /// The tolerance \a tol is the minimum squared length of the cross-product
  /// of the two plane normals. It is also compared to the squared distance of
  /// the base point of \a other away from \a this plane when considering whether
  /// the planes are coincident.
  VTKM_EXEC_CONT
  bool Intersect(const Plane<CoordType>& other,
                 Ray<CoordType, 3, true>& ray,
                 bool& coincident,
                 CoordType tol2 = static_cast<CoordType>(1e-6f)) const;
};

/// Represent a sphere of the given \a Dimension.
/// If a constructor is given an invalid specification, then
/// the Radius of the resulting sphere will be -1.
template <typename CoordType = vtkm::FloatDefault, int Dim = 3>
struct Sphere
{
  static constexpr int Dimension = Dim;
  using Vector = vtkm::Vec<CoordType, Dim>;
  Vector Center;
  CoordType Radius;

  /// Construct a default sphere (unit radius at the origin).
  VTKM_EXEC_CONT Sphere();

  /// Construct a sphere from a center point and radius.
  VTKM_EXEC_CONT Sphere(const Vector& center, CoordType radius);

  /// Return true if the sphere is valid (i.e., has a strictly positive radius).
  VTKM_EXEC_CONT
  bool IsValid() const { return this->Radius > 0.f; }

  /// Return whether the point lies strictly inside the sphere.
  VTKM_EXEC_CONT
  bool Contains(const Vector& point, CoordType tol2 = 0.f) const;

  /// Classify a point as inside (-1), on (0), or outside (+1) of the sphere.
  ///
  /// The tolerance \a tol2 is the maximum allowable difference in squared
  /// magnitude between the squared radius and the squared distance between
  /// the \a point and Center.
  VTKM_EXEC_CONT
  int Classify(const Vector& point, CoordType tol2 = 0.f) const;
};

// -----------------------------------------------------------------------------
// Synonyms
//
// These "using" statements aim to make it easier to use the templated
// structs above when working with a particular dimension and/or the
// default floating-point type.

/// Lines are two-sided rays:
template <typename CoordType, int Dim = 3>
using Line = Ray<CoordType, Dim, true>;

// Shortcuts for 2D and 3D rays, lines, and line segments:
template <typename CoordType>
using Ray2 = Ray<CoordType, 2>;
template <typename CoordType>
using Ray3 = Ray<CoordType, 3>;
template <typename CoordType>
using Line2 = Line<CoordType, 2>;
template <typename CoordType>
using Line3 = Line<CoordType, 3>;
template <typename CoordType>
using LineSegment2 = LineSegment<CoordType, 2>;
template <typename CoordType>
using LineSegment3 = LineSegment<CoordType, 3>;

/// Circle is an alias for a 2-Dimensional sphere.
template <typename T>
using Circle = Sphere<T, 2>;

// Aliases for d-dimensional spheres.
template <typename T>
using Sphere2 = Sphere<T, 2>;
template <typename T>
using Sphere3 = Sphere<T, 3>;

// Shortcuts for default floating-point types
using Ray2d = Ray2<vtkm::FloatDefault>;
using Ray3d = Ray3<vtkm::FloatDefault>;
using Line2d = Line2<vtkm::FloatDefault>;
using Line3d = Line3<vtkm::FloatDefault>;
using LineSegment2d = LineSegment2<vtkm::FloatDefault>;
using LineSegment3d = LineSegment3<vtkm::FloatDefault>;
using Plane3d = Plane<vtkm::FloatDefault>;
using Circle2d = Circle<vtkm::FloatDefault>;
using Sphere2d = Sphere2<vtkm::FloatDefault>;
using Sphere3d = Sphere3<vtkm::FloatDefault>;

// -----------------------------------------------------------------------------
// Construction techniques
//
// These are free functions that create instances of geometric structs by taking
// in data that is not identical to the state of the struct and converting it
// into state for the struct.

/// Construct a plane from a point plus one of: a line, a ray, or a line segment.
///
/// The plane returned will contain the point and the line/ray/segment.
/// The plane normal will be the cross product of the line/ray/segment's direction
/// and the vector from the line/ray/segment's origin to the given \a point.
/// If the \a point is collinear with the line/ray/line-segment, an invalid
/// plane will be returned.
template <typename CoordType, bool IsTwoSided>
VTKM_EXEC_CONT vtkm::Plane<CoordType> make_PlaneFromPointAndLine(
  const vtkm::Vec<CoordType, 3>& point,
  const vtkm::Ray<CoordType, 3, IsTwoSided>& ray,
  CoordType tol2 = static_cast<CoordType>(1e-8f));

template <typename CoordType>
VTKM_EXEC_CONT vtkm::Plane<CoordType> make_PlaneFromPointAndLineSegment(
  const vtkm::Vec<CoordType, 3>& point,
  const vtkm::LineSegment3<CoordType>& segment,
  CoordType tol2 = static_cast<CoordType>(1e-8f));

/// Construct a circle from 3 points.
template <typename CoordType>
VTKM_EXEC_CONT vtkm::Circle<CoordType> make_CircleFrom3Points(
  const typename vtkm::Vec<CoordType, 2>& p0,
  const typename vtkm::Vec<CoordType, 2>& p1,
  const typename vtkm::Vec<CoordType, 2>& p2,
  CoordType tol = static_cast<CoordType>(1e-6f));

/// Construct a sphere from 4 points.
template <typename CoordType>
VTKM_EXEC_CONT vtkm::Sphere<CoordType, 3> make_SphereFrom4Points(
  const vtkm::Vec<CoordType, 3>& a0,
  const vtkm::Vec<CoordType, 3>& a1,
  const vtkm::Vec<CoordType, 3>& a2,
  const vtkm::Vec<CoordType, 3>& a3,
  CoordType tol = static_cast<CoordType>(1e-6f));

} // namespace vtkm

#include "vtkm/Geometry.hxx"

#endif // vtk_m_Geometry_h