Go back to Richel Bilderbeek's homepage.

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

 

 

 

 

 

(C++) Plane

 

Technical facts

 

 

 

 

 

 

./CppPlane/CppPlane.pri

 

INCLUDEPATH += \
    ../../Classes/CppPlane

SOURCES += \
    ../../Classes/CppPlane/plane.cpp \
    ../../Classes/CppPlane/planez.cpp \
    ../../Classes/CppPlane/planex.cpp \
    ../../Classes/CppPlane/planey.cpp \
    ../../Classes/CppPlane/plane_test.cpp \
    ../../Classes/CppPlane/planez_test.cpp \
    ../../Classes/CppPlane/planey_test.cpp \
    ../../Classes/CppPlane/planex_test.cpp

HEADERS  += \
    ../../Classes/CppPlane/plane.h \
    ../../Classes/CppPlane/planez.h \
    ../../Classes/CppPlane/planex.h \
    ../../Classes/CppPlane/planey.h

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

 

 

 

 

 

./CppPlane/plane.h

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.htm
//---------------------------------------------------------------------------
#ifndef RIBI_PLANE_H
#define RIBI_PLANE_H

#include <vector>
#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 <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/make_shared.hpp>
#include <boost/shared_ptr.hpp>
#ifndef _WIN32
#include <boost/geometry/geometries/polygon.hpp>
#endif

#include "apfloat.h"
#pragma GCC diagnostic pop

namespace ribi {

struct PlaneX;
struct PlaneY;
struct PlaneZ;

///Any 3D plane, even a single point
///Can be constructed from its equation and at least three 3D points
//A plane stores its 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
struct Plane
{
  typedef apfloat Double;
  typedef boost::geometry::model::d2::point_xy<apfloat> Coordinat2D;
  typedef boost::geometry::model::point<apfloat,3,boost::geometry::cs::cartesian> Coordinat3D;
  typedef std::vector<Coordinat2D> Coordinats2D;
  typedef std::vector<Coordinat3D> Coordinats3D;
  typedef std::vector<Double> Doubles;

  ///Construct a Plane from three points
  /*

   |
   |   /
   +--+
   |\/|
   |/\|
--+--+---------
  /|
/ |

  */
  explicit Plane(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  ) noexcept;

  ///Get the 2D projection of these 3D points,
  ///Assumes these are in a Plane
  /*

  A: (0,0,1)                  A: (0,0)
  B: (1,0,0)                  B: (SQRT(2),0)
  C: (1,1,0)                  C: (SQRT(2),SQRT(2))

  |    /
  |   /                       |
  A-----C                     |  C
  |\/  /      -> becomes ->   | /|
  |/\ /                       |/ |
  +--B---                     A--B-----

  */
  Coordinats2D CalcProjection(const Coordinats3D& points) const;

  ///If the Plane can be expressed as X = A*Y + B*Z + C, return the X
  Double CalcX(const Double& y, const Double& z) const;

  ///If the Plane can be expressed as Y = A*X + B*Z + C, return the Y
  Double CalcY(const Double& y, const Double& z) const;

  ///If the Plane can be expressed as Z = A*X + B*Y + C, return the Z
  Double CalcZ(const Double& y, const Double& z) const;

  ///Can the Plane be expressed as X = A*Y + B*Z + C ?
  bool CanCalcX() const noexcept;

  ///Can the Plane be expressed as Y = A*X + B*Z + C ?
  bool CanCalcY() const noexcept;

  ///Can the Plane be expressed as Z = A*X + B*Y + C ?
  bool CanCalcZ() const noexcept;

  static boost::shared_ptr<PlaneX> CreatePlaneX(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  ) noexcept;

  static boost::shared_ptr<PlaneY> CreatePlaneY(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  ) noexcept;

  static boost::shared_ptr<PlaneZ> CreatePlaneZ(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  ) noexcept;

  ///Calculates the error between plane and coordinat
  Double CalcError(const Coordinat3D& coordinat) const noexcept;

  ///Calculates the maximum allowed error for that coordinat for it to be in the plane
  Double CalcMaxError(const Coordinat3D& coordinat) const noexcept;

  ///If the Plane can be expressed as X = A*Y + B*Z + C, return the coefficients
  Doubles GetCoefficientsX() const;

  ///If the Plane can be expressed as Y = A*X + B*Z + C, return the coefficients
  Doubles GetCoefficientsY() const;

  ///If the Plane can be expressed as Z = A*X + B*Y + C, return the coefficients
  Doubles GetCoefficientsZ() const;

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

  ///Checks if the coordinat is in the plane
  bool IsInPlane(const Coordinat3D& coordinat) const noexcept;


  private:
  ~Plane() noexcept {}

  ///A non-horizontal plane; a plane that can be expressed as 'X(Y,Z) = A*Y + B*Z + C'
  const boost::shared_ptr<const PlaneX> m_plane_x;

  ///A non-horizontal plane; a plane that can be expressed as 'Y(X,Z) = A*X + B*Z + C'
  const boost::shared_ptr<const PlaneY> m_plane_y;

  ///A non-vertical plane; a plane that can be expressed as 'Z(X,Y) = A*X + B*Y + C'
  const boost::shared_ptr<const PlaneZ> m_plane_z;

  const Coordinats3D m_points;

  static boost::shared_ptr<PlaneX> CreatePlaneX(const Doubles& coefficients_x) noexcept;
  static boost::shared_ptr<PlaneY> CreatePlaneY(const Doubles& coefficients_y) noexcept;
  static boost::shared_ptr<PlaneZ> CreatePlaneZ(const Doubles& coefficients_z) noexcept;

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

  friend void boost::checked_delete<>(      Plane*);
  friend void boost::checked_delete<>(const Plane*);
  friend struct std::default_delete<      Plane>;
  friend struct std::default_delete<const Plane>;
  friend class boost::detail::sp_ms_deleter<      Plane>;
  friend class boost::detail::sp_ms_deleter<const Plane>;
  friend std::ostream& operator<<(std::ostream& os, const Plane& plane) noexcept;
};

std::ostream& operator<<(std::ostream& os, const Plane& plane) noexcept;

} //~namespace ribi

#endif // RIBI_PLANE_H

 

 

 

 

 

./CppPlane/plane.cpp

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.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 "plane.h"

#include <cassert>

#include "container.h"
#include "geometry.h"
#include "planex.h"
#include "planey.h"
#include "planez.h"
#include "trace.h"
#pragma GCC diagnostic pop

ribi::Plane::Plane(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) noexcept
: m_plane_x(CreatePlaneX(p1,p2,p3)),
  m_plane_y(CreatePlaneY(p1,p2,p3)),
  m_plane_z(CreatePlaneZ(p1,p2,p3)),
  m_points( {p1,p2,p3} )
{
  #ifndef NDEBUG
  Test();
  
  assert(Geometry().IsEqual3d(m_points[0],p1));
  assert(Geometry().IsEqual3d(m_points[1],p2));
  assert(Geometry().IsEqual3d(m_points[2],p3));
  #endif

  if (m_plane_z)
  {
    try
    {
      m_plane_z->GetFunctionA();
      m_plane_z->GetFunctionB();
      m_plane_z->GetFunctionC();
    }
    catch (...)
    {
      assert(!"Should not get here");
    }
  }

}

apfloat ribi::Plane::CalcError(const Coordinat3D& coordinat) const noexcept
{
  const bool verbose{false};
  const apfloat x = boost::geometry::get<0>(coordinat);
  const apfloat y = boost::geometry::get<1>(coordinat);
  const apfloat z = boost::geometry::get<2>(coordinat);

  apfloat min_error = std::numeric_limits<double>::max();

  if (verbose)
  {
    std::stringstream s;
    s << std::setprecision(99) << (*this);
    TRACE(s.str());
    if (CanCalcX()) { TRACE(Container().ToStr(GetCoefficientsX())); }
    if (CanCalcY()) { TRACE(Container().ToStr(GetCoefficientsY())); }
    if (CanCalcZ()) { TRACE(Container().ToStr(GetCoefficientsZ())); }
  }

  //Absolute member function
  if (CanCalcX())
  {
    const auto expected = x;
    const apfloat calculated = CalcX(y,z);
    const apfloat error = abs(calculated - expected);
    min_error = std::min(error,min_error);
  }
  if (CanCalcY())
  {
    const auto expected = y;
    const auto calculated = CalcY(x,z);
    const auto error = abs(calculated - expected);
    min_error = std::min(error,min_error);
  }
  if (CanCalcZ())
  {
    const auto expected = z;
    const auto calculated = CalcZ(x,y);
    const auto error = abs(calculated - expected);
    min_error = std::min(error,min_error);
  }
  assert(min_error >= 0.0);
  return min_error;
}

apfloat ribi::Plane::CalcMaxError(const Coordinat3D& coordinat) const noexcept
{
  apfloat max_error{std::numeric_limits<double>::denorm_min()};
  if (CanCalcX())
  {
    max_error = std::max(max_error,m_plane_x->CalcMaxError(coordinat));
  }
  if (CanCalcY())
  {
    max_error = std::max(max_error,m_plane_y->CalcMaxError(coordinat));
  }
  if (CanCalcZ())
  {
    max_error = std::max(max_error,m_plane_z->CalcMaxError(coordinat));
  }
  return max_error;
}


ribi::Plane::Coordinats2D ribi::Plane::CalcProjection(
  const Coordinats3D& points
) const
{
  if (!m_plane_x && !m_plane_y && !m_plane_z)
  {
    throw std::logic_error("Plane::CalcProjection: cannot express any plane");
  }
  try { if (m_plane_x) { return m_plane_x->CalcProjection(points); }}
  catch (std::logic_error&) { /* OK, try next plane */ }

  try { if (m_plane_y) { return m_plane_y->CalcProjection(points); }}
  catch (std::logic_error&) { /* OK, try next plane */ }

  try { if (m_plane_z) { return m_plane_z->CalcProjection(points); }}
  catch (std::logic_error&) { /* OK, try next plane */ }

  TRACE("ERROR");
  TRACE("INITIAL POINTS");
  for (const auto& point: m_points)
  {
    std::stringstream s;
    s
      << "("
      << boost::geometry::get<0>(point) << ","
      << boost::geometry::get<1>(point) << ","
      << boost::geometry::get<2>(point)
      << ")"
    ;
    TRACE(s.str());
  }

  TRACE(points.size());
  {
    if (m_plane_x)
    {
      try { TRACE(*m_plane_x); } catch(std::logic_error&) { TRACE("Failed m_plane_x"); }
      try { m_plane_x->CalcProjection(points); } catch (std::logic_error&) { TRACE("Failed m_plane_x->CalcProjection"); }
    }
  }
  {
    if (m_plane_y)
    {
      try { TRACE(*m_plane_y); } catch(std::logic_error&) { TRACE("Failed m_plane_y"); }
      try { m_plane_y->CalcProjection(points); } catch (std::logic_error&) { TRACE("Failed m_plane_y->CalcProjection"); }
    }
  }
  {
    if (m_plane_z)
    {
      try { TRACE(*m_plane_z); } catch(std::logic_error&) { TRACE("Failed m_plane_z"); }
      try { m_plane_z->CalcProjection(points); } catch (std::logic_error&) { TRACE("Failed m_plane_z->CalcProjection"); }
    }
  }
  for (const auto& point: points)
  {
    std::stringstream s;
    s
      << "("
      << boost::geometry::get<0>(point) << ","
      << boost::geometry::get<1>(point) << ","
      << boost::geometry::get<2>(point)
      << ")"
    ;
    TRACE(s.str());
  }
  assert(!"Should not get here");
  throw std::logic_error("Plane::CalcProjection: unexpected behavior");
}

ribi::Plane::Double ribi::Plane::CalcX(const Double& y, const Double& z) const
{
  if (!CanCalcX())
  {
    throw std::logic_error("Plane::CalcX: cannot express the plane as 'X = A*Y + B*Z + C'");
  }
  return m_plane_x->CalcX(y,z);
}

ribi::Plane::Double ribi::Plane::CalcY(const ribi::Plane::Double& x, const ribi::Plane::Double& z) const
{
  if (!CanCalcY())
  {
    throw std::logic_error("Plane::CalcY: cannot express the plane as 'Y = A*X + B*Y + C'");
  }
  return m_plane_y->CalcY(x,z);
}

ribi::Plane::Double ribi::Plane::CalcZ(const ribi::Plane::Double& x, const ribi::Plane::Double& y) const
{
  if (!CanCalcZ())
  {
    throw std::logic_error("Plane::CalcZ: cannot express the plane as 'Z = A*X + B*Y + C'");
  }
  return m_plane_z->CalcZ(x,y);
}

bool ribi::Plane::CanCalcX() const noexcept
{
  return m_plane_x.get();
}
bool ribi::Plane::CanCalcY() const noexcept
{
  return m_plane_y.get();
}

bool ribi::Plane::CanCalcZ() const noexcept
{
  return m_plane_z.get();
}

boost::shared_ptr<ribi::PlaneX> ribi::Plane::CreatePlaneX(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) noexcept
{
  const bool verbose{false};
  try
  {
    const boost::shared_ptr<PlaneX> p(
      boost::make_shared<PlaneX>(p1,p2,p3)
    );
    assert(p);
    return p;
  }
  catch (std::exception& e)
  {
    if (verbose) { TRACE(e.what()); }
    return boost::shared_ptr<PlaneX>();
  }
}

boost::shared_ptr<ribi::PlaneY> ribi::Plane::CreatePlaneY(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) noexcept
{
  try
  {
    const boost::shared_ptr<PlaneY> p
      = boost::make_shared<PlaneY>(p1,p2,p3);
    assert(p);
    return p;
  }
  catch (std::exception&)
  {
    return boost::shared_ptr<PlaneY>();
  }
}

boost::shared_ptr<ribi::PlaneZ> ribi::Plane::CreatePlaneZ(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) noexcept
{
  try
  {
    const boost::shared_ptr<PlaneZ> p
      = boost::make_shared<PlaneZ>(p1,p2,p3);
    assert(p);
    return p;
  }
  catch (std::exception&)
  {
    return boost::shared_ptr<PlaneZ>();
  }
}

std::vector<apfloat> ribi::Plane::GetCoefficientsX() const
{
  if (!m_plane_x)
  {
    throw std::logic_error("Plane::GetCoefficientsX: cannot express the plane as 'X = A*Y + B*Z + C'");
  }
  return m_plane_x->GetCoefficients();
}

std::vector<apfloat> ribi::Plane::GetCoefficientsY() const
{
  if (!m_plane_y)
  {
    throw std::logic_error("Plane::GetCoefficientsY: cannot express the plane as 'Y = A*X + B*Z + C'");
  }
  return m_plane_y->GetCoefficients();
}

std::vector<apfloat> ribi::Plane::GetCoefficientsZ() const
{
  if (!m_plane_z)
  {
    throw std::logic_error("Plane::GetCoefficientsZ: cannot express the plane as 'Z = A*X + B*Y + C'");
  }
  return m_plane_z->GetCoefficients();
}


std::string ribi::Plane::GetVersion() noexcept
{
  return "1.9";
}

std::vector<std::string> ribi::Plane::GetVersionHistory() noexcept
{
  return {
    "2014-03-07: version 1.0: initial version",
    "2014-03-10: version 1.1: allow vertical planes",
    "2014-03-13: version 1.2: bug fixed",
    "2014-04-01: version 1.3: use of std::unique_ptr",
    "2014-06-13: version 1.4: added operator<<, ToStr calls operator<<, shortened time to compile",
    "2014-06-16: version 1.5: improved detection of planes that can be expressed in less than three dimensions"
    "2014-07-03: version 1.6: use of apfloat, improved accuracy",
    "2014-07-10: version 1.7: use of apfloat only",
    "2014-07-15: version 1.8: multiple bugfixes",
    "2014-08-02: version 1.9: use of stubs, to speed up testing"
  };
}

bool ribi::Plane::IsInPlane(const Coordinat3D& coordinat) const noexcept
{
  //return CalcError(coordinat) <= CalcMaxError(coordinat);;
  assert(m_plane_x || m_plane_y || m_plane_z);
  const bool is_in_plane {
       (m_plane_x && m_plane_x->IsInPlane(coordinat))
    || (m_plane_y && m_plane_y->IsInPlane(coordinat))
    || (m_plane_z && m_plane_z->IsInPlane(coordinat))
  };
  #ifndef NDEBUG
  const bool has_error_below_max = CalcError(coordinat) <= CalcMaxError(coordinat);
  if (is_in_plane != has_error_below_max)
  {
    TRACE("ERROR");
    TRACE(is_in_plane);
    TRACE(has_error_below_max);
    TRACE(CalcError(coordinat));
    TRACE(CalcMaxError(coordinat));
    TRACE("BREAK");
  }
  assert(is_in_plane == has_error_below_max);
  #endif
  return is_in_plane;
}

std::ostream& ribi::operator<<(std::ostream& os, const Plane& plane) noexcept
{
  os << '(';
  const auto n_points = static_cast<int>(plane.m_points.size());
  for (/* const */ auto i=0; i!=n_points; ++i)
  {
    os << Geometry().ToStr(plane.m_points[i]);
    os << (i != n_points - 1 ? ',' : ')');
  }
  os << ',';
  if (plane.m_plane_x)
  {
    try { os << (*plane.m_plane_x); }
    catch (std::exception&) { os << "divnull"; }
  }
  else
  {
    os << "null";
  }
  os << ',';
  if (plane.m_plane_y)
  {
    try { os << (*plane.m_plane_y); }
    catch (std::exception&) { os << "divnull"; }
  }
  else
  {
    os << "null";
  }
  os << ',';
  if (plane.m_plane_z)
  {
    try { os << (*plane.m_plane_z); }
    catch (std::exception&) { os << "divnull"; }
  }
  else
  {
    os << "null";
  }
  return os;
}

 

 

 

 

 

./CppPlane/plane_test.cpp

 

#include "plane.h"

#include <cassert>

#include <boost/limits.hpp>
#include <boost/make_shared.hpp>

#include "container.h"
#include "geometry.h"
#include "planex.h"
#include "planey.h"
#include "planez.h"
#include "testtimer.h"
#include "trace.h"

#ifndef NDEBUG
void ribi::Plane::Test() noexcept
{
  {
    static bool is_tested{false};
    if (is_tested) return;
    is_tested = true;
  }
  Container();
  Geometry();
  { const auto planex = boost::make_shared<PlaneX>(); assert(planex); }
  { const auto planey = boost::make_shared<PlaneY>(); assert(planey); }
  { const auto planez = boost::make_shared<PlaneZ>(); assert(planez); }
  const TestTimer test_timer(__func__,__FILE__,1.0);
  const bool verbose{false};
  const bool show_warning{false};
  using boost::geometry::get;

  //Sorted by difficulty
  const auto series = PlaneZ::GetTestSeries();

  if (verbose) TRACE("Plane that can be expressed in all three forms");
  {
    const Coordinat3D p1( 1.0, 2.0,3.0);
    const Coordinat3D p2( 4.0, 6.0,9.0);
    const Coordinat3D p3(12.0,11.0,9.0);
    const Plane p(p1,p2,p3);
    assert(p.CanCalcX());
    assert(p.CanCalcY());
    assert(p.CanCalcZ());
    assert(
      !p.CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      ).empty()
    );

  }
  if (verbose) TRACE("Plane X = 2");
  {
    const Coordinat3D p1(2.0, 2.0,3.0);
    const Coordinat3D p2(2.0, 6.0,9.0);
    const Coordinat3D p3(2.0,11.0,9.0);
    const Plane p(p1,p2,p3);
    assert(p.CanCalcX());
    assert(!p.CanCalcY());
    assert(!p.CanCalcZ());
    assert(
      !p.CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      ).empty()
    );
  }
  if (verbose) TRACE("Plane X = 123");
  {
    const Coordinat3D p1(123.0, 2.0,3.0);
    const Coordinat3D p2(123.0, 6.0,9.0);
    const Coordinat3D p3(123.0,11.0,9.0);
    const Plane p(p1,p2,p3);
    assert(p.CanCalcX());
    assert(!p.CanCalcY());
    assert(!p.CanCalcZ());
    assert(
      !p.CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      ).empty()
    );
  }
  if (verbose) TRACE("Plane Y = 3");
  {
    const Coordinat3D p1( 2.0, 3.0, 5.0);
    const Coordinat3D p2( 7.0, 3.0, 9.0);
    const Coordinat3D p3(11.0,3.0,13.0);
    const Plane p(p1,p2,p3);

    assert(p.CanCalcY());
    assert(!p.CanCalcX());
    assert(!p.CanCalcZ());
    assert(
      !p.CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      ).empty()
    );
  }
  if (verbose) TRACE("Plane Z = 5");
  {
    const Coordinat3D p1( 2.0, 3.0,5.0);
    const Coordinat3D p2( 7.0,11.0,5.0);
    const Coordinat3D p3(13.0,17.0,5.0);
    const Plane p(p1,p2,p3);
    assert(p.CanCalcZ());
    assert(!p.CanCalcX());
    assert(!p.CanCalcY());
    assert(
      !p.CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      ).empty()
    );
  }

  //IsInPlane for Z=0 plane
  if (verbose) TRACE("CanCalcZ and IsInPlane, Z = 0 plane, from 1.0 coordinat");
  {
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,1.0,0.0);
    const Coordinat3D p3(1.0,0.0,0.0);
    const Plane p(p1,p2,p3);
    assert(!p.CanCalcX());
    assert(!p.CanCalcY());
    assert( p.CanCalcZ());

    assert( p.IsInPlane(Coordinat3D( 0.0, 0.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0, 1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,0.0)));
  }
  if (verbose) TRACE("CanCalcZ and IsInPlane, Z = 0 plane, from smallest possible coordinat");
  {
    const double i = std::numeric_limits<double>::denorm_min();
    assert(i > 0.0);
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,  i,0.0);
    const Coordinat3D p3(  i,0.0,0.0);
    const Plane p(p1,p2,p3);
    assert(!p.CanCalcX());
    assert(!p.CanCalcY());
    assert( p.CanCalcZ());

    assert( p.IsInPlane(Coordinat3D( 0.0, 0.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0, 1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,0.0)));

  }
  if (verbose) TRACE("CanCalcZ, Z = 0 plane, from biggest possible coordinat");
  {
    const double i = std::numeric_limits<double>::max();
    assert(i > 0.0);
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,  i,0.0);
    const Coordinat3D p3(  i,0.0,0.0);
    const Plane p(p1,p2,p3);
    assert(!p.CanCalcX());
    assert(!p.CanCalcY());
    assert( p.CanCalcZ());

    assert( p.IsInPlane(Coordinat3D( 0.0, 0.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0, 1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,0.0)));
  }
  if (verbose) TRACE("CanCalcZ, Z = 0 plane, zooming in");
  {
    for (const double i:series)
    {
      if (i == 0.0) continue;
      assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
      const Coordinat3D p1(0.0,0.0,0.0);
      const Coordinat3D p2(0.0,  i,0.0);
      const Coordinat3D p3(  i,0.0,0.0);
      const Plane p(p1,p2,p3);
      assert(!p.CanCalcX());
      assert(!p.CanCalcY());
      assert( p.CanCalcZ());

      for (const double j:series)
      {
        assert(p.IsInPlane(Coordinat3D(0.0,0.0,0.0)));
        assert(p.IsInPlane(Coordinat3D(  j,  j,0.0)));
        assert(p.IsInPlane(Coordinat3D(  j, -j,0.0)));
        assert(p.IsInPlane(Coordinat3D( -j,  j,0.0)));
        assert(p.IsInPlane(Coordinat3D( -j, -j,0.0)));
      }
    }
  }
  /*

    |    /#/##########
    |   B#/###########
    |  /#/############
    | /#/#############
    |/#/##############
    A-------C--------- Z = z
    |/
  --O----------------- Z = 0
   /|


  */
  if (verbose) TRACE("CanCalcZ, Z = z plane, varying height");
  {
    for (const double z:series)
    {
      const double i = 1.0;
      const Coordinat3D p1(0.0,0.0,z);
      const Coordinat3D p2(0.0,  i,z);
      const Coordinat3D p3(  i,0.0,z);
      const Plane p(p1,p2,p3);
      assert(!p.CanCalcX());
      assert(!p.CanCalcY());
      if (!p.CanCalcZ())
      {
        TRACE("ERROR");
        TRACE(z);
        TRACE(i);
      }

      assert( p.CanCalcZ());

      for (const double j:series)
      {
        if (show_warning && !p.IsInPlane(Coordinat3D(j,j,z)))
        {
          std::stringstream s;
          s << "Warning: coordinat " << Geometry().ToStr(Coordinat3D(j,j,z))
            << " is determined not to be in a Plane that was created from points "
            << Geometry().ToStr(p1) << ", "
            << Geometry().ToStr(p2) << " and "
            << Geometry().ToStr(p3) << ". Error: "
            << p.CalcError(Coordinat3D(j,j,z))
            << ", max error: "
            << p.CalcMaxError(Coordinat3D(j,j,z))
            << " ("
            << (p.CalcError(Coordinat3D(j,j,z)) / p.CalcMaxError(Coordinat3D(j,j,z)))
            << "x)"
          ;
          TRACE(s.str());
        }
        if (!p.IsInPlane(Coordinat3D(j,j,z)))
        {
          const auto max_error = p.CalcMaxError(Coordinat3D(j,j,z));
          const auto error = p.CalcError(Coordinat3D(j,j,z));
          if (error == 0.0) continue;
          TRACE(max_error);
          TRACE(error);
          if (abs(1.0 - (max_error / error)) < 0.01)
          {
            //Allow another percent of freedom
            continue;
          }

          TRACE("ERROR");
          TRACE(z);
          TRACE(i);
          TRACE(j);
          TRACE(p.CalcError(Coordinat3D(j,j,z)));
          TRACE(p.CalcMaxError(Coordinat3D(j,j,z)));
          TRACE(p);

          TRACE("AGAIN");
          TRACE(z);
          TRACE(i);
          TRACE(j);
          TRACE(p.CalcError(Coordinat3D(j,j,z)));
          TRACE(p.CalcMaxError(Coordinat3D(j,j,z)));
          TRACE(p);

        }
        assert(p.IsInPlane(Coordinat3D(  j,  j,z)));
        assert(p.IsInPlane(Coordinat3D(  j, -j,z)));
        assert(p.IsInPlane(Coordinat3D( -j,  j,z)));
        assert(p.IsInPlane(Coordinat3D( -j, -j,z)));
      }
    }
  }
  if (verbose) TRACE("CanCalcZ, Z = z plane, zooming in");
  {
    const double min = 10e+8;
    const double max = 10e-8;
    //for (const double z:series)
    for (double z=min; z<max; z*=10.0)
    {
      //for (const double i:series)
      for (double i=min; i<max; i*=10.0)
      {
        if (i == 0.0) continue;
        assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
        const Coordinat3D p1(0.0,0.0,z);
        const Coordinat3D p2(0.0,  i,z);
        const Coordinat3D p3(  i,0.0,z);
        const Plane p(p1,p2,p3);
        assert(!p.CanCalcX());
        assert(!p.CanCalcY());
        if (!p.CanCalcZ())
        {
          TRACE("ERROR");
          TRACE(z);
          TRACE(i);
        }

        assert( p.CanCalcZ());

        for (const double j:series)
        {
          if (!p.IsInPlane(Coordinat3D(j,j,z)))
          {
            std::stringstream s;
            s << "Warning: coordinat " << Geometry().ToStr(Coordinat3D(j,j,z))
              << " is determined not to be in a Plane that was created from points "
              << Geometry().ToStr(p1) << ", "
              << Geometry().ToStr(p2) << " and "
              << Geometry().ToStr(p3) << ". Error: "
              << p.CalcError(Coordinat3D(j,j,z))
              << ", max error: "
              << p.CalcMaxError(Coordinat3D(j,j,z))
              << " ("
              << (p.CalcError(Coordinat3D(j,j,z)) / p.CalcMaxError(Coordinat3D(j,j,z)))
              << "x)"
            ;
            TRACE(s.str());
            //continue;
          }
          if (!p.IsInPlane(Coordinat3D(j,j,z)))
          {
            TRACE("ERROR");
            TRACE(z);
            TRACE(i);
            TRACE(j);
            TRACE(p.CalcError(Coordinat3D(j,j,z)));
            TRACE(p.CalcMaxError(Coordinat3D(j,j,z)));
            TRACE(p);

            TRACE("AGAIN");
            TRACE(z);
            TRACE(i);
            TRACE(j);
            TRACE(p.CalcError(Coordinat3D(j,j,z)));
            TRACE(p.CalcMaxError(Coordinat3D(j,j,z)));
            TRACE(p);

          }
          assert(p.IsInPlane(Coordinat3D(  j,  j,z)));
          assert(p.IsInPlane(Coordinat3D(  j, -j,z)));
          assert(p.IsInPlane(Coordinat3D( -j,  j,z)));
          assert(p.IsInPlane(Coordinat3D( -j, -j,z)));
        }
      }
    }
  }
  //IsInPlane for X=0 plane
  /*

     |####/
     B###C
     |##/
     |#/
     |/
  ---A------
    /|
   /#|
  /##|

  */

  if (verbose) TRACE("IsInPlane, X = 0 plane, from 1.0 coordinats");
  {
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,0.0,1.0);
    const Coordinat3D p3(0.0,1.0,0.0);
    const Plane p(p1,p2,p3);
    assert( p.CanCalcX());
    assert(!p.CanCalcY());
    assert(!p.CanCalcZ());

    for (double i = std::numeric_limits<double>::denorm_min();
      i < 1.0e-8; //std::numeric_limits<double>::max();
      i *= 10.0)
    {
      assert(p.IsInPlane(Coordinat3D(0.0, i, i)));
      assert(p.IsInPlane(Coordinat3D(0.0, i,-i)));
      assert(p.IsInPlane(Coordinat3D(0.0,-i, i)));
      assert(p.IsInPlane(Coordinat3D(0.0,-i,-i)));
    }
  }
  if (verbose) TRACE("IsInPlane, X = 0 plane, from smallest possible coordinats");
  {
    const double i = std::numeric_limits<double>::denorm_min();
    assert(i > 0.0);
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,0.0,i);
    const Coordinat3D p3(0.0,i,0.0);
    const Plane p(p1,p2,p3);
    assert( p.CanCalcX());
    assert(!p.CanCalcY());
    assert(!p.CanCalcZ());
  }
  if (verbose) TRACE("IsInPlane, X = 0 plane, zooming in, #223");
  {
    for (const double i:series)
    {
      if (i == 0.0) continue;
      assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
      const Coordinat3D p1(0.0,0.0,0.0);
      const Coordinat3D p2(0.0,0.0,  i);
      const Coordinat3D p3(0.0,  i,0.0);
      const Plane p(p1,p2,p3);

      for (const double j:series)
      {
        if (!p.IsInPlane(Coordinat3D(0.0, j, j)))
        {
          TRACE(i);
          TRACE(j);
          TRACE(p.CalcError(Coordinat3D(0.0, j, j)));
          TRACE(p.CalcMaxError(Coordinat3D(0.0, j, j)));
          if (p.CanCalcX()) { TRACE(Container().ToStr(p.GetCoefficientsX())); }
          if (p.CanCalcY()) { TRACE(Container().ToStr(p.GetCoefficientsY())); }
          if (p.CanCalcZ()) { TRACE(Container().ToStr(p.GetCoefficientsZ())); }
        }
        assert(p.IsInPlane(Coordinat3D(0.0, j, j)));
        assert(p.IsInPlane(Coordinat3D(0.0, j,-j)));
        assert(p.IsInPlane(Coordinat3D(0.0,-j, j)));
        assert(p.IsInPlane(Coordinat3D(0.0,-j,-j)));
      }
    }
  }



  //IsInPlane for Y=0 plane
  /*

     |####/
     B###/#
     |##/##
     |#/###
     |/####
  ---A---C-
    /|
   / |
  /  |

  */

  if (verbose) TRACE("IsInPlane, Y = 0 plane, from 1.0 coordinats");
  {
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,0.0,1.0);
    const Coordinat3D p3(1.0,0.0,0.0);
    const Plane p(p1,p2,p3);
    assert(!p.CanCalcX());
    assert( p.CanCalcY());
    assert(!p.CanCalcZ());

    for (double i = std::numeric_limits<double>::denorm_min();
      i < 1.0e8 /*std::numeric_limits<double>::max()*/;
      i *= 10.0
    )
    {
      assert(p.IsInPlane(Coordinat3D( i,0.0, i)));
      assert(p.IsInPlane(Coordinat3D( i,0.0,-i)));
      assert(p.IsInPlane(Coordinat3D(-i,0.0, i)));
      assert(p.IsInPlane(Coordinat3D(-i,0.0,-i)));
    }
  }
  if (verbose) TRACE("IsInPlane, Y = 0 plane, from smallest possible coordinats");
  {
    const double i = std::numeric_limits<double>::denorm_min();
    assert(i > 0.0);
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,0.0,  i);
    const Coordinat3D p3(  i,0.0,0.0);
    const Plane p(p1,p2,p3);
    assert(!p.CanCalcX());
    assert( p.CanCalcY());
    assert(!p.CanCalcZ());
  }
  if (verbose) TRACE("IsInPlane, Y = 0 plane, zooming in, #223");
  {
    for (const double i:series)
    {
      if (i == 0.0) continue;
      assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
      const Coordinat3D p1(0.0,0.0,0.0);
      const Coordinat3D p2(0.0,0.0,  i);
      const Coordinat3D p3(  i,0.0,0.0);
      const Plane p(p1,p2,p3);

      for (const double j:series)
      {
        if (!p.IsInPlane(Coordinat3D(j,0.0, j)))
        {
          TRACE(i);
          TRACE(j);
          TRACE(p.CalcError(Coordinat3D(j,0.0, j)));
          TRACE(p.CalcMaxError(Coordinat3D(j,0.0,j)));
          if (p.CanCalcX()) { TRACE(Container().ToStr(p.GetCoefficientsX())); }
          if (p.CanCalcY()) { TRACE(Container().ToStr(p.GetCoefficientsY())); }
          if (p.CanCalcZ()) { TRACE(Container().ToStr(p.GetCoefficientsZ())); }
        }
        assert(p.IsInPlane(Coordinat3D( j,0.0, j)));
        assert(p.IsInPlane(Coordinat3D( j,0.0,-j)));
        assert(p.IsInPlane(Coordinat3D(-j,0.0, j)));
        assert(p.IsInPlane(Coordinat3D(-j,0.0,-j)));
      }
    }
  }
  //Create plane with different slopes
  /*
  |          /
  |         /
  |        /         __--D
  |       /      __--   /|
  |      /   __--      / |
  |     /__--         /  |
  |    B             /   +
  |   /         __--C   /
  |  /      __--    |  /
  | /   __--        | /
  |/__--            |/
  A-----------------+---------------

  A = p1 = Origin, fixed
  B = p2 = At y axis, fixed
  C = p3 = Above X axis, dependent on slope
  D = p4 = To be calculated
  */
  if (verbose) TRACE("IsInPlane, Slope, Slope in Z direction");
  {
    for (double slope = 1.0; slope > 1.0e-8; slope /= 10.0)
    {
      //if (verbose) { TRACE(slope); }
      const double slope_less = slope * 0.999999;
      const double slope_more = slope * 1.000001;

      const Coordinat3D p1(0.0,0.0,0.0);
      const Coordinat3D p2(0.0,1.0,0.0);
      const Coordinat3D p3(1.0,0.0,slope);
      const Plane p(p1,p2,p3);
      assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,slope)));

      if (slope_less < slope) //Not always true, when slope is very small
      {
        if (p.IsInPlane(Coordinat3D(1.0, 1.0,slope_less)))
        {
          TRACE("ERROR");
          TRACE(p.CalcError(Coordinat3D(1.0, 1.0,slope_less)));
          TRACE(p.CalcMaxError(Coordinat3D(1.0, 1.0,slope_less)));
        }
        assert(!p.IsInPlane(Coordinat3D( 1.0, 1.0,slope_less)));
      }
      if (slope_more > slope) //Not always true, when slope is very big
      {
        if (p.IsInPlane(Coordinat3D( 1.0, 1.0,slope_more)))
        {
          TRACE("ERROR");
          TRACE(p.CalcError(Coordinat3D( 1.0, 1.0,slope_more)));
        }
        assert(!p.IsInPlane(Coordinat3D( 1.0, 1.0,slope_more)));
      }
    }
  }
  //Create plane with different slopes
  /*
  |          /
  |         /
  |        /__-D
  |     __--   |
  | __-- /     |
  B-    /      |
  |    /       |
  |   /        |
  |  /      __-C
  | /   __--  /
  |/__--     /
  A---------+-----------------------

  A = p1 = Origin, fixed
  B = p2 = At y axis, fixed
  C = p3 = Above X axis, dependent on slope
  D = p4 = To be calculated
  */
  if (verbose) TRACE("IsInPlane, Slope, vertical plane");
  {
    typedef std::pair<double,double> Pair;
    for (const Pair co:
      {
        std::make_pair(    1.0,    1.0),
        std::make_pair(   10.0,   10.0),
        std::make_pair(  100.0,  100.0),
        std::make_pair( 1000.0, 1000.0),
        std::make_pair(10000.0,10000.0),
        std::make_pair(1.0,1.000),
        std::make_pair(1.0,0.100),
        std::make_pair(1.0,0.010),
        std::make_pair(1.0,0.001),
      }
    )
    {
      const Coordinat3D p1(0.0,0.0,0.0);
      const Coordinat3D p2(0.0,0.0,1.0);
      const Coordinat3D p3(co.first,co.second,0.0);
      const Plane p(p1,p2,p3);
      assert( p.IsInPlane(Coordinat3D(co.first,co.second,1.0)));
    }
  }
  if (verbose) TRACE("CalcProjection, from a crash in the program");
  {
    const Coordinat3D p1( 1.0,-0.0,0.0);
    const Coordinat3D p2(-1.0, 0.0,0.0);
    const Coordinat3D p3( 1.0,-0.0,1.0);
    const Coordinat3D p4(-1.0, 0.0,1.0);
    const Plane p(p1,p2,p3);

    assert(p.CanCalcY());
    assert(!p.CanCalcX());
    assert(!p.CanCalcZ());
    assert(!p.CalcProjection( { p1,p2,p3,p4 } ).empty());
  }
  if (verbose) TRACE("IsInPlane, from #218");
  {
    const double x1 = -0.5500000000000004884981308350688777863979339599609375;
    const double y1 = 2.000000000000000444089209850062616169452667236328125;
    const double z1 = 0; //left out the '.0' intentionally
    const double x2 = -3.78623595505618038004058689693920314311981201171875;
    const double y2 = 2; //left out the '.0' intentionally
    const double z2 = 0; //left out the '.0' intentionally
    const double x3 = -0.5500000000000004884981308350688777863979339599609375;
    const double y3 = 2.000000000000000444089209850062616169452667236328125;
    const double z3 = 10; //left out the '.0' intentionally
    const double x4 = -3.78623595505618038004058689693920314311981201171875;
    const double y4 = 2; //left out the '.0' intentionally
    const double z4 = 10; //left out the '.0' intentionally

    const Coordinat3D p1(x1,y1,z1);
    const Coordinat3D p2(x2,y2,z2);
    const Coordinat3D p3(x3,y3,z3);
    const Coordinat3D p4(x4,y4,z4);
    const Plane p(p1,p2,p3);
    try
    {
      const bool must_be_true = p.IsInPlane(p4);
      assert(must_be_true);
    }
    catch (std::exception&)
    {
      TRACE("ERROR");
      std::stringstream s;
      s
        << std::setprecision(99)
        << p << '\n'
        << Container().ToStr(p.GetCoefficientsX()) << '\n'
        << Container().ToStr(p.GetCoefficientsY()) << '\n'
        << Container().ToStr(p.GetCoefficientsZ()) << '\n'
      ;
      TRACE(s.str());
      TRACE("BREAK");
    }
    assert(p.IsInPlane(p4));
  }

  if (verbose) TRACE("IsInPlane, crashes with Plane v1.5");
  {
    const double x1 = 0.0004035051226622692510832834944523028752882964909076690673828125;
    const double y1 = 0.00023296416881187433805568132161312178141088224947452545166015625;
    const double z1 = 0; //left out the '.0' intentionally

    const double x2 = 0.000403505141811931846741734464245610070065595209598541259765625;
    const double y2 = 0.00023296414405748076185791173298156309101614169776439666748046875;
    const double z2 = 0; //left out the '.0' intentionally

    const double x3 = 0.0004035051226622692510832834944523028752882964909076690673828125;
    const double y3 = 0.00023296416881187433805568132161312178141088224947452545166015625;
    const double z3 = 0.00025000000000000000520417042793042128323577344417572021484375;

    const double x4 = 0.000403505141811931846741734464245610070065595209598541259765625;
    const double y4 = 0.00023296414405748076185791173298156309101614169776439666748046875;
    const double z4 = 0.00025000000000000000520417042793042128323577344417572021484375;

    {
      const apfloat x1_as_apfloat(x1);
      const double  x1_as_double(Geometry().ToDoubleSafe(x1_as_apfloat));
      assert(x1 == x1_as_double);
    }

    const Coordinat3D p1(x1,y1,z1);
    const Coordinat3D p2(x2,y2,z2);
    const Coordinat3D p3(x3,y3,z3);
    const Coordinat3D p4(x4,y4,z4);
    const Plane p(p1,p2,p3);

    if (!p.IsInPlane(p4))
    {
      TRACE("ERROR");
      TRACE(p.CalcError(p4));
      TRACE(p.CalcMaxError(p4));

      TRACE("BREAK");
    }

    assert(p.IsInPlane(p4));
  }
  #ifdef NOT_TODAY_20140714
  if (verbose) TRACE("IsInPlane, crashes with Plane v1.6");
  {
    // TRACE '"ERROR"' line 392 in file '..\..\Classes\CppTriangleMesh\trianglemeshcellscreator.cpp': 'ERROR'
    // TRACE 's.str()' line 399 in file '..\..\Classes\CppTriangleMesh\trianglemeshcellscreator.cpp': '4
    // (-5,5,0) (index: 659) (-5,-0.9999999999999997779553950749686919152736663818359375,0) (index: 666) (-5,5,1) (index: 700) (-5,-0.9999999999999997779553950749686919152736663818359375,1) (index: 707) '
    // TRACE '"BREAK"' line 400 in file '..\..\Classes\CppTriangleMesh\trianglemeshcellscreator.cpp': 'BREAK'
    // Assertion failed!
    const double x1 = -5; //left out the '.0' intentionally
    const double y1 = 5; //left out the '.0' intentionally
    const double z1 = 0; //left out the '.0' intentionally

    const double x2 = -5; //left out the '.0' intentionally
    const double y2 = -0.9999999999999997779553950749686919152736663818359375;
    const double z2 = 0; //left out the '.0' intentionally

    const double x3 = -5; //left out the '.0' intentionally
    const double y3 = 5; //left out the '.0' intentionally
    const double z3 = 1; //left out the '.0' intentionally

    const double x4 = -5; //left out the '.0' intentionally
    const double y4 = -0.9999999999999997779553950749686919152736663818359375;
    const double z4 = 1; //left out the '.0' intentionally


    const Coordinat3D p1(x1,y1,z1);
    const Coordinat3D p2(x2,y2,z2);
    const Coordinat3D p3(x3,y3,z3);
    const Coordinat3D p4(x4,y4,z4);
    const Plane p(p1,p2,p3);
    try
    {
      assert(p.IsInPlane(p4));
    }
    catch (std::exception&)
    {
      TRACE("ERROR");
      std::stringstream s;
      s
        << std::setprecision(99)
        << p << '\n'
        << Container().ToStr(p.GetCoefficientsX()) << '\n'
        << Container().ToStr(p.GetCoefficientsY()) << '\n'
        << Container().ToStr(p.GetCoefficientsZ()) << '\n'
      ;
      TRACE(s.str());
      TRACE("BREAK");
    }
    assert(p.IsInPlane(p4));
  }
  #endif
}
#endif

 

 

 

 

 

./CppPlane/planex.h

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.htm
//---------------------------------------------------------------------------
#ifndef RIBI_PLANEX_H
#define RIBI_PLANEX_H

#include <vector>
#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 <boost/geometry/geometries/point_xy.hpp>
#include <boost/make_shared.hpp>
//#include "planez.h"
#include "apfloat.h"
#pragma GCC diagnostic pop

namespace ribi {

struct PlaneZ;

///A 3D plane that can have its X expressed as a function of Y and Z.
///Can be constructed from its equation and at least three 3D points
///A plane stores its coefficients in the following form:
/// A.z + B.x + C.y = D
///Converting this to x being a function of y and z:
/// -B.x =  A  .z + C  .y - D
///    x = -A/B.z - C/B.y + D/B
///where A,B,C and D can be obtained by GetCoefficients
///Easier is to express X as
///  x = Ay + Bz + C
///where A,B,C can be obtained by GetFunctionA/B/C
struct PlaneX
{
  typedef apfloat Double;
  typedef boost::geometry::model::d2::point_xy<apfloat> Coordinat2D;
  typedef boost::geometry::model::point<apfloat,3,boost::geometry::cs::cartesian> Coordinat3D;
  typedef std::vector<Double> Doubles;
  typedef std::vector<Coordinat2D> Coordinats2D;
  typedef std::vector<Coordinat3D> Coordinats3D;

  ///Create plane X = 0.0
  PlaneX() noexcept;

  ///Construct from three points
  explicit PlaneX(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  );

  Double CalcError(const Coordinat3D& coordinat) const noexcept;

  Double CalcMaxError(const Coordinat3D& coordinat) const noexcept;

  ///Get the 2D projection of these 3D points,
  /*

  A: (0,0,1)                  A: (0,0)
  B: (1,0,0)                  B: (SQRT(2),0)
  C: (1,1,0)                  C: (SQRT(2),SQRT(2))

  |    /
  |   /                       |
  A-----C                     |  C
  |\/  /      -> becomes ->   | /|
  |/\ /                       |/ |
  +--B---                     A--B-----

  */
  Coordinats2D CalcProjection(const Coordinats3D& points) const;

  ///Throws when cannot calculate X, which is when the plane is horizontal
  Double CalcX(const Double& y, const Double& z) const;

  Doubles GetCoefficients() const noexcept;

  ///x = Ay + Bz + C
  ///Will throw if A cannot be calculated
  Double GetFunctionA() const;

  ///x = Ay + Bz + C
  ///Will throw if B cannot be calculated
  Double GetFunctionB() const;

  ///x = Ay + Bz + C
  ///Will throw if C cannot be calculated
  Double GetFunctionC() const;

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

  ///Checks if the coordinat is in the plane
  bool IsInPlane(const Coordinat3D& coordinat) const noexcept;

  private:
  ~PlaneX() noexcept;

  ///A PlaneX is actually a PlaneZ used with its coordinats rotated from (X,Y,Z) to (Z,Y,Y)
  const std::unique_ptr<PlaneZ> m_plane_z;

  ///Calculates m_min_error per GetFunctionC()
  static Double CalcMinErrorPerC() noexcept;

  ///Will throw if plane cannot be created
  static std::unique_ptr<PlaneZ> Create(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  );

  static Doubles Rotate(const Doubles& coefficients) noexcept;

  ///Rotates the X,Y and Z value of a Coordinat
  static Coordinat3D Rotate(const Coordinat3D& point) noexcept;

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

  ///Convert the PlaneX to a x(y,z), e.g 'x=(2*y) + (3*z) + 5' (spaces exactly as shown)
  std::string ToFunction() const;

  friend void boost::checked_delete<>(      PlaneX*);
  friend void boost::checked_delete<>(const PlaneX*);
  friend struct std::default_delete<      PlaneX>;
  friend struct std::default_delete<const PlaneX>;
  friend class boost::detail::sp_ms_deleter<      PlaneX>;
  friend class boost::detail::sp_ms_deleter<const PlaneX>;
  friend std::ostream& operator<<(std::ostream& os,const PlaneX& planex);
};

std::ostream& operator<<(std::ostream& os,const PlaneX& planex);

} //~namespace ribi

#endif // RIBI_PLANEX_H

 

 

 

 

 

./CppPlane/planex.cpp

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.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 "planex.h"

#include <cassert>

#include "container.h"
#include "geometry.h"
#include "planez.h"
#include "trace.h"
#pragma GCC diagnostic pop

///Create plane X = 0.0
ribi::PlaneX::PlaneX() noexcept
  : PlaneX(
    Coordinat3D(0.0,0.0,0.0),
    Coordinat3D(0.0,1.0,0.0),
    Coordinat3D(0.0,0.0,1.0)
  )
{
  #ifndef NDEBUG
  Test();
  #endif
}


ribi::PlaneX::PlaneX(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) : m_plane_z{Create(p1,p2,p3)}
{
  #ifndef NDEBUG
  Test();
  #endif
}

ribi::PlaneX::~PlaneX() noexcept
{
  //Nothing to do
}

apfloat ribi::PlaneX::CalcError(const Coordinat3D& coordinat) const noexcept
{
  const apfloat x = boost::geometry::get<0>(coordinat);
  const apfloat y = boost::geometry::get<1>(coordinat);
  const apfloat z = boost::geometry::get<2>(coordinat);
  const auto expected = x;
  const apfloat calculated = CalcX(y,z);
  const apfloat error = abs(calculated - expected);
  return error;
}

ribi::PlaneX::Double ribi::PlaneX::CalcMinErrorPerC() noexcept
{
  //min_error_per_c will be about 0.000000001
  //stub_value increases this jut a little, by a 0.000001%
  const double stub_value = 0.000000001 * 1.00000001;
  #define USE_STUB
  #ifdef USE_STUB
  return stub_value;
  #else //USE_STUB
  //PlaneX calculates its own tolerance for errors, by measuring it
  static Double min_error_per_c = 0.0;
  if (min_error_per_c > 0.0) return min_error_per_c;

  //const double low = std::numeric_limits<double>::denorm_min();
  //const double high = std::numeric_limits<double>::max();
  const double low  = 1.0e-16;
  const double high = 1.0e+16;
  const double min_x = low;
  const double max_x = high;
  const double min_y = low;
  const double max_y = high;
  const double min_z = low;
  const double max_z = high;
  const apfloat zero(0.0);

  for (double z = min_z; z < max_z; z*=10.0)
  {
    for (double y = min_y; y < max_y; y*=10.0)
    {
      for (double x = min_x; x < max_x; x*=10.0)
      //const double x = y;
      {
        const Coordinat3D p1(  x,0.0,0.0);
        const Coordinat3D p2(  x,  y,0.0);
        const Coordinat3D p3(  x,0.0,z);
        const PlaneX p(p1,p2,p3);
        for (const auto& p4: { p1, p2, p3 } )
        {
          const auto error = p.CalcError(p4);
          const auto error_per_c = error / p.GetFunctionC();
          assert(error_per_c >= zero);
          //TRACE(apfloat(min_error_per_c) / p.GetFunctionC());
          if (error_per_c > min_error_per_c)
          {
            min_error_per_c = error_per_c;
            //TRACE(min_error_per_c);
            //TRACE(x);
            //TRACE(y);
            //TRACE(z);
            //TRACE(p.GetFunctionC());
            //TRACE(apfloat(min_error_per_c) / p.GetFunctionC());
            //std::stringstream s;
            //s << Geometry().ToStr(p4) << " " << min_error;
            //TRACE(s.str());
          }
        }
      }
    }
    //TRACE(min_error_per_c);
  }
  //TRACE(min_error_per_c); //0.000000001e0
  assert(min_error_per_c > zero);
  assert(min_error_per_c < stub_value);
  assert(min_error_per_c > 0.99 * stub_value);
  return min_error_per_c;
  #endif //USE_STUB
}

apfloat ribi::PlaneX::CalcMaxError(const Coordinat3D& /*coordinat*/) const noexcept
{
  assert(CalcMinErrorPerC() > apfloat(0.0));
  const auto max_error = abs(CalcMinErrorPerC() * GetFunctionC());
  assert(max_error >= apfloat(0.0));
  return max_error;
  /*
  //const apfloat x = boost::geometry::get<0>(coordinat);
  const apfloat y = boost::geometry::get<1>(coordinat);
  const apfloat z = boost::geometry::get<2>(coordinat);
  const auto coefficients = GetCoefficients();
  const auto a = coefficients[0];
  const auto b = coefficients[1];
  const auto c = coefficients[2];
  const apfloat e = boost::numeric::bounds<double>::smallest();
  assert(e > 0.0);
  //    x = -A/B.z - C/B.y + D/B
  //If B is zero, the slope in X and Y cannot be calculated
  if (b.sign())
  {
    const auto rc_y = -c / b;
    const auto rc_z = -a / b;
    const auto max_error_x = 0.0;
    const auto max_error_z = abs(e * rc_z * z) + 0.0;
    const auto max_error_y = abs(e * rc_y * y) + 0.0;
    const auto max_error = max_error_x + max_error_y + max_error_z;
    assert(max_error >= 0.0);
    return max_error;
  }
  assert(e > 0.0);
  return e;
  */
}

ribi::PlaneX::Coordinats2D ribi::PlaneX::CalcProjection(
  const Coordinats3D& points
) const
{
  assert(m_plane_z);
  auto v(points);
  for(auto& i: v) { i = Rotate(i); }
  try
  {
    return m_plane_z->CalcProjection(v);
  }
  catch (std::logic_error&)
  {
    throw std::logic_error("PlaneX::CalcProjection: cannot calculate projection");
  }
}

ribi::PlaneX::Double ribi::PlaneX::CalcX(const Double& y, const Double& z) const
{
  assert(m_plane_z);
  try
  {
    return m_plane_z->CalcZ(y,z);
  }
  catch (std::logic_error&)
  {
    throw std::logic_error("ribi::PlaneX::CalcX: cannot calculate X of a horizontal plane");
  }
}

std::unique_ptr<ribi::PlaneZ> ribi::PlaneX::Create(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
)
{
  std::unique_ptr<PlaneZ> p(
    new PlaneZ(
      Rotate(p1),
      Rotate(p2),
      Rotate(p3)
    )
  );
  assert(p);
  return p;
}

std::vector<apfloat> ribi::PlaneX::GetCoefficients() const noexcept
{
  const auto v(m_plane_z->GetCoefficients());
  assert(v.size() == 4);
  return { v[2],v[0],v[1],v[3] };
}

apfloat ribi::PlaneX::GetFunctionA() const
{
  assert(m_plane_z);
  return m_plane_z->GetFunctionA();
}

apfloat ribi::PlaneX::GetFunctionB() const
{
  assert(m_plane_z);
  return m_plane_z->GetFunctionB();
}

apfloat ribi::PlaneX::GetFunctionC() const
{
  assert(m_plane_z);
  return m_plane_z->GetFunctionC();
}

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

std::vector<std::string> ribi::PlaneX::GetVersionHistory() const noexcept
{
  return {
    "2014-03-10: version 1.0: initial version, split off from PlaneZ",
    "2014-03-13: version 1.1: bug fixed",
    "2014-04-01: version 1.2: use of std::unique_ptr",
    "2014-06-13: version 1.3: shortened time to compile, allow obtaining the constants in function 'x = Ay + Bz + C'",
    "2014-07-03: version 1.4: use of apfloat",
    "2014-07-09: version 1.5: use double in interface only",
    "2014-07-10: version 1.6: use of apfloat only"
  };
}

bool ribi::PlaneX::IsInPlane(const Coordinat3D& coordinat) const noexcept
{
  try
  {
    const apfloat error = CalcError(coordinat);
    const apfloat max_error = CalcMaxError(coordinat);
    return error <= max_error;
  }
  catch (std::exception& e)
  {
    TRACE("ERROR");
    TRACE(e.what());
    assert(!"Should not get here");
    throw;
  }
}

std::vector<apfloat> ribi::PlaneX::Rotate(const Doubles& coefficients) noexcept
{
  assert(coefficients.size() == 4);
  return
  {
    coefficients[2],
    coefficients[0],
    coefficients[1],
    coefficients[3]
  };
}

ribi::PlaneX::Coordinat3D ribi::PlaneX::Rotate(
  const Coordinat3D& point
) noexcept
{
  return Coordinat3D(
    boost::geometry::get<1>(point),
    boost::geometry::get<2>(point),
    boost::geometry::get<0>(point)
  );
}



std::string ribi::PlaneX::ToFunction() const
{
  std::stringstream s;
  s << (*this);
  return s.str();
}

std::ostream& ribi::operator<<(std::ostream& os,const PlaneX& planex)
{

  assert(planex.m_plane_z);
  try
  {
    os
      << "x=("
      << planex.m_plane_z->GetFunctionA() << "*y) + ("
      << planex.m_plane_z->GetFunctionB() << "*z) + "
      << planex.m_plane_z->GetFunctionC()
    ;
  }
  catch (std::logic_error&)
  {
    throw std::logic_error("ribi::PlaneX::operator<<: cannot display function of a horizontal plane");
  }
  return os;
}

 

 

 

 

 

./CppPlane/planex_test.cpp

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.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 "planex.h"

#include <cassert>

#include "container.h"
#include "geometry.h"
#include "planez.h"
#include "testtimer.h"
#include "trace.h"
#pragma GCC diagnostic pop

#ifndef NDEBUG
void ribi::PlaneX::Test() noexcept
{
  {
    static bool is_tested { false };
    if (is_tested) return;
    is_tested = true;
  }
  const TestTimer test_timer(__func__,__FILE__,1.0);
  using boost::geometry::get;

  const bool verbose{false};
  if (verbose) TRACE("Default construction");
  {
    const PlaneX p;
    assert(!p.ToFunction().empty());
    assert(!p.GetCoefficients().empty());
  }
  /*
  {


    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 };
    PlaneX p(
      Point(p1_x,p1_y,p1_z),
      Point(p2_x,p2_y,p2_z),
      Point(p3_x,p3_y,p3_z)
    );
    const auto t(p.GetCoefficients());
    assert(t.size() == 4);
    const double a { t[0] };
    const double b { t[1] };
    const double c { t[2] };
    const double d { t[3] };
    const double a_expected {  30.0 };
    const double b_expected { -48.0 };
    const double c_expected {  17.0 };
    const double d_expected { -15.0 };
    assert(abs(a - a_expected) < 0.001);
    assert(abs(b - b_expected) < 0.001);
    assert(abs(c - c_expected) < 0.001);
    assert(abs(d - d_expected) < 0.001);
    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'
      ;
    }
    assert(abs(d - d_p1_expected) < 0.001);
    assert(abs(d - d_p2_expected) < 0.001);
    assert(abs(d - d_p3_expected) < 0.001);
  }
  */
  //CalcPlaneX
  /*
  {
    //CalcPlaneX 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


    const Point3D p1(1.0,1.0,10.0);
    const Point3D p2(1.0,2.0,13.0);
    const Point3D p3(2.0,1.0,12.0);
    PlaneX p(p1,p2,p3);
    const auto t(p.GetCoefficients());
    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(abs(a - a_expected) < 0.001);
    assert(abs(b - b_expected) < 0.001);
    assert(abs(c - c_expected) < 0.001);
    assert(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(abs(d - d_p1_expected) < 0.001);
    assert(abs(d - d_p2_expected) < 0.001);
    assert(abs(d - d_p3_expected) < 0.001);
    TRACE(p.ToFunction());
  }
  */
  if (verbose) TRACE("CalcX, diagonal plane");
  {


    const Coordinat3D p1(1.0,2.0,3.0);
    const Coordinat3D p2(2.0,5.0,8.0);
    const Coordinat3D p3(3.0,7.0,11.0);
    PlaneX p(p1,p2,p3);
    assert( abs(p.CalcX(2.0, 3.0)- 1.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcX(5.0, 8.0)- 2.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcX(7.0,11.0)-3.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
  if (verbose) TRACE("CalcX, vertical plane X = 2.0");
  /*

    |####/
    |###/
    |##/
    |#/
    |/
---+--------
   /|
  / |
/  |

  */
  {


    const Coordinat3D p1(2.0, 3.0, 5.0);
    const Coordinat3D p2(2.0, 7.0,11.0);
    const Coordinat3D p3(2.0,13.0,17.0);
    PlaneX p(p1,p2,p3);
    assert( abs(p.CalcX(1.0,2.0)-2.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcX(3.0,5.0)-2.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcX(7.0,9.0)-2.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
  if (verbose) TRACE("ToFunction, 3 points and 4 points");
  {
    std::function<double(double,double)> f {
      [](const double y, const double z)
      {
        return (2.0 * y) + (3.0 * z) + 5.0;
      }
    };


    const double y1 { 2.0 };
    const double z1 { 3.0 };
    const double y2 { 5.0 };
    const double z2 { 7.0 };
    const double y3 { 11.0 };
    const double z3 { 13.0 };
    const double y4 { 17.0 };
    const double z4 { 29.0 };
    const Coordinat3D p1(f(y1,z1),y1,z1);
    const Coordinat3D p2(f(y2,z2),y2,z2);
    const Coordinat3D p3(f(y3,z3),y3,z3);
    const PlaneX a(p1,p2,p3);
    //assert(a.ToFunction() == "x=(2*y) + (3*z) + 5");
    assert(!a.ToFunction().empty());
    const Coordinat3D p4(f(y4,z4),y4,z4);
    assert(a.ToFunction() == PlaneX(p1,p2,p4).ToFunction());
    assert(a.ToFunction() == PlaneX(p1,p3,p4).ToFunction());
    assert(a.ToFunction() == PlaneX(p1,p4,p3).ToFunction());
    assert(a.ToFunction() == PlaneX(p2,p1,p3).ToFunction());
    assert(a.ToFunction() == PlaneX(p2,p1,p4).ToFunction());
    assert(a.ToFunction() == PlaneX(p2,p3,p1).ToFunction());
    assert(a.ToFunction() == PlaneX(p2,p3,p4).ToFunction());
    assert(a.ToFunction() == PlaneX(p2,p4,p1).ToFunction());
    assert(a.ToFunction() == PlaneX(p2,p4,p3).ToFunction());
    assert(a.ToFunction() == PlaneX(p3,p1,p2).ToFunction());
    assert(a.ToFunction() == PlaneX(p3,p1,p4).ToFunction());
    assert(a.ToFunction() == PlaneX(p3,p2,p1).ToFunction());
    assert(a.ToFunction() == PlaneX(p3,p2,p4).ToFunction());
    assert(a.ToFunction() == PlaneX(p3,p4,p1).ToFunction());
    assert(a.ToFunction() == PlaneX(p3,p4,p2).ToFunction());
    assert(a.ToFunction() == PlaneX(p4,p1,p2).ToFunction());
    assert(a.ToFunction() == PlaneX(p4,p1,p3).ToFunction());
    assert(a.ToFunction() == PlaneX(p4,p2,p1).ToFunction());
    assert(a.ToFunction() == PlaneX(p4,p2,p3).ToFunction());
    assert(a.ToFunction() == PlaneX(p4,p3,p1).ToFunction());
    assert(a.ToFunction() == PlaneX(p4,p3,p2).ToFunction());
  }

  if (verbose) TRACE("IsInPlane, X = 1, zooming to smallest three points to determine a plane, point above origin");
  {
    for (double i = 1.0; i > 1.0e-8; i/=10.0)
    {
      const Coordinat3D p1(1.0,0.0,0.0);
      const Coordinat3D p2(1.0,  i,0.0);
      const Coordinat3D p3(1.0,0.0,  i);
      const PlaneX p(p1,p2,p3);
      if (verbose)
      {
        TRACE("----------------------------");
        TRACE(i);
        TRACE(p.CalcMaxError(p1));
        TRACE(p.CalcError(p1));
        TRACE(p.CalcMaxError(p2));
        TRACE(p.CalcError(p2));
        TRACE(p.CalcMaxError(p3));
        TRACE(p.CalcError(p3));
        TRACE(p.GetFunctionA());
        TRACE(p.GetFunctionB());
        TRACE(p.GetFunctionC());
        TRACE(p.GetCoefficients()[0]);
        TRACE(p.GetCoefficients()[1]);
        TRACE(p.GetCoefficients()[2]);
        TRACE(p.GetCoefficients()[3]);
        TRACE(std::numeric_limits<double>::epsilon());
        TRACE(std::sqrt(std::numeric_limits<double>::epsilon()));
        TRACE(std::numeric_limits<double>::denorm_min());
      }
      assert(p.IsInPlane(p1));
      assert(p.IsInPlane(p2));
      assert(p.IsInPlane(p3));
    }
  }

  if (verbose) TRACE("IsInPlane, for bug #228");
  {
    //Values copied literally, so no .0 on purpose
    const Coordinat3D p1(-5,-5,0);
    const Coordinat3D p2(-5,-0.999999999999999880,0);
    const Coordinat3D p3(-5,-5,10);
    const Coordinat3D p4(-5,-0.999999999999999880,10);
    const PlaneX p(p1,p2,p3);
    if (verbose)
    {
      TRACE("----------------------------");
      TRACE(p.CalcMaxError(p1));
      TRACE(p.CalcError(p1));
      TRACE(p.CalcMaxError(p2));
      TRACE(p.CalcError(p2));
      TRACE(p.CalcMaxError(p3));
      TRACE(p.CalcError(p3));
      TRACE(p.CalcMaxError(p4));
      TRACE(p.CalcError(p4));
      TRACE(p.GetFunctionA());
      TRACE(p.GetFunctionB());
      TRACE(p.GetFunctionC());
      TRACE(p.GetCoefficients()[0]);
      TRACE(p.GetCoefficients()[1]);
      TRACE(p.GetCoefficients()[2]);
      TRACE(p.GetCoefficients()[3]);
    }
    assert(p.IsInPlane(p1));
    assert(p.IsInPlane(p2));
    assert(p.IsInPlane(p3));
    assert(p.IsInPlane(p4));
  }
  if (verbose) TRACE("GetProjection");
  {
    /*

    A: (0,0,1)                  A: ( 0,1)
    B: (1,0,0)                  B: ( 0,0)
    C: (1,1,0)                  C: (SQRT(2),0)

    |    /
    |   /                          |
    A-----C                        A
    |\/  /      -> becomes ->      |\
    |/\ /                          | \
  --+--B---                    ----B--C--
   /|                              |
  / |                              |
    */
    const std::vector<Coordinat2D> v {
      PlaneX().CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      )
    };
    assert(v.size() == 3);
    assert(abs(get<0>(v[0]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[0]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[1]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[1]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[2]) - std::sqrt(2.0) ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[2]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace

  }
  if (verbose) TRACE("GetProjection, for X = 2 plane");
  {
    /*

    A: (2,0,0)                  A: (0,0)
    B: (2,1,0)                  B: (1,0)
    C: (2,0,1)                  C: (0,1)

    |######/                 |
    C#####/                  C
    |\###/                   |\
    |#\#/                    | \
    |##B       -> becomes -> |  \
    |#/                      |   \
    |/                       |    \
---A--------             ---A-----B---
   /|                        |
  / |                        |
/  |                        |

    */
    const std::vector<Coordinat2D> v {
      PlaneX(
        Coordinat3D(0.0+2.0,0.0,0.0),
        Coordinat3D(0.0+2.0,1.0,0.0),
        Coordinat3D(0.0+2.0,0.0,1.0)

      ).CalcProjection(
        {
          Coordinat3D(0.0+2.0,0.0,0.0),
          Coordinat3D(0.0+2.0,1.0,0.0),
          Coordinat3D(0.0+2.0,0.0,1.0)
        }
      )
    };
    assert(v.size() == 3);
    assert(abs(get<0>(v[0]) - 0.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[0]) - 0.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[1]) - 1.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[1]) - 0.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[2]) - 0.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[2]) - 1.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
}
#endif

 

 

 

 

 

./CppPlane/planey.h

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.htm
//---------------------------------------------------------------------------
#ifndef RIBI_PLANEY_H
#define RIBI_PLANEY_H

#include <vector>
#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 <boost/geometry/geometries/point_xy.hpp>
#include <boost/make_shared.hpp>
//#include "planez.h"

#include "apfloat.h"
#pragma GCC diagnostic pop

namespace ribi {

struct PlaneZ;

///A 3D plane that can have its X expressed as a function of Y and Z.
///Can be constructed from its equation and at least three 3D points
///A plane stores its coefficients in the following form:
/// A.z + B.x + C.y = D
///Converting this to y being a function of x and z:
/// -C.y =  A  .z + B  .x - D
///    y = -A/C.z - B/C.x + D/C
///Where A,B,C and D can be obtained by using GetCoefficients
///Easier might be to express Y as:
///y = Ax + Bz + C
///Where A,B and C can be obtained by using GetFunctionA/B/C

struct PlaneY
{
  typedef boost::geometry::model::d2::point_xy<apfloat> Coordinat2D;
  typedef boost::geometry::model::point<apfloat,3,boost::geometry::cs::cartesian> Coordinat3D;
  typedef std::vector<Coordinat2D> Coordinats2D;
  typedef std::vector<Coordinat3D> Coordinats3D;
  typedef apfloat Double;
  typedef std::vector<Double> Doubles;

  ///Create plane Y = 0.0
  PlaneY() noexcept;

  ///Construct from three points
  explicit PlaneY(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  );

  Double CalcError(const Coordinat3D& coordinat) const noexcept;

  Double CalcMaxError(const Coordinat3D& coordinat) const noexcept;

  ///Get the 2D projection of these 3D points,
  /*

  A: (0,0,1)                  A: (0,0)
  B: (1,0,0)                  B: (SQRT(2),0)
  C: (1,1,0)                  C: (SQRT(2),SQRT(2))

  |    /
  |   /                       |
  A-----C                     |  C
  |\/  /      -> becomes ->   | /|
  |/\ /                       |/ |
  +--B---                     A--B-----

  */
  Coordinats2D CalcProjection(const Coordinats3D& points) const;

  ///Throws when cannot calculate Y, which is when the plane is horizontal
  Double CalcY(const Double& y, const Double& z) const;

  ///y = Ax + Bz + C
  ///Will throw if A cannot be calculated
  Double GetFunctionA() const;

  ///y = Ax + Bz + C
  ///Will throw if B cannot be calculated
  Double GetFunctionB() const;

  ///y = Ax + Bz + C
  ///Will throw if C cannot be calculated
  Double GetFunctionC() const;

  Doubles GetCoefficients() const noexcept;

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

  ///Checks if the coordinat is in the plane
  bool IsInPlane(const Coordinat3D& coordinat) const noexcept;

  private:
  ~PlaneY() noexcept;

  ///A PlaneY is actually a PlaneZ used with its coordinats rotated from (X,Y,Z) to (Z,Y,Y)
  const std::unique_ptr<PlaneZ> m_plane_z;

  ///Calculates m_min_error per GetFunctionC()
  static Double CalcMinErrorPerC() noexcept;

  static std::unique_ptr<PlaneZ> Create(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  );

  static Doubles Rotate(const Doubles& coefficients) noexcept;

  static Coordinat3D Rotate(const Coordinat3D& point) noexcept;

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

  ///Convert the PlaneY to a y(x,z), e.g 'y=(2*x) + (3*z) + 5' (spaces exactly as shown)
  std::string ToFunction() const;

  friend void boost::checked_delete<>(      PlaneY*);
  friend void boost::checked_delete<>(const PlaneY*);
  friend struct std::default_delete<      PlaneY>;
  friend struct std::default_delete<const PlaneY>;
  friend class boost::detail::sp_ms_deleter<      PlaneY>;
  friend class boost::detail::sp_ms_deleter<const PlaneY>;
  friend std::ostream& operator<<(std::ostream& os,const PlaneY& planey);
};

std::ostream& operator<<(std::ostream& os,const PlaneY& planey);

} //~namespace ribi

#endif // RIBI_PLANEY_H

 

 

 

 

 

./CppPlane/planey.cpp

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.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 "planey.h"

#include <cassert>

#include "geometry.h"
#include "planez.h"
#include "trace.h"
#pragma GCC diagnostic pop

ribi::PlaneY::PlaneY() noexcept
  : PlaneY(
    Coordinat3D(0.0,0.0,0.0),
    Coordinat3D(1.0,0.0,0.0),
    Coordinat3D(0.0,0.0,1.0)
  )
{
  #ifndef NDEBUG
  Test();
  #endif
}

ribi::PlaneY::PlaneY(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
)
  : m_plane_z{Create(p1,p2,p3)}
{
  #ifndef NDEBUG
  Test();
  #endif
}

ribi::PlaneY::~PlaneY()
{
  //OK
}

apfloat ribi::PlaneY::CalcError(const Coordinat3D& coordinat) const noexcept
{
  const apfloat x = boost::geometry::get<0>(coordinat);
  const apfloat y = boost::geometry::get<1>(coordinat);
  const apfloat z = boost::geometry::get<2>(coordinat);
  const auto expected = y;
  const auto calculated = CalcY(x,z);
  const auto error = abs(calculated - expected);
  return error;
}

ribi::PlaneY::Double ribi::PlaneY::CalcMinErrorPerC() noexcept
{
  //min_error_per_c will be about 0.000000001
  //stub_value increases this jut a little, by a 0.000001%
  const double stub_value = 0.000000001 * 1.00000001;
  #define USE_STUB
  #ifdef USE_STUB
  return stub_value;
  #else //USE_STUB
  //PlaneX calculates its own tolerance for errors, by measuring it
  static Double min_error_per_c = 0.0;
  if (min_error_per_c > 0.0) return min_error_per_c;

  //const double low = std::numeric_limits<double>::denorm_min();
  //const double high = std::numeric_limits<double>::max();
  const double low  = 1.0e-16;
  const double high = 1.0e+16;
  const double min_x = low;
  const double max_x = high;
  const double min_y = low;
  const double max_y = high;
  const double min_z = low;
  const double max_z = high;
  const apfloat zero(0.0);

  for (double z = min_z; z < max_z; z*=10.0)
  {
    for (double y = min_y; y < max_y; y*=10.0)
    {
      for (double x = min_x; x < max_x; x*=10.0)
      {
        const Coordinat3D p1(0.0,  y,0.0);
        const Coordinat3D p2(  x,  y,0.0);
        const Coordinat3D p3(0.0,  y,z);
        const PlaneY p(p1,p2,p3);
        for (const auto& p4: { p1, p2, p3 } )
        {
          const auto error = p.CalcError(p4);
          const auto error_per_c = error / p.GetFunctionC();
          assert(error_per_c >= zero);
          //TRACE(apfloat(min_error_per_c) / p.GetFunctionC());
          if (error_per_c > min_error_per_c)
          {
            min_error_per_c = error_per_c;
            //TRACE(min_error_per_c);
            //TRACE(x);
            //TRACE(y);
            //TRACE(z);
            //TRACE(p.GetFunctionC());
            //TRACE(apfloat(min_error_per_c) / p.GetFunctionC());
            //std::stringstream s;
            //s << Geometry().ToStr(p4) << " " << min_error;
            //TRACE(s.str());
          }
        }
      }
    }
    //TRACE(min_error_per_c);
  }
  //TRACE(min_error_per_c); //Output: TRACE 'min_error_per_c' line 127 in file '..\..\Classes\CppPlane\planey.cpp': '0.000000001e0'
  assert(min_error_per_c > zero);
  assert(min_error_per_c < stub_value);
  assert(min_error_per_c > 0.99 * stub_value);
  return min_error_per_c;
  #endif // USE_STUB
}

apfloat ribi::PlaneY::CalcMaxError(const Coordinat3D& /*coordinat*/) const noexcept
{
  assert(CalcMinErrorPerC() > 0.0);
  const auto max_error = abs(CalcMinErrorPerC() * GetFunctionC());
  assert(max_error >= 0.0);
  return max_error;
  /*
  const apfloat x = boost::geometry::get<0>(coordinat);
  //const apfloat y = boost::geometry::get<1>(coordinat);
  const apfloat z = boost::geometry::get<2>(coordinat);
  const auto coefficients = GetCoefficients();
  const auto a = coefficients[0];
  const auto b = coefficients[1];
  const auto c = coefficients[2];
  const apfloat e = boost::numeric::bounds<double>::smallest();
  assert(e > 0.0);
  //    y = -A/C.z - B/C.x + D/C
  //If C is zero, the slope in X and Y cannot be calculated
  if (c.sign())
  {
    const auto rc_x = -b / c;
    const auto rc_z = -a / c;
    const auto max_error_x = abs(e * rc_x * x) + 0.0;
    const auto max_error_y = 0.0;
    const auto max_error_z = abs(e * rc_z * z) + 0.0;
    const auto max_error = max_error_x + max_error_y + max_error_z;
    assert(max_error >= 0.0);
    return max_error;
  }
  assert(e > 0.0);
  return e;
  */
}

ribi::PlaneY::Coordinats2D ribi::PlaneY::CalcProjection(
  const Coordinats3D& points
) const
{
  auto v(points);
  for(auto& i: v) { i = Rotate(i); }
  try
  {
    return m_plane_z->CalcProjection(v);
  }
  catch (std::logic_error&)
  {
    throw std::logic_error("PlaneY::CalcProjection: cannot calculate projection");
  }
}

ribi::PlaneY::Double ribi::PlaneY::CalcY(const Double& x, const Double& z) const
{
  try
  {
    return m_plane_z->CalcZ(x,z);
  }
  catch (std::logic_error&)
  {
    throw std::logic_error("ribi::PlaneY::CalcY: cannot calculate Y of a horizontal plane");
  }
}

std::unique_ptr<ribi::PlaneZ> ribi::PlaneY::Create(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
)
{
  std::unique_ptr<PlaneZ> p(
    new PlaneZ(Rotate(p1), Rotate(p2), Rotate(p3))
  );
  assert(p);
  return p;
}

std::vector<apfloat> ribi::PlaneY::GetCoefficients() const noexcept
{
  const auto v(m_plane_z->GetCoefficients());
  assert(v.size() == 4);
  return { v[1],v[2],v[0],v[3] };
}

ribi::PlaneY::Double ribi::PlaneY::GetFunctionA() const
{
  assert(m_plane_z);
  return m_plane_z->GetFunctionA();
}

ribi::PlaneY::Double ribi::PlaneY::GetFunctionB() const
{
  assert(m_plane_z);
  return m_plane_z->GetFunctionB();
}

ribi::PlaneY::Double ribi::PlaneY::GetFunctionC() const
{
  assert(m_plane_z);
  return m_plane_z->GetFunctionC();
}

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

std::vector<std::string> ribi::PlaneY::GetVersionHistory() const noexcept
{
  return {
    "2014-03-10: version 1.0: initial version, split off from PlaneX",
    "2014-03-13: version 1.1: bug fixed",
    "2014-04-01: version 1.2: use of std::unique_ptr",
    "2014-06-13: version 1.3: shortened time to compile",
    "2014-07-03: version 1.4: use of apfloat",
    "2014-07-09: version 1.5: use double in interface only",
    "2014-07-10: version 1.6: use of apfloat only"
  };
}

bool ribi::PlaneY::IsInPlane(const Coordinat3D& coordinat) const noexcept
{
  try
  {
    const apfloat error = CalcError(coordinat);
    const apfloat max_error = CalcMaxError(coordinat);
    return error <= max_error;
  }
  catch (std::exception& e)
  {
    TRACE("ERROR");
    TRACE(e.what());
    assert(!"Should not get here");
    throw;
  }
}

ribi::PlaneY::Doubles ribi::PlaneY::Rotate(const Doubles& coefficients) noexcept
{
  assert(coefficients.size() == 4);
  return
  {
    coefficients[1],
    coefficients[2],
    coefficients[0],
    coefficients[3]
  };
}

ribi::PlaneY::Coordinat3D ribi::PlaneY::Rotate(
  const Coordinat3D& point
) noexcept
{
  //The 0-2-1 order is confirmed by doing a projection of a triangle on the Y=0 plane
  //on a Y=0 plane
  return Coordinat3D(
    boost::geometry::get<0>(point),
    boost::geometry::get<2>(point),
    boost::geometry::get<1>(point)
  );
}

std::string ribi::PlaneY::ToFunction() const
{
  std::stringstream s;
  s << (*this);
  return s.str();
}

std::ostream& ribi::operator<<(std::ostream& os,const PlaneY& planey)
{
  assert(planey.m_plane_z);
  try
  {
    os
      << "y=("
      << planey.m_plane_z->GetFunctionA() << "*x) + ("
      << planey.m_plane_z->GetFunctionB() << "*z) + "
      << planey.m_plane_z->GetFunctionC()
    ;
  }
  catch (std::logic_error&)
  {
    const std::string error
      = "ribi::PlaneY::ToFunction: cannot calculate X of a horizontal plane";
    throw std::logic_error(error.c_str());
  }
  return os;
}

 

 

 

 

 

./CppPlane/planey_test.cpp

 

#include "planey.h"

#include <cassert>

#include "geometry.h"
#include "planez.h"
#include "testtimer.h"
#include "trace.h"

#ifndef NDEBUG
void ribi::PlaneY::Test() noexcept
{
  {
    static bool is_tested { false };
    if (is_tested) return;
    is_tested = true;
  }
  const TestTimer test_timer(__func__,__FILE__,1.0);
  using boost::geometry::get;

  const bool verbose{false};
  if (verbose) TRACE("Default construction");
  {
    const PlaneY p;
    assert(!p.ToFunction().empty());
    assert(!p.GetCoefficients().empty());
  }
  if (verbose) TRACE("GetProjection, for plane Y = 0.0, points on plane");
  {
    /*

    A: (0,0,1)                  A: (0,1)
    B: (0,0,0)                  B: (0,0)
    C: (1,0,0)                  C: (1,0)

    |    /
    |   /                       |
    A  /                        A
    |\/         -> becomes ->   |\
    |/\                         | \
  --B--C---                   --B--C-----
   /|                           |
  / |                           |

    */
    const Coordinat3D a(0.0,0.0,1.0);
    const Coordinat3D b(0.0,0.0,0.0);
    const Coordinat3D c(1.0,0.0,0.0);
    assert(PlaneY().ToFunction() == PlaneY(a,b,c).ToFunction()
      && "The three points are on the Y = 0 plane");
    const std::vector<Coordinat2D> v { PlaneY().CalcProjection( { a,b,c } ) };
    assert(v.size() == 3);
    assert(abs(get<0>(v[0]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[0]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[1]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[1]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[2]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[2]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
  /*
  {


    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 };
    PlaneY p(
      Point(p1_x,p1_y,p1_z),
      Point(p2_x,p2_y,p2_z),
      Point(p3_x,p3_y,p3_z)
    );
    const auto t(p.GetCoefficients());
    assert(t.size() == 4);
    const double a { t[0] };
    const double b { t[1] };
    const double c { t[2] };
    const double d { t[3] };
    const double a_expected {  30.0 };
    const double b_expected { -48.0 };
    const double c_expected {  17.0 };
    const double d_expected { -15.0 };
    assert(abs(a - a_expected) < 0.001);
    assert(abs(b - b_expected) < 0.001);
    assert(abs(c - c_expected) < 0.001);
    assert(abs(d - d_expected) < 0.001);
    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'
      ;

    }
    assert(abs(d - d_p1_expected) < 0.001);
    assert(abs(d - d_p2_expected) < 0.001);
    assert(abs(d - d_p3_expected) < 0.001);
  }
  */
  if (verbose) TRACE("CalcPlaneY");
  /*
  {
    //CalcPlaneY 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


    const Point3D p1(1.0,1.0,10.0);
    const Point3D p2(1.0,2.0,13.0);
    const Point3D p3(2.0,1.0,12.0);
    PlaneY p(p1,p2,p3);
    const auto t(p.GetCoefficients());
    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(abs(a - a_expected) < 0.001);
    assert(abs(b - b_expected) < 0.001);
    assert(abs(c - c_expected) < 0.001);
    assert(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(abs(d - d_p1_expected) < 0.001);
    assert(abs(d - d_p2_expected) < 0.001);
    assert(abs(d - d_p3_expected) < 0.001);
    TRACE(p.ToFunction());
  }
  */
  if (verbose) TRACE("CalcY, diagonal plane");
  {


    const Coordinat3D p1(1.0,2.0,3.0);
    const Coordinat3D p2(2.0,5.0,8.0);
    const Coordinat3D p3(3.0,7.0,11.0);
    PlaneY p(p1,p2,p3);
    assert( abs(p.CalcY(1.0, 3.0)- 2.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcY(2.0, 8.0)- 5.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcY(3.0,11.0)- 7.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
  if (verbose) TRACE("CalcY, vertical plane Y = 3.0");
  /*

    |####/
    |###/#
    |##/##
    |#/###
    |/####
---+-----
   /|
  / |
/  |

  */
  {


    const Coordinat3D p1( 2.0,3.0, 5.0);
    const Coordinat3D p2( 7.0,3.0,11.0);
    const Coordinat3D p3(13.0,3.0,17.0);
    PlaneY p(p1,p2,p3);
    assert( abs(p.CalcY(1.0,2.0)-3.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcY(3.0,5.0)-3.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcY(7.0,9.0)-3.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }

  if (verbose) TRACE("IsInPlane, Y = 1, zooming to smallest three points to determine a plane, point above origin");
  {
    for (double i = 1.0;
      i > 1.0e-8; //i > 0.0;
      i/=10.0
    )
    {
      const Coordinat3D p1(0.0,1.0,0.0);
      const Coordinat3D p2(  i,1.0,0.0);
      const Coordinat3D p3(0.0,1.0,  i);
      const PlaneY p(p1,p2,p3);
      if (verbose)
      {
        TRACE("----------------------------");
        TRACE(i);
        TRACE(p.CalcMaxError(p1));
        TRACE(p.CalcError(p1));
        TRACE(p.CalcMaxError(p2));
        TRACE(p.CalcError(p2));
        TRACE(p.CalcMaxError(p3));
        TRACE(p.CalcError(p3));
        TRACE(p.GetFunctionA());
        TRACE(p.GetFunctionB());
        TRACE(p.GetFunctionC());
        TRACE(p.GetCoefficients()[0]);
        TRACE(p.GetCoefficients()[1]);
        TRACE(p.GetCoefficients()[2]);
        TRACE(p.GetCoefficients()[3]);
        TRACE(std::numeric_limits<double>::epsilon());
        TRACE(std::sqrt(std::numeric_limits<double>::epsilon()));
        TRACE(std::numeric_limits<double>::denorm_min());
      }
      assert(p.IsInPlane(p1));
      assert(p.IsInPlane(p2));
      assert(p.IsInPlane(p3));
    }
  }
  if (verbose) TRACE("IsPlane");
  {
    const Coordinat3D p1(-3.64472,-0.25,0.0);
    const Coordinat3D p2(-4.52988,-0.25,0.0);
    const Coordinat3D p3(-3.64472,-0.25,10.0);
    const Coordinat3D p4(-4.52988,-0.25,10.0);
    const PlaneY p(p1,p2,p3);
    if (verbose)
    {
      TRACE("----------------------------");
      TRACE(Geometry().ToStrSafe(p.CalcMaxError(p1)));
      TRACE(Geometry().ToStrSafe(p.CalcError(p1)));
      TRACE(Geometry().ToStrSafe(p.CalcMaxError(p2)));
      TRACE(Geometry().ToStrSafe(p.CalcError(p2)));
      TRACE(Geometry().ToStrSafe(p.CalcMaxError(p3)));
      TRACE(Geometry().ToStrSafe(p.CalcError(p3)));
      TRACE(Geometry().ToStrSafe(p.CalcMaxError(p4)));
      TRACE(Geometry().ToStrSafe(p.CalcError(p4)));
      TRACE(Geometry().ToStrSafe(p.GetFunctionA()));
      TRACE(Geometry().ToStrSafe(p.GetFunctionB()));
      TRACE(Geometry().ToStrSafe(p.GetFunctionC()));
      TRACE(Geometry().ToStrSafe(p.GetCoefficients()[0]));
      TRACE(Geometry().ToStrSafe(p.GetCoefficients()[1]));
      TRACE(Geometry().ToStrSafe(p.GetCoefficients()[2]));
      TRACE(Geometry().ToStrSafe(p.GetCoefficients()[3]));
    }
    assert(p.IsInPlane(p1));
    assert(p.IsInPlane(p2));
    assert(p.IsInPlane(p3));
    assert(p.IsInPlane(p4));
    const Coordinats3D v{p1,p2,p3,p4};
    assert(Geometry().IsPlane(v));
  }
  if (verbose) TRACE("ToFunction, 3 points and 4 points");
  {
    std::function<double(double,double)> f {
      [](const double x, const double z)
      {
        return (2.0 * x) + (3.0 * z) + 5.0;
      }
    };


    const double x1 { 2.0 };
    const double z1 { 3.0 };
    const double x2 { 5.0 };
    const double z2 { 7.0 };
    const double x3 { 11.0 };
    const double z3 { 13.0 };
    const double x4 { 17.0 };
    const double z4 { 29.0 };
    const Coordinat3D p1(x1,f(x1,z1),z1);
    const Coordinat3D p2(x2,f(x2,z2),z2);
    const Coordinat3D p3(x3,f(x3,z3),z3);
    const PlaneY a(p1,p2,p3);
    //assert(a.ToFunction() == "y=(2*x) + (3*z) + 5");
    assert(!a.ToFunction().empty());
    const Coordinat3D p4(x4,f(x4,z4),z4);
    assert(a.ToFunction() == PlaneY(p1,p2,p4).ToFunction());
    assert(a.ToFunction() == PlaneY(p1,p3,p4).ToFunction());
    assert(a.ToFunction() == PlaneY(p1,p4,p3).ToFunction());
    assert(a.ToFunction() == PlaneY(p2,p1,p3).ToFunction());
    assert(a.ToFunction() == PlaneY(p2,p1,p4).ToFunction());
    assert(a.ToFunction() == PlaneY(p2,p3,p1).ToFunction());
    assert(a.ToFunction() == PlaneY(p2,p3,p4).ToFunction());
    assert(a.ToFunction() == PlaneY(p2,p4,p1).ToFunction());
    assert(a.ToFunction() == PlaneY(p2,p4,p3).ToFunction());
    assert(a.ToFunction() == PlaneY(p3,p1,p2).ToFunction());
    assert(a.ToFunction() == PlaneY(p3,p1,p4).ToFunction());
    assert(a.ToFunction() == PlaneY(p3,p2,p1).ToFunction());
    assert(a.ToFunction() == PlaneY(p3,p2,p4).ToFunction());
    assert(a.ToFunction() == PlaneY(p3,p4,p1).ToFunction());
    assert(a.ToFunction() == PlaneY(p3,p4,p2).ToFunction());
    assert(a.ToFunction() == PlaneY(p4,p1,p2).ToFunction());
    assert(a.ToFunction() == PlaneY(p4,p1,p3).ToFunction());
    assert(a.ToFunction() == PlaneY(p4,p2,p1).ToFunction());
    assert(a.ToFunction() == PlaneY(p4,p2,p3).ToFunction());
    assert(a.ToFunction() == PlaneY(p4,p3,p1).ToFunction());
    assert(a.ToFunction() == PlaneY(p4,p3,p2).ToFunction());
  }


  if (verbose) TRACE("GetProjection, for plane Y = 0.0, points not on plane");
  {
    /*

    A: (0,0,1)                  A: (0,1)
    B: (1,0,0)                  B: (1,0)
    C: (1,1,0)                  C: (SQRT(2),0)

    |    /
    |   /                       |
    A-----C                     A
    |\/  /      -> becomes ->   |\
    |/\ /                       | \
  --+--B---                   --B--C-----
   /|                           |
  / |                           |

    */
    const std::vector<Coordinat2D> v {
      PlaneY().CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      )
    };
    assert(v.size() == 3);
    assert(abs(get<0>(v[0]) -           0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[0]) -           1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[1]) -           1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[1]) -           0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[2]) - std::sqrt(2.0)) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[2]) -           0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
}
#endif

 

 

 

 

 

./CppPlane/planez.h

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.htm
//---------------------------------------------------------------------------
#ifndef RIBI_PLANEZ_H
#define RIBI_PLANEZ_H

#include <vector>
#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 <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/make_shared.hpp>
#ifndef _WIN32
#include <boost/geometry/geometries/polygon.hpp>
#endif
#include "apfloat.h"
#pragma GCC diagnostic pop

namespace ribi {

///A 3D plane that can have its Z expressed as a function of X and Y;
///A 3D plane that encompasses all X and Y coordinats;
///Examples are:
/// Z = X + Y
/// Z = X
/// Z = Y
/// Z = 0
///Can be constructed from its equation and at least three 3D points
//A plane stores its 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
//
//
struct PlaneZ
{
  typedef boost::geometry::model::d2::point_xy<apfloat> Coordinat2D;
  typedef boost::geometry::model::point<apfloat,3,boost::geometry::cs::cartesian> Coordinat3D;
  typedef std::vector<Coordinat2D> Coordinats2D;
  typedef std::vector<Coordinat3D> Coordinats3D;
  typedef apfloat Double;
  typedef std::vector<Double> Doubles;

  ///Create plane Z = 0.0
  PlaneZ() noexcept;


  ///Construct from three points
  explicit PlaneZ(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  );

  Double CalcError(const Coordinat3D& coordinat) const noexcept;

  ///Get the 2D projection of these 3D points,
  ///Assumes these are in a plane
  /*

  A: (0,0,1)                  A: (0,0)
  B: (1,0,0)                  B: (SQRT(2),0)
  C: (1,1,0)                  C: (SQRT(2),SQRT(2))

  |    /
  |   /                       |
  A-----C                     |  C
  |\/  /      -> becomes ->   | /|
  |/\ /                       |/ |
  +--B---                     A--B-----

  */
  Coordinats2D CalcProjection(const Coordinats3D& points) const;

  ///Calculates the maximum allowed error for that coordinat for it to be in the plane
  Double CalcMaxError(const Coordinat3D& coordinat) const noexcept;

  ///Throws when cannot calculate Z, which is when the plane is vertical
  Double CalcZ(const Double& x, const Double& y) const;

  const Doubles& GetCoefficients() const noexcept { return m_coefficients; }

  ///This plane has equation 'z = Ax + By + C'
  ///Will throw if A cannot be calculated
  Double GetFunctionA() const;

  ///This plane has equation 'z = Ax + By + C'
  ///Will throw if B cannot be calculated
  Double GetFunctionB() const;

  ///This plane has equation 'z = Ax + By + C'
  ///Will throw if C cannot be calculated
  Double GetFunctionC() const;

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

  ///Checks if the coordinat is in the plane
  bool IsInPlane(const Coordinat3D& coordinat) const noexcept;

  ///Obtain a testing series of doubles (to be used as coordinat elements)
  ///in increasing order of difficulty
  static std::vector<double> GetTestSeries() noexcept;

  private:
  ///Construct from its coefficients
  PlaneZ(const Doubles& coefficients);

  ~PlaneZ() noexcept;

  //m_coefficients.size == 4
  const Doubles m_coefficients;

  ///Calculates m_min_error per GetFunctionC()
  static Double CalcMinErrorPerC() noexcept;

  static Doubles CalcPlaneZ(
    const Coordinat3D& p1,
    const Coordinat3D& p2,
    const Coordinat3D& p3
  ) noexcept;

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

  ///Convert the Plane to function z(x,y), e.g
  ///'z=(2*x) + (3*y) + 5' (spaces exactly as shown)
  ///Where 2,3 and 5 can be obtained with GetFunctionA,GetFunctionB and GetFunctionC
  ///respectively
  std::string ToFunction() const;

  friend void boost::checked_delete<>(      PlaneZ*);
  friend void boost::checked_delete<>(const PlaneZ*);
  friend struct std::default_delete<      PlaneZ>;
  friend struct std::default_delete<const PlaneZ>;
  friend class boost::detail::sp_ms_deleter<      PlaneZ>;
  friend class boost::detail::sp_ms_deleter<const PlaneZ>;
};

std::ostream& operator<<(std::ostream& os,const PlaneZ& planez);

} //~namespace ribi

#endif // RIBI_PLANEZ_H

 

 

 

 

 

./CppPlane/planez.cpp

 

//---------------------------------------------------------------------------
/*
Plane, 3D plane class
Copyright (C) 2013-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/CppPlane.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 "planez.h"

#include <cassert>

#include "container.h"
#include "geometry.h"
#include "trace.h"
#pragma GCC diagnostic pop

ribi::PlaneZ::PlaneZ() noexcept
  : PlaneZ(
    Coordinat3D(0.0,0.0,0.0),
    Coordinat3D(1.0,0.0,0.0),
    Coordinat3D(0.0,1.0,0.0)
  )
{
  #ifndef NDEBUG
  Test();
  #endif
}

ribi::PlaneZ::PlaneZ(
  const Doubles& coefficients
)
  : m_coefficients(coefficients)
{
  #ifndef NDEBUG
  Test();
  #endif
  const bool verbose{false};
  assert(GetCoefficients().size() == 4);

  if (m_coefficients[2] == 0.0)
  {
    if (verbose) { TRACE(Container().ToStr(m_coefficients)); }
    throw std::logic_error("PlaneZ (from coeffients) cannot be expressed in less than 3D space");
  }

  try
  {
    //TRACE(GetFunctionA());
    //TRACE(GetFunctionB());
    //TRACE(GetFunctionC());
    assert(GetFunctionA() == 0.0 || GetFunctionA() != 0.0);
    assert(GetFunctionB() == 0.0 || GetFunctionB() != 0.0);
    assert(GetFunctionC() == 0.0 || GetFunctionC() != 0.0);
  }
  catch (...)
  {
    assert(!"Should not get here");
  }
}

ribi::PlaneZ::PlaneZ(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) : PlaneZ(CalcPlaneZ(p1,p2,p3))
{
  #ifndef NDEBUG
  Test();
  #endif
  const bool verbose{false};
  assert(GetCoefficients().size() == 4);

  if (m_coefficients[2] == 0.0)
  {
    if (verbose) { TRACE(Container().ToStr(m_coefficients)); }
    throw std::logic_error("Plane (from points) that can be expressed in less than 3D space");
  }
}

ribi::PlaneZ::~PlaneZ() noexcept
{
  //OK
}

apfloat ribi::PlaneZ::CalcError(const Coordinat3D& coordinat) const noexcept
{
  const apfloat x = boost::geometry::get<0>(coordinat);
  const apfloat y = boost::geometry::get<1>(coordinat);
  const apfloat z = boost::geometry::get<2>(coordinat);
  const auto expected = z;
  const auto calculated = CalcZ(x,y);
  const auto error = abs(calculated - expected);
  return error;
}

ribi::PlaneZ::Double ribi::PlaneZ::CalcMinErrorPerC() noexcept
{
  //min_error_per_c will be about 0.000000001
  //stub_value increases this jut a little, by a 0.000001%
  const double stub_value = 0.000000001 * 1.00000001;
  #define USE_STUB
  #ifdef USE_STUB
  return stub_value;
  #else //USE_STUB
  //PlaneZ calculates its own tolerance for errors, by measuring it
  static Double min_error_per_c = 0.0;
  if (min_error_per_c > 0.0) return min_error_per_c;

  //const double low = std::numeric_limits<double>::denorm_min();
  //const double high = std::numeric_limits<double>::max();
  const double low  = 1.0e-16;
  const double high = 1.0e+16;
  const double min_x = low;
  const double max_x = high;
  const double min_y = low;
  const double max_y = high;
  const double min_z = low;
  const double max_z = high;
  const apfloat zero(0.0);

  for (double z = min_z; z < max_z; z*=10.0)
  {
    for (double y = min_y; y < max_y; y*=10.0)
    {
      for (double x = min_x; x < max_x; x*=10.0)
      {
        const Coordinat3D p1(0.0,0.0,z);
        const Coordinat3D p2(0.0,  y,z);
        const Coordinat3D p3(  x,0.0,z);
        const PlaneZ p(p1,p2,p3);
        for (const auto& p4: { p1,p2,p3 } )
        {
          const auto error = p.CalcError(p4);
          const auto error_per_c = error / p.GetFunctionC();
          assert(error_per_c >= zero);
          if (error_per_c > min_error_per_c)
          {
            min_error_per_c = error_per_c;
            //TRACE(min_error_per_c);
            //TRACE(x);
            //TRACE(y);
            //TRACE(z);
            //TRACE(p.GetFunctionC());
            //TRACE(apfloat(min_error_per_c) / p.GetFunctionC());
            //std::stringstream s;
            //s << Geometry().ToStr(p4) << " " << min_error;
            //TRACE(s.str());
          }
        }
      }
    }
    //TRACE(min_error_per_c);
  }
  //TRACE(min_error_per_c);
  assert(min_error_per_c > zero);
  assert(min_error_per_c < stub_value);
  assert(min_error_per_c > 0.99 * stub_value);
  return min_error_per_c;
  #endif // USE_STUB
}

apfloat ribi::PlaneZ::CalcMaxError(const Coordinat3D& /*coordinat*/) const noexcept
{
  assert(CalcMinErrorPerC() > apfloat(0.0));
  const auto max_error = abs(CalcMinErrorPerC() * GetFunctionC());
  assert(max_error >= apfloat(0.0));
  return max_error;
  /*
  const bool verbose{false};
  const apfloat x = boost::geometry::get<0>(coordinat);
  const apfloat y = boost::geometry::get<1>(coordinat);
  //const apfloat z = boost::geometry::get<2>(coordinat);
  const auto coefficients = GetCoefficients();
  const auto a = coefficients[0];
  const auto b = coefficients[1];
  const auto c = coefficients[2];
  //const apfloat e = boost::numeric::bounds<double>::smallest();
  const apfloat e = CalcMinErrorPerC(); //std::sqrt(std::numeric_limits<double>::epsilon());
  assert(e > 0.0);

  if (verbose)
  {
    TRACE(Geometry().ToStr(coordinat));
    TRACE(x);
    TRACE(y);
    TRACE(a);
    TRACE(b);
    TRACE(c);
  }

  // z = -A/C.x - B/C.y + D/C = (-A.x - B.y + D) / C
  //If C is zero, the slope in X and Y cannot be calculated
  if (c.sign())
  {
    const auto rc_x = -a / c;
    const auto rc_y = -b / c;
    const auto max_error_x = abs(e * rc_x * x) + 0.0;
    const auto max_error_y = abs(e * rc_y * y) + 0.0;
    const auto max_error_z = e;
    const auto max_error = max_error_x + max_error_y + max_error_z;

    if (verbose)
    {
      TRACE(rc_x);
      TRACE(rc_y);
      TRACE(max_error_x);
      TRACE(max_error_y);
      TRACE(max_error);
    }
    assert(max_error > 0.0);
    return max_error;
  }
  assert(e > 0.0);
  return e;
  */
}

std::vector<apfloat> ribi::PlaneZ::CalcPlaneZ(
  const Coordinat3D& p1,
  const Coordinat3D& p2,
  const Coordinat3D& p3
) noexcept
{
  
  const auto v(
    Geometry().CalcPlane(
      p1,
      p2,
      p3
    )
  );
  assert(v.size() == 4);
  return v;
}

ribi::PlaneZ::Coordinats2D ribi::PlaneZ::CalcProjection(
  const Coordinats3D& points
) const
{
  assert(points.size() >= 3);
  const apfloat x_origin = 0.0;
  const apfloat y_origin = 0.0;
  const apfloat z_origin = CalcZ(x_origin,y_origin);

  Coordinats2D v;
  for (const auto& point: points)
  {
    const Double x(boost::geometry::get<0>(point));
    const Double y(boost::geometry::get<1>(point));
    const Double z(boost::geometry::get<2>(point));
    const Double dx =
      sqrt( //Apfloat does not add the std::
          ((x - x_origin) * (x - x_origin))
        + ((z - z_origin) * (z - z_origin))
      ) * (x - x_origin)
    ;
    const Double dy =
      sqrt( //Apfloat does not add the std::
          ((y - y_origin) * (y - y_origin))
        + ((z - z_origin) * (z - z_origin))
      ) * (y - y_origin)
    ;

    Coordinat2D point_xy(dx,dy);
    v.push_back(point_xy);
  }
  return v;
}

ribi::PlaneZ::Double ribi::PlaneZ::CalcZ(const Double& x, const Double& y) const
{
  // z = -A/C.x - B/C.y + D/C = (-A.x - B.y + D) / C
  const bool verbose{false};
  const auto a = m_coefficients[0];
  const auto b = m_coefficients[1];
  const auto c = m_coefficients[2];
  const auto d = m_coefficients[3];
  if (c == 0.0)
  {
    throw std::logic_error("ribi::PlaneZ::CalcZ: cannot calculate Z of a vertical plane");
  }
  const apfloat term1 = -a*x;
  const apfloat term2 = -b*y;
  const apfloat numerator = term1 + term2 + d;
  if (verbose)
  {
    TRACE(numerator);
    TRACE(c);
  }
  const apfloat result = numerator / c;
  return result;
}

apfloat ribi::PlaneZ::GetFunctionA() const
{
  const bool verbose{false};
  const auto coeff_a = m_coefficients[0];
  const auto coeff_c = m_coefficients[2];
  static const apfloat zero(0.0);
  assert(coeff_c != zero);
  if (verbose)
  {
    TRACE(coeff_c);
  }
  const auto a = -coeff_a/coeff_c;
  return a;
}

apfloat ribi::PlaneZ::GetFunctionB() const
{
  const auto coeff_b = m_coefficients[1];
  const auto coeff_c = m_coefficients[2];
  static const apfloat zero(0.0);
  assert(coeff_c != zero);
  const auto b = -coeff_b/coeff_c;
  return b;
}

apfloat ribi::PlaneZ::GetFunctionC() const
{
  const auto coeff_c = m_coefficients[2];
  const auto coeff_d = m_coefficients[3];
  static const apfloat zero(0.0);
  assert(coeff_c != zero);

  try
  {
    const auto c = coeff_d/coeff_c;
    return c;
  }
  catch (...)
  {
    assert(!"Should not get here");
    throw;
  }
}

std::vector<double> ribi::PlaneZ::GetTestSeries() noexcept
{
  return
  {
     1.0,
    -1.0,
     std::numeric_limits<double>::epsilon(),
    -std::numeric_limits<double>::epsilon(),
     1.e8,
    -1.e8,
    0.0
    // std::numeric_limits<double>::denorm_min(),
    //-std::numeric_limits<double>::denorm_min(),
    // 1.e64,
    //-1.e64,
    //std::numeric_limits<double>::min(),
    //std::numeric_limits<double>::max()
  };
}

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

std::vector<std::string> ribi::PlaneZ::GetVersionHistory() const noexcept
{
  return {
    "2014-03-10: version 1.0: initial version, split off from Plane",
    "2014-03-10: version 1.1: bug fixed, only occurred at debugging",
    "2014-03-13: version 1.2: bug fixed",
    "2014-04-01: version 1.3: use of std::unique_ptr",
    "2014-07-03: version 1.4: use of apfloat",
    "2014-07-09: version 1.5: use double in interface only"
    "2014-07-10: version 1.6: use of apfloat only"
  };
}

bool ribi::PlaneZ::IsInPlane(const Coordinat3D& coordinat) const noexcept
{
  try
  {
    const apfloat error = CalcError(coordinat);
    const apfloat max_error = CalcMaxError(coordinat);
    return error <= max_error;
  }
  catch (std::exception& e)
  {
    TRACE("ERROR");
    TRACE(e.what());
    assert(!"Should not get here");
    throw;
  }
}

std::string ribi::PlaneZ::ToFunction() const
{
  std::stringstream s;
  s << (*this);
  return s.str();
}

std::ostream& ribi::operator<<(std::ostream& os, const PlaneZ& planez)
{
  try
  {
    os
      << "z=("
      << planez.GetFunctionA() << "*x) + ("
      << planez.GetFunctionB() << "*y) + "
      << planez.GetFunctionC()
    ;
  }
  catch (std::logic_error&)
  {
    const std::string error
      = "ribi::PlaneZ::ToFunction: cannot calculate Z of a vertical plane";
    throw std::logic_error(error.c_str());
  }
  return os;
}

 

 

 

 

 

./CppPlane/planez_test.cpp

 

#include "planez.h"

#include <cassert>

#include "container.h"
#include "geometry.h"
#include "testtimer.h"
#include "trace.h"

#ifndef NDEBUG
void ribi::PlaneZ::Test() noexcept
{
  {
    static bool is_tested{false};
    if (is_tested) return;
    is_tested = true;
  }
  const TestTimer test_timer(__func__,__FILE__,1.0);
  using boost::geometry::get;
  const bool verbose{false};
  const bool show_warning{false};
  const auto series = GetTestSeries();

  if (verbose) { TRACE("Default construction"); }
  {
    const PlaneZ p;
    assert(!p.ToFunction().empty());
    assert(!p.GetCoefficients().empty());
  }

  if (verbose) TRACE("PlaneZ, Z = 5");
  {
    const Coordinat3D p1( 2.0, 3.0,5.0);
    const Coordinat3D p2( 7.0,11.0,5.0);
    const Coordinat3D p3(13.0,17.0,5.0);
    const PlaneZ p(p1,p2,p3);
    assert(
      !p.CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      ).empty()
    );
  }

  //IsInPlane for Z=0 plane
  if (verbose) TRACE("PlaneZ, preparation for Plane::CanCalcZ and Plane::IsInPlane, Z = 0 plane, from 1.0 coordinat");
  {
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,1.0,0.0);
    const Coordinat3D p3(1.0,0.0,0.0);
    const PlaneZ p(p1,p2,p3);
    assert( p.IsInPlane(Coordinat3D( 0.0, 0.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0, 1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,0.0)));
  }
  if (verbose) TRACE("PlaneZ, preparation for Plane::CanCalcZ and Plane::IsInPlane, Z = 0 plane, from smallest possible coordinat");
  {
    const double i = std::numeric_limits<double>::denorm_min();
    assert(i > 0.0);
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,  i,0.0);
    const Coordinat3D p3(  i,0.0,0.0);
    const PlaneZ p(p1,p2,p3);
    assert( p.IsInPlane(Coordinat3D( 0.0, 0.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0, 1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,0.0)));

  }
  if (verbose) TRACE("PlaneZ, preparation for Plane::CanCalcZ and Plane::IsInPlane, Z = 0 plane, from biggest possible coordinat");
  {
    const double i = std::numeric_limits<double>::max();
    assert(i > 0.0);
    const Coordinat3D p1(0.0,0.0,0.0);
    const Coordinat3D p2(0.0,  i,0.0);
    const Coordinat3D p3(  i,0.0,0.0);
    const PlaneZ p(p1,p2,p3);
    assert( p.IsInPlane(Coordinat3D( 0.0, 0.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0,-1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D(-1.0, 1.0,0.0)));
    assert( p.IsInPlane(Coordinat3D( 1.0, 1.0,0.0)));
  }
  if (verbose) TRACE("PlaneZ, preparation for Plane::CanCalcZ and Plane::IsInPlane, Z = 0 plane, zooming in");
  {
    for (const double i:series)
    {
      if (i == 0.0) continue;
      assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
      const Coordinat3D p1(0.0,0.0,0.0);
      const Coordinat3D p2(0.0,  i,0.0);
      const Coordinat3D p3(  i,0.0,0.0);
      const PlaneZ p(p1,p2,p3);
      for (const double j:series)
      {
        assert(p.IsInPlane(Coordinat3D(0.0,0.0,0.0)));
        assert(p.IsInPlane(Coordinat3D(  j,  j,0.0)));
        assert(p.IsInPlane(Coordinat3D(  j, -j,0.0)));
        assert(p.IsInPlane(Coordinat3D( -j,  j,0.0)));
        assert(p.IsInPlane(Coordinat3D( -j, -j,0.0)));
      }
    }
  }
  if (verbose) TRACE("IsInPlane, Z = 1, zooming to smallest three points to determine a plane, point above origin");
  {
    for (double i = 1.0;
      i > 1.0e-8; //i > 0.0;
      i/=10.0
    )
    {
      const Coordinat3D p1(0.0,0.0,1.0);
      const Coordinat3D p2(0.0,  i,1.0);
      const Coordinat3D p3(  i,0.0,1.0);
      const PlaneZ p(p1,p2,p3);
      if (verbose)
      {
        TRACE("----------------------------");
        TRACE(i);
        TRACE(p.CalcMaxError(p1));
        TRACE(p.CalcError(p1));
        TRACE(p.GetFunctionA());
        TRACE(p.GetFunctionB());
        TRACE(p.GetFunctionC());
        TRACE(p.GetCoefficients()[0]);
        TRACE(p.GetCoefficients()[1]);
        TRACE(p.GetCoefficients()[2]);
        TRACE(p.GetCoefficients()[3]);
        TRACE(std::numeric_limits<double>::epsilon());
        TRACE(std::sqrt(std::numeric_limits<double>::epsilon()));
        TRACE(std::numeric_limits<double>::denorm_min());
      }
      assert(p.IsInPlane(p1));
      assert(p.IsInPlane(p2));
      assert(p.IsInPlane(p3));
    }
  }
  if (verbose) TRACE("IsInPlane, Z = 1, zooming to smallest three points to determine a plane, point above origin");
  {
    const double min = 1.0e-8;
    const double max = 1.0e+8;
    for (double z = min; z < max; z*=10.0)
    {
      for (double i = min; i < max; i*=10.0)
      {
        const Coordinat3D p1(0.0,0.0,z);
        const Coordinat3D p2(0.0,  i,z);
        const Coordinat3D p3(  i,0.0,z);
        assert(i != 0.0);
        const PlaneZ p(p1,p2,p3);
        if ( (!p.IsInPlane(p1) || !p.IsInPlane(p2) || !p.IsInPlane(p3)))
        {
          std::stringstream s;
          s << "Warning: coordinats " << Geometry().ToStr(p1)
            << ", " << Geometry().ToStr(p2)
            << " and " << Geometry().ToStr(p3)
            << " are determined not to be in a PlaneZ that was created from points"
          ;
          if (abs(1.0 - (p.CalcMaxError(p1) / p.CalcError(p1))) < 0.01)
          {
            //Allow another percent of freedom
            continue;
          }
          if (!show_warning) continue;

          TRACE(s.str());
          assert(p.IsInPlane(p1) == p.IsInPlane(p2) && p.IsInPlane(p2) == p.IsInPlane(p3));
          TRACE(p);
          TRACE(z);
          TRACE(i);
          TRACE(p.CalcMaxError(p1));
          TRACE(p.CalcError(p1));
          TRACE(p.CalcMaxError(p1) / p.CalcError(p1));
          TRACE(p.CalcMaxError(p2));
          TRACE(p.CalcError(p2));
          TRACE(p.CalcMaxError(p2) / p.CalcError(p2));
          TRACE(p.CalcMaxError(p3));
          TRACE(p.CalcError(p3));
          TRACE(p.CalcMaxError(p3) / p.CalcError(p3));
          TRACE(p.GetFunctionA());
          TRACE(p.GetFunctionB());
          TRACE(p.GetFunctionC());
        }
        assert(p.IsInPlane(p1));
        assert(p.IsInPlane(p2));
        assert(p.IsInPlane(p3));
      }
    }
  }
  if (verbose) TRACE("CanCalcZ, Z = 1.0 plane, zooming in");
  {
    for (const double i:series)
    {
      if (i == 0.0) continue;
      assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
      const Coordinat3D p1(0.0,0.0,1.0);
      const Coordinat3D p2(0.0,  i,1.0);
      const Coordinat3D p3(  i,0.0,1.0);
      const PlaneZ p(p1,p2,p3);
      for (const double j:series)
      {
        if (!p.IsInPlane(Coordinat3D(0.0,0.0,1.0)))
        {
          TRACE(i);
          TRACE(p.CalcMaxError(Coordinat3D(0.0,0.0,1.0)));
          TRACE(p.CalcError(Coordinat3D(0.0,0.0,1.0)));
        }
        assert(p.IsInPlane(Coordinat3D(0.0,0.0,1.0)));
        assert(p.IsInPlane(Coordinat3D(  j,  j,1.0)));
        assert(p.IsInPlane(Coordinat3D(  j, -j,1.0)));
        assert(p.IsInPlane(Coordinat3D( -j,  j,1.0)));
        assert(p.IsInPlane(Coordinat3D( -j, -j,1.0)));
      }
    }
  }

  /*

    |    /#/##########
    |   B#/###########
    |  /#/############
    | /#/#############
    |/#/##############
    A-------C--------- Z = z
    |/
  --O----------------- Z = 0
   /|


  */
  if (verbose) TRACE("PlaneZ, preparation for Plane::CanCalcZ and Plane::IsInPlane, Z = z plane, zooming in");
  {
    //The height of the plane
    for (const double z:series)
    {
      //The distance from the origin, will be used by the two construction points
      for (const double i:series)
      {
        if (i == 0.0) continue;
        assert(i != 0.0 && "Cannot express plane when all its coordinats are at origin");
        const Coordinat3D p1(0.0,0.0,z);
        const Coordinat3D p2(0.0,  i,z);
        const Coordinat3D p3(  i,0.0,z);
        const PlaneZ p(p1,p2,p3);
        //The distance (actually, half the Manhattan distance) from the origin,
        //will be used by points tested to be in this plane
        for (const double j:series)
        {
          if (!p.IsInPlane(Coordinat3D(j,j,z)))
          {
            if (show_warning)
            {
              std::stringstream s;
              s << "Warning: coordinat " << Geometry().ToStr(Coordinat3D(j,j,z))
                << " is determined not to be in a PlaneZ that was created from points "
                << Geometry().ToStr(p1) << ", "
                << Geometry().ToStr(p2) << " and "
                << Geometry().ToStr(p3) << "."
              ;
              TRACE(s.str());
            }
            continue;
          }
          if (!p.IsInPlane(Coordinat3D(j,j,z)))
          {
            TRACE("ERROR");
            TRACE(z);
            TRACE(i);
            TRACE(j);
            TRACE(p.CalcError(Coordinat3D(j,j,z)));
            TRACE(p.CalcMaxError(Coordinat3D(j,j,z)));
            TRACE(p);

            TRACE("AGAIN");
            TRACE(z);
            TRACE(i);
            TRACE(j);
            TRACE(p.CalcError(Coordinat3D(j,j,z)));
            TRACE(p.CalcMaxError(Coordinat3D(j,j,z)));
            TRACE(p);

          }
          assert(p.IsInPlane(Coordinat3D(  j,  j,z)));
          assert(p.IsInPlane(Coordinat3D(  j, -j,z)));
          assert(p.IsInPlane(Coordinat3D( -j,  j,z)));
          assert(p.IsInPlane(Coordinat3D( -j, -j,z)));
        }
      }
    }
  }

  if (verbose) { TRACE("Check formulas"); }
  {
    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};
    PlaneZ p(
      Coordinat3D(p1_x,p1_y,p1_z),
      Coordinat3D(p2_x,p2_y,p2_z),
      Coordinat3D(p3_x,p3_y,p3_z)
    );
    const auto t(Geometry().ToDouble(p.GetCoefficients()));
    assert(t.size() == 4);
    const auto a = t[0];
    const auto b = t[1];
    const auto c = t[2];
    const auto d = t[3];
    const auto a_expected =  30.0;
    const auto b_expected = -48.0;
    const auto c_expected =  17.0;
    const auto d_expected = -15.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 auto d_p1_expected = (a * p1_x) + (b * p1_y) + (c * p1_z);
    const auto d_p2_expected = (a * p2_x) + (b * p2_y) + (c * p2_z);
    const auto d_p3_expected = (a * p3_x) + (b * p3_y) + (c * p3_z);
    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("CalcPlaneZ"); }
  {
    //CalcPlaneZ 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
    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 PlaneZ p(p1,p2,p3);
    const auto t(Geometry().ToDouble(p.GetCoefficients()));
    const auto a = t[0];
    const auto b = t[1];
    const auto c = t[2];
    const auto d = t[3];
    const auto a_expected = -2.0;
    const auto b_expected = -3.0;
    const auto c_expected =  1.0;
    const auto 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 auto d_p1_expected = (a * 1.0) + (b * 1.0) + (c * 10.0);
    const auto d_p2_expected = (a * 1.0) + (b * 2.0) + (c * 13.0);
    const auto 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("CalcZ, diagonal plane"); }
  {
    const Coordinat3D p1(1.0,2.0,3.0);
    const Coordinat3D p2(2.0,5.0,8.0);
    const Coordinat3D p3(3.0,7.0,11.0);
    const PlaneZ p(p1,p2,p3);
    assert(abs(p.CalcZ(1.0,2.0)- 3.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(p.CalcZ(2.0,5.0)- 8.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(p.CalcZ(3.0,7.0)-11.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
  if (verbose) { TRACE("CalcZ, horizontal plane Z = 5.0"); }
  /*

    |    /
    |   /#
    |  /##
    | /###
    |/####
---+-----
   /|
  / |
/  |

  */
  {


    const Coordinat3D p1( 2.0, 3.0,5.0);
    const Coordinat3D p2( 7.0,11.0,5.0);
    const Coordinat3D p3(13.0,17.0,5.0);
    const PlaneZ p(p1,p2,p3);
    assert( abs(p.CalcZ(1.0,2.0)-5.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcZ(3.0,5.0)-5.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert( abs(p.CalcZ(7.0,9.0)-5.0) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
  if (verbose) { TRACE("ToFunction, 3 points and 4 points"); }
  {
    std::function<double(double,double)> f {
      [](const double x, const double y)
      {
        return (2.0 * x) + (3.0 * y) + 5.0;
      }
    };


    const double x1 =  2.0;
    const double y1 =  3.0;
    const double x2 =  5.0;
    const double y2 =  7.0;
    const double x3 = 11.0;
    const double y3 = 13.0;
    const double x4 = 17.0;
    const double y4 = 29.0;
    const Coordinat3D p1(x1,y1,f(x1,y1));
    const Coordinat3D p2(x2,y2,f(x2,y2));
    const Coordinat3D p3(x3,y3,f(x3,y3));
    const PlaneZ a(p1,p2,p3);
    //assert(a.ToFunction() == "z=(2*x) + (3*y) + 5");
    assert(!a.ToFunction().empty());

    const Coordinat3D p4(x4,y4,f(x4,y4));
    assert(a.ToFunction() == PlaneZ(p1,p2,p4).ToFunction());
    assert(a.ToFunction() == PlaneZ(p1,p3,p4).ToFunction());
    assert(a.ToFunction() == PlaneZ(p1,p4,p3).ToFunction());
    assert(a.ToFunction() == PlaneZ(p2,p1,p3).ToFunction());
    assert(a.ToFunction() == PlaneZ(p2,p1,p4).ToFunction());
    assert(a.ToFunction() == PlaneZ(p2,p3,p1).ToFunction());
    assert(a.ToFunction() == PlaneZ(p2,p3,p4).ToFunction());
    assert(a.ToFunction() == PlaneZ(p2,p4,p1).ToFunction());
    assert(a.ToFunction() == PlaneZ(p2,p4,p3).ToFunction());
    assert(a.ToFunction() == PlaneZ(p3,p1,p2).ToFunction());
    assert(a.ToFunction() == PlaneZ(p3,p1,p4).ToFunction());
    assert(a.ToFunction() == PlaneZ(p3,p2,p1).ToFunction());
    assert(a.ToFunction() == PlaneZ(p3,p2,p4).ToFunction());
    assert(a.ToFunction() == PlaneZ(p3,p4,p1).ToFunction());
    assert(a.ToFunction() == PlaneZ(p3,p4,p2).ToFunction());
    assert(a.ToFunction() == PlaneZ(p4,p1,p2).ToFunction());
    assert(a.ToFunction() == PlaneZ(p4,p1,p3).ToFunction());
    assert(a.ToFunction() == PlaneZ(p4,p2,p1).ToFunction());
    assert(a.ToFunction() == PlaneZ(p4,p2,p3).ToFunction());
    assert(a.ToFunction() == PlaneZ(p4,p3,p1).ToFunction());
    assert(a.ToFunction() == PlaneZ(p4,p3,p2).ToFunction());

  }
  if (verbose) { TRACE("GetProjection, for Z = 0 plane"); }
  {
    /*

    A: (0,0,1)                  A: (0,0)
    B: (1,0,0)                  B: (1,0)
    C: (1,1,0)                  C: (1,1)

    |    /
    |   /                       |
    A-----C                     |  C
    |\/  /      -> becomes ->   | /|
    |/\ /                       |/ |
    +--B---                     A--B-----

    */
    const std::vector<Coordinat2D> v {
      PlaneZ().CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0),
          Coordinat3D(1.0,0.0,0.0),
          Coordinat3D(1.0,1.0,0.0)
        }
      )
    };
    assert(v.size() == 3);
    assert(abs(get<0>(v[0]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[0]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[1]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[1]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[2]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[2]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace

  }
  if (verbose) { TRACE("CalcProjection, for Z = 2 plane"); }
  {
    /*

    A: (0,0,1+2)                  A: (0,0)
    B: (1,0,0+2)                  B: (1,0)
    C: (1,1,0+2)                  C: (1,1)

    |    /
    |   /                       |
    A-----C                     |  C
    |\/  /      -> becomes ->   | /|
    |/\ /                       |/ |
    +--B---                     A--B-----

    */
    const std::vector<Coordinat2D> v {
      PlaneZ(
        Coordinat3D(0.0,0.0,0.0+2.0),
        Coordinat3D(0.0,1.0,0.0+2.0),
        Coordinat3D(1.0,0.0,0.0+2.0)
      ).CalcProjection(
        {
          Coordinat3D(0.0,0.0,1.0+2.0),
          Coordinat3D(1.0,0.0,0.0+2.0),
          Coordinat3D(1.0,1.0,0.0+2.0)
        }
      )
    };
    assert(v.size() == 3);
    assert(abs(get<0>(v[0]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[0]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[1]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[1]) - 0.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<0>(v[2]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
    assert(abs(get<1>(v[2]) - 1.0 ) < 0.001); //no std:: , as apfloat puts abs in the global namespace
  }
}
#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