Go back to Richel Bilderbeek's homepage.

Go back to Richel Bilderbeek's C++ page.

 

 

 

 

 

(C++) Geometry

 

STLQt CreatorLubuntu

 

CppGeometry is a folder with my geometry functions.

 

For the Boost.Geometry library, go to the page about Boost.Geometry.

 

 

 

 

 

CppGeometry code snippets

 

 

 

 

 

 

Technical facts

 

 

 

 

 

 

./CppGeometry/CppGeometry.pri

 

INCLUDEPATH += \
    ../../Classes/CppGeometry

SOURCES += \
    ../../Classes/CppGeometry/geometry.cpp \
    ../../Classes/CppGeometry/geometry_test.cpp \
    ../../Classes/CppGeometry/geometry_calccrossproduct.cpp \
    ../../Classes/CppGeometry/geometry_calccenter.cpp \
    ../../Classes/CppGeometry/geometry_getangle.cpp \
    ../../Classes/CppGeometry/geometry_is_clockwise.cpp \
    ../../Classes/CppGeometry/geometry_is_counter_clockwise.cpp

HEADERS  += \
    ../../Classes/CppGeometry/geometry.h

OTHER_FILES += \
    ../../Classes/CppGeometry/Licence.txt

 

 

 

 

 

./CppGeometry/geometry.h

 

//---------------------------------------------------------------------------
/*
Geometry, class with geometry functions
Copyright (C) 2014-2015 Richel Bilderbeek

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppGeometry.htm
//---------------------------------------------------------------------------
#ifndef RIBI_GEOMETRY_H
#define RIBI_GEOMETRY_H

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
#pragma GCC diagnostic ignored "-Wunused-variable"
#include <cassert>
#include <functional>
#include <string>
#include <tuple>
#include <vector>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#ifndef _WIN32
#include <boost/geometry/geometries/polygon.hpp>
#endif

#include "apfloat.h"

struct QPen;
struct QRect;
struct QRectF;
struct QPoint;
struct QPointF;
#pragma GCC diagnostic pop

namespace ribi {

///Geometry functions, working with Boost.Geometry
struct Geometry
{
  typedef boost::geometry::model::d2::point_xy<double> Coordinat2D;
  typedef boost::geometry::model::d2::point_xy<apfloat> ApCoordinat2D;
  typedef boost::geometry::model::point<double,3,boost::geometry::cs::cartesian> Coordinat3D;
  typedef boost::geometry::model::point<apfloat,3,boost::geometry::cs::cartesian> ApCoordinat3D;

  typedef boost::geometry::model::linestring<Coordinat2D> Linestring;
  typedef boost::geometry::model::polygon<Coordinat2D> Polygon;
  typedef boost::geometry::model::box<Coordinat2D> Rect;

  typedef apfloat Apfloat;
  typedef double Double;
  typedef std::vector<ApCoordinat2D> ApCoordinats2D;
  typedef std::vector<ApCoordinat3D> ApCoordinats3D;
  typedef std::vector<Apfloat> Apfloats;
  typedef std::vector<Coordinat2D> Coordinats2D;
  typedef std::vector<Coordinat3D> Coordinats3D;
  typedef std::vector<Double> Doubles;
  typedef std::vector<Linestring> Linestrings;
  typedef std::vector<Polygon> Polygons;

  typedef std::pair<Polygons,Linestrings> Shapes;

  Geometry();

  ///apfloat its atan2 does not correspond to std::atan2
  ///Atan2 makes it correspond
  Apfloat Atan2(const Apfloat& dx, const Apfloat& b) const noexcept;

  Coordinat2D CalcCenter(const Coordinats2D& v) const noexcept;
  Coordinat3D CalcCenter(const Coordinats3D& v) const noexcept;
  ApCoordinat2D CalcCenter(const ApCoordinats2D& v) const noexcept;
  ApCoordinat3D CalcCenter(const ApCoordinats3D& v) const noexcept;

  Coordinat3D CalcCrossProduct(const Coordinat3D& a,const Coordinat3D& b) const noexcept;
  ApCoordinat3D CalcCrossProduct(const ApCoordinat3D& a,const ApCoordinat3D& b) const noexcept;


  double CalcDotProduct(
    const Coordinat3D& a,
    const Coordinat3D& b
  ) const noexcept;

  Apfloat CalcDotProduct(
    const ApCoordinat3D& a,
    const ApCoordinat3D& b
  ) const noexcept;

  Coordinat3D CalcNormal(
    const Coordinat3D& a,
    const Coordinat3D& b,
    const Coordinat3D& c
  ) const noexcept;

  ApCoordinat3D CalcNormal(
    const ApCoordinat3D& a,
    const ApCoordinat3D& b,
    const ApCoordinat3D& c
  ) const noexcept;

  ///Determine the plane that goes through these three points
  ///Returns a std::vector of size 4
  ///CalcPlane return the coefficients in the following form:
  /// A.x + B.y + C.z = D
  ///Converting this to z being a function of x and y:
  /// -C.z = A.x + B.y - D
  /// z = -A/C.x - B/C.y + D/C
  std::vector<double> CalcPlane(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  ) const noexcept;

  std::vector<apfloat> CalcPlane(
    const ApCoordinat3D& p1,
    const ApCoordinat3D& p2,
    const ApCoordinat3D& p3
  ) const noexcept;

    Coordinats2D CalcProjection(const   Coordinats3D& v) const;
  ApCoordinats2D CalcProjection(const ApCoordinats3D& v) const;

  Coordinat2D Coordinat2DToBoostGeometryPointXy(
    const Coordinat2D& c
  ) const noexcept;

  std::vector<Coordinat2D> Coordinats2DToBoostGeometryPointsXy(
    const std::vector<Coordinat2D>& v
  ) const noexcept;

  Coordinat3D CreatePoint(
    const double x,
    const double y,
    const double z
  ) const noexcept;

  Rect CreateRect(
    const double left,
    const double top,
    const double width,
    const double height
  ) const noexcept;

  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
    CreateHouseShape() const noexcept;

  Polygon CreateShapeHeart(const double scale = 1.0) const noexcept;
  Polygon CreateShapeHouse(const double scale = 1.0) const noexcept;
  Polygon CreateShapePolygon(const int n, const double rotation = 0.0, const double scale = 0.0) const noexcept;
  Polygon CreateShapeTriangle(const double scale = 1.0) const noexcept;

  ///Functor for X-Y-Z ordering
  std::function<bool(const Coordinat2D& lhs, const Coordinat2D& rhs)> Equals2d() const noexcept;
  std::function<bool(const Coordinat3D& lhs, const Coordinat3D& rhs)> Equals() const noexcept;
  std::function<bool(const ApCoordinat3D& lhs, const ApCoordinat3D& rhs)> Equals3d() const noexcept;

  ///Take the floating point modulus
  double Fmod(const double x, const double mod) const noexcept;
  Apfloat Fmod(const Apfloat& x, const Apfloat& mod) const noexcept;

  ///Obtain the angle in radians between two deltas, using a screen geometry and following a clock
  ///12 o'clock is 0.0 * pi
  /// 3 o'clock is 0.5 * pi
  /// 6 o'clock is 1.0 * pi
  /// 9 o'clock is 1.5 * pi

  /*
   Y          Y
   |    (11)  |  (1)
-2|          |
   |          |
-1| (10)     |      (2)
   |          |
  0+----------0--------X
   |          |
+1| (8)      |      (4)
   |          |
+2|          |
   |     (7)  |  (5)
   +----------+--------X
        - - -   + + +
        3 2 1 0 1 2 3

  Approximate coordinat for a point for every hour, with the approximate angle
   1: ( 1,-2) :  1/6 * pi
   2: ( 2,-1) :  2/6 * pi
   3: ( 3, 0) :  3/6 * pi
   4: ( 2, 1) :  4/6 * pi
   5: ( 1, 2) :  5/6 * pi
   6: ( 0, 3) :  6/6 * pi
   7: (-1, 2) :  7/6 * pi
   8: (-2, 1) :  8/6 * pi
   9: (-3, 0) :  9/6 * pi
  10: (-2,-1) : 10/6 * pi
  11: (-1,-2) : 11/6 * pi
  12: ( 0,-3) : 12/6 * pi

  */
  //From www.richelbilderbeek.nl/CppGetAngle.htm
  double GetAngleClockScreen(const double dx, const double dy) const noexcept; //The original GetAngle
  Apfloat GetAngleClockScreen(const Apfloat& dx, const Apfloat& dy) const noexcept;
  double GetAngleClockScreen(const Coordinat2D& p) const noexcept;
  Apfloat GetAngleClockScreen(const ApCoordinat2D& p) const noexcept;

  ///Obtain the angle in radians between two deltas, using a Cartesian plane and following a clock
  ///12 o'clock is 0.0 * pi
  /// 3 o'clock is 0.5 * pi
  /// 6 o'clock is 1.0 * pi
  /// 9 o'clock is 1.5 * pi

  /*
   Y          Y
   |    (11)  |  (1)
+2|          |
   |          |
+1| (10)     |      (2)
   |          |
  0+----------0--------X
   |          |
-1| (8)      |      (4)
   |          |
-2|          |
   |     (7)  |  (5)
   +----------+--------X
        - - -   + + +
        3 2 1 0 1 2 3

  Approximate coordinat for a point for every hour, with the approximate angle
   1: ( 1, 2) :  1/6 * pi
   2: ( 2, 1) :  2/6 * pi
   3: ( 3, 0) :  3/6 * pi
   4: ( 2,-1) :  4/6 * pi
   5: ( 1,-2) :  5/6 * pi
   6: ( 0,-3) :  6/6 * pi
   7: (-1,-2) :  7/6 * pi
   8: (-2,-1) :  8/6 * pi
   9: (-3, 0) :  9/6 * pi
  10: (-2, 1) : 10/6 * pi
  11: (-1, 2) : 11/6 * pi
  12: ( 0, 3) : 12/6 * pi

  */
  //From www.richelbilderbeek.nl/CppGetAngle.htm
  double GetAngleClockCartesian(const double dx, const double dy) const noexcept;
  Apfloat GetAngleClockCartesian(const Apfloat& dx, const Apfloat& dy) const noexcept;
  double GetAngleClockCartesian(const Coordinat2D& p) const noexcept;
  Apfloat GetAngleClockCartesian(const ApCoordinat2D& p) const noexcept;


  template <class T>
  T GetBottom(const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T>>& r) const noexcept
  {
    using boost::geometry::get;
    using boost::geometry::min_corner;
    using boost::geometry::max_corner;
    const auto bottom = get<max_corner,1>(r);
    #ifndef NDEBUG
    const auto top = get<min_corner,1>(r);
    assert(top <= bottom);
    #endif
    return bottom;
  }

  //From www.richelbilderbeek.nl/CppGetDistance.htm
  double GetDistance(const double dx, const double dy) const noexcept;
  double GetDistance(const Coordinat2D& a, const Coordinat2D& b) const noexcept;
  double GetDistance(const double x1, const double y1, const double x2, const double y2) const noexcept { return GetDistance(x2-x1,y2-y1); }

  template <class T>
  T GetHeight(const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T>>& r) const noexcept
  {
    using boost::geometry::get;
    using boost::geometry::min_corner;
    using boost::geometry::max_corner;
    const auto top   (get<min_corner,1>(r));
    const auto bottom(get<max_corner,1>(r));
    assert(top <= bottom);
    return bottom - top;
  }


  template <class T>
  T GetLeft(const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T>>& r) const noexcept
  {
    using boost::geometry::get;
    using boost::geometry::min_corner;
    using boost::geometry::max_corner;
    const auto left (get<min_corner,0>(r));
    #ifndef NDEBUG
    const auto right(get<max_corner,0>(r));
    assert(left <= right);
    #endif
    return left;
  }

  template <class T>
  T GetRight(const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T>>& r) const noexcept
  {
    using boost::geometry::get;
    using boost::geometry::min_corner;
    using boost::geometry::max_corner;
    const auto right(get<max_corner,0>(r));
    #ifndef NDEBUG
    const auto left (get<min_corner,0>(r));
    assert(left <= right);
    #endif
    return right;
  }

  template <class T>
  T GetTop(const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T>>& r) const noexcept
  {
    using boost::geometry::get;
    using boost::geometry::min_corner;
    using boost::geometry::max_corner;
    const auto top   (get<min_corner,1>(r));
    #ifndef NDEBUG
    const auto bottom(get<max_corner,1>(r));
    assert(top <= bottom);
    #endif
    return top;
  }

  std::string GetVersion() const noexcept;
  std::vector<std::string> GetVersionHistory() const noexcept;

  template <class T>
  T GetWidth(const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T>>& r) const noexcept
  {
    using boost::geometry::get;
    using boost::geometry::min_corner;
    using boost::geometry::max_corner;
    const auto left (get<min_corner,0>(r));
    const auto right(get<max_corner,0>(r));
    assert(left <= right);
    return right - left;
  }
  ///Are two angles ordered clockwise
  ///12 o'clock is 0.0 * pi
  /// 3 o'clock is 0.5 * pi
  /// 6 o'clock is 1.0 * pi
  /// 9 o'clock is 1.5 * pi
  ///
  ///Yes: 0.0 * pi and 0.5 * pi
  ///No : 0.5 * pi and 0.0 * pi
  bool IsClockwise(const double angle_a_radians, const double angle_b_radians) const noexcept;
  bool IsClockwise(const Apfloat& angle_a_radians, const Apfloat& angle_b_radians) const noexcept;
  bool IsClockwise(const std::vector<double>& angles_radians) const noexcept;
  bool IsClockwise(const Apfloats& angles_radians) const noexcept;

  ///Are the points ordered clockwise in the XYZ plane, seen from the observer?
  /*

  Uses the following geometry

                |
  (-2, 1) = D   +   A = (2,1)
                |
           -+-+-O-+-+-
                |
  (-2,-1) = C   +   B = (2,-1)
                |

  */
  bool IsClockwiseCartesian(const   Coordinats3D& points, const   Coordinat3D& observer) const noexcept;
  bool IsClockwiseCartesian(const ApCoordinats3D& points, const ApCoordinat3D& observer) const noexcept;
  ///Are the points ordered clockwise in the XY plane seen from above
  /// (e.g. from coordinat {0,0,1} )
  //bool IsClockwiseHorizontal(const std::vector<Coordinat3D>& points) const noexcept;
  //bool IsClockwiseHorizontal(const std::vector<Coordinat2D>& points) const noexcept;
  bool IsClockwiseCartesianHorizontal(const std::vector<  Coordinat2D>& v) const noexcept;
  bool IsClockwiseCartesianHorizontal(const std::vector<ApCoordinat2D>& v) const noexcept;
  bool IsClockwiseCartesianHorizontal(const std::vector<  Coordinat3D>& v) const noexcept;
  bool IsClockwiseCartesianHorizontal(const std::vector<ApCoordinat3D>& v) const noexcept;

  ///Are two angles ordered counter-clockwise
  ///12 o'clock is 0.0 * pi
  /// 3 o'clock is 0.5 * pi
  /// 6 o'clock is 1.0 * pi
  /// 9 o'clock is 1.5 * pi
  ///
  ///Yes: 0.5 * pi and 0.0 * pi
  ///No : 0.0 * pi and 0.5 * pi
  bool IsCounterClockwise(const double angle_a_radians, const double angle_b_radians) const noexcept;
  bool IsCounterClockwise(const Apfloat& angle_a_radians, const Apfloat& angle_b_radians) const noexcept;
  bool IsCounterClockwise(const Doubles& angles_radians) const noexcept;
  bool IsCounterClockwise(const Apfloats& angles_radians) const noexcept;

  ///Are the points ordered counterclockwise in the XYZ plane, seen from the observer?
  /*

  Uses the following geometry

                |
  (-2, 1) = D   +   A = (2,1)
                |
           -+-+-O-+-+-
                |
  (-2,-1) = C   +   B = (2,-1)
                |

  */
  bool IsCounterClockwiseCartesian(const   Coordinats3D& points, const   Coordinat3D& observer) const noexcept;
  bool IsCounterClockwiseCartesian(const ApCoordinats3D& points, const ApCoordinat3D& observer) const noexcept;

  ///Are the points ordered counter-clockwise in the XY plane seen from above
  /// (e.g. from coordinat {0,0,1} )
  //bool IsCounterClockwiseHorizontal(const std::vector<Coordinat3D>& points) const noexcept;
  bool IsCounterClockwiseCartesianHorizontal(const std::vector<  Coordinat2D>& v) const noexcept;
  bool IsCounterClockwiseCartesianHorizontal(const std::vector<ApCoordinat2D>& v) const noexcept;
  bool IsCounterClockwiseCartesianHorizontal(const std::vector<  Coordinat3D>& v) const noexcept;
  bool IsCounterClockwiseCartesianHorizontal(const std::vector<ApCoordinat3D>& v) const noexcept;



  ///Creates a polygon from the points and checks if the polygon is convex
  /*

   A---------B      A---------B
   E        /        \   D   /
    \      /          \ / \ /
     D----C            E   C

     Convex           Concave

  */
  bool IsConvex(boost::geometry::model::polygon<Coordinat2D> polygon) const noexcept;
  bool IsConvex(const   Coordinats2D& points) const noexcept;
  bool IsConvex(const   Coordinats3D& points) const noexcept;
  bool IsConvex(const ApCoordinats3D& points) const noexcept;

  bool IsEqual2d(const Coordinat2D& lhs, const Coordinat2D& rhs) const noexcept { return Equals2d()(lhs,rhs); }
  bool IsEqual(const Coordinat3D& lhs, const Coordinat3D& rhs) const noexcept { return Equals()(lhs,rhs); }
  bool IsEqual3d(const ApCoordinat3D& lhs, const ApCoordinat3D& rhs) const noexcept { return Equals3d()(lhs,rhs); }

  ///Determines if these coordinats are in a plane
  bool IsPlane(const std::vector<Coordinat3D>& v) const noexcept;
  bool IsPlane(const std::vector<ApCoordinat3D>& v) const noexcept;

  ///Functor for X-Y ordering
  std::function<bool(const ribi::Geometry::Coordinat2D& lhs, const ribi::Geometry::Coordinat2D& rhs)> Order2dByX() const noexcept;

  ///Functor for X-Y-Z ordering
  std::function<bool(const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)> OrderByX() const noexcept;

  ///Functor for Z-Y-X ordering
  std::function<bool(const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)> OrderByZ() const noexcept;

  Linestring Rescale(
    const Linestring& polygon,
    const double scale,
    const double scale_origin_x = 0.0,
    const double scale_origin_y = 0.0
  ) const noexcept;

  Linestrings Rescale(
    const Linestrings& linestrings,
    const double scale,
    const double scale_origin_x = 0.0,
    const double scale_origin_y = 0.0
  ) const noexcept;

  Polygon Rescale(
    const Polygon& polygon,
    const double scale,
    const double scale_origin_x = 0.0,
    const double scale_origin_y = 0.0
  ) const noexcept;

  Polygons Rescale(
    const Polygons& polygons,
    const double scale,
    const double scale_origin_x = 0.0,
    const double scale_origin_y = 0.0
  ) const noexcept;

  Shapes Rescale(
    const Shapes& shapes,
    const double scale,
    const double scale_origin_x = 0.0,
    const double scale_origin_y = 0.0
  ) const noexcept;

  Apfloats ToApfloat(const Doubles& v) const noexcept;
  ApCoordinat2D ToApfloat(const Coordinat2D& p) const noexcept;
  ApCoordinat3D ToApfloat(const Coordinat3D& p) const noexcept;
  ApCoordinats2D ToApfloat(const Coordinats2D& p) const noexcept;
  ApCoordinats3D ToApfloat(const Coordinats3D& p) const noexcept;

  ///Will throw a boost::bad_lexical_cast if it fails
  double ToDouble(const apfloat& a) const;
  std::vector<double> ToDouble(const std::vector<apfloat>& a) const;

  ///Will convert the apfloat to its smallest value, zero or its heighest value
  double ToDoubleSafe(const apfloat& a) const noexcept;
  Coordinat2D ToDoubleSafe(const ApCoordinat2D& a) const noexcept;
  Coordinat3D ToDoubleSafe(const ApCoordinat3D& a) const noexcept;
  Coordinats2D ToDoubleSafe(const ApCoordinats2D& a) const noexcept;

  ///Convert a linestring to a closed polygon. If the linestring
  ///is open, close it
  /*

    +-+      +-+
    |    ->  | |
    +-+      +-+

    +-+      +-+
    | |  ->  | |
    +-+      +-+

  */
  Polygon ToPolygon(const Linestring& linestring) const noexcept;


  std::string ToStr(const Apfloats& p) const noexcept;
  std::string ToStr(const   Coordinat2D& p) const noexcept;
  std::string ToStr(const   Coordinats2D& p) const noexcept;
  std::string ToStr(const ApCoordinat2D& p) const noexcept;
  std::string ToStr(const ApCoordinats2D& p) const noexcept;
  std::string ToStr(const   Coordinat3D& p) const noexcept;
  std::string ToStr(const ApCoordinat3D& p) const noexcept;
  std::string ToStr(const ApCoordinats3D& p) const noexcept;
  std::string ToStr(const Linestring& polygon) const noexcept;
  std::string ToStr(const Polygon& polygon) const noexcept;
  std::string ToStr(const QPen& pen) noexcept;
  std::string ToStr(const QPoint& rect) noexcept;
  std::string ToStr(const QPointF& rect) noexcept;
  std::string ToStr(const QRect& rect) noexcept;
  std::string ToStr(const QRectF& rect) noexcept;

  ///Uses ToDoubleSafe(const apfloat&)
  std::string ToStrSafe(const Apfloat& f) const noexcept;

  //Create a complete SVG file contents
  std::string ToSvg(
    const Polygons& polygons,
    const double stroke_width = 1.0
  ) const noexcept;

  std::string ToSvg(
    const Polygons& polygons,
    const Linestrings& linestrings,
    const double stroke_width = 1.0
  ) const noexcept;

  std::string ToSvg(
    const Shapes& shapes,
    const double stroke_width = 1.0
  ) const noexcept { return ToSvg(shapes.first,shapes.second,stroke_width); }

  std::string ToWkt(
    const Polygons& polygons
  ) const noexcept { return ToWkt(polygons, {}); }

  std::string ToWkt(
    const Polygons& polygons,
    const Linestrings& linestrings
  ) const noexcept;

  std::string ToWkt(
    const Shapes& shapes
  ) const noexcept { return ToWkt(shapes.first,shapes.second); }


  Polygon Translate(
    const Polygon& shape,
    const double dx,
    const double dy
  ) const noexcept;


  //Throws if WKT does not convert to a linestring
  Linestring WktToLinestring(const std::string& wkt) const;

  //Throws if WKT does not convert to a polygon
  Polygon WktToPolygon(const std::string& wkt) const;
  Shapes WktToShapes(const std::string& wkt) const;
  Shapes WktToShapes(const std::vector<std::string>& wkt) const;
  std::string WktToSvg(const std::string& wkt, const double svg_stroke_width) const;
  std::string WktToSvg(const std::vector<std::string>& wkt, const double svg_stroke_width) const;

  private:

  #ifndef NDEBUG
  static void Test() noexcept;
  #endif

  ///Replaces what static_cast<int>(a) would do
  int ToInt(const Apfloat& a) const noexcept;

  //Create the line this geometry is in an SVG file
  std::string ToSvgStr(
    const std::vector<Polygon>& polygons,
    const double stroke_width = 1.0
  ) const noexcept;

  std::string ToSvgStr(
    const Polygon& polygon,
    const double stroke_width = 1.0
  ) const noexcept;

  std::string ToSvgStr(
    const Linestrings& linestrings,
    const double stroke_width = 1.0
  ) const noexcept;

  std::string ToSvgStr(
    const Linestring& linestring,
    const double stroke_width = 1.0
  ) const noexcept;

};

ribi::Geometry::Coordinat2D operator-(
  const ribi::Geometry::Coordinat2D& a,
  const ribi::Geometry::Coordinat2D& b
) noexcept;

ribi::Geometry::ApCoordinat2D operator-(
  const ribi::Geometry::ApCoordinat2D& a,
  const ribi::Geometry::ApCoordinat2D& b
) noexcept;

ribi::Geometry::Coordinat3D operator-(
  const ribi::Geometry::Coordinat3D& a,
  const ribi::Geometry::Coordinat3D& b
) noexcept;

ribi::Geometry::ApCoordinat3D operator-(
  const ribi::Geometry::ApCoordinat3D& a,
  const ribi::Geometry::ApCoordinat3D& b
) noexcept;



std::ostream& operator<<(std::ostream& os, const Geometry::Coordinat2D& p) noexcept;
std::ostream& operator<<(std::ostream& os, const Geometry::Coordinat3D& p) noexcept;
std::ostream& operator<<(std::ostream& os, const Geometry::Linestring& linestring) noexcept;
std::ostream& operator<<(std::ostream& os, const Geometry::Polygon& p) noexcept;

std::ostream& operator<<(std::ostream& os,const QPen& pen) noexcept;
std::ostream& operator<<(std::ostream& os,const QPoint& rect) noexcept;
std::ostream& operator<<(std::ostream& os,const QPointF& rect) noexcept;
std::ostream& operator<<(std::ostream& os,const QRect& rect) noexcept;
std::ostream& operator<<(std::ostream& os,const QRectF& rect) noexcept;


} //~namespace ribi

#endif // RIBI_GEOMETRY_H

 

 

 

 

 

./CppGeometry/geometry.cpp

 

//---------------------------------------------------------------------------
/*
Geometry, class with geometry functions
Copyright (C) 2014-2015 Richel Bilderbeek

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppGeometry.htm
//---------------------------------------------------------------------------
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
#pragma GCC diagnostic ignored "-Wunused-but-set-parameter"
#include "geometry.h"

#include <cassert>
#include <cmath>

#include <boost/math/constants/constants.hpp>
#include <boost/xpressive/xpressive.hpp>
#ifndef WIN32
#include <boost/geometry/geometries/linestring.hpp>
#endif

#include "apcplx.h"

#include <QPen>
#include <QPoint>
#include <QRect>

#include "ribi_regex.h"
#include "plane.h"
#include "trace.h"

#pragma GCC diagnostic pop

ribi::Geometry::Geometry()
{
  #ifndef NDEBUG
  Test();
  #endif
}

ribi::Geometry::Apfloat ribi::Geometry::Atan2(const Apfloat& dx, const Apfloat& dy) const noexcept
{
  const apfloat pi(boost::math::constants::pi<double>());
  return atan2(dy,dx); //apfloats's atan has reversed arguments
}

ribi::Geometry::Apfloat ribi::Geometry::CalcDotProduct(
  const ApCoordinat3D& a,
  const ApCoordinat3D& b
) const noexcept
{
  using boost::geometry::get;
  return
      (get<0>(a) * get<0>(b))
    + (get<1>(a) * get<1>(b))
    + (get<2>(a) * get<2>(b))
  ;
}

double ribi::Geometry::CalcDotProduct(
  const Coordinat3D& a,
  const Coordinat3D& b
) const noexcept
{
  return ToDoubleSafe(
    CalcDotProduct(
      ToApfloat(a),
      ToApfloat(b)
    )
  );
}

ribi::Geometry::ApCoordinat3D ribi::Geometry::CalcNormal(
  const ApCoordinat3D& a,
  const ApCoordinat3D& b,
  const ApCoordinat3D& c
) const noexcept
{
  const auto u = c - a;
  const auto v = b - a;
  return CalcCrossProduct(u,v);
}

ribi::Geometry::Coordinat3D ribi::Geometry::CalcNormal(
  const Coordinat3D& a,
  const Coordinat3D& b,
  const Coordinat3D& c
) const noexcept
{
  return ToDoubleSafe(CalcNormal(ToApfloat(a),ToApfloat(b),ToApfloat(c)));
}

std::vector<double> ribi::Geometry::CalcPlane(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) const noexcept
{
  const std::vector<apfloat> v
    = CalcPlane(ToApfloat(p1),ToApfloat(p2),ToApfloat(p3)
  );
  return ToDouble(v);
}

std::vector<apfloat> ribi::Geometry::CalcPlane(
  const ApCoordinat3D& p1,
  const ApCoordinat3D& p2,
  const ApCoordinat3D& p3
) const noexcept
{
  using boost::geometry::get;
  const auto v1 = p3 - p1;
  const auto v2 = p2 - p1;

  const auto cross = CalcCrossProduct(v1,v2);

  const auto a = get<0>(cross);
  const auto b = get<1>(cross);
  const auto c = get<2>(cross);

  const auto x = get<0>(p1);
  const auto y = get<1>(p1);
  const auto z = get<2>(p1);

  const auto term1 = a * x;
  const auto term2 = b * y;
  const auto term3 = c * z;
  const auto d = term1 + term2 + term3;

  return { a,b,c,d };
}

ribi::Geometry::ApCoordinats2D ribi::Geometry::CalcProjection(
  const ribi::Geometry::ApCoordinats3D& points) const
{
  assert(points.size() >= 3);
  assert(IsPlane(points));
  const std::unique_ptr<Plane> plane(new Plane(points[0],points[1],points[2]));
  assert(plane);

  #ifndef NDEBUG
  const bool verbose{false};
  if (verbose)
  {
    TRACE(plane->CanCalcX());
    TRACE(plane->CanCalcY());
    TRACE(plane->CanCalcZ());
  }
  #endif

  return plane->CalcProjection(points);
}

ribi::Geometry::Coordinats2D ribi::Geometry::CalcProjection(const Coordinats3D& v) const
{
  return ToDoubleSafe(CalcProjection(ToApfloat(v)));
}

boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
  ribi::Geometry::CreateHouseShape() const noexcept
{
  /*

    3-
     |
    2-    0
     |   / \
    1-  4   1
     |  |   |
    0-  3---2
     |
     +--|---|---|
        0   1   2

  */
  typedef boost::geometry::model::d2::point_xy<double> Coordinat2D;
  //typedef boost::geometry::model::polygon<Coordinat2D> Polygon;


  const std::vector<Coordinat2D> points {
    {0.5, 2.0}, //0
    {1.0, 1.0}, //1
    {1.0, 0.0}, //2
    {0.0, 0.0}, //3
    {0.0, 1.0}  //4
  };

  boost::geometry::model::polygon<Coordinat2D> house;
  boost::geometry::append(house, points);
  return house;
}

ribi::Geometry::Coordinat3D ribi::Geometry::CreatePoint(
  const double x,
  const double y,
  const double z
) const noexcept
{
  const auto c(Coordinat3D(x,y,z));
  return c;
}

ribi::Geometry::Rect ribi::Geometry::CreateRect(
  const double left,
  const double top,
  const double width,
  const double height
) const noexcept
{
  boost::geometry::model::box<Coordinat2D> rect(
    Coordinat2D(left        , top         ),
    Coordinat2D(left + width, top + height)
  );
  return rect;
}

boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
  ribi::Geometry::CreateShapeHeart(const double scale) const noexcept
{
  const std::vector<boost::geometry::model::d2::point_xy<double>> points {
    { 0.0 * scale, 1.0 * scale}, //0
    { 1.0 * scale, 2.0 * scale}, //1
    { 2.0 * scale, 1.0 * scale}, //2
    { 2.0 * scale, 0.0 * scale}, //3
    { 0.0 * scale,-2.0 * scale}, //4
    {-2.0 * scale, 0.0 * scale}, //5
    {-2.0 * scale, 1.0 * scale}, //6
    {-1.0 * scale, 2.0 * scale}  //7
  };
  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>> v;
  boost::geometry::append(v, points);
  return v;
}

boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
  ribi::Geometry::CreateShapeHouse(const double scale) const noexcept
{
  const std::vector<boost::geometry::model::d2::point_xy<double>> points {
    { 0.0 * scale, 2.0 * scale}, //0
    { 1.0 * scale, 1.0 * scale}, //1
    { 1.0 * scale,-1.0 * scale}, //2
    {-1.0 * scale,-1.0 * scale}, //3
    {-1.0 * scale, 1.0 * scale}  //4
  };
  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>> v;
  boost::geometry::append(v, points);
  return v;
}

boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
  ribi::Geometry::CreateShapePolygon(
  const int n,
  const double rotation,
  const double scale
) const noexcept
{
  assert(n >= 3 && "A polygon has at least three edges");
  const double tau { boost::math::constants::two_pi<double>() };
  std::vector<boost::geometry::model::d2::point_xy<double>> points;
  for (int i=0; i!=n; ++i)
  {
    const double angle { tau * static_cast<double>(i) / static_cast<double>(n) };
    const double x {  std::sin(angle + rotation) };
    const double y { -std::cos(angle + rotation) };
    boost::geometry::model::d2::point_xy<double> point(x * scale, y * scale);
    points.push_back(point);
  }

  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>> v;
  boost::geometry::append(v, points);
  return v;
}

boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
  ribi::Geometry::CreateShapeTriangle(const double scale) const noexcept
{
  const std::vector<boost::geometry::model::d2::point_xy<double>> points {
    { 0.0 * scale, 0.0 * scale}, //0
    { 1.0 * scale, 0.0 * scale}, //1
    { 0.0 * scale, 1.0 * scale}  //2
  };
  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>> v;
  boost::geometry::append(v, points);
  return v;
}

double ribi::Geometry::Fmod(const double x, const double mod) const noexcept
{
  #ifndef NDEBUG
  if (std::isnan(x)
    || std::isnan(mod)
    || mod <= 0.0
  )
  {
    TRACE("ERROR");
    TRACE(std::isnan(x));
    TRACE(std::isnan(mod));
    TRACE(x);
    TRACE("BREAK");
  }
  #endif

  assert(!std::isnan(x) && "Fmod its input value cannot be NaN");
  assert(!std::isnan(mod) && "Fmod its input mod cannot be NaN");
  assert(mod != 0.0 && "Cannot do a modulus of zero");
  assert(mod >= 0.0 && "Cannot do a modulus by a negative value");

  if (x < 0.0)
  {
    const double new_x {
      x + (mod * (1 + static_cast<int>(-x / mod)))
    };
    assert(new_x >= 0.0);
    return Fmod(new_x,mod);
  }
  double result = x - (mod * static_cast<int>(x / mod));

  assert(!std::isnan(result) && "Fmod its result cannot be NaN");
  assert(result >= 0.0 && "Fmod must return a value between zero and mod");
  assert(result < mod  && "Fmod must return a value between zero and mod");
  return result;
}

ribi::Geometry::Apfloat ribi::Geometry::Fmod(const Apfloat& x, const Apfloat& mod) const noexcept
{
  if (x < 0.0)
  {
    const Apfloat new_x(
      x + (mod * (1 + ToInt(-x / mod)))
    );
    assert(new_x >= 0.0);
    return Fmod(new_x,mod);
  }
  const Apfloat result = x - (mod * Apfloat(ToInt(x / mod)));

  if (result < 0.0) //Happens due to rounding errors
  {
    assert(result > -0.000000001);
    return 0.0;
  }
  if (result >= mod)
  {
    assert(result - mod < Apfloat(0.000000001));
    const Apfloat new_result = mod - Apfloat(std::numeric_limits<double>::epsilon());
    assert(new_result < mod);
    return new_result;
  }

  #ifndef NDEBUG
  if(result < 0.0 || result >= mod)
  {
    TRACE("ERROR");
    TRACE(ToStrSafe(result));
    TRACE(ToStrSafe(mod));
    TRACE("BREAK");
  }
  #endif
  assert(result >= 0.0 && "Fmod must return a value between zero and mod");
  assert(result < mod  && "Fmod must return a value between zero and mod");
  return result;
}



double ribi::Geometry::GetDistance(const double dx, const double dy) const noexcept
{
  //return GetDistance(
  //  boost::geometry::model::d2::point_xy<double>(0.0,0.0),
  //  boost::geometry::model::d2::point_xy<double>(dx,dy)
  //);
  return std::sqrt( (dx * dx) + (dy * dy) );
}

double ribi::Geometry::GetDistance(
  const Coordinat2D& a,
  const Coordinat2D& b
  ) const noexcept
{
  //return GetDistance(a.x() - b.x(),a.y() - b.y());
  return boost::geometry::distance(a,b);
}



std::string ribi::Geometry::GetVersion() const noexcept
{
  return "1.6";
}

std::vector<std::string> ribi::Geometry::GetVersionHistory() const noexcept
{
  return {
    "2014-02-25: version 1.0: initial version",
    "2014-03-11: version 1.1: removed use of custom coordinat classes, use Boost.Geometry instead"
    "2014-03-25: version 1.2: fixed bug in IsConvex",
    "2014-05-01: version 1.3: minor improments and cleanup",
    "2014-07-15: version 1.4: fixed bug in IsClockwise",
    "2014-07-17: version 1.5: fixed bug in IsCounterClockwise",
    "2014-07-18: version 1.6: renamed some member functions to explicitly mention the coordinat system (Cartesian or screen) used"
  };
}



bool ribi::Geometry::IsConvex(Polygon polygon) const noexcept
{
  //assert(boost::geometry::num_points(polygon) == points.size()
  //  && "Points and polygon have the same number of points");
  const bool verbose{false};
  boost::geometry::correct(polygon);
  Polygon hull;
  boost::geometry::convex_hull(polygon, hull);
  const bool is_convex_first = boost::geometry::equals(polygon,hull);
  if (verbose) { TRACE(is_convex_first); }
  if (is_convex_first)
  {
    return true;
  }
  //Correct for bug https://svn.boost.org/trac/boost/ticket/9873
  //A polygon is convex when it does not use the edge with the possibly longest distance
  const std::vector<Coordinat2D> v = polygon.outer();
  //Find longest distance
  double max_dist = 0.0;
  const int sz = static_cast<int>(v.size());
  for (int i=0; i!=sz; ++i)
  {
    for (int j=i+1; j!=sz; ++j)
    {
      const double dist = GetDistance(v[i],v[j]);
      if (dist > max_dist) { max_dist = dist; }
    }
  }

  for (int i=1; i!=sz; ++i)
  {
    const double dist = GetDistance(v[i-1],v[i]);
    if (dist == max_dist) return false;
  }
  return true;
}

bool ribi::Geometry::IsConvex(const Coordinats2D& points) const noexcept
{
  Polygon polygon;
  for (const auto& point: points)
  {
    boost::geometry::append(polygon,point);
  };
  assert(boost::geometry::num_points(polygon) == points.size());
  const bool is_convex = IsConvex(polygon);
  return is_convex;
}

bool ribi::Geometry::IsConvex(const ApCoordinats3D& points) const noexcept
{
  const bool verbose{false};

  #ifndef NDEBUG
  assert(points.size() >= 3);
  if (points.size() == 3)
  {
    //Three different points are always convex
    assert(!boost::geometry::equals(points[0],points[1]));
    assert(!boost::geometry::equals(points[0],points[2]));
    assert(!boost::geometry::equals(points[1],points[2]));
    return true;
  }
  assert(points.size() == 4);
  if(!IsPlane(points))
  {
    TRACE("ERROR");
    TRACE(points.size());
    for (const auto& point: points) { TRACE(Geometry().ToStr(point)); }
    TRACE("BREAK");
  }
  assert(IsPlane(points));
  #endif // NDEBUG
  if (verbose)
  {
    std::stringstream s;
    s << "{";
    for (const auto& point3d: points)
    {
      s << ToStr(point3d) << ",";
    }
    std::string po_str(s.str());
    po_str[po_str.size() - 1] = '}';
    TRACE(po_str);
  }
  //Use the first three points for a Plane
  for (const std::vector<int> v:
    {
      std::vector<int>( {0,1,2} ),
      std::vector<int>( {0,1,3} ),
      std::vector<int>( {0,2,1} ),
      std::vector<int>( {0,2,3} ),
      std::vector<int>( {0,3,1} ),
      std::vector<int>( {0,3,2} ),
      std::vector<int>( {1,0,2} ),
      std::vector<int>( {1,0,3} ),
      std::vector<int>( {1,2,0} ),
      std::vector<int>( {1,2,3} ),
      std::vector<int>( {1,3,0} ),
      std::vector<int>( {1,3,2} ),
      std::vector<int>( {2,0,1} ),
      std::vector<int>( {2,0,3} ),
      std::vector<int>( {2,1,0} ),
      std::vector<int>( {2,1,3} ),
      std::vector<int>( {2,3,0} ),
      std::vector<int>( {2,3,1} ),
      std::vector<int>( {3,0,1} ),
      std::vector<int>( {3,0,2} ),
      std::vector<int>( {3,1,0} ),
      std::vector<int>( {3,1,2} ),
      std::vector<int>( {3,2,0} ),
      std::vector<int>( {3,2,1} )
    }
  )
  {
    const boost::shared_ptr<Plane> plane(
      new Plane(
        points[v[0]],
        points[v[1]],
        points[v[2]]
      )
    );
    assert(plane);

    #ifndef NDEBUG
    if (verbose)
    {
      TRACE(*plane);
      TRACE(plane->CanCalcX());
      TRACE(plane->CanCalcY());
      TRACE(plane->CanCalcZ());
    }
    #endif

    const ApCoordinats2D apcoordinats2d = plane->CalcProjection(points);
    if (verbose)
    {
      TRACE(ToStr(apcoordinats2d));
    }
    const Coordinats2D coordinats2d = ToDoubleSafe(apcoordinats2d);
    if (verbose)
    {
      TRACE(ToStr(coordinats2d));
    }
    if (IsConvex(coordinats2d)) return true;
  }
  return false;
}

bool ribi::Geometry::IsConvex(const Coordinats3D& points) const noexcept
{
  return IsConvex(ToApfloat(points));
}



std::function<bool(const ribi::Geometry::Coordinat2D& lhs, const ribi::Geometry::Coordinat2D& rhs)>
  ribi::Geometry::Equals2d() const noexcept
{
  return [](const ribi::Geometry::Coordinat2D& lhs, const ribi::Geometry::Coordinat2D& rhs)
  {
    using boost::geometry::get;
    return
      get<0>(lhs) == get<0>(rhs)
      && get<1>(lhs) == get<1>(rhs)
    ;
  };
}

std::function<bool(const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)>
  ribi::Geometry::Equals() const noexcept
{
  return [](const Coordinat3D& lhs, const Coordinat3D& rhs)
  {
    using boost::geometry::get;
    return
      get<0>(lhs) == get<0>(rhs)
      && get<1>(lhs) == get<1>(rhs)
      && get<2>(lhs) == get<2>(rhs)
    ;
  };
}

std::function<bool(const ribi::Geometry::ApCoordinat3D& lhs, const ribi::Geometry::ApCoordinat3D& rhs)>
  ribi::Geometry::Equals3d() const noexcept
{
  return [](const ApCoordinat3D& lhs, const ApCoordinat3D& rhs)
  {
    using boost::geometry::get;
    return
      get<0>(lhs) == get<0>(rhs)
      && get<1>(lhs) == get<1>(rhs)
      && get<2>(lhs) == get<2>(rhs)
    ;
  };
}

bool ribi::Geometry::IsPlane(const std::vector<ApCoordinat3D>& v) const noexcept
{
  const bool verbose{false};
  using boost::geometry::get;

  if (v.size() < 3) return false;
  if (v.size() == 3) return true;
  #ifndef NDEBUG
  if (v.size() > 4)
  {
    TRACE("ERROR");
    TRACE(v.size());
    TRACE("BREAK");
  }
  #endif
  assert(v.size() == 4);

  try
  {
    const std::unique_ptr<Plane> plane(new Plane(v[0],v[1],v[2]));
    assert(plane);
    assert(plane->IsInPlane(v[3]) == (plane->CalcError(v[3]) <= plane->CalcMaxError(v[3])));
    return plane->IsInPlane(v[3]);
  }
  catch (std::logic_error& e)
  {
    if (verbose)
    {
      std::stringstream s;
      s << "Geometry::IsPlane: not in plane, as plane cannot be constructed ('"
        << e.what() << "')"
      ;
      TRACE(s.str());
    }
    return false;
  }
}

bool ribi::Geometry::IsPlane(const std::vector<Coordinat3D>& v) const noexcept
{
  return IsPlane(ToApfloat(v));
}


std::function<bool(const ribi::Geometry::Coordinat2D& lhs, const ribi::Geometry::Coordinat2D& rhs)>
  ribi::Geometry::Order2dByX() const noexcept
{
  return [](const ribi::Geometry::Coordinat2D& lhs, const ribi::Geometry::Coordinat2D& rhs)
  {
    using boost::geometry::get;
    if (get<0>(lhs) < get<0>(rhs)) return true;
    if (get<0>(lhs) > get<0>(rhs)) return false;
    return get<1>(lhs) < get<1>(rhs);
  };
}

std::function<bool(const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)>
  ribi::Geometry::OrderByX() const noexcept
{
  return [](const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)
  {
    using boost::geometry::get;
    if (get<0>(lhs) < get<0>(rhs)) return true;
    if (get<0>(lhs) > get<0>(rhs)) return false;
    if (get<1>(lhs) < get<1>(rhs)) return true;
    if (get<1>(lhs) > get<1>(rhs)) return false;
    return get<2>(lhs) < get<2>(rhs);
  };
}

std::function<bool(const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)>
  ribi::Geometry::OrderByZ() const noexcept
{
  return [](const ribi::Geometry::Coordinat3D& lhs, const ribi::Geometry::Coordinat3D& rhs)
  {
    using boost::geometry::get;
    if (get<2>(lhs) < get<2>(rhs)) return true;
    if (get<2>(lhs) > get<2>(rhs)) return false;
    if (get<1>(lhs) < get<1>(rhs)) return true;
    if (get<1>(lhs) > get<1>(rhs)) return false;
    return get<0>(lhs) < get<0>(rhs);
  };
}

ribi::Geometry::Linestring ribi::Geometry::Rescale(
  const Linestring& linestring,
  const double scale,
  const double scale_origin_x,
  const double scale_origin_y
) const noexcept
{
  std::vector<Coordinat2D> v;


  for (const auto& point:linestring)
  {
    const double x = point.x();
    const double dx = x - scale_origin_x;
    const double new_x = scale_origin_x + (scale * dx);

    const double y = point.y();
    const double dy = y - scale_origin_y;
    const double new_y = scale_origin_y + (scale * dy);

    v.push_back(Coordinat2D(new_x,new_y));
  }

  Linestring new_linestring;
  boost::geometry::append(new_linestring,v);
  return new_linestring;
}


ribi::Geometry::Linestrings ribi::Geometry::Rescale(
  const Linestrings& linestrings,
  const double scale,
  const double scale_origin_x,
  const double scale_origin_y
) const noexcept
{
  Linestrings new_linestrings;
  for (const auto& linestring: linestrings)
  {
    new_linestrings.push_back(
      Rescale(linestring,scale,scale_origin_x,scale_origin_y)
    );
  }
  return new_linestrings;
}

ribi::Geometry::Polygon ribi::Geometry::Rescale(
  const Polygon& polygon,
  const double scale,
  const double scale_origin_x,
  const double scale_origin_y
) const noexcept
{
  typedef boost::geometry::model::d2::point_xy<double> Coordinat2D;
  typedef boost::geometry::model::polygon<Coordinat2D> Polygon;

  boost::geometry::model::ring<Coordinat2D> points;
  boost::geometry::convert(polygon,points);

  for (auto& point:points)
  {
    const double x = point.x();
    const double dx = x - scale_origin_x;
    const double new_x = scale_origin_x + (scale * dx);

    const double y = point.y();
    const double dy = y - scale_origin_y;
    const double new_y = scale_origin_y + (scale * dy);

    point = Coordinat2D(new_x,new_y);
  }

  Polygon new_polygon;
  boost::geometry::append(new_polygon, points);
  return new_polygon;
}

ribi::Geometry::Polygons ribi::Geometry::Rescale(
  const Polygons& polygons,
  const double scale,
  const double scale_origin_x,
  const double scale_origin_y
) const noexcept
{
  Polygons new_polygons;
  for (const auto& polygon: polygons)
  {
    new_polygons.push_back(
      Rescale(polygon,scale,scale_origin_x,scale_origin_y)
    );
  }
  return new_polygons;
}

ribi::Geometry::Shapes ribi::Geometry::Rescale(
  const Shapes& shapes,
  const double scale,
  const double scale_origin_x,
  const double scale_origin_y
) const noexcept
{
  Shapes new_shapes;
  new_shapes.first
    = Rescale(shapes.first,scale,scale_origin_x,scale_origin_y);
  new_shapes.second
    = Rescale(shapes.second,scale,scale_origin_x,scale_origin_y);
  return new_shapes;
}

ribi::Geometry::Apfloats ribi::Geometry::ToApfloat(const Doubles& v) const noexcept
{
  Apfloats w;
  std::transform(std::begin(v),std::end(v),std::back_inserter(w),
    [](const double d) { return Apfloat(d); }
  );
  return w;
}

ribi::Geometry::ApCoordinat2D ribi::Geometry::ToApfloat(const Coordinat2D& p) const noexcept
{
  return ApCoordinat2D(
    boost::geometry::get<0>(p),
    boost::geometry::get<1>(p)
  );
}

ribi::Geometry::ApCoordinat3D ribi::Geometry::ToApfloat(const Coordinat3D& p) const noexcept
{
  return ApCoordinat3D(
    boost::geometry::get<0>(p),
    boost::geometry::get<1>(p),
    boost::geometry::get<2>(p)
  );
}

ribi::Geometry::ApCoordinats2D ribi::Geometry::ToApfloat(const Coordinats2D& v) const noexcept
{
  ApCoordinats2D w;
  std::transform(std::begin(v),std::end(v),std::back_inserter(w),
    [this](const Coordinat2D& c) { return ToApfloat(c); }
  );
  return w;
}

ribi::Geometry::ApCoordinats3D ribi::Geometry::ToApfloat(const Coordinats3D& v) const noexcept
{
  ApCoordinats3D w;
  std::transform(std::begin(v),std::end(v),std::back_inserter(w),
    [this](const Coordinat3D& c) { return ToApfloat(c); }
  );
  return w;
}


double ribi::Geometry::ToDouble(const apfloat& a) const
{
  std::stringstream s;
  s << std::fixed << std::setprecision(99) << a;
  //s << pretty << a;
  const double x = boost::lexical_cast<double>(s.str());
  return x;
}

std::vector<double> ribi::Geometry::ToDouble(const std::vector<apfloat>& v) const
{
  std::vector<double> w;
  std::transform(std::begin(v),std::end(v),std::back_inserter(w),
    [this](const apfloat& f) { return ToDouble(f); }
  );
  assert(v.size() == w.size());
  return w;
}

double ribi::Geometry::ToDoubleSafe(const apfloat& a) const noexcept
{
  try
  {
    return ToDouble(a);
  }
  catch (boost::bad_lexical_cast&)
  {
    if (a < apfloat(std::numeric_limits<double>::lowest())) { return std::numeric_limits<double>::lowest(); }
    if (a > apfloat(std::numeric_limits<double>::max())) { return std::numeric_limits<double>::max(); }
    if (a > apfloat(0.0))
    {
       if (a < apfloat(0.5) * apfloat(std::numeric_limits<double>::denorm_min())) { return 0.0; }
       return std::numeric_limits<double>::denorm_min();
    }
    if (a < apfloat(0.0))
    {
       if (a > apfloat(-0.5) * apfloat(std::numeric_limits<double>::denorm_min())) { return 0.0; }
       const double r = -std::numeric_limits<double>::denorm_min();
       assert(r != 0.0);
       return r;
    }
    TRACE(a);
    assert(!"Should not get here");
    return 0.0;
  }
}

ribi::Geometry::Coordinat2D ribi::Geometry::ToDoubleSafe(const ApCoordinat2D& a) const noexcept
{
  return Coordinat2D(
    ToDoubleSafe(boost::geometry::get<0>(a)),
    ToDoubleSafe(boost::geometry::get<1>(a))
  );
}

ribi::Geometry::Coordinat3D ribi::Geometry::ToDoubleSafe(const ApCoordinat3D& a) const noexcept
{
  return Coordinat3D(
    ToDoubleSafe(boost::geometry::get<0>(a)),
    ToDoubleSafe(boost::geometry::get<1>(a)),
    ToDoubleSafe(boost::geometry::get<2>(a))
  );
}

ribi::Geometry::Coordinats2D ribi::Geometry::ToDoubleSafe(const ribi::Geometry::ApCoordinats2D& v) const noexcept
{
  Coordinats2D w;
  std::transform(std::begin(v),std::end(v),std::back_inserter(w),
    [this](const ApCoordinat2D& c) { return ToDoubleSafe(c); }
  );
  return w;
}

int ribi::Geometry::ToInt(const Apfloat& a) const noexcept
{
  return static_cast<int>(ToDoubleSafe(a));
}

ribi::Geometry::Polygon ribi::Geometry::ToPolygon(const Linestring& linestring) const noexcept
{
  std::vector<Coordinat2D> v;
  for (const auto& point: linestring) { v.push_back(point); }

  Polygon polygon;
  if (v.empty()) return polygon;
  if (boost::geometry::equals(v.front(),v.back())) { v.pop_back(); }

  boost::geometry::append(polygon,v);
  return polygon;
}

std::string ribi::Geometry::ToStrSafe(const apfloat& f) const noexcept
{
  try
  {
    std::stringstream s;
    s << pretty << f;
    return s.str();
  }
  catch (boost::bad_lexical_cast&)
  {
    std::stringstream s;
    s << ToDoubleSafe(f);
    return s.str();
  }
}

std::string ribi::Geometry::ToStr(const Apfloats& v) const noexcept
{
  std::stringstream s;
  for (const auto& i:v) { s << ToStrSafe(i) << ','; }
  std::string t = s.str();
  if (!t.empty()) t.pop_back();
  return t;
}

std::string ribi::Geometry::ToStr(const ApCoordinat2D& p) const noexcept
{
  std::stringstream s;
  s << '('
    << ToStrSafe(boost::geometry::get<0>(p)) << ','
    << ToStrSafe(boost::geometry::get<1>(p))
    << ')'
  ;
  return s.str();
}

std::string ribi::Geometry::ToStr(const ApCoordinats2D& v) const noexcept
{
  std::stringstream s;
  for (const auto& i:v) { s << ToStr(i) << '-'; }
  std::string t = s.str();
  if (!t.empty()) t.pop_back();
  return t;
}

std::string ribi::Geometry::ToStr(const ApCoordinat3D& p) const noexcept
{
  std::stringstream s;
  s << '('
    << ToStrSafe(boost::geometry::get<0>(p)) << ','
    << ToStrSafe(boost::geometry::get<1>(p)) << ','
    << ToStrSafe(boost::geometry::get<2>(p))
    << ')'
  ;
  return s.str();
}

std::string ribi::Geometry::ToStr(const ApCoordinats3D& v) const noexcept
{
  std::stringstream s;
  for (const auto& i:v) { s << ToStr(i) << '-'; }
  std::string t = s.str();
  if (!t.empty()) t.pop_back();
  return t;
}

std::string ribi::Geometry::ToStr(const Coordinat2D& p) const noexcept
{
  std::stringstream s;
  s << p;
  return s.str();
}

std::string ribi::Geometry::ToStr(const Coordinat3D& p) const noexcept
{
  std::stringstream s;
  s << p;
  return s.str();
}

std::string ribi::Geometry::ToStr(const Coordinats2D& v) const noexcept
{
  std::stringstream s;
  for (const auto& i:v) { s << ToStr(i) << '-'; }
  std::string t = s.str();
  if (!t.empty()) t.pop_back();
  return t;
}

std::string ribi::Geometry::ToStr(const Linestring& linestring) const noexcept
{
  std::stringstream s;
  s << linestring;
  return s.str();
}

std::string ribi::Geometry::ToStr(const Polygon& polygon) const noexcept
{
  std::stringstream s;
  s << polygon;
  return s.str();
}

std::string ribi::Geometry::ToSvg(
  const Polygons& polygons,
  const double stroke_width
) const noexcept
{
  assert(stroke_width > 0.0);
  std::stringstream s;
  s
    << std::setprecision(99)
    << R"*(<svg xmlns="http://www.w3.org/2000/svg" version="1.1">)*"
    << '\n'
    << ToSvgStr(polygons,stroke_width) << '\n'
    << R"*(</svg>)*"
  ;
  return s.str();
}

std::string ribi::Geometry::ToSvg(
  const Polygons& polygons,
  const Linestrings& linestrings,
  const double stroke_width
) const noexcept
{
  assert(stroke_width > 0.0);
  std::stringstream s;
  s
    << std::setprecision(99)
    << R"*(<svg xmlns="http://www.w3.org/2000/svg" version="1.1">)*"
    << '\n'
    << ToSvgStr(polygons,stroke_width) << '\n'
    << ToSvgStr(linestrings,stroke_width) << '\n'
    << R"*(</svg>)*"
  ;
  return s.str();
}

std::string ribi::Geometry::ToSvgStr(
  const Polygons& polygons,
  const double stroke_width
) const noexcept
{
  assert(stroke_width > 0.0);
  std::stringstream s;
  for (const Polygon& polygon: polygons)
  {
    s << ToSvgStr(polygon,stroke_width) << '\n';
  }
  return s.str();
}

std::string ribi::Geometry::ToSvgStr(
  const Linestrings& linestrings,
  const double stroke_width
) const noexcept
{
  assert(stroke_width > 0.0);
  std::stringstream s;
  for (const Linestring& linestring: linestrings)
  {
    s << ToSvgStr(linestring,stroke_width) << '\n';
  }
  return s.str();
}

std::string ribi::Geometry::ToSvgStr(
  const Polygon& polygon,
  const double stroke_width
) const noexcept
{
  if (stroke_width <= 0.0) { return ""; }
  assert(stroke_width > 0.0);
  const std::vector<Coordinat2D> points = polygon.outer();
  const int n_points = static_cast<int>(points.size());
  if (n_points < 3) { return ""; }
  assert(n_points >= 3 && "A polygon has at least three points");
  //Move to first point
  std::stringstream s;
  s
    << std::setprecision(99)
    <<  R"*(  <path d="M )*"
    << points[0].x()
    << " "
    << points[0].y()
  ;
  //Draw lines to others
  s << " L ";
  for (int i=1; i!=n_points; ++i)
  {
    s
      << points[i].x()
      << " "
      << points[i].y()
      << " "
    ;

    //No trailing comma
    if (i != n_points - 1)
    {
      s << ",";
    }
  }
  s
    << R"*(z" stroke="black" fill="none" stroke-width=")*"
    << stroke_width
    << R"*("/>)*"
  ;
  return s.str();
}

std::string ribi::Geometry::ToSvgStr(
  const Linestring& linestring,
  const double stroke_width
) const noexcept
{
  if (stroke_width <= 0.0) { return ""; }
  assert(stroke_width > 0.0);
  const std::vector<Coordinat2D> points = linestring;
  const int n_points = static_cast<int>(points.size());
  if (n_points < 2) { return ""; }
  assert(n_points >= 2 && "A linestring has at least two points");
  //Move to first point
  std::stringstream s;
  s
    << std::setprecision(99)
    <<  R"*(  <path d="M )*"
    << points[0].x()
    << " "
    << points[0].y()
  ;
  //Draw lines to others
  s << " L ";
  for (int i=1; i!=n_points; ++i)
  {
    s
      << points[i].x()
      << " "
      << points[i].y()
      << " "
    ;

    //No trailing comma
    if (i != n_points - 1) { s << ","; }
  }
  s
    << R"*(" stroke="black" fill="none" stroke-width=")*"
    << stroke_width
    << R"*("/>)*"
  ;
  return s.str();
}

std::string ribi::Geometry::ToWkt(
  const Polygons& polygons,
  const Linestrings& linestrings
) const noexcept
{
  std::string s;
  for (const auto& polygon:polygons)
  {
    const auto w = boost::geometry::wkt<>(polygon);
    std::stringstream t;
    t << w << ',';
    s += t.str();
  }
  for (const auto& linestring:linestrings)
  {
    const auto w = boost::geometry::wkt<>(linestring);
    std::stringstream t;
    t << w << ',';
    s += t.str();
  }
  if (!s.empty()) s.pop_back();
  return s;
}

boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>
  ribi::Geometry::Translate(
    const boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>& shape,
    const double dx,
    const double dy
  ) const noexcept
{
  typedef boost::geometry::model::d2::point_xy<double> Coordinat2D;
  typedef boost::geometry::model::polygon<Coordinat2D> Polygon;

  boost::geometry::model::ring<Coordinat2D> points;
  boost::geometry::convert(shape,points);

  for (auto& point:points)
  {
    point = Coordinat2D(point.x() + dx, point.y() + dy);
  }

  Polygon new_shape;
  boost::geometry::append(new_shape, points);
  return new_shape;
}

ribi::Geometry::Linestring ribi::Geometry::WktToLinestring(const std::string& wkt) const
{
  Linestring linestring;
  boost::geometry::read_wkt(wkt,linestring);
  return linestring;
}

ribi::Geometry::Polygon ribi::Geometry::WktToPolygon(const std::string& wkt) const
{
  Polygon polygon;
  boost::geometry::read_wkt(wkt,polygon);
  return polygon;
}

ribi::Geometry::Shapes ribi::Geometry::WktToShapes(const std::string& wkt) const
{
  std::vector<std::string> v;
  const std::string regex_str = Regex().GetRegexShapes();
  for (const auto& s: Regex().GetRegexMatches(wkt,regex_str))
  {
    v.push_back(s);
  }
  return WktToShapes(v);
}

ribi::Geometry::Shapes ribi::Geometry::WktToShapes(const std::vector<std::string>& wkt) const
{
  Polygons polygons;
  Linestrings linestrings;
  for (const auto& s: wkt)
  {
    try
    {
      const auto polygon = WktToPolygon(s);
      //boost::geometry::read_wkt(s,polygon);
      polygons.push_back(polygon);
    }
    catch (std::exception&)
    {
      //OK
    }

    try
    {
      const auto linestring = WktToLinestring(s);
      //Linestring linestring;
      //boost::geometry::read_wkt(s,linestring);
      linestrings.push_back(linestring);
    }
    catch (std::exception&)
    {
      //OK
    }
  }
  const auto shapes = std::make_pair(polygons,linestrings);
  return shapes;
}


std::string ribi::Geometry::WktToSvg(const std::string& wkt, const double svg_stroke_width) const
{
  return ToSvg(WktToShapes(wkt),svg_stroke_width);
}

std::string ribi::Geometry::WktToSvg(const std::vector<std::string>& wkt, const double svg_stroke_width) const
{
  return ToSvg(WktToShapes(wkt),svg_stroke_width);
}

ribi::Geometry::Coordinat2D ribi::operator-(
  const ribi::Geometry::Coordinat2D& a,
  const ribi::Geometry::Coordinat2D& b
) noexcept
{
  return std::remove_const<std::remove_reference<decltype(a)>::type>::type(
    boost::geometry::get<0>(a) - boost::geometry::get<0>(b),
    boost::geometry::get<1>(a) - boost::geometry::get<1>(b)
  );
}

ribi::Geometry::ApCoordinat2D ribi::operator-(
  const ribi::Geometry::ApCoordinat2D& a,
  const ribi::Geometry::ApCoordinat2D& b
) noexcept
{
  return std::remove_const<std::remove_reference<decltype(a)>::type>::type(
    boost::geometry::get<0>(a) - boost::geometry::get<0>(b),
    boost::geometry::get<1>(a) - boost::geometry::get<1>(b)
  );
}

ribi::Geometry::Coordinat3D ribi::operator-(
    const ribi::Geometry::Coordinat3D& a,
    const ribi::Geometry::Coordinat3D& b
  ) noexcept
{
  return std::remove_const<std::remove_reference<decltype(a)>::type>::type(
    boost::geometry::get<0>(a) - boost::geometry::get<0>(b),
    boost::geometry::get<1>(a) - boost::geometry::get<1>(b),
    boost::geometry::get<2>(a) - boost::geometry::get<2>(b)
  );
}

ribi::Geometry::ApCoordinat3D ribi::operator-(
    const ribi::Geometry::ApCoordinat3D& a,
    const ribi::Geometry::ApCoordinat3D& b
  ) noexcept
{
  return std::remove_const<std::remove_reference<decltype(a)>::type>::type(
    boost::geometry::get<0>(a) - boost::geometry::get<0>(b),
    boost::geometry::get<1>(a) - boost::geometry::get<1>(b),
    boost::geometry::get<2>(a) - boost::geometry::get<2>(b)
  );
}

std::ostream& ribi::operator<<(std::ostream& os, const Geometry::Coordinat2D& p) noexcept
{
  using boost::geometry::get;
  os << '(' << get<0>(p) << ',' << get<1>(p) << ')';
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os, const Geometry::Coordinat3D& p) noexcept
{
  using boost::geometry::get;
  os << '(' << get<0>(p) << ',' << get<1>(p) << ',' << get<2>(p) << ')';
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os, const Geometry::Linestring& p) noexcept
{
  const auto points = p;
  const int n_points = static_cast<int>(points.size());
  for (int i=0; i!=n_points; ++i)
  {
    os << points[i];
    if (i != n_points-1) { os << ','; }
  }
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os, const Geometry::Polygon& p) noexcept
{
  const auto points = p.outer();
  const int n_points = static_cast<int>(points.size());
  for (int i=0; i!=n_points; ++i)
  {
    os << points[i];
    if (i != n_points-1) { os << ','; }
  }
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os,const QPen& pen) noexcept
{
  os << pen.color() << " (" << pen.widthF() << ')';
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os,const QPoint& point) noexcept
{
  os << '(' << point.x() << ',' << point.y() << ')';
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os,const QPointF& point) noexcept
{
  os << '(' << point.x() << ',' << point.y() << ')';
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os,const QRect& rect) noexcept
{
  os
    << '(' << rect.left() << ',' << rect.top() << ')'
    << ','
    << '(' << rect.width() << 'x' << rect.height() << ')'
  ;
  return os;
}

std::ostream& ribi::operator<<(std::ostream& os,const QRectF& rect) noexcept
{
  os
    << '(' << rect.left() << ',' << rect.top() << ')'
    << ','
    << '(' << rect.width() << 'x' << rect.height() << ')'
  ;
  return os;

}

 

 

 

 

 

./CppGeometry/geometry_calccenter.cpp

 

#include "geometry.h"

ribi::Geometry::ApCoordinat2D ribi::Geometry::CalcCenter(const ApCoordinats2D& points) const noexcept
{
  ApCoordinat2D sum(0.0,0.0);
  for (const auto& point: points)
  {
    sum.x(sum.x() + point.x());
    sum.y(sum.y() + point.y());
  }
  const Apfloat n(static_cast<double>(points.size()));
  const ApCoordinat2D center(
    sum.x() / n,
    sum.y() / n
  );
  return center;

}

ribi::Geometry::ApCoordinat3D ribi::Geometry::CalcCenter(const ApCoordinats3D& points) const noexcept
{
  Apfloat sum_x = 0.0;
  Apfloat sum_y = 0.0;
  Apfloat sum_z = 0.0;
  for (const auto& point: points)
  {
    sum_x += boost::geometry::get<0>(point);
    sum_y += boost::geometry::get<1>(point);
    sum_z += boost::geometry::get<2>(point);
  }
  const Apfloat n(static_cast<double>(points.size()));
  const ApCoordinat3D center(
    sum_x / n,
    sum_y / n,
    sum_z / n
  );
  return center;
}

 

 

 

 

 

./CppGeometry/geometry_calccrossproduct.cpp

 

#include "geometry.h"

template <class T>
boost::geometry::model::point<T,3,boost::geometry::cs::cartesian>
CalcCrossProductImpl(
  const boost::geometry::model::point<T,3,boost::geometry::cs::cartesian>& a,
  const boost::geometry::model::point<T,3,boost::geometry::cs::cartesian>& b)
{
  using boost::geometry::get;
  const boost::geometry::model::point<T,3,boost::geometry::cs::cartesian> result{
    (get<1>(a) * get<2>(b)) - (get<2>(a) * get<1>(b)),
    (get<2>(a) * get<0>(b)) - (get<0>(a) * get<2>(b)),
    (get<0>(a) * get<1>(b)) - (get<1>(a) * get<0>(b))
  };

  return result;
}

ribi::Geometry::Coordinat3D ribi::Geometry::CalcCrossProduct(
  const Coordinat3D& a,
  const Coordinat3D& b
) const noexcept
{
  return CalcCrossProductImpl<double>(a,b);
  /*
  using boost::geometry::get;
  return Coordinat3D(
    (get<1>(a) * get<2>(b)) - (get<2>(a) * get<1>(b)),
    (get<2>(a) * get<0>(b)) - (get<0>(a) * get<2>(b)),
    (get<0>(a) * get<1>(b)) - (get<1>(a) * get<0>(b))
  );
  */
}

ribi::Geometry::ApCoordinat3D ribi::Geometry::CalcCrossProduct(
  const ApCoordinat3D& a,
  const ApCoordinat3D& b
) const noexcept
{
  return CalcCrossProductImpl<Apfloat>(a,b);
  /*
  using boost::geometry::get;
  return ApCoordinat3D(
    (get<1>(a) * get<2>(b)) - (get<2>(a) * get<1>(b)),
    (get<2>(a) * get<0>(b)) - (get<0>(a) * get<2>(b)),
    (get<0>(a) * get<1>(b)) - (get<1>(a) * get<0>(b))
  );
  */
}

 

 

 

 

 

./CppGeometry/geometry_getangle.cpp

 

#include "geometry.h"

#include "apcplx.h"

template <class T> T Atan2Impl(const T& a, const T& b) noexcept;
template <> double Atan2Impl(const double& a, const double& b) noexcept { return std::atan2(a,b); }
template <> apfloat Atan2Impl(const apfloat& a, const apfloat& b) noexcept { return ribi::Geometry().Atan2(a,b); }

template <class T>
T GetAngleClockCartesianImpl(const T& dx, const T& dy) noexcept
{
  const T pi = boost::math::constants::pi<double>();
  const T tau = boost::math::constants::two_pi<double>();
  const T angle = ribi::Geometry().Fmod(pi - Atan2Impl(dx,-dy),tau);
  assert(angle >= 0.0 && "GetAngle must return a value between zero and two pi");
  assert(angle < tau  && "GetAngle must return a value between zero and two pi");
  return angle;
}

double ribi::Geometry::GetAngleClockCartesian(const double dx, const double dy) const noexcept
{
  return GetAngleClockCartesianImpl(dx,dy);
}

ribi::Geometry::Apfloat ribi::Geometry::GetAngleClockCartesian(const Apfloat& dx, const Apfloat& dy) const noexcept
{
  return GetAngleClockCartesianImpl(dx,dy);
}

double ribi::Geometry::GetAngleClockCartesian(const Coordinat2D& p) const noexcept
{
  return GetAngleClockCartesianImpl(boost::geometry::get<0>(p),boost::geometry::get<1>(p));
}

ribi::Geometry::Apfloat ribi::Geometry::GetAngleClockCartesian(const ApCoordinat2D& p) const noexcept
{
  return GetAngleClockCartesianImpl(boost::geometry::get<0>(p),boost::geometry::get<1>(p));
}


template <class T>
T GetAngleClockScreenImpl(const T& dx, const T& dy) noexcept
{
  const T pi = boost::math::constants::pi<double>();
  const T tau = boost::math::constants::two_pi<double>();
  const T angle = ribi::Geometry().Fmod(pi - Atan2Impl(dx,dy),tau);
  assert(angle >= 0.0 && "GetAngle must return a value between zero and two pi");
  assert(angle < tau  && "GetAngle must return a value between zero and two pi");
  return angle;
}

double ribi::Geometry::GetAngleClockScreen(const double dx, const double dy) const noexcept
{
  //const double pi = boost::math::constants::pi<double>();
  //const double tau = boost::math::constants::two_pi<double>();
  //const double angle = Fmod(pi - (std::atan2(dx,dy)),tau);
  //assert(angle >= 0.0 && "GetAngle must return a value between zero and two pi");
  //assert(angle < tau  && "GetAngle must return a value between zero and two pi");
  //return angle;
  //return ToDoubleSafe(GetAngle(Apfloat(dx),Apfloat(dy)));
  return GetAngleClockScreenImpl(dx,dy);
}

ribi::Geometry::Apfloat ribi::Geometry::GetAngleClockScreen(const Apfloat& dx, const Apfloat& dy) const noexcept
{
  return GetAngleClockScreenImpl(dx,dy);
  //return GetAngle(ApCoordinat2D(dx,dy));
}

double ribi::Geometry::GetAngleClockScreen(const Coordinat2D& p) const noexcept
{
  return GetAngleClockScreenImpl(boost::geometry::get<0>(p),boost::geometry::get<1>(p));

  //return ToDoubleSafe(GetAngle(ToApfloat(p)));
}

ribi::Geometry::Apfloat ribi::Geometry::GetAngleClockScreen(const ApCoordinat2D& p) const noexcept
{
  return GetAngleClockScreenImpl(boost::geometry::get<0>(p),boost::geometry::get<1>(p));
  /*
  const Apfloat pi = boost::math::constants::pi<double>();
  const Apfloat tau = boost::math::constants::two_pi<double>();
  const Apfloat dx = boost::geometry::get<0>(p);
  const Apfloat dy = boost::geometry::get<1>(p);
  const Apfloat part = ::atan2(dx,dy);
  const Apfloat angle = Fmod(pi - ::atan2(dx,dy),tau);
  assert(angle >= 0.0 && "GetAngle must return a value between zero and two pi");
  assert(angle < tau  && "GetAngle must return a value between zero and two pi");
  return angle;
  */
}

 

 

 

 

 

./CppGeometry/geometry_is_clockwise.cpp

 

#include "geometry.h"

#include <cassert>

#include "plane.h"
#include "trace.h"

bool ribi::Geometry::IsClockwise(const Apfloat& a, const Apfloat&b) const noexcept
{
  const Apfloat pi  = boost::math::constants::pi<double>();
  const Apfloat tau = boost::math::constants::two_pi<double>();
  const Apfloat c { Fmod(a,tau) };
  const Apfloat d { Fmod(b,tau) };
  assert(c >= 0);
  assert(c < tau);
  assert(d >= 0);
  assert(d < tau);
  if (d - c >= 0.0 && d - c < pi) return true;
  if (d + tau - c >= 0.0 && d + tau - c < pi) return true;
  return false;
}

bool ribi::Geometry::IsClockwise(const double a, const double b) const noexcept
{
  return IsClockwise(Apfloat(a),Apfloat(b));
}

bool ribi::Geometry::IsClockwise(const Apfloats& angles) const noexcept
{
  const int sz = static_cast<int>(angles.size());
  assert(sz >= 2 && "Need at least two angles to determine if these are clockwise");
  for (int i=0; i!=sz-1; ++i)
  {
    if (!IsClockwise(angles[i],angles[i+1])) return false;
  }

  //Calculate cumulative clockwise distance
  const Apfloat tau{boost::math::constants::two_pi<double>()};
  Apfloat sum{0.0};
  for (int i=0; i!=sz-1; ++i)
  {
    const Apfloat diff{angles[i+1] - angles[i]};
    if (diff > 0.0)
    {
      sum += diff;
    }
    else
    {
      assert(diff + tau > 0.0);
      sum += (diff + tau);
    }
  }
  return sum < tau;
}

bool ribi::Geometry::IsClockwise(const Doubles& angles) const noexcept
{
  return IsClockwise(ToApfloat(angles));
}

bool ribi::Geometry::IsClockwiseCartesian(
  const ApCoordinats3D& points,
  const ApCoordinat3D& observer
) const noexcept
{
  const bool verbose{false};
  const int n_points = static_cast<int>(points.size());
  assert(n_points == 3 || n_points == 4);
  if (n_points == 3)
  {
    const ApCoordinat3D a{points[0]};
    const ApCoordinat3D b{points[1]};
    const ApCoordinat3D c{points[2]};
    const ApCoordinat3D normal{CalcNormal(a,b,c)};
    const Apfloat direction{CalcDotProduct(normal,a - observer)};
    const bool is_clockwise { direction < Apfloat(0.0) }; //Difference between CW ('<') and CCW ('>')
    if (verbose)
    {
      TRACE(ToStr(a));
      TRACE(ToStr(b));
      TRACE(ToStr(c));
      TRACE(ToStr(normal));
      TRACE(ToStrSafe(direction));
      TRACE(is_clockwise);
    }
    return is_clockwise;
  }
  else
  {
    assert(n_points == 4);
    //See if the points in the projection are in the same direction
    assert(Geometry().IsPlane(points));
    const std::unique_ptr<Plane> plane(new Plane(points[0],points[1],points[2]));
    assert(plane);
    const auto v(
        plane->CalcProjection(
          {
            points[0],
            points[1],
            points[2],
            points[3]
          }
        )
    );
    //If the points are messed up, they cannot be clockwise
    if (!IsClockwiseCartesianHorizontal(v) && !IsCounterClockwiseCartesianHorizontal(v)) return false;
    //The neatly orderder points have the same winding as the first three
    std::remove_const<std::remove_reference<decltype(points)>::type>::type a;
    std::copy(points.begin() + 0,points.begin() + 3,std::back_inserter(a));
    return IsClockwiseCartesian(a,observer);
  }
}

bool ribi::Geometry::IsClockwiseCartesian(
  const Coordinats3D& points,
  const Coordinat3D& observer
) const noexcept
{
  return IsClockwiseCartesian(ToApfloat(points),ToApfloat(observer));
}

bool ribi::Geometry::IsClockwiseCartesianHorizontal(const ApCoordinats3D& points) const noexcept
{
  const bool verbose{false};
  using boost::geometry::get;
  const auto center = CalcCenter(points);
  Apfloats angles;
  angles.reserve(points.size());
  std::transform(points.begin(),points.end(),
    std::back_inserter(angles),
    [this,center](const ApCoordinat3D& coordinat)
    {
      return GetAngleClockCartesian(
        get<0>(coordinat) - get<0>(center),
        get<1>(coordinat) - get<1>(center)
      );
    }
  );
  const bool is_clockwise = IsClockwise(angles);
  if (verbose)
  {
    TRACE(ToStr(center));
    for (const auto point: points) { TRACE(ToStr(point)); }
    for (const auto angle: angles)
    {
      const Apfloat pi{boost::math::constants::pi<double>()};
      std::stringstream s;
      s << ToStrSafe(angle / pi) << " pi";
      const std::string angle_str{s.str()};
      TRACE(angle_str);
    }
    TRACE(is_clockwise);
  }
  return is_clockwise;
}

bool ribi::Geometry::IsClockwiseCartesianHorizontal(const Coordinats3D& points) const noexcept
{
  return IsClockwiseCartesianHorizontal(ToApfloat(points));
}

bool ribi::Geometry::IsClockwiseCartesianHorizontal(const ApCoordinats2D& points) const noexcept
{
  //Points are determined from their center
  const auto center = CalcCenter(points);
  Apfloats angles;
  std::transform(points.begin(),points.end(),
    std::back_inserter(angles),
    [this,center](const ApCoordinat2D& coordinat)
    {
      return GetAngleClockCartesian(
        coordinat.x() - center.x(),
        coordinat.y() - center.y()
      );
    }
  );
  return IsClockwise(angles);
}

 

 

 

 

 

./CppGeometry/geometry_is_counter_clockwise.cpp

 

#include "geometry.h"

#include <cassert>

#include "plane.h"
#include "trace.h"

bool ribi::Geometry::IsCounterClockwise(const Apfloat& a, const Apfloat& b) const noexcept
{
  return !IsClockwise(a,b);
}

bool ribi::Geometry::IsCounterClockwise(const double a, const double b) const noexcept
{
  return IsCounterClockwise(Apfloat(a),Apfloat(b));
}

bool ribi::Geometry::IsCounterClockwise(const Apfloats& angles) const noexcept
{
  const int sz = static_cast<int>(angles.size());
  assert(sz >= 2 && "Need at least two angles to determine if these are counter-clockwise");
  for (int i=0; i!=sz-1; ++i)
  {
    if (!IsCounterClockwise(angles[i],angles[i+1])) return false;
  }

  //Calculate cumulative clockwise distance
  const Apfloat tau{boost::math::constants::two_pi<double>()};
  Apfloat sum{0.0};
  for (int i=0; i!=sz-1; ++i)
  {
    const Apfloat diff{angles[i] - angles[i+1]}; //Order reversed compared to IsClockwise
    if (diff > 0.0)
    {
      sum += diff;
    }
    else
    {
      assert(diff + tau > 0.0);
      sum += (diff + tau);
    }
  }
  return sum < tau;
}

bool ribi::Geometry::IsCounterClockwise(const Doubles& angles) const noexcept
{
  return IsCounterClockwise(ToApfloat(angles));
}

bool ribi::Geometry::IsCounterClockwiseCartesian(
  const ApCoordinats3D& points,
  const ApCoordinat3D& observer
) const noexcept
{
  const bool verbose{false};
  const int n_points{static_cast<int>(points.size())};
  assert(n_points == 3 || n_points == 4);
  if (n_points == 3)
  {
    const auto a = points[0];
    const auto b = points[1];
    const auto c = points[2];
    const auto normal = CalcNormal(a,b,c);
    const Apfloat direction{CalcDotProduct(normal,a - observer)};
    const bool is_counter_clockwise{direction > Apfloat(0.0)}; //Difference between CW ('<') and CCW ('>')
    if (verbose)
    {
      TRACE(ToStr(a));
      TRACE(ToStr(b));
      TRACE(ToStr(c));
      TRACE(ToStr(normal));
      TRACE(ToStrSafe(direction));
      TRACE(is_counter_clockwise);
    }
    return is_counter_clockwise;
  }
  else
  {
    assert(n_points == 4);
    //See if the points in the projection are in the same direction
    assert(Geometry().IsPlane(points));
    const std::unique_ptr<Plane> plane(new Plane(points[0],points[1],points[2]));
    assert(plane);
    const auto v =
      plane->CalcProjection(
        {
          points[0],
          points[1],
          points[2],
          points[3]
        }
      )
    ;
    #ifndef NDEBUG
    if (verbose)
    {
      for (const auto& p: v) TRACE(ToStr(p));
      TRACE(IsClockwiseCartesianHorizontal(v));
      TRACE(IsCounterClockwiseCartesianHorizontal(v));
    }
    #endif
    //If the points are messed up, they cannot be clockwise
    if (!IsClockwiseCartesianHorizontal(v) && !IsCounterClockwiseCartesianHorizontal(v)) return false;
    //The neatly orderder point have the same winding as the first three
    std::remove_const<std::remove_reference<decltype(points)>::type>::type a;
    std::copy(points.begin() + 0,points.begin() + 3,std::back_inserter(a));
    return IsCounterClockwiseCartesian(a,observer);
  }
}

bool ribi::Geometry::IsCounterClockwiseCartesian(
  const Coordinats3D& points,
  const Coordinat3D& observer
) const noexcept
{
  return IsCounterClockwiseCartesian(ToApfloat(points),ToApfloat(observer));
}

bool ribi::Geometry::IsCounterClockwiseCartesianHorizontal(const ApCoordinats2D& points) const noexcept
{
  //Points are determined from their center
  const auto center(CalcCenter(points));
  Apfloats angles;
  angles.reserve(points.size());
  std::transform(points.begin(),points.end(),
    std::back_inserter(angles),
    [this,center](const ApCoordinat2D& coordinat)
    {
      return GetAngleClockCartesian(
        coordinat - center
      );
    }
  );
  return IsCounterClockwise(angles);
}

bool ribi::Geometry::IsCounterClockwiseCartesianHorizontal(const Coordinats2D& points) const noexcept
{
  return IsCounterClockwiseCartesianHorizontal(ToApfloat(points));
}

bool ribi::Geometry::IsCounterClockwiseCartesianHorizontal(const ApCoordinats3D& points3d) const noexcept
{
  ApCoordinats2D points2d;
  for (const auto& point3d: points3d)
  {
    points2d.push_back(
      {
        boost::geometry::get<0>(point3d),
        boost::geometry::get<1>(point3d)
      }
    );
  }
  return IsCounterClockwiseCartesianHorizontal(points2d);
}

bool ribi::Geometry::IsCounterClockwiseCartesianHorizontal(const Coordinats3D& points) const noexcept
{
  return IsCounterClockwiseCartesianHorizontal(ToApfloat(points));
}

 

 

 

 

 

./CppGeometry/geometry_test.cpp

 

#include "geometry.h"

#include <cassert>

#include "apcplx.h" //apfloat's atan2

#include "plane.h"
#include "ribi_regex.h"
#include "testtimer.h"
#include "trace.h"

#ifndef NDEBUG
void ribi::Geometry::Test() noexcept
{
  {
    static bool is_tested{false};
    if (is_tested) return;
    is_tested = true;
  }
  {
    const boost::shared_ptr<Plane> plane{
      new Plane(
        Plane::Coordinat3D(0.0,0.0,0.0),
        Plane::Coordinat3D(0.0,1.0,0.0),
        Plane::Coordinat3D(1.0,0.0,0.0)
      )
    };
  }
  ::ribi::Regex();
  const TestTimer test_timer(__func__,__FILE__,1.1);
  const bool verbose{false};
  const double pi = boost::math::constants::pi<double>();
  const Geometry g;
  if (verbose) { TRACE("CalcPlane #2"); }
  {
    //const TestTimer test_timer(boost::lexical_cast<std::string>(__LINE__),__FILE__,1.0);
    using boost::geometry::cs::cartesian;
    const double p1_x =  1.0;
    const double p1_y =  2.0;
    const double p1_z =  3.0;
    const double p2_x =  4.0;
    const double p2_y =  6.0;
    const double p2_z =  9.0;
    const double p3_x = 12.0;
    const double p3_y = 11.0;
    const double p3_z =  9.0;
    const Coordinat3D p1(p1_x,p1_y,p1_z);
    const Coordinat3D p2(p2_x,p2_y,p2_z);
    const Coordinat3D p3(p3_x,p3_y,p3_z);
    const auto t(g.CalcPlane(p1,p2,p3));
    const double a = t[0];
    const double b = t[1];
    const double c = t[2];
    const double d = t[3];
    const double d_p1_expected = (a * p1_x) + (b * p1_y) + (c * p1_z);
    const double d_p2_expected = (a * p2_x) + (b * p2_y) + (c * p2_z);
    const double d_p3_expected = (a * p3_x) + (b * p3_y) + (c * p3_z);
    const bool verbose{false};
    if (verbose)
    {
      std::clog
        << "(a * x) + (b * y) + (c * z) = d" << '\n'
        << "(" << a << " * x) + (" << b << " * y) + (" << c << " * z) = " << d << '\n'
        << "(" << a << " * " << p1_x << ") + (" << b << " * " << p1_y << ") + (" << c << " * " << p1_z << ") = " << d << '\n'
        << "(" << (a * p1_x) << ") + (" << (b * p1_y) << ") + (" << (c * p1_z) << ") = " << d << '\n'
        << "(" << a << " * " << p2_x << ") + (" << b << " * " << p2_y << ") + (" << c << " * " << p2_z << ") = " << d << '\n'
        << "(" << (a * p2_x) << ") + (" << (b * p2_y) << ") + (" << (c * p2_z) << ") = " << d << '\n'
        << "(" << a << " * " << p3_x << ") + (" << b << " * " << p3_y << ") + (" << c << " * " << p3_z << ") = " << d << '\n'
        << "(" << (a * p3_x) << ") + (" << (b * p3_y) << ") + (" << (c * p3_z) << ") = " << d << '\n'
      ;
      /* Screen output

      (a * x) + (b * y) + (c * z) = d
      (30 * x) + (-48 * y) + (17 * z) = -15
      (30 * 1) + (-48 * 2) + (17 * 3) = -15
      (30) + (-96) + (51) = -15
      (30 * 4) + (-48 * 6) + (17 * 9) = -15
      (120) + (-288) + (153) = -15
      (30 * 12) + (-48 * 11) + (17 * 9) = -15
      (360) + (-528) + (153) = -15

      */
    }
    assert(std::abs(d - d_p1_expected) < 0.001);
    assert(std::abs(d - d_p2_expected) < 0.001);
    assert(std::abs(d - d_p3_expected) < 0.001);
  }
  if (verbose) { TRACE("CalcPlane #2"); }
  {
    //CalcPlane return the coefficients in the following form:
    // A.x + B.y + C.z = D
    //Converting this to z being a function of x and y:
    // -C.z = A.x + B.y - D
    // z = -A/C.x - B/C.y + D/C
    //In this test, use the formula:
    //  z = (2.0 * x) + (3.0 * y) + (5.0)
    //Coefficients must then become:
    //  -A/C = 2.0
    //  -B/C = 3.0
    //   D/C = 5.0
    //Coefficients are, when setting C to 1.0:
    //  -A = 2.0 => A = -2.0
    //  -B = 3.0 => B = -3.0
    //   C = 1.0
    //   D = 5.0
    using boost::geometry::model::point;
    using boost::geometry::cs::cartesian;
    const Coordinat3D p1(1.0,1.0,10.0);
    const Coordinat3D p2(1.0,2.0,13.0);
    const Coordinat3D p3(2.0,1.0,12.0);
    const auto t(g.CalcPlane(p1,p2,p3));
    const double a = t[0];
    const double b = t[1];
    const double c = t[2];
    const double d = t[3];
    const double a_expected = -2.0;
    const double b_expected = -3.0;
    const double c_expected =  1.0;
    const double d_expected =  5.0;
    assert(std::abs(a - a_expected) < 0.001);
    assert(std::abs(b - b_expected) < 0.001);
    assert(std::abs(c - c_expected) < 0.001);
    assert(std::abs(d - d_expected) < 0.001);
    const double d_p1_expected = (a * 1.0) + (b * 1.0) + (c * 10.0);
    const double d_p2_expected = (a * 1.0) + (b * 2.0) + (c * 13.0);
    const double d_p3_expected = (a * 2.0) + (b * 1.0) + (c * 12.0);
    assert(std::abs(d - d_p1_expected) < 0.001);
    assert(std::abs(d - d_p2_expected) < 0.001);
    assert(std::abs(d - d_p3_expected) < 0.001);

  }
  if (verbose) { TRACE("Fmod"); }
  {
    const double expected_min = 1.0 - 0.00001;
    const double expected_max = 1.0 + 0.00001;
    assert(g.Fmod(3.0,2.0) > expected_min && g.Fmod(3.0,2.0) < expected_max);
    assert(g.Fmod(13.0,2.0) > expected_min && g.Fmod(13.0,2.0) < expected_max);
    assert(g.Fmod(-1.0,2.0) > expected_min && g.Fmod(-1.0,2.0) < expected_max);
    assert(g.Fmod(-3.0,2.0) > expected_min && g.Fmod(-3.0,2.0) < expected_max);
    assert(g.Fmod(-13.0,2.0) > expected_min && g.Fmod(-13.0,2.0) < expected_max);
  }
  if (verbose) { TRACE("std::atan2 versus apfloat's atan2"); }
  {
    for (double dx = -1.0; dx < 1.01; dx += 1.0)
    {
      for (double dy = -1.0; dy < 1.01; dy += 1.0)
      {
        if (dx == 0.0 && dy == 0.0) continue;
        const auto a = std::atan2(dx,dy);
        const auto b = atan2(Apfloat(dx),Apfloat(dy));
        const auto c = g.Atan2(Apfloat(dx),Apfloat(dy));
        const auto error = abs(a - c);
        assert(error < 0.01); //apfloat does not use namespace std::
      }
    }
  }

  if (verbose) { TRACE("GetAngleClockCartesian, double"); }
  {
    {
      const double angle =  g.GetAngleClockCartesian(0.0, 1.0); //North
      const double expected = 0.0 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(1.0, 1.0); //North-East
      const double expected = 0.25 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(1.0,0.0); //East
      const double expected = 0.5 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(1.0,-1.0); //South-East
      const double expected = 0.75 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(0.0,-1.0); //South
      const double expected = 1.0 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(-1.0,-1.0); //South-West
      const double expected = 1.25 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(-1.0,0.0); //West
      const double expected = 1.5 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockCartesian(-1.0,1.0); //North-West
      const double expected = 1.75 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
  }
  if (verbose) { TRACE("GetAngleClockCartesian, Apfloat"); }
  {
    const auto angle =  g.GetAngleClockCartesian(Apfloat(0.0),Apfloat(1.0)); //North
    const auto expected = 0.0 * pi;
    const auto error = abs(angle-expected); //apfloat does not use namespace std::
    assert(error < 0.01);
  }
  if (verbose) { TRACE("GetAngleClockCartesian, Apfloat, test eight directions"); }
  {
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(1.0),Apfloat(1.0)); //North-East
      const auto expected = 0.25 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(1.0),Apfloat(0.0)); //East
      const auto expected = 0.5 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(1.0),Apfloat(-1.0)); //South-East
      const auto expected = 0.75 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(0.0),Apfloat(-1.0)); //South
      const auto expected = 1.0 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(-1.0),Apfloat(-1.0)); //South-West
      const auto expected = 1.25 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(-1.0),Apfloat(0.0)); //West
      const auto expected = 1.5 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockCartesian(Apfloat(-1.0),Apfloat(1.0)); //North-West
      const auto expected = 1.75 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
  }
  if (verbose) { TRACE("GetAngleClockScreen, double, test eight directions"); }
  {
    {
      const double angle =  g.GetAngleClockScreen(0.0,-1.0); //North
      const double expected = 0.0 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(1.0,-1.0); //North-East
      const double expected = 0.25 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(1.0,0.0); //East
      const double expected = 0.5 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(1.0,1.0); //South-East
      const double expected = 0.75 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(0.0,1.0); //South
      const double expected = 1.0 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(-1.0,1.0); //South-West
      const double expected = 1.25 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(-1.0,0.0); //West
      const double expected = 1.5 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
    {
      const double angle =  g.GetAngleClockScreen(-1.0,-1.0); //North-West
      const double expected = 1.75 * pi;
      assert(std::abs(angle-expected) < 0.01);
    }
  }
  if (verbose) { TRACE("GetAngleClockScreen, Apfloat"); }
  {
    const auto angle =  g.GetAngleClockScreen(Apfloat(0.0),Apfloat(-1.0)); //North
    const auto expected = 0.0 * pi;
    const auto error = abs(angle-expected); //apfloat does not use namespace std::
    assert(error < 0.01);
  }
  if (verbose) { TRACE("GetAngleClockScreen, Apfloat, eight directions"); }
  {
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(1.0),Apfloat(-1.0)); //North-East
      const auto expected = 0.25 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(1.0),Apfloat(0.0)); //East
      const auto expected = 0.5 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(1.0),Apfloat(1.0)); //South-East
      const auto expected = 0.75 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(0.0),Apfloat(1.0)); //South
      const auto expected = 1.0 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(-1.0),Apfloat(1.0)); //South-West
      const auto expected = 1.25 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(-1.0),Apfloat(0.0)); //West
      const auto expected = 1.5 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
    {
      const auto angle =  g.GetAngleClockScreen(Apfloat(-1.0),Apfloat(-1.0)); //North-West
      const auto expected = 1.75 * pi;
      assert(abs(angle-expected) < 0.01); //apfloat does not use namespace std::
    }
  }
  if (verbose) { TRACE("GetDistance"); }
  {
    {
      const double distance = g.GetDistance(3.0,4.0);
      const double expected = 5.0;
      assert(std::abs(distance-expected) < 0.01);
    }
    {
      const double distance = g.GetDistance(-3.0,4.0);
      const double expected = 5.0;
      assert(std::abs(distance-expected) < 0.01);
    }
    {
      const double distance = g.GetDistance(3.0,-4.0);
      const double expected = 5.0;
      assert(std::abs(distance-expected) < 0.01);
    }
    {
      const double distance = g.GetDistance(-3.0,-4.0);
      const double expected = 5.0;
      assert(std::abs(distance-expected) < 0.01);
    }
  }
  if (verbose) TRACE("double -> apfloat -> double conversion");
  {
    for (const double x:
      {
        1.0,
        0.0,
        boost::math::constants::pi<double>(),
        boost::numeric::bounds<double>::smallest(),
        boost::numeric::bounds<double>::lowest(),
        boost::numeric::bounds<double>::highest()
      }
    )
    {
      const apfloat y = x;
      const double z = g.ToDouble(y);
      assert(x == z);
    }
  }
  if (verbose) TRACE("IsClockWise of two angles, doubles");
  {
    assert( g.IsClockwise(-2.5 * pi,-2.0 * pi));
    assert(!g.IsClockwise(-2.0 * pi,-2.5 * pi));
    assert( g.IsClockwise(-1.5 * pi,-1.0 * pi));
    assert(!g.IsClockwise(-1.0 * pi,-1.5 * pi));
    assert( g.IsClockwise(-0.5 * pi, 0.0 * pi));
    assert(!g.IsClockwise( 0.0 * pi,-0.5 * pi));
    assert( g.IsClockwise( 0.0 * pi, 0.5 * pi));
    assert(!g.IsClockwise( 0.5 * pi, 0.0 * pi));
    assert( g.IsClockwise( 0.5 * pi, 1.0 * pi));
    assert(!g.IsClockwise( 1.0 * pi, 0.5 * pi));
    assert( g.IsClockwise( 1.5 * pi, 2.0 * pi));
    assert(!g.IsClockwise( 2.0 * pi, 1.5 * pi));
    assert( g.IsClockwise( 2.5 * pi, 3.0 * pi));
    assert(!g.IsClockwise( 3.0 * pi, 2.5 * pi));
  }
  if (verbose) TRACE("IsCounterClockWise of two angles, doubles");
  {
    assert(!g.IsCounterClockwise(-2.5 * pi,-2.0 * pi));
    assert( g.IsCounterClockwise(-2.0 * pi,-2.5 * pi));
    assert(!g.IsCounterClockwise(-1.5 * pi,-1.0 * pi));
    assert( g.IsCounterClockwise(-1.0 * pi,-1.5 * pi));
    assert(!g.IsCounterClockwise(-0.5 * pi, 0.0 * pi));
    assert( g.IsCounterClockwise( 0.0 * pi,-0.5 * pi));
    assert(!g.IsCounterClockwise( 0.0 * pi, 0.5 * pi));
    assert( g.IsCounterClockwise( 0.5 * pi, 0.0 * pi));
    assert(!g.IsCounterClockwise( 0.5 * pi, 1.0 * pi));
    assert( g.IsCounterClockwise( 1.0 * pi, 0.5 * pi));
    assert(!g.IsCounterClockwise( 1.5 * pi, 2.0 * pi));
    assert( g.IsCounterClockwise( 2.0 * pi, 1.5 * pi));
    assert(!g.IsCounterClockwise( 2.5 * pi, 3.0 * pi));
    assert( g.IsCounterClockwise( 3.0 * pi, 2.5 * pi));
  }
  if (verbose) TRACE("IsClockWise of two, vector");
  {
    assert( g.IsClockwise(Doubles( {-2.5 * pi,-2.0 * pi } )));
    assert(!g.IsClockwise(Doubles( {-2.0 * pi,-2.5 * pi } )));
    assert( g.IsClockwise(Doubles( {-1.5 * pi,-1.0 * pi } )));
    assert(!g.IsClockwise(Doubles( {-1.0 * pi,-1.5 * pi } )));
    assert( g.IsClockwise(Doubles( {-0.5 * pi, 0.0 * pi } )));
    assert(!g.IsClockwise(Doubles( { 0.0 * pi,-0.5 * pi } )));
    assert( g.IsClockwise(Doubles( { 0.0 * pi, 0.5 * pi } )));
    assert(!g.IsClockwise(Doubles( { 0.5 * pi, 0.0 * pi } )));
    assert( g.IsClockwise(Doubles( { 0.5 * pi, 1.0 * pi } )));
    assert(!g.IsClockwise(Doubles( { 1.0 * pi, 0.5 * pi } )));
    assert( g.IsClockwise(Doubles( { 1.5 * pi, 2.0 * pi } )));
    assert(!g.IsClockwise(Doubles( { 2.0 * pi, 1.5 * pi } )));
    assert( g.IsClockwise(Doubles( { 2.5 * pi, 3.0 * pi } )));
    assert(!g.IsClockwise(Doubles( { 3.0 * pi, 2.5 * pi } )));
  }
  if (verbose) TRACE("IsCounterClockWise of two, vector");
  {
    assert(!g.IsCounterClockwise(Doubles( {-2.5 * pi,-2.0 * pi } )));
    assert( g.IsCounterClockwise(Doubles( {-2.0 * pi,-2.5 * pi } )));
    assert(!g.IsCounterClockwise(Doubles( {-1.5 * pi,-1.0 * pi } )));
    assert( g.IsCounterClockwise(Doubles( {-1.0 * pi,-1.5 * pi } )));
    assert(!g.IsCounterClockwise(Doubles( {-0.5 * pi, 0.0 * pi } )));
    assert( g.IsCounterClockwise(Doubles( { 0.0 * pi,-0.5 * pi } )));
    assert(!g.IsCounterClockwise(Doubles( { 0.0 * pi, 0.5 * pi } )));
    assert( g.IsCounterClockwise(Doubles( { 0.5 * pi, 0.0 * pi } )));
    assert(!g.IsCounterClockwise(Doubles( { 0.5 * pi, 1.0 * pi } )));
    assert( g.IsCounterClockwise(Doubles( { 1.0 * pi, 0.5 * pi } )));
    assert(!g.IsCounterClockwise(Doubles( { 1.5 * pi, 2.0 * pi } )));
    assert( g.IsCounterClockwise(Doubles( { 2.0 * pi, 1.5 * pi } )));
    assert(!g.IsCounterClockwise(Doubles( { 2.5 * pi, 3.0 * pi } )));
    assert( g.IsCounterClockwise(Doubles( { 3.0 * pi, 2.5 * pi } )));
  }
  if (verbose) TRACE("IsClockWise of three, vector");
  {

    assert( g.IsClockwise(Doubles( {-2.5 * pi,-2.0 * pi,-1.8 * pi} ))); //CW
    assert(!g.IsClockwise(Doubles( {-1.8 * pi,-2.0 * pi,-2.5 * pi} ))); //CCW
    assert(!g.IsClockwise(Doubles( {-2.0 * pi,-2.5 * pi,-1.8 * pi} ))); //Mess

    assert( g.IsClockwise(Doubles( {-1.5 * pi,-1.0 * pi,-0.8 * pi} ))); //CW
    assert(!g.IsClockwise(Doubles( {-0.8 * pi,-1.0 * pi,-1.5 * pi} ))); //CCW
    assert(!g.IsClockwise(Doubles( {-1.0 * pi,-1.5 * pi,-0.8 * pi} ))); //Mess

    assert( g.IsClockwise(Doubles( {-0.5 * pi, 0.0 * pi, 0.3 * pi } ))); //CW
    assert(!g.IsClockwise(Doubles( { 0.3 * pi, 0.0 * pi,-0.5 * pi } ))); //CCW
    assert(!g.IsClockwise(Doubles( { 0.0 * pi,-0.5 * pi, 0.3 * pi } ))); //Mess

    assert( g.IsClockwise(Doubles( { 0.0 * pi, 0.5 * pi, 0.8 * pi } ))); //CW
    assert(!g.IsClockwise(Doubles( { 0.8 * pi, 0.5 * pi, 0.0 * pi } ))); //CCW
    assert(!g.IsClockwise(Doubles( { 0.5 * pi, 0.0 * pi, 0.8 * pi } ))); //Mess

    assert( g.IsClockwise(Doubles( { 0.5 * pi, 1.0 * pi, 1.3 * pi } ))); //CW
    assert(!g.IsClockwise(Doubles( { 1.3 * pi, 1.0 * pi, 0.5 * pi } ))); //CCW
    assert(!g.IsClockwise(Doubles( { 1.0 * pi, 0.5 * pi, 1.3 * pi } ))); //Mess

    assert( g.IsClockwise(Doubles( { 1.5 * pi, 2.0 * pi, 2.3 * pi } ))); //CW
    assert(!g.IsClockwise(Doubles( { 2.3 * pi, 2.0 * pi, 1.5 * pi } ))); //CCW
    assert(!g.IsClockwise(Doubles( { 2.0 * pi, 1.5 * pi, 2.3 * pi } ))); //Mess

    assert( g.IsClockwise(Doubles( { 2.5 * pi, 3.0 * pi, 3.3 * pi } ))); //CW
    assert(!g.IsClockwise(Doubles( { 3.3 * pi, 3.0 * pi, 2.5 * pi } ))); //CCW
    assert(!g.IsClockwise(Doubles( { 3.0 * pi, 2.5 * pi, 3.3 * pi } ))); //Mess
  }
  if (verbose) TRACE("IsCounterClockWise of three, vector");
  {

    assert(!g.IsCounterClockwise(Doubles( {-2.5 * pi,-2.0 * pi,-1.8 * pi} ))); //CW
    assert( g.IsCounterClockwise(Doubles( {-1.8 * pi,-2.0 * pi,-2.5 * pi} ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( {-2.0 * pi,-2.5 * pi,-1.8 * pi} ))); //Mess

    assert(!g.IsCounterClockwise(Doubles( {-1.5 * pi,-1.0 * pi,-0.8 * pi} ))); //CW
    assert( g.IsCounterClockwise(Doubles( {-0.8 * pi,-1.0 * pi,-1.5 * pi} ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( {-1.0 * pi,-1.5 * pi,-0.8 * pi} ))); //Mess

    assert(!g.IsCounterClockwise(Doubles( {-0.5 * pi, 0.0 * pi, 0.3 * pi } ))); //CW
    assert( g.IsCounterClockwise(Doubles( { 0.3 * pi, 0.0 * pi,-0.5 * pi } ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( { 0.0 * pi,-0.5 * pi, 0.3 * pi } ))); //Mess

    assert(!g.IsCounterClockwise(Doubles( { 0.0 * pi, 0.5 * pi, 0.8 * pi } ))); //CW
    assert( g.IsCounterClockwise(Doubles( { 0.8 * pi, 0.5 * pi, 0.0 * pi } ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( { 0.5 * pi, 0.0 * pi, 0.8 * pi } ))); //Mess

    assert(!g.IsCounterClockwise(Doubles( { 0.5 * pi, 1.0 * pi, 1.3 * pi } ))); //CW
    assert( g.IsCounterClockwise(Doubles( { 1.3 * pi, 1.0 * pi, 0.5 * pi } ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( { 1.0 * pi, 0.5 * pi, 1.3 * pi } ))); //Mess

    assert(!g.IsCounterClockwise(Doubles( { 1.5 * pi, 2.0 * pi, 2.3 * pi } ))); //CW
    assert( g.IsCounterClockwise(Doubles( { 2.3 * pi, 2.0 * pi, 1.5 * pi } ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( { 2.0 * pi, 1.5 * pi, 2.3 * pi } ))); //Mess

    assert(!g.IsCounterClockwise(Doubles( { 2.5 * pi, 3.0 * pi, 3.3 * pi } ))); //CW
    assert( g.IsCounterClockwise(Doubles( { 3.3 * pi, 3.0 * pi, 2.5 * pi } ))); //CCW
    assert(!g.IsCounterClockwise(Doubles( { 3.0 * pi, 2.5 * pi, 3.3 * pi } ))); //Mess
  }
  if (verbose) TRACE("IsClockWise of four, vector");
  {

    //Correct order
    assert( g.IsClockwise(Doubles( {-2.5 * pi,-2.0 * pi,-1.5 * pi,-1.0 * pi} )));
    assert( g.IsClockwise(Doubles( {-2.0 * pi,-1.5 * pi,-1.0 * pi,-0.5 * pi} )));
    assert( g.IsClockwise(Doubles( {-1.5 * pi,-1.0 * pi,-0.5 * pi, 0.0 * pi} )));
    assert( g.IsClockwise(Doubles( {-1.0 * pi,-0.5 * pi, 0.0 * pi, 0.5 * pi} )));
    assert( g.IsClockwise(Doubles( {-0.5 * pi, 0.0 * pi, 0.5 * pi, 1.0 * pi} )));
    assert( g.IsClockwise(Doubles( { 0.0 * pi, 0.5 * pi, 1.0 * pi, 1.5 * pi} )));
    assert( g.IsClockwise(Doubles( { 0.5 * pi, 1.0 * pi, 1.5 * pi, 2.0 * pi} )));
    assert( g.IsClockwise(Doubles( { 1.0 * pi, 1.5 * pi, 2.0 * pi, 2.5 * pi} )));
    assert( g.IsClockwise(Doubles( { 1.5 * pi, 2.0 * pi, 2.5 * pi, 3.0 * pi} )));
    assert( g.IsClockwise(Doubles( { 2.0 * pi, 2.5 * pi, 3.0 * pi, 3.5 * pi} )));
    assert( g.IsClockwise(Doubles( { 2.5 * pi, 3.0 * pi, 3.5 * pi, 4.0 * pi} )));

    //Swap [0] and [1]
    assert(!g.IsClockwise(Doubles( {-2.0 * pi,-2.5 * pi,-1.5 * pi,-1.0 * pi} )));
    assert(!g.IsClockwise(Doubles( {-1.5 * pi,-2.0 * pi,-1.0 * pi,-0.5 * pi} )));
    assert(!g.IsClockwise(Doubles( {-1.0 * pi,-1.5 * pi,-0.5 * pi, 0.0 * pi} )));
    assert(!g.IsClockwise(Doubles( {-0.5 * pi,-1.0 * pi, 0.0 * pi, 0.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 0.0 * pi,-0.5 * pi, 0.5 * pi, 1.0 * pi} )));
    assert(!g.IsClockwise(Doubles( { 0.5 * pi, 0.0 * pi, 1.0 * pi, 1.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 1.0 * pi, 0.5 * pi, 1.5 * pi, 2.0 * pi} )));
    assert(!g.IsClockwise(Doubles( { 1.5 * pi, 1.0 * pi, 2.0 * pi, 2.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 2.0 * pi, 1.5 * pi, 2.5 * pi, 3.0 * pi} )));
    assert(!g.IsClockwise(Doubles( { 2.5 * pi, 2.0 * pi, 3.0 * pi, 3.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 3.0 * pi, 2.5 * pi, 3.5 * pi, 4.0 * pi} )));

    //Swap [2] and [3]
    assert(!g.IsClockwise(Doubles( {-2.5 * pi,-2.0 * pi,-1.0 * pi,-1.5 * pi} )));
    assert(!g.IsClockwise(Doubles( {-2.0 * pi,-1.5 * pi,-0.5 * pi,-1.0 * pi} )));
    assert(!g.IsClockwise(Doubles( {-1.5 * pi,-1.0 * pi, 0.0 * pi,-0.5 * pi} )));
    assert(!g.IsClockwise(Doubles( {-1.0 * pi,-0.5 * pi, 0.5 * pi, 0.0 * pi} )));
    assert(!g.IsClockwise(Doubles( {-0.5 * pi, 0.0 * pi, 1.0 * pi, 0.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 0.0 * pi, 0.5 * pi, 1.5 * pi, 1.0 * pi} )));
    assert(!g.IsClockwise(Doubles( { 0.5 * pi, 1.0 * pi, 2.0 * pi, 1.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 1.0 * pi, 1.5 * pi, 2.5 * pi, 2.0 * pi} )));
    assert(!g.IsClockwise(Doubles( { 1.5 * pi, 2.0 * pi, 3.0 * pi, 2.5 * pi} )));
    assert(!g.IsClockwise(Doubles( { 2.0 * pi, 2.5 * pi, 3.5 * pi, 3.0 * pi} )));
    assert(!g.IsClockwise(Doubles( { 2.5 * pi, 3.0 * pi, 4.0 * pi, 3.5 * pi} )));
  }
  if (verbose) TRACE("IsCounterClockWise of four, vector");
  {

    //Clockwise order
    assert(!g.IsCounterClockwise(Doubles( {-2.5 * pi,-2.0 * pi,-1.5 * pi,-1.0 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( {-2.0 * pi,-1.5 * pi,-1.0 * pi,-0.5 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( {-1.5 * pi,-1.0 * pi,-0.5 * pi, 0.0 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( {-1.0 * pi,-0.5 * pi, 0.0 * pi, 0.5 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( {-0.5 * pi, 0.0 * pi, 0.5 * pi, 1.0 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( { 0.0 * pi, 0.5 * pi, 1.0 * pi, 1.5 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( { 0.5 * pi, 1.0 * pi, 1.5 * pi, 2.0 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( { 1.0 * pi, 1.5 * pi, 2.0 * pi, 2.5 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( { 1.5 * pi, 2.0 * pi, 2.5 * pi, 3.0 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( { 2.0 * pi, 2.5 * pi, 3.0 * pi, 3.5 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( { 2.5 * pi, 3.0 * pi, 3.5 * pi, 4.0 * pi} )));
  }

  if (verbose) TRACE("IsClockWise of four, vector, going round more than once");
  /*
       A
       | _-D
     _-+-
    C  |
     \  |
      \ |
       \|
        B

  */
  {

    assert(!g.IsClockwise(Doubles( {0.0 * pi,0.9 * pi,1.8 * pi,2.7 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( {0.0 * pi,0.9 * pi,1.8 * pi,2.7 * pi} )));
  }
  if (verbose) TRACE("IsCounterClockWise of four, vector, going round more than once");
  /*
       D
       | _-A
     _-+-
    B  |
     \  |
      \ |
       \|
        C

  */
  {

    assert(!g.IsCounterClockwise(Doubles( {2.0 * pi,1.1 * pi,0.2 * pi,-0.7 * pi} )));
    assert(!g.IsCounterClockwise(Doubles( {2.7 * pi,1.8 * pi,0.9 * pi,0.0 * pi} )));
  }
  if (verbose) { TRACE("IsClockwise, 3D, from #228, test for positive, three of four points"); }
  {
    /*
    (-3.78624,2,10)
    (-3.78624,2,0)
    (-0.55,2,0)
    Removed: (-0.55,2,10)
    Observer: (-2.1871,3.74169,5)
    */

    const Coordinats3D points {
      Coordinat3D(-3.78624,2,10),
      Coordinat3D(-3.78624,2,0),
      Coordinat3D(-0.55,2,0)
    };
    const Coordinat3D observer{-2.1871,3.74169,5};
    assert(g.IsClockwiseCartesian(points,observer));
  }

  if (verbose) { TRACE("IsCounterClockwise, 3D, from #228, test for positive, three of four points"); }
  {
    /*
    (-0.55,2,0)
    (-3.78624,2,0)
    (-3.78624,2,10)
    Removed: (-0.55,2,10)
    Observer: (-2.1871,3.74169,5)
    */

    const Coordinats3D points {
      Coordinat3D(-0.55,2,0),
      Coordinat3D(-3.78624,2,0),
      Coordinat3D(-3.78624,2,10)
    };
    const Coordinat3D observer{-2.1871,3.74169,5};
    assert(g.IsCounterClockwiseCartesian(points,observer));
  }
  if (verbose) { TRACE("IsCounterClockwise, 3D, from #228, test for positive, all four points"); }
  {
    /*
    (-0.55,2,0)
    (-3.78624,2,0)
    (-3.78624,2,10)
    (-0.55,2,10)
    Observer: (-2.1871,3.74169,5)
    */

    const Coordinats3D points {
      Coordinat3D(-0.55,2,0),
      Coordinat3D(-3.78624,2,0),
      Coordinat3D(-3.78624,2,10),
      Coordinat3D(-0.55,2,10)
    };
    const Coordinat3D observer{-2.1871,3.74169,5};
    assert(g.IsCounterClockwiseCartesian(points,observer));
  }
  if (verbose) TRACE("IsConvex, two dimensions");
  {
    /* Polygons used:

    0123456789012    0123456789012
   0+------------   0+------------
   1|               1|
   2| A---------B   2| A---------B
   3| E        /    3|  \   D   /
   4|  \      /     4|   \ / \ /
   5|   D----C      5|    E   C
    |                |

       Convex           Concave

    */


    //Convex shape
    {
      const std::vector<Coordinat2D> points {
        { 2.0, 2.0}, //A
        {12.0, 2.0}, //B
        { 9.0, 5.0}, //C
        { 4.0, 5.0}, //D
        { 2.0, 3.0}  //E
      };

      Polygon polygon;
      boost::geometry::append(polygon, points);
      assert(g.IsConvex(polygon));
    }
    //Concave shape
    {
      const std::vector<Coordinat2D> points {
        { 2.0, 2.0}, //A
        {12.0, 2.0}, //B
        { 9.0, 5.0}, //C
        { 7.0, 3.0}, //D
        { 5.0, 5.0}  //E
      };

      Polygon polygon;
      boost::geometry::append(polygon, points);
      assert(!g.IsConvex(polygon));
    }
  }
  if (verbose) TRACE("Convex shape, 2D, from error #1");
  {

    /*

       __--3
     2-
     |
     |
     |
     |  __-0
     1--


    */


    const std::vector<Coordinat2D> points {
      {17.4971,33.0765},
      {9.28854,29.5636},
      {9.28854,40.6764},
      {17.4971,44.4009}
    };
    assert(g.IsConvex(points));
  }
  if (verbose) TRACE("Convex shape, 2D, from error #2");
  {
    //From 3D points:
    /*
    // (2.35114,3.23607,5) (index: 644)'
    // (1.17557,2.35781,5) (index: 658)'
    // (2.35114,3.23607,6) (index: 668)'
    // (1.17557,2.35781,6) (index: 682)'
    */
    /*

       __--2
     3-   /
         /
        /
       /
      / __-0
     1--

    */


    const std::vector<Coordinat2D> points {
      {17.497 ,33.0765},
      { 9.2885,29.5636},
      {17.497 ,44.4009},
      { 9.2885,40.6764}
    };
    assert(!g.IsConvex(points));
  }


  if (verbose) TRACE("Convex shape, 2D, from error #3, point 0");
  {
    /*
        __--2
     1--    |
     |      |
     |      |
     |      |
     |      |
     |      3
     0

    */


    const std::vector<Coordinat2D> points {
      { 9.2885,29.5636},
      { 9.2885,40.6764},
      {17.497 ,44.4009},
      {17.497 ,33.0765}
    };
    assert(g.IsConvex(points));
  }

  if (verbose) TRACE("Convex shape, 2D, from error #3, point 5");
  {
    /*
            3
     0      |
     |      |
     |      |
     |      |
     |      |
     |  __--2
     1--

    */


    const std::vector<Coordinat2D> points {
      { 9.2885,40.6764},
      { 9.2885,29.5636},
      {17.497 ,33.0765},
      {17.497 ,44.4009}
    };
    assert(g.IsConvex(points));
  }

  if (verbose) TRACE("Convex shape, 2D, from error #3, point 5");
  {
    /*

    2-3
    |
    1-0

    */


    const std::vector<Coordinat2D> points {
      {15.9835,631.923},
      {8.24075,628.579},
      {8.24075,679.58 },
      {15.9835,682.926}
    };
    assert(g.IsConvex(points));
  }


  if (verbose) TRACE("Convex shape, 3D, points in Y=0 plane");
  {

    if (verbose) TRACE("Convex shape, 3D, points in Y=0 plane, Z shape");
    {
      /*
                  |####/
        3---2     |###/#
           /      +##+##
          /       |#/###
         /        |/####
        1---0   --+--+--
                 /|
      */
      const ApCoordinats3D points {
        ApCoordinat3D(1.0,0.0,0.0),
        ApCoordinat3D(0.0,0.0,0.0),
        ApCoordinat3D(1.0,0.0,1.0),
        ApCoordinat3D(0.0,0.0,1.0)
      };

      assert(!g.IsConvex(points) && "This is an hourglass shape, so it is not convex");
    }
    if (verbose) TRACE("Convex shape, 3D, points in Y=0 plane, C shape");
    {
      /*

        2---3 Z=1.0
        |
        |
        |
        1---0 Z=0.0

      */
      const ApCoordinats3D points {
        ApCoordinat3D(1.0,0.0,0.0),
        ApCoordinat3D(0.0,0.0,0.0),
        ApCoordinat3D(0.0,0.0,1.0),
        ApCoordinat3D(1.0,0.0,1.0)
      };
      assert(g.IsConvex(points) && "This is a corrected hourglass shape, so it is convex");
    }
  }
  if (verbose) TRACE("Convex shape, 3D, points in X=Y plane");
  {

    {
      /*
                  |    |/     _-
                  |  _-2    _-
        3---2     |_- /   _-
           /      3  /  _-
          /       | / _0
         /        |/_-
        1---0   --1-----
                 /|
      */
      const ApCoordinats3D points {
        ApCoordinat3D(1.0,1.0,0.0),
        ApCoordinat3D(0.0,0.0,0.0),
        ApCoordinat3D(1.0,1.0,1.0),
        ApCoordinat3D(0.0,0.0,1.0)
      };
      assert(!g.IsConvex(points) && "This is an hourglass shape, so it is not convex");
    }
    //Convex shape, 3D, points in X=Y plane
    {
      /*

                  |    |/     _-
                  |  _-3    _-
        2---3     |_- /   _-
        |         2  /  _-
        |         | / _0
        |         |/_-
        1---0   --1-----
                 /|

      */
      const ApCoordinats3D points {
        ApCoordinat3D(1.0,1.0,0.0),
        ApCoordinat3D(0.0,0.0,0.0),
        ApCoordinat3D(0.0,0.0,1.0),
        ApCoordinat3D(1.0,1.0,1.0)
      };
      assert(g.IsConvex(points) && "This is a corrected hourglass shape, so it is convex");
    }
  }
  if (verbose) TRACE("Convex shape, 3D, points in 2X=Y plane, not at origin");
  {

    {
      const ApCoordinats3D points {
        ApCoordinat3D(2.0,4.0,0.0),
        ApCoordinat3D(1.0,1.0,0.0),
        ApCoordinat3D(2.0,4.0,1.0),
        ApCoordinat3D(1.0,1.0,1.0)
      };
      assert(!g.IsConvex(points) && "This is an hourglass shape, so it is not convex");
    }
    {
      const ApCoordinats3D points {
        ApCoordinat3D(2.0,2.0,0.0),
        ApCoordinat3D(1.0,1.0,0.0),
        ApCoordinat3D(1.0,1.0,1.0),
        ApCoordinat3D(2.0,2.0,1.0)
      };
      assert(g.IsConvex(points) && "This is a corrected hourglass shape, so it is convex");
    }
  }
  if (verbose) TRACE("Convex shape, 3D, points in 2X=Y plane, above Z=0");
  {
    {
      const ApCoordinats3D points {
        ApCoordinat3D(2.0,4.0,1.0),
        ApCoordinat3D(1.0,1.0,1.0),
        ApCoordinat3D(2.0,4.0,2.0),
        ApCoordinat3D(1.0,1.0,2.0)
      };
      assert(!g.IsConvex(points) && "This is an hourglass shape, so it is not convex");
    }
    {
      const ApCoordinats3D points {
        ApCoordinat3D(2.0,2.0,1.0),
        ApCoordinat3D(1.0,1.0,1.0),
        ApCoordinat3D(1.0,1.0,2.0),
        ApCoordinat3D(2.0,2.0,2.0)
      };
      assert(g.IsConvex(points) && "This is a corrected hourglass shape, so it is convex");
    }
  }
  if (verbose) TRACE("Convex shape, 3D, from error");
  {

    {
      /*

        3---2 Z=6.0
           /
          /
         /
        1---0 Z=5.0

      */
      const ApCoordinats3D points {
        ApCoordinat3D(2.35114,3.23607,5.0),
        ApCoordinat3D(1.17557,2.35781,5.0),
        ApCoordinat3D(2.35114,3.23607,6.0),
        ApCoordinat3D(1.17557,2.35781,6.0)
      };
      assert(!g.IsConvex(points) && "This is an hourglass shape, so it is not convex");
    }
    {
      /*

        2---3 Z=6.0
        |
        |
        |
        1---0 Z=5.0

      */
      const ApCoordinats3D points {
        ApCoordinat3D(2.35114,3.23607,5.0),
        ApCoordinat3D(1.17557,2.35781,5.0),
        ApCoordinat3D(1.17557,2.35781,6.0),
        ApCoordinat3D(2.35114,3.23607,6.0)
      };
      assert(g.IsConvex(points) && "This is a corrected hourglass shape, so it is convex");
    }
  }

  if (verbose) TRACE("IsClockwise in three dimensions, points in XY plane");
  {
    /*

       Y
     3 |
       |
     2 | C-B
       | |/
     1 | A       (seen from above, Z = 1)
       |
     0 +------X
       0 1 2 3

    */


    const ApCoordinat3D a(1.0, 1.0, 1.0);
    const ApCoordinat3D b(2.0, 2.0, 1.0);
    const ApCoordinat3D c(1.0, 2.0, 1.0);
    const ApCoordinats3D coordinats { a,b,c };
    for (int i = 0; i!=3; ++i)
    {
      for (int j=i; j!=3; ++j)
      {
        const double x = static_cast<double>(i-1) * 5.0;
        const double y = static_cast<double>(j-1) * 5.0;
        const ApCoordinat3D above(x, y, 2.0);
        const ApCoordinat3D below(x, y,-1.0);
        assert(!g.IsClockwiseCartesian(coordinats,above));
        assert( g.IsClockwiseCartesian(coordinats,below));
        assert( g.IsCounterClockwiseCartesian(coordinats,above));
        assert(!g.IsCounterClockwiseCartesian(coordinats,below));
      }
    }
  }
  if (verbose) TRACE("IsClockwise in three dimensions, points in XY plane");
  {
    /*

       Y
     3 |
       |
     2 | B-C
       | |/
     1 | A       (seen from above, Z = 1)
       |
     0 +------X
       0 1 2 3

    */

    const ApCoordinat3D a(1.0, 1.0, 1.0);
    const ApCoordinat3D b(1.0, 2.0, 1.0);
    const ApCoordinat3D c(2.0, 2.0, 1.0);
    const ApCoordinats3D coordinats { a,b,c };
    for (int i = 0; i!=3; ++i)
    {
      for (int j=i; j!=3; ++j)
      {
        const double x = static_cast<double>(i-1) * 5.0;
        const double y = static_cast<double>(j-1) * 5.0;
        const ApCoordinat3D above(x, y, 2.0);
        const ApCoordinat3D below(x, y,-1.0);
        assert( g.IsClockwiseCartesian(coordinats,above));
        assert(!g.IsClockwiseCartesian(coordinats,below));
        assert(!g.IsCounterClockwiseCartesian(coordinats,above));
        assert( g.IsCounterClockwiseCartesian(coordinats,below));
      }
    }
  }


  if (verbose) TRACE("CalcCenter");
  {
    /*

   3.5 Y
     3 |
   2.5 |
     2 | D-C
   1.5 | | |
     1 | A-B     (seen from above, Z = 1)
   0.5 |
     0 +------X
       0 1 2 3

    */

    const ApCoordinat3D a(1.0, 1.0, 1.0);
    const ApCoordinat3D b(2.0, 1.0, 1.0);
    const ApCoordinat3D c(2.0, 2.0, 1.0);
    const ApCoordinat3D d(1.0, 2.0, 1.0);
    const ApCoordinats3D coordinats { a,b,c,d };
    const auto center = g.CalcCenter(coordinats);
    const auto center_x = boost::geometry::get<0>(center);
    const auto center_y = boost::geometry::get<1>(center);
    const auto center_z = boost::geometry::get<2>(center);
    assert(center_x > 1.49 && center_x < 1.51);
    assert(center_y > 1.49 && center_y < 1.51);
    assert(center_z > 0.99 && center_z < 1.01);
  }

  if (verbose) TRACE("IsClockwise, clockwise, in three dimensions, four points in XY plane");
  {
    /*

   3.5 Y
     3 |
   2.5 |
     2 | D-C
   1.5 | | |
     1 | A-B     (seen from above, Z = 1)
   0.5 |
     0 +------X
       0 1 2 3

    */

    const ApCoordinat3D a(1.0, 1.0, 1.0);
    const ApCoordinat3D b(2.0, 1.0, 1.0);
    const ApCoordinat3D c(2.0, 2.0, 1.0);
    const ApCoordinat3D d(1.0, 2.0, 1.0);
    const ApCoordinats3D coordinats { a,b,c,d };
    for (int i = 1; i!=3; ++i)
    {
      for (int j=i; j!=3; ++j)
      {
        const double x = static_cast<double>(i-1) * 5.0;
        const double y = static_cast<double>(j-1) * 5.0;
        const ApCoordinat3D above(x, y, 2.0);
        const ApCoordinat3D below(x, y,-1.0);
        assert(!g.IsClockwiseCartesianHorizontal(coordinats));
        assert( g.IsCounterClockwiseCartesianHorizontal(coordinats));
        assert(!g.IsClockwiseCartesian(coordinats,above));
        assert( g.IsClockwiseCartesian(coordinats,below));
        assert( g.IsCounterClockwiseCartesian(coordinats,above));
        assert(!g.IsCounterClockwiseCartesian(coordinats,below));
      }
    }
  }
  if (verbose) TRACE("IsClockwise, counter-clockwise, in three dimensions, four points in XY plane");
  {
    /*
   3.5 Y
     3 |
   2.5 |
     2 | B-C
   1.5 | | |
     1 | A-D     (seen from above, Z = 1)
   0.5 |
     0 +------X
       0 1 2 3

    */

    const ApCoordinat3D a(1.0, 1.0, 1.0);
    const ApCoordinat3D b(1.0, 2.0, 1.0);
    const ApCoordinat3D c(2.0, 2.0, 1.0);
    const ApCoordinat3D d(2.0, 1.0, 1.0);
    const ApCoordinats3D coordinats { a,b,c,d };
    for (int i = 1; i!=3; ++i)
    {
      for (int j=i; j!=3; ++j)
      {
        const double x = static_cast<double>(i-1) * 5.0;
        const double y = static_cast<double>(j-1) * 5.0;
        const ApCoordinat3D above(x, y, 2.0);
        const ApCoordinat3D below(x, y,-1.0);
        assert( g.IsClockwiseCartesian(coordinats,above));
        assert(!g.IsClockwiseCartesian(coordinats,below));
        assert(!g.IsCounterClockwiseCartesian(coordinats,above));
        assert( g.IsCounterClockwiseCartesian(coordinats,below));
      }
    }
  }
  if (verbose) TRACE("IsClockwise, indetermined direction, in three dimensions, four points in XY plane");
  {
    /*
   3.5 Y
     3 |
   2.5 |
     2 | C-B
   1.5 |  x
     1 | A-D     (seen from above, Z = 1)
   0.5 |
     0 +------X
       0 1 2 3

    */

    const ApCoordinat3D a(1.0, 1.0, 1.0);
    const ApCoordinat3D b(2.0, 2.0, 1.0);
    const ApCoordinat3D c(1.0, 2.0, 1.0);
    const ApCoordinat3D d(2.0, 1.0, 1.0);
    const ApCoordinats3D coordinats { a,b,c,d };
    for (int i = 1; i!=3; ++i)
    {
      for (int j=i; j!=3; ++j)
      {
        const Apfloat x = static_cast<double>(i-1) * 5.0;
        const Apfloat y = static_cast<double>(j-1) * 5.0;
        const ApCoordinat3D above(x, y, 2.0);
        const ApCoordinat3D below(x, y,-1.0);
        assert(!g.IsClockwiseCartesian(coordinats,above));
        assert(!g.IsClockwiseCartesian(coordinats,below));
        assert(!g.IsCounterClockwiseCartesian(coordinats,above));
        assert(!g.IsCounterClockwiseCartesian(coordinats,below));
      }
    }
  }
  //IsPlane
  {

    const ApCoordinat3D a(-3.64472,-0.25,0.0);
    const ApCoordinat3D b(-4.52988,-0.25,0.0);
    const ApCoordinat3D c(-3.64472,-0.25,10.0);
    const ApCoordinat3D d(-4.52988,-0.25,10.0);
    const ApCoordinats3D v{a,b,c,d};
    assert(Geometry().IsPlane(v));
  }
  if (verbose) TRACE("Translate");
  {

    //Translate (0.5,1.0) to origin
    const auto house = g.CreateHouseShape();
    const auto translated_house = g.Translate(house,-0.5,-1.0);
    boost::geometry::model::ring<Coordinat2D> translated_points;
    boost::geometry::convert(translated_house,translated_points);
    const std::vector<Coordinat2D> translated_points_expected {
      { 0.0, 1.0}, //0
      { 0.5, 0.0}, //1
      { 0.5,-1.0}, //2
      {-0.5,-1.0}, //3
      {-0.5, 0.0}  //4
    };
    assert(
      std::equal(
        translated_points.begin(),translated_points.end(),
        translated_points_expected.begin(),
        [](const Coordinat2D& a,const Coordinat2D& b)
        {
          return boost::geometry::equals(a,b);
        }
      )
      && "Points must be translated as expected"
    );
  }
  if (verbose) TRACE("Scale up twice, from origin");
  {
    //Scale twice up, from origin
    const auto house = g.CreateHouseShape();
    const auto rescaled_house = g.Rescale(house,2.0);
    boost::geometry::model::ring<Coordinat2D> rescaled_points;
    boost::geometry::convert(rescaled_house,rescaled_points);
    const std::vector<Coordinat2D> rescaled_points_expected {
      {1.0, 4.0}, //0
      {2.0, 2.0}, //1
      {2.0, 0.0}, //2
      {0.0, 0.0}, //3
      {0.0, 2.0}  //4
    };
    assert(
      std::equal(
        rescaled_points.begin(),rescaled_points.end(),
        rescaled_points_expected.begin(),
        [](const Coordinat2D& a,const Coordinat2D& b)
        {
          //std::cout << "(" << a.x() << "," << a.y() << ")-("
          //  << b.x() << "," << b.y() << ")" << std::endl;
          return boost::geometry::equals(a,b);
        }
      )
      && "Points must be rescaled as expected"
    );
  }
  if (verbose) TRACE("Scale up twice, from non-origin");
  {
    //Scale up, from center at (0.5,1.0)
    /*

      3-
       |
      2-    0
       |   / \
      1-  4 * 1
       |  |   |
      0-  3---2
       |
       +--|---|---|
          0   1   2

    */

    const auto house = g.CreateHouseShape();
    const auto rescaled_house = g.Rescale(house,2.0,0.5,1.0);
    boost::geometry::model::ring<Coordinat2D> rescaled_points;
    boost::geometry::convert(rescaled_house,rescaled_points);
    const std::vector<Coordinat2D> rescaled_points_expected {
      { 0.5, 3.0}, //0
      { 1.5, 1.0}, //1
      { 1.5,-1.0}, //2
      {-0.5,-1.0}, //3
      {-0.5, 1.0}  //4
    };
    assert(
      std::equal(
        rescaled_points.begin(),rescaled_points.end(),
        rescaled_points_expected.begin(),
        [](const Coordinat2D& a,const Coordinat2D& b)
        {
          //std::cout << "(" << a.x() << "," << a.y() << ")-("
          //  << b.x() << "," << b.y() << ")" << std::endl;
          return boost::geometry::equals(a,b);
        }
      )
      && "Points must be rescaled as expected"
    );
  }
  if (verbose) { TRACE("WktToSvg"); }
  {

    const std::string s
      = g.WktToSvg("POLYGON((0 0,0 1,1 0)),LINESTRING(0 0,0 1,1 0)",1.0)
    ;
    assert(!s.empty());
  }

  if (verbose) { TRACE("GetLineLineIntersections"); }
  #ifdef TODO_RICHEL
  {
    typedef boost::geometry::model::d2::point_xy<double> Point;
    typedef boost::geometry::model::linestring<Point> Line;
    typedef boost::geometry::model::box<Point> Rect;
    //Assume line segment (0,0)-(2,2) intersects line segment (2,0)-(0,2) at point (1,1)
    {
      const Line line1 = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(2.0,2.0) } ));
      const Line line2 = CreateLine( std::vector<Point>( { Point(2.0,0.0), Point(0.0,2.0) } ));
      assert( GetLineLineIntersections(line1,line2).size() == 1);
      assert( fuzzy_equal_to()(GetLineLineIntersections(line1,line2)[0].x(),1.0) );
      assert( fuzzy_equal_to()(GetLineLineIntersections(line1,line2)[0].y(),1.0) );
    }
    //Lines are finite, however, as the line segment
    //(0,0)-(0.2,0.2) does not intersect line segment (2,0)-(0,2) at point (1,1)
    {
      const Line line1 = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(0.2,0.2) } ));
      const Line line2 = CreateLine( std::vector<Point>( { Point(2.0,0.0), Point(0.0,2.0) } ));
      assert( GetLineLineIntersections(line1,line2).size() == 0);
    }
    //Assume line segment (0,0)-(2,2) intersects with rectangle (1,0)-(3,9) at point (1,1)
    {
      const Line line = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(2.0,2.0) } ));
      const Rect rect(Point(1.0,0.0), Point(3.0,3.9));
      GetLineRectIntersections(line,rect);
      const std::vector<Point> v = GetLineRectIntersections(line,rect);
      assert(v.size() == 1);
      assert( fuzzy_equal_to()(v[0].x(),1.0) );
      assert( fuzzy_equal_to()(v[0].y(),1.0) );
    }
    //Assume line segment (0,0)-(2,2) intersects with rectangle (0,1)-(3,9) at point (1,1)
    {
      const Line line = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(2.0,2.0) } ));
      const Rect rect(Point(0.0,1.0), Point(3.0,9.0));
      GetLineRectIntersections(line,rect);
      const std::vector<Point> v = GetLineRectIntersections(line,rect);
      assert(v.size() == 1);
      assert( fuzzy_equal_to()(v[0].x(),1.0) );
      assert( fuzzy_equal_to()(v[0].y(),1.0) );
    }
    //Assume line segment (0,0)-(2,2) intersects with rectangle (1,1)-(3,3) at point (1,1) once
    {
      const Line line = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(2.0,2.0) } ));
      const Rect rect(Point(1.0,1.0), Point(3.0,3.0));
      GetLineRectIntersections(line,rect);
      const std::vector<Point> v = GetLineRectIntersections(line,rect);
      assert(v.size() == 1);
      assert( fuzzy_equal_to()(v[0].x(),1.0) );
      assert( fuzzy_equal_to()(v[0].y(),1.0) );
    }

    //Assume line segment (0,0)-(4,4) intersects with rectangle (1,0)-(3,9) at points (1,1) and (3,3)
    {
      const Line line = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(4.0,4.0) } ));
      const Rect rect(Point(1.0,0.0), Point(3.0,3.9));
      GetLineRectIntersections(line,rect);
      const std::vector<Point> v = GetLineRectIntersections(line,rect);
      assert(v.size() == 2);
      assert( fuzzy_equal_to()(v[0].x(),1.0) );
      assert( fuzzy_equal_to()(v[0].y(),1.0) );
      assert( fuzzy_equal_to()(v[1].x(),3.0) );
      assert( fuzzy_equal_to()(v[1].y(),3.0) );
    }
    //Assume line segment (0,0)-(4,4) intersects with rectangle (0,1)-(3,9) at points (1,1) and (3,3)
    {
      const Line line = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(4.0,4.0) } ));
      const Rect rect(Point(0.0,1.0), Point(3.0,9.0));
      GetLineRectIntersections(line,rect);
      const std::vector<Point> v = GetLineRectIntersections(line,rect);
      assert(v.size() == 2);
      assert( fuzzy_equal_to()(v[0].x(),1.0) );
      assert( fuzzy_equal_to()(v[0].y(),1.0) );
      assert( fuzzy_equal_to()(v[1].x(),3.0) );
      assert( fuzzy_equal_to()(v[1].y(),3.0) );
    }
    //Assume line segment (0,0)-(4,4) intersects with rectangle (1,1)-(3,3) at points (1,1) and (3,3)
    {
      const Line line = CreateLine( std::vector<Point>( { Point(0.0,0.0), Point(4.0,4.0) } ));
      const Rect rect(Point(1.0,1.0), Point(3.0,3.0));
      GetLineRectIntersections(line,rect);
      const std::vector<Point> v = GetLineRectIntersections(line,rect);
      assert(v.size() == 2);
      assert( fuzzy_equal_to()(v[0].x(),1.0) );
      assert( fuzzy_equal_to()(v[0].y(),1.0) );
      assert( fuzzy_equal_to()(v[1].x(),3.0) );
      assert( fuzzy_equal_to()(v[1].y(),3.0) );
    }
  }
  #endif //#ifdef TODO_RICHEL
  //WktToLinestring: open linestring to polygon
  {

    const auto p = g.WktToLinestring("LINESTRING(0 0 0 1 1 1 1 0)");
    assert(p.size() == 4);
  }
  //WktToPolygon: open linestring to polygon
  {

    const auto p = g.WktToPolygon("POLYGON((0 0 0 1 1 1 1 0))");
    assert(p.outer().size() == 4);
  }
  //ToPolygon: open linestring to polygon
  {
    /*

    |        |
    + 2-3    + 2-3
    | |   -> | | |
    + 1-0    + 1-0
    |        |
    +-+-+-   +-+-+-

    */

    const auto l = g.WktToLinestring("LINESTRING(2 1 1 1 1 2 2 2)");
    const auto p = g.ToPolygon(l);
    assert(p.outer().size() == 4);
  }
  //ToPolygon: closed linestring to polygon
  {
    /*

    |        |
    + 2-3    + 2-3
    | | | -> | | |
    + 1-0    + 1-0
    |        |
    +-+-+-   +-+-+-

    */

    const auto l = g.WktToLinestring("LINESTRING(2 1 1 1 1 2 2 2 2 1)");
    const auto p = g.ToPolygon(l);
    assert(p.outer().size() == 4);
  }
}
#endif

 

 

 

 

 

Go back to Richel Bilderbeek's C++ page.

Go back to Richel Bilderbeek's homepage.

 

Valid XHTML 1.0 Strict

This page has been created by the tool CodeToHtml