/************************************************************************** This file is part of the C++ library "homog2d", dedicated to handle 2D lines and points, see https://github.com/skramm/homog2d Author & Copyright 2019-2024 Sebastien Kramm Contact: firstname.lastname@univ-rouen.fr Licence: MPL v2 This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/. **************************************************************************/ /** \file homog2d.hpp \brief single header file, implements some 2D homogeneous stuff. See https://github.com/skramm/homog2d */ #ifndef HG_HOMOG2D_HPP #define HG_HOMOG2D_HPP #define _USE_MATH_DEFINES #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // required for uint8_t #include // required for std::unique_ptr #ifdef HOMOG2D_USE_SVG_IMPORT #define HOMOG2D_ENABLE_VRTP // #include // why was that needed in the first place? #include "tinyxml2.h" #endif #ifdef HOMOG2D_ENABLE_VRTP #include #endif #ifdef HOMOG2D_USE_EIGEN #include #endif #ifdef HOMOG2D_USE_TTMATH #include #ifndef HOMOG2D_INUMTYPE #define HOMOG2D_INUMTYPE ttmath::Big<2,2> #endif #define homog2d_abs ttmath::Abs #define homog2d_sin ttmath::Sin #define homog2d_cos ttmath::Cos #define homog2d_acos ttmath::ACos #define homog2d_asin ttmath::ASin #define homog2d_atan ttmath::ATan #define homog2d_sqrt ttmath::Sqrt #else #define homog2d_abs std::abs #define homog2d_sin std::sin #define homog2d_cos std::cos #define homog2d_acos std::acos #define homog2d_asin std::asin #define homog2d_atan std::atan #define homog2d_sqrt std::sqrt #endif #ifdef HOMOG2D_USE_OPENCV #include "opencv2/imgproc.hpp" #include "opencv2/highgui.hpp" #endif #ifdef HOMOG2D_USE_BOOSTGEOM #include #endif #ifdef _MSC_VER #define HOMOG2D_PRETTY_FUNCTION __FUNCSIG__ #else #define HOMOG2D_PRETTY_FUNCTION __PRETTY_FUNCTION__ #endif #ifdef HOMOG2D_DEBUGMODE #define HOMOG2D_START std::cout << "START: line:" << __LINE__ \ << " func=\n" << HOMOG2D_PRETTY_FUNCTION << "()\n" // #define HOMOG2D_START std::cout << __FUNCTION__ << "()\n" #else #define HOMOG2D_START #endif #ifdef HOMOG2D_DEBUGMODE #define HOMOG2D_LOG(a) \ std::cout << '-' << __FUNCTION__ << "(), line " << __LINE__ << ": " \ << std::scientific << std::setprecision(10) \ << a << std::endl; #else #define HOMOG2D_LOG(a) {;} #endif #define HOMOG2D_ASSERT( a ) \ if( !(a) ) \ { \ std::cerr << "ASSERT FAILURE, line=" << __LINE__ << std::endl; \ std::exit(1); \ } /// Assert debug macro, used internally if \c HOMOG2D_DEBUGMODE is defined #define HOMOG2D_DEBUG_ASSERT(a,b) \ { \ if( (a) == false ) \ { \ std::cerr << "Homog2d assert failure, version:" << HOMOG2D_VERSION \ << ", line:" << __LINE__ << "\n -details: " << b << '\n'; \ std::cout << "homog2d: internal failure, please check stderr and report this on https://github.com/skramm/homog2d/issues\n"; \ std::exit(1); \ } \ } #define HOMOG2D_CHECK_ROW_COL \ if( r > 2 ) \ throw std::runtime_error( "Error: invalid row value: r=" + std::to_string(r) ); \ if( c > 2 ) \ throw std::runtime_error( "Error: invalid col value: r=" + std::to_string(r) ) #ifdef HOMOG2D_USE_TTMATH # define HOMOG2D_CHECK_IS_NUMBER(T) #else # define HOMOG2D_CHECK_IS_NUMBER(T) \ static_assert( (std::is_arithmetic::value && !std::is_same::value), "Type of value must be numerical" ) #endif #ifndef HOMOG2D_NOWARNINGS #define HOMOG2D_LOG_WARNING( a ) \ std::cerr << "homog2d warning (" << ++err::warningCount() << "), line " << __LINE__ << "\n msg=" << a << "\n"; #else #define HOMOG2D_LOG_WARNING #endif /* \todo 20230212 ttmath support: this definition does not work, I don't know why !!! see namespace \ref trait \verbatim #define HOMOG2D_CHECK_IS_NUMBER(T) \ static_assert( \ ((std::is_arithmetic::value && !std::is_same::value) || trait::IsBigNumType::value), \ "Type of value must be numerical" ) \endverbatim */ /// Internal type used for numerical computations, possible values: \c double, long double #if !defined(HOMOG2D_INUMTYPE) #define HOMOG2D_INUMTYPE double #endif #if !defined(HOMOG2D_BIND_X) #define HOMOG2D_BIND_X x #endif #if !defined(HOMOG2D_BIND_Y) #define HOMOG2D_BIND_Y y #endif /// Error throw wrapper macro #define HOMOG2D_THROW_ERROR_1( msg ) \ { \ std::ostringstream oss; \ oss << "homog2d: line " << __LINE__ << ", function:" << __FUNCTION__ << "(): " \ << msg << "\n -full function name: " << HOMOG2D_PRETTY_FUNCTION \ << "\n -Error count=" << ++err::errorCount(); \ throw std::runtime_error( oss.str() ); \ } /// Error throw wrapper macro, first arg is the function name #define HOMOG2D_THROW_ERROR_2( func, msg ) \ { \ std::ostringstream oss; \ oss << "homog2d: line " << __LINE__ << ", function:" << func << "(): " \ << msg << "\n -full function name: " << HOMOG2D_PRETTY_FUNCTION \ << "\n -Error count=" << ++err::errorCount(); \ throw std::runtime_error( oss.str() ); \ } /////////////////////////////////////// // Default values for thresholds #ifndef HOMOG2D_THR_ZERO_DIST #define HOMOG2D_THR_ZERO_DIST 1E-10 #endif #ifndef HOMOG2D_THR_ZERO_ORTHO_DIST #define HOMOG2D_THR_ZERO_ORTHO_DIST 1E-14 #endif // default value: 1 thousand of a radian (tan = 0.001 too) #ifndef HOMOG2D_THR_ZERO_ANGLE #define HOMOG2D_THR_ZERO_ANGLE 0.001 #endif #ifndef HOMOG2D_THR_ZERO_DENOM #define HOMOG2D_THR_ZERO_DENOM 1E-10 #endif #ifndef HOMOG2D_THR_ZERO_DETER #define HOMOG2D_THR_ZERO_DETER 1E-15 #endif /////////////////////////////////////// /// Max number of iterations for "Point inside Polygon" algorithm. /// May be adjusted, see manual #ifndef HOMOG2D_MAXITER_PIP #define HOMOG2D_MAXITER_PIP 5 #endif #define HOMOG2D_VERSION "2.11.2" // some MS environments seem to lack Pi definition, even if _USE_MATH_DEFINES is defined #ifndef M_PI #define M_PI 3.14159265358979323846 #endif namespace h2d { /// Holds static counters, for runtime errors and warnings namespace err { /// Used to count the errors /** This is used in the HOMOG2D_THROW_ERROR_1 macros. Some user code could catch the exceptions, thus this will enable the counting of errors */ inline size_t& errorCount() { static size_t c; return c; } /// Used in macro HOMOG2D_LOG_WARNING inline size_t& warningCount() { static size_t c; return c; } } //namespace err /// Holds the types needed for policy based design namespace typ { struct IsLine {}; struct IsPoint {}; struct IsHomogr {}; struct IsEpipmat {}; struct IsClosed {}; struct IsOpen {}; } // namespace typ namespace detail { template class Matrix_; /// Helper class for Root (Point/Line) type, used as a trick to allow partial specialization of member functions template struct BaseHelper {}; /// Helper class for used to get the underlying floating-point type, see Dtype and Common::dtype() template struct DataFpType {}; /// Helper class for PolylineBase, used as a trick to allow partial specialization of member functions template struct PlHelper {}; #ifdef HOMOG2D_FUTURE_STUFF /// Helper class for Matrix type template struct HelperMat {}; template<> struct HelperMat { using M_OtherType = IsEpipmat; }; template<> struct HelperMat { using M_OtherType = IsHomogr; }; #endif template struct HelperPL; template<> struct HelperPL { using OtherType = typ::IsLine; }; template<> struct HelperPL { using OtherType = typ::IsPoint; }; /// A trick used in static_assert, so it aborts only if function is instanciated template struct AlwaysFalse { enum { value = false }; }; } // namespace detail /// Holds base classes, not part of API namespace base { template class PolylineBase; template class LPBase; } template class Hmatrix_; template using Homogr_ = Hmatrix_; #ifdef HOMOG2D_FUTURE_STUFF template using Epipmat_ = Hmatrix_; #endif template class Segment_; template class Circle_; template class FRect_; template class Ellipse_; namespace img { struct SvgImage; // forward declaration } template using Point2d_ = base::LPBase; template using Line2d_ = base::LPBase; template using PointPair1_ = std::pair,Point2d_>; template using PointPair2_ = std::pair,Point2d_>; template using CPolyline_ = base::PolylineBase; template using OPolyline_ = base::PolylineBase; #ifdef HOMOG2D_ENABLE_VRTP /// A variant type, holding all possible types. Used to achieve runtime polymorphism /** See https://github.com/skramm/homog2d/blob/master/docs/homog2d_manual.md#section_rtp */ template using CommonType_ = std::variant< Segment_, Point2d_, Line2d_, Circle_, Ellipse_, FRect_, CPolyline_, OPolyline_ >; #endif // HOMOG2D_ENABLE_VRTP //------------------------------------------------------------------ /// Holds drawing related code, independent of back-end library namespace img { /// Color type , see DrawParams struct Color { uint8_t r = 80; uint8_t g = 80; uint8_t b = 80; Color( uint8_t rr, uint8_t gg, uint8_t bb ): r(rr),g(gg),b(bb) {} Color() = default; friend std::ostream& operator << ( std::ostream& f, const Color& c ) { f << "Color:" << (int)c.r << '-' << (int)c.g << '-' << (int)c.b; return f; } }; /// Helper function, will generate a vector of \c nb random RGB colors /** - RGB values will be between \c minval and \c minval+coeff */ inline std::vector genRandomColors( size_t nb, int minval=20, int maxval=250 ) { if( maxval<=minval ) HOMOG2D_THROW_ERROR_1( "Illegal values for minval and maxval" ); std::vector vcol( nb ); std::srand( std::time(nullptr) ); for( size_t i=0; i friend class h2d::Circle_; template friend class h2d::Segment_; template friend class h2d::FRect_; template friend class h2d::Ellipse_; template friend class h2d::base::PolylineBase; template friend class h2d::base::LPBase; /// Inner struct, holds the values. Needed so we can assign a default value as static member /// \todo 20240329 maybe we can merge parameters _ptDelta and _pointSize into a single one? struct Dp_values { Color _color; int _lineThickness = 1; int _pointSize = 4; int _lineType = 1; /// if OpenCv: 1 for cv::LINE_AA, 2 for cv::LINE_8 uint8_t _ptDelta = 5; ///< pixels, used for drawing points PtStyle _ptStyle = PtStyle::Plus; bool _enhancePoint = false; ///< to draw selected points bool _showPoints = false; ///< show the points (useful only for Segment_ and Polyline_) bool _showIndex = false; ///< show the index as number int _fontSize = 20; ///< font size for drawText() std::string _attrString; ///< added attributes (SVG only) /// Returns the point style following the current one PtStyle nextPointStyle() const { if( _ptStyle == PtStyle::Dot ) return PtStyle::Plus; auto curr = static_cast(_ptStyle); return static_cast(curr+1); } }; friend std::ostream& operator << ( std::ostream& f, const DrawParams& dp ) { f << "-" << dp._dpValues._color << "\n-line width=" << dp._dpValues._lineThickness << "\n-pointSize=" << dp._dpValues._pointSize << "\n-showPoints=" << dp._dpValues._showPoints << "\n-fontSize=" << dp._dpValues._fontSize << '\n'; return f; } public: Dp_values _dpValues; private: static Dp_values& p_getDefault() { static Dp_values s_defValue; return s_defValue; } public: DrawParams() { _dpValues = p_getDefault(); } void setDefault() { p_getDefault() = this->_dpValues; } static void resetDefault() { p_getDefault() = Dp_values(); } DrawParams& setPointStyle( PtStyle ps ) { if( (int)ps > (int)PtStyle::Dot ) throw std::runtime_error( "Error: invalid value for point style"); _dpValues._ptStyle = ps; return *this; } DrawParams& setPointSize( uint8_t ps ) { _dpValues._pointSize = ps; _dpValues._ptDelta = ps; return *this; } DrawParams& setThickness( uint8_t t ) { _dpValues._lineThickness = t; return *this; } DrawParams& setColor( uint8_t r, uint8_t g, uint8_t b ) { _dpValues._color = Color{r,g,b}; return *this; } DrawParams& setColor( Color col ) { _dpValues._color = col; return *this; } DrawParams& selectPoint() { _dpValues._enhancePoint = true; return *this; } /// Set or unset the drawing of points (useful only for Segment_ and Polyline_) DrawParams& showPoints( bool b=true ) { _dpValues._showPoints = b; return *this; } /// Set font size for drawText() DrawParams& setFontSize( int value /* pixels */ ) { assert( value > 1 ); _dpValues._fontSize = value; return *this; } /// Set or unset the drawing of points (useful only for Segment_ and Polyline_) DrawParams& showIndex( bool b=true ) { _dpValues._showIndex = b; return *this; } /// Add some specific SVG attributes (ignored for Opencv renderings) /** \sa getAttrString() */ DrawParams& setAttrString( std::string attr ) { _dpValues._attrString = attr; return *this; } Color color() const { return _dpValues._color; } #ifdef HOMOG2D_USE_OPENCV cv::Scalar cvColor() const { return cv::Scalar( _dpValues._color.b, _dpValues._color.g, _dpValues._color.r ); } #endif // HOMOG2D_USE_OPENCV private: /// Checks if the user-given SVG attribute string (with \ref setAttrString() ) holds the fill="none" mention. /** This is because to have no filling, the object needs to have the fill="none" attribute, so the default behavior is to always add that attribute.
But then if the user wants some filling, then we would have both fill="none" and fill="somecolor", and that would render the svg invalid.
So the drawing code checks if user has added some filling, and if so, does not add the fill="none" attribute. */ bool holdsFill() const { if( !_dpValues._attrString.empty() ) if( _dpValues._attrString.find("fill=") != std::string::npos ) return true; return false; } /// \sa setAttrString() std::string getAttrString() const { if( _dpValues._attrString.empty() ) return std::string(); return _dpValues._attrString + ' '; } std::string getSvgRgbColor() const { std::ostringstream oss; oss << "rgb(" << (int)_dpValues._color.r << ',' << (int)_dpValues._color.g << ',' << (int)_dpValues._color.b << ')'; return oss.str(); } }; // class DrawParams //------------------------------------------------------------------ /// Opaque data structure, will hold the image type, depending on back-end library. /// This type is the one used in all the drawing functions. /** At present the two allowed types are cv::Mat (external Opencv library, requires the symbol HOMOG2D_USE_OPENCV to be defined) or SvgImage (no dependency) */ template class Image { private: T _realImg; size_t _width = 500; size_t _height = 500; public: Image() = default; Image( T& m ): _realImg(m) {} /// Returns a reference on the underlying image T& getReal() { return _realImg; } /// Returns a const reference on the underlying image const T& getReal() const { return _realImg; } std::pair size() const { return std::make_pair( _width, _height ); } /// That constructor is the default, shouln't be instanciated, see specializations Image( size_t, size_t ) { assert(0); // static_assert( detail::AlwaysFalse::value, "no concrete implementation available" ); // static_assert( std::false_type, "no concrete implementation available" ); } void setSize( size_t width, size_t height ); template void setSize( const std::pair& ); void write( std::string ) const // will be specialized { assert(0); } int cols() const { return _width; } int rows() const { return _height; } void clear( Color c=Color(255,255,255) ) { clear(c.r,c.g,c.b); } void clear( uint8_t, uint8_t, uint8_t ) { assert(0); } void clear( uint8_t ) { assert(0); } void drawText( std::string, Point2d_, img::DrawParams dp=img::DrawParams() ); template void draw( const U& object, img::DrawParams dp=img::DrawParams() ); template void draw( const std::pair& p_objects, img::DrawParams dp=img::DrawParams() ); #ifdef HOMOG2D_USE_OPENCV /// Show image on window \c wname (not available for SVG !) void show( std::string wname ) { cv::imshow( wname, _realImg ); } private: void p_setSize( size_t width, size_t height ) { _width = width; _height = height; _realImg.create( (int)height, (int)width, CV_8UC3 ); clear(); } #endif }; // class Image #ifdef HOMOG2D_USE_OPENCV template <> inline Image::Image( size_t width, size_t height ) { p_setSize( width, height ); } template <> inline void Image::setSize( size_t width, size_t height ) { p_setSize( width, height ); } #endif template void Image::setSize( size_t width, size_t height ) { _width = width; _height = height; } template template void Image::setSize( const std::pair& pa ) { setSize( pa.first, pa.second ); } template <> inline Image::Image( size_t width, size_t height ) { setSize( width, height ); } template <> inline void Image::write( std::string fname ) const { std::ofstream file( fname ); if( !file.is_open() ) { HOMOG2D_THROW_ERROR_1( "unable to open output file '" + fname + "'" ); } file << "\n" << "\n"; file << _realImg._svgString.str(); file << "\n"; } template <> inline void Image::clear( uint8_t, uint8_t, uint8_t ) { _realImg._svgString.str(""); _realImg._svgString.clear(); } #ifdef HOMOG2D_USE_OPENCV template <> inline void Image::clear( uint8_t r, uint8_t g, uint8_t b ) { _realImg = cv::Scalar(b,g,r); } template <> inline void Image::clear( uint8_t col ) { _realImg = cv::Scalar(col,col,col); } template <> inline void Image::write( std::string fname ) const { cv::imwrite( fname, _realImg ); } #endif // HOMOG2D_USE_OPENCV template template void Image::draw( const U& object, img::DrawParams dp ) { object.draw( *this, dp ); } template template void Image::draw( const std::pair& pairp, img::DrawParams dp ) { pairp.first.draw( *this, dp ); pairp.second.draw( *this, dp ); } } // namespace img ///////////////////////////////////////////////////////////////////////////// // SECTION - PUBLIC ENUM DECLARATIONS ///////////////////////////////////////////////////////////////////////////// /// Used in base::PolylineBase_::rotate() member function enum class Rotate: int8_t { CCW, ///< Counter ClockWise rotation CW, ///< ClockWise rotation Full, ///< 180° rotation VMirror, ///< vertical symmetry HMirror ///< horizontal symmetry }; enum class CardDir: int8_t { Bottom,Top, Left, Right }; /// Used in Line2d::getValue() and getOrthogonalLine() enum class GivenCoord: uint8_t { X, Y }; /// Used in line constructor, to instanciate a H or V line, see base::LPBase( LineDir, T ) enum class LineDir: uint8_t { H, V }; /// Type of Root object, see rtp::Root::type(). /// Maybe printed out with getString() enum class Type: uint8_t { Line2d, Point2d, Segment, FRect, Circle, Ellipse, OPolyline, CPolyline }; /// Type of underlying floating point, see LPBase::dtype(). /// Maybe printed out with getString() enum class Dtype: uint8_t { Float, Double, LongDouble, Other, // ? #ifdef HOMOG2D_USE_TTMATH Ttmath ///< only if HOMOG2D_USE_TTMATH is defined, see manual #endif }; /// Returns stringified version of \ref type() inline const char* getString( Type t ) { const char* s=0; switch( t ) { case Type::Line2d: s="Line2d"; break; case Type::Point2d: s="Point2d"; break; case Type::Segment: s="Segment"; break; case Type::FRect: s="FRect"; break; case Type::Circle: s="Circle"; break; case Type::Ellipse: s="Ellipse"; break; case Type::OPolyline: s="OPolyline"; break; case Type::CPolyline: s="CPolyline"; break; assert(0); } return s; } #ifdef HOMOG2D_ENABLE_VRTP //------------------------------------------------------------------ /// Holds functors, used to manage runtime polymorphism using \c std::variant namespace fct { //------------------------------------------------------------------ /// A functor to get the type of an object in a std::variant, call with std::visit() /** \sa CommonType_ \sa type() */ struct TypeFunct { template Type operator ()(const T& a) { return a.type(); } }; struct DTypeFunct { template Dtype operator ()(const T& a) { return a.dtype(); } }; //------------------------------------------------------------------ /// A functor to get the length of an object in a std::variant, call with std::visit() /** \sa CommonType_ \sa length() */ struct LengthFunct { template HOMOG2D_INUMTYPE operator ()(const T& a) { return a.length(); } }; //------------------------------------------------------------------ /// A functor to get the area of an object in a std::variant, call with std::visit() /** \sa CommonType_ \sa area() */ struct AreaFunct { template HOMOG2D_INUMTYPE operator ()(const T& a) { return a.area(); } }; //------------------------------------------------------------------ /// A functor used to apply a homography matrix to an object template class TransformFunct { public: TransformFunct( const Homogr_& h ): _h(h) {} template CommonType_ operator ()(const T& a) { return CommonType_{_h * a}; } private: const Homogr_& _h; }; //------------------------------------------------------------------ /// A functor used to draw objects. To use with std::variant and std::visit() template struct DrawFunct { DrawFunct( img::Image& img, img::DrawParams dp=img::DrawParams() ): _img(img), _drawParams(dp) {} img::Image& _img; const img::DrawParams _drawParams; template void operator ()(const T& a) { a.draw( _img, _drawParams ); } }; //------------------------------------------------------------------ /// Convert std::variant object into the underlying type /** source: https://stackoverflow.com/a/72955535/193789 */ template struct VariantUnwrapper { const std::variant& var; template operator T() { return std::get(var); } }; #if __cplusplus < 202002L /// Fix for the above VariantUnwrapper for C++17 /** (may be removed when we switch to C++20) \sa VariantUnwrapper */ template VariantUnwrapper( const std::variant& ) -> VariantUnwrapper; #endif } // namespace fct #endif // HOMOG2D_ENABLE_VRTP inline const char* getString( Dtype t ) { const char* s=0; switch( t ) { case Dtype::Float: s="Float"; break; case Dtype::Double: s="Double"; break; case Dtype::LongDouble: s="LongDouble"; break; case Dtype::Other: s="Other"; break; #ifdef HOMOG2D_USE_TTMATH case Dtype::Ttmath: s="ttmath"; break; #endif assert(0); } return s; } //------------------------------------------------------------------ /// Holds private stuff namespace priv { /// Implementation of dsize(), returns nb of bits of mantissa and exponent (default implementation) template inline std::pair impl_dsize( const detail::DataFpType& ) { return std::make_pair( std::numeric_limits::digits, sizeof(T)*8-std::numeric_limits::digits-1 ); } inline Dtype impl_dtype( const detail::DataFpType& ) { return Dtype::Float; } inline Dtype impl_dtype( const detail::DataFpType& ) { return Dtype::Double; } inline Dtype impl_dtype( const detail::DataFpType& ) { return Dtype::LongDouble; } template inline Dtype impl_dtype( const detail::DataFpType& ) { return Dtype::Other; } #ifdef HOMOG2D_USE_TTMATH /// Implementation for ttmath types template inline Dtype impl_dtype( const detail::DataFpType>& ) { return Dtype::Ttmath; } /// Implementation for ttmath types template inline std::pair impl_dsize( const detail::DataFpType>& ) { return std::make_pair( M*sizeof(size_t)*8, E*sizeof(size_t)*8 ); } #endif } // namespace priv //------------------------------------------------------------------ /// Holds threshold values and api to access these namespace thr { static HOMOG2D_INUMTYPE& nullDistance() { static HOMOG2D_INUMTYPE s_zeroDistance = HOMOG2D_THR_ZERO_DIST; return s_zeroDistance; } static HOMOG2D_INUMTYPE& nullOrthogDistance() { static HOMOG2D_INUMTYPE s_zeroOrthoDistance = HOMOG2D_THR_ZERO_ORTHO_DIST; return s_zeroOrthoDistance; } static HOMOG2D_INUMTYPE& nullAngleValue() { static HOMOG2D_INUMTYPE s_zeroAngleValue = HOMOG2D_THR_ZERO_ANGLE; return s_zeroAngleValue; } static HOMOG2D_INUMTYPE& nullDenom() { static HOMOG2D_INUMTYPE _zeroDenom = HOMOG2D_THR_ZERO_DENOM; return _zeroDenom; } static HOMOG2D_INUMTYPE& nullDeter() { static HOMOG2D_INUMTYPE _zeroDeter = HOMOG2D_THR_ZERO_DETER; return _zeroDeter; } /// This one is used for the Welzl minimum enclosing circle static bool& doNotCheckRadius() { static bool _doNotCheckRadius = false; return _doNotCheckRadius; } /// Helper function, could be needed inline void printThresholds( std::ostream& f ) { f << "homog2d: current threshold values:" << "\n -nullDistance()=" << nullDistance() << "\n -nullOrthogDistance()=" << nullOrthogDistance() << "\n -nullAngleValue()=" << nullAngleValue() << "\n -nullDenom()=" << nullDenom() << "\n -nullDeter()=" << nullDeter() << '\n'; } } // namespace thr // forward declaration namespace base { template auto operator << ( std::ostream&, const h2d::base::LPBase& ) -> std::ostream&; template auto operator << ( std::ostream&, const h2d::base::PolylineBase& ) -> std::ostream&; } // forward declaration, related to https://github.com/skramm/homog2d/issues/2 template Line2d_ operator * ( const Homogr_&, const Line2d_& ); namespace detail { // forward declaration template void product( Matrix_&, const Matrix_&, const Matrix_& ); //------------------------------------------------------------------ template bool shareCommonCoord( const Point2d_& p1, const Point2d_& p2 ) { if( homog2d_abs( p1.getX() - p2.getX() ) < thr::nullOrthogDistance() || homog2d_abs( p1.getY() - p2.getY() ) < thr::nullOrthogDistance() ) return true; return false; } //------------------------------------------------------------------ /// Private free function, get top-left and bottom-right points from two arbitrary points /** Throws if one of the coordinates is equal to the other (x1=x2 or y1=y2)*/ template PointPair1_ getCorrectPoints( const Point2d_& p0, const Point2d_& p1 ) { #ifndef HOMOG2D_NOCHECKS if( shareCommonCoord( p0, p1 ) ) HOMOG2D_THROW_ERROR_1( "a coordinate of the 2 points is identical, does not define a rectangle:\n p0=" << p0 << " p1=" << p1 ); #endif Point2d_ p00( std::min(p0.getX(), p1.getX()), std::min(p0.getY(), p1.getY()) ); Point2d_ p11( std::max(p0.getX(), p1.getX()), std::max(p0.getY(), p1.getY()) ); return std::make_pair( p00, p11 ); } template using matrix_t = std::array,3>; //------------------------------------------------------------------ /// Common templated class for all the geometric primitives template class Common { public: /// Get numerical data type as a Dtype value, can be stringified with h2d::getString(Dtype) /// \sa h2d::dtype( const T& ) Dtype dtype() const { return priv::impl_dtype( detail::DataFpType() ); } /// Get data size expressed as number of bits for, respectively, mantissa and exponent. /// \sa h2d::dsize(const T&) std::pair dsize() const { return priv::impl_dsize( detail::DataFpType() ); } /// This function is a fallback for all sub-classes that do not provide such a method. /** It is necessary in a runtime polymorphism context, as we would have build failures if a given type disallows providing such a method (for example, when trying to check if some object is inside an open polyline, which makes no sense). */ template constexpr bool isInside( const Common& ) const { HOMOG2D_START; return false; } #if 0 /// Fallback for classes not implementing this /// \todo 20240327: this class is inherited by HMatrix... on which the concept of bounding box makes no sense! /// Fix this. template FRect_ getBB() const { HOMOG2D_THROW_ERROR_1( "unable to compute BB for object of type " << getString(this->type) ); return FRect_(); // to avoid a compile warning } #endif size_t size() const { return 0; } }; // class Common } // namespace detail //------------------------------------------------------------------ #ifdef HOMOG2D_ENABLE_PRTP /// Holds pointer-based runtime polymorphism stuff namespace rtp { /// Non-templated root class, to achieve dynamic (runtime) polymorphism /** Only exists if symbol \c HOMOG2D_ENABLE_PRTP is defined, see build options. */ class Root { public: virtual void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const = 0; #ifdef HOMOG2D_USE_OPENCV virtual void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const = 0; #endif virtual HOMOG2D_INUMTYPE length() const = 0; virtual HOMOG2D_INUMTYPE area() const = 0; virtual Type type() const = 0; friend std::ostream& operator << ( std::ostream& f, const Root& p ); virtual ~Root() {} }; } // namespace rtp #endif namespace detail { //------------------------------------------------------------------ /// A simple wrapper over a 3x3 matrix, provides root functionalities /** Homogeneous (thus the 'mutable' attribute). \todo 20240326: we might need to add another level of inheritance. This class inherits \c Common, which is designed to be inherited geometric primitives and as such holds member function that cannot be used on a matrix ! (example: isInside() )
So either we remove the latter function and find a way to put it somewhere else, either we create another intermediate class. */ template class Matrix_: public Common { template friend class Matrix_; template friend void product( base::LPBase&, const detail::Matrix_&, const base::LPBase& ); template friend void product( Matrix_&, const Matrix_&, const Matrix_& ); protected: mutable matrix_t _mdata; mutable bool _isNormalized = false; public: /// Constructor Matrix_() { p_fillZero(); } /// Copy-Constructor template Matrix_( const Matrix_& other ) { for( int i=0; i<3; i++ ) for( int j=0; j<3; j++ ) _mdata[i][j] = other._mdata[i][j]; _isNormalized = other._isNormalized; } matrix_t& getRaw() { return _mdata; } const matrix_t& getRaw() const { return _mdata; } /* template void set( size_t r, size_t c, T v ) { #ifndef HOMOG2D_NOCHECKS HOMOG2D_CHECK_ROW_COL; #endif _mdata[r][c] = v; }*/ const FPT& value( size_t r, size_t c ) const { #ifndef HOMOG2D_NOCHECKS HOMOG2D_CHECK_ROW_COL; #endif return _mdata[r][c]; } FPT& value( size_t r, size_t c ) { #ifndef HOMOG2D_NOCHECKS HOMOG2D_CHECK_ROW_COL; #endif return _mdata[r][c]; } /// Return determinant of matrix /** See https://en.wikipedia.org/wiki/Determinant */ HOMOG2D_INUMTYPE determ() const { auto det = _mdata[0][0] * p_det2x2( {1,1, 1,2, 2,1, 2,2} ); det -= _mdata[0][1] * p_det2x2( {1,0, 1,2, 2,0, 2,2} ); det += _mdata[0][2] * p_det2x2( {1,0, 1,1, 2,0, 2,1} ); return det; } /// Transpose and return matrix Matrix_& transpose() { matrix_t out; for( int i=0; i<3; i++ ) for( int j=0; j<3; j++ ) out[i][j] = _mdata[j][i]; _mdata = out; _isNormalized = false; return *this; } /// Inverse matrix Matrix_& inverse() { auto det = determ(); if( homog2d_abs(det) < thr::nullDeter() ) HOMOG2D_THROW_ERROR_1( "matrix is not invertible, det=" << std::scientific << homog2d_abs(det) ); auto adjugate = p_adjugate(); p_divideAll(adjugate, det); _mdata = adjugate._mdata; _isNormalized = false; return *this; } bool isNormalized() const { return _isNormalized; } protected: void p_normalize( int r, int c ) const { #ifndef HOMOG2D_NOCHECKS if( std::fabs(_mdata[r][c]) < thr::nullDenom() ) HOMOG2D_THROW_ERROR_1( "Unable to normalize matrix, value at (" << r << ',' << c << ") less than " << thr::nullDenom() ); #endif p_divideBy( r, c ); if( std::signbit( _mdata[r][c] ) ) for( auto& li: _mdata ) for( auto& e: li ) e = -e; _isNormalized = true; } /// Divide all elements by the value at (r,c), used for normalization. /** No need to check value, done by caller */ void p_divideBy( size_t r, size_t c ) const { // assert( std::fabs( _mdata[r][c] ) > 1000*std::numeric_limits::epsilon() ); for( auto& li: _mdata ) for( auto& e: li ) e /= _mdata[r][c]; } void p_fillZero() { for( auto& li: _mdata ) for( auto& e: li ) e = 0.; } void p_fillEye() { p_fillZero(); _mdata[0][0] = 1.; _mdata[1][1] = 1.; // "eye" matrix => unit transformation _mdata[2][2] = 1.; } template void p_fillWith( const T& in ) { for( auto i=0; i<3; i++ ) for( auto j=0; j<3; j++ ) _mdata[i][j] = in[i][j]; _isNormalized = false; } /// Divide all elements of \c mat by \c value template void p_divideAll( detail::Matrix_& mat, FPT2 value ) const { for( int i=0; i<3; i++ ) for( int j=0; j<3; j++ ) mat._mdata[i][j] /= value; _isNormalized = false; } template friend Matrix_ operator * ( const Matrix_&, const Matrix_& ); private: HOMOG2D_INUMTYPE p_det2x2( const std::vector& v ) const { auto det = static_cast( _mdata[v[0]][v[1]] ) * _mdata[v[6]][v[7]]; det -= static_cast( _mdata[v[2]][v[3]] ) * _mdata[v[4]][v[5]]; return det; } /// Computes adjugate matrix, see https://en.wikipedia.org/wiki/Adjugate_matrix#3_%C3%97_3_generic_matrix detail::Matrix_ p_adjugate() const { detail::Matrix_ mat_out; matrix_t& out = mat_out._mdata; out[ 0 ][ 0 ] = p_det2x2( {1,1, 1,2, 2,1, 2,2} ); out[ 0 ][ 1 ] = -p_det2x2( {0,1, 0,2, 2,1, 2,2} ); out[ 0 ][ 2 ] = p_det2x2( {0,1, 0,2, 1,1, 1,2} ); out[ 1 ][ 0 ] = -p_det2x2( {1,0, 1,2, 2,0, 2,2} ); out[ 1 ][ 1 ] = p_det2x2( {0,0, 0,2, 2,0, 2,2} ); out[ 1 ][ 2 ] = -p_det2x2( {0,0, 0,2, 1,0, 1,2} ); out[ 2 ][ 0 ] = p_det2x2( {1,0, 1,1, 2,0, 2,1} ); out[ 2 ][ 1 ] = -p_det2x2( {0,0, 0,1, 2,0, 2,1} ); out[ 2 ][ 2 ] = p_det2x2( {0,0, 0,1, 1,0, 1,1} ); return mat_out; } friend std::ostream& operator << ( std::ostream& f, const Matrix_& h ) { for( const auto& li: h._mdata ) { f << "| "; for( const auto& e: li ) f << std::setw(6) << e << ' '; f << " |\n"; } return f; } }; // class Matrix_ template void product( base::LPBase&, const detail::Matrix_&, const base::LPBase& ); template Matrix_ operator * ( const Matrix_& h1, const Matrix_& h2 ) { // HOMOG2D_START; Matrix_ out; product( out, h1, h2 ); return out; } } // namespace detail namespace trait { template struct IsBigNumType : std::false_type {}; #ifdef HOMOG2D_USE_TTMATH template struct IsBigNumType> : std::true_type {}; #endif } // namespace trait //------------------------------------------------------------------ /// A 2D homography, defining a planar transformation /** To define an affine or rigid transformation, you can use: - setRotation() - setTranslation() - setScale() To add an affine or rigid transformation to the current one, you can use: - addRotation() - addTranslation() - addScale() To return to unit transformation, use init() Implemented as a 3x3 matrix Templated by Floating-Point Type (FPT) and by type M (typ::IsEpipmat or typ::IsHomogr) */ template class Hmatrix_ : public detail::Matrix_ { template friend class LPBase; template friend Line2d_ operator * ( const Homogr_&, const Line2d_& ); template friend Point2d_ operator * ( const Homogr_&, const Point2d_& ); template friend base::LPBase::OtherType,V> operator * ( const Hmatrix_& h, const base::LPBase& in ); public: /// \name Constructors ///@{ /// Default constructor, initialize to unit transformation Hmatrix_() { init(); } /// Constructor, set homography to a rotation matrix of angle \c val template explicit Hmatrix_( T val ) { HOMOG2D_CHECK_IS_NUMBER(T); init(); setRotation( val ); } /// Constructor, set homography to a translation matrix ( see Hmatrix_( T ) ) template explicit Hmatrix_( T1 tx, T2 ty ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); init(); setTranslation( tx, ty ); } /// Constructor, used to fill with a vector of vector matrix /** \warning - Input matrix \b must be 3 x 3, but type can be anything that can be copied to \c double - no checking is done on validity of matrix as an homography. Thus some assert can get triggered elsewhere. */ template Hmatrix_( const std::vector>& in ) { #ifndef HOMOG2D_NOCHECKS HOMOG2D_CHECK_IS_NUMBER(T); if( in.size() != 3 ) HOMOG2D_THROW_ERROR_1( "Invalid line size for input: " << in.size() ); for( auto li: in ) if( li.size() != 3 ) HOMOG2D_THROW_ERROR_1( "Invalid column size for input: " << li.size() ); #endif detail::Matrix_::p_fillWith( in ); normalize(); } /// Constructor, used to fill with a std::array template Hmatrix_( const std::array,3>& in ) { HOMOG2D_CHECK_IS_NUMBER(T); detail::Matrix_::p_fillWith( in ); normalize(); } /// Copy-constructor Hmatrix_( const Hmatrix_& other ) : detail::Matrix_( other) , _hasChanged ( true ) , _hmt ( nullptr ) { detail::Matrix_::getRaw() = other.getRaw(); } #ifdef HOMOG2D_USE_OPENCV /// Constructor used to initialise with a cv::Mat, call the assignment operator Hmatrix_( const cv::Mat& mat ) { *this = mat; } #endif ///@} /// Assignment operator Hmatrix_& operator = ( const Hmatrix_& other ) { if( this != &other ) detail::Matrix_::getRaw() = other.getRaw(); _hasChanged = true; return *this; } /// Inverse matrix Hmatrix_& inverse() { detail::Matrix_::inverse(); normalize(); return *this; } #if 0 /// Setter \warning No normalization is done, as this can be done /// several times to store values, we therefore must not normalize in between template void set( size_t r, ///< row size_t c, ///< col T v ///< value ) { #ifndef HOMOG2D_NOCHECKS HOMOG2D_CHECK_ROW_COL; #endif _data[r][c] = v; _isNormalized = false; _hasChanged = true; } /// Getter FPT get( size_t r, size_t c ) const { #ifndef HOMOG2D_NOCHECKS HOMOG2D_CHECK_ROW_COL; #endif return _data[r][c]; } detail::Matrix_& getMat() { return static_cast>(*this); } const detail::Matrix_& getMat() const { return static_cast>(*this); } #endif // 0 void init() { impl_mat_init0( detail::BaseHelper() ); } /// \name Adding/assigning a transformation ///@{ /// Adds a translation \c tx,ty to the matrix template Hmatrix_& addTranslation( T tx, T ty ) { HOMOG2D_CHECK_IS_NUMBER(T); Hmatrix_ out; out.setTranslation( tx, ty ); *this = out * *this; return *this; } /// Sets the matrix as a translation \c tx,ty template Hmatrix_& setTranslation( T1 tx, T2 ty ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); init(); auto& mat = detail::Matrix_::_mdata; mat[0][2] = tx; mat[1][2] = ty; detail::Matrix_::_isNormalized = true; _hasChanged = true; return *this; } /// Adds a rotation with an angle \c theta (radians) to the matrix template Hmatrix_& addRotation( T theta ) { HOMOG2D_CHECK_IS_NUMBER(T); Hmatrix_ out; out.setRotation( theta ); *this = out * *this; return *this; } /// Sets the matrix as a rotation with an angle \c theta (radians) template Hmatrix_& setRotation( T theta ) { HOMOG2D_CHECK_IS_NUMBER(T); auto& mat = detail::Matrix_::_mdata; init(); mat[0][0] = mat[1][1] = std::cos(theta); mat[1][0] = std::sin(theta); mat[0][1] = -mat[1][0]; detail::Matrix_::_isNormalized = true; _hasChanged = true; return *this; } /// Adds the same scale factor to the matrix template Hmatrix_& addScale( T k ) { HOMOG2D_CHECK_IS_NUMBER(T); return this->addScale( k, k ); } /// Adds a scale factor to the matrix template Hmatrix_& addScale( T1 kx, T2 ky ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); Hmatrix_ out; out.setScale( kx, ky ); *this = out * *this; _hasChanged = true; return *this; } /// Sets the matrix as a scaling transformation (same on two axis) template Hmatrix_& setScale( T k ) { HOMOG2D_CHECK_IS_NUMBER(T); return setScale( k, k ); } /// Sets the matrix as a scaling transformation template Hmatrix_& setScale( T1 kx, T2 ky ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); init(); auto& mat = detail::Matrix_::_mdata; mat[0][0] = kx; mat[1][1] = ky; detail::Matrix_::_isNormalized = true; _hasChanged = true; return *this; } ///@} template void applyTo( T& ) const; #ifdef HOMOG2D_USE_OPENCV void copyTo( cv::Mat&, int type=CV_64F ) const; Hmatrix_& operator = ( const cv::Mat& ); #endif /// Homography normalisation void normalize() const { detail::Matrix_::p_normalize(2,2); _hasChanged = true; } void buildFrom4Points( const std::vector>&, const std::vector>&, int method=1 ); /// Matrix multiplication, call the base class product friend Hmatrix_ operator * ( const Hmatrix_& h1, const Hmatrix_& h2 ) { Hmatrix_ out; detail::product( out, static_cast>(h1), static_cast>(h2) ) ; out.normalize(); out._hasChanged = true; return out; } /// Comparison operator. Does normalization if required /** This does an absolute comparison of all matrix elements, one by one, and if one differs more than the threshold, it will return false */ bool operator == ( const Hmatrix_& h ) const { auto& data = detail::Matrix_::_mdata; if( !detail::Matrix_::isNormalized() ) normalize(); if( !h.isNormalized() ) h.normalize(); for( int i=0; i<3; i++ ) for( int j=0; j<3; j++ ) if( std::fabs( static_cast( data[i][j] ) - h.value(i,j) ) >= thr::nullDeter() ) return false; return true; } /// Comparison operator. Does normalization if required bool operator != ( const Hmatrix_& h ) const { return !(*this == h); } ////////////////////////// // PRIVATE FUNCTIONS // ////////////////////////// private: #ifdef HOMOG2D_FUTURE_STUFF /// Implementation for epipolar matrices: initialize to aligned axis void impl_mat_init0( const detail::BaseHelper& ) { _data.fillZero(); _data[2][1] = 1.; _data[1][2] = 1.; _isNormalized = true; } #endif /// Implementation for homographies: initialize to unit transformation void impl_mat_init0( const detail::BaseHelper& ) { detail::Matrix_::p_fillEye(); detail::Matrix_::_isNormalized = true; } ////////////////////////// // DATA SECTION // ////////////////////////// private: mutable bool _hasChanged = true; mutable std::unique_ptr> _hmt; ///< used to store \f$ H^{-1} \f$, but only if required friend std::ostream& operator << ( std::ostream& f, const Hmatrix_& h ) { f << "Hmatrix:\n" << static_cast&>(h); return f; } }; // class Hmatrix_ //------------------------------------------------------------------ /// Holds traits classes namespace trait { /// Traits class, used in generic draw() function /** \todo 20240504 this is only used in the draw() free function but actually not really needed: as that function just calls the member function, we can just let the build fail is no suitable member function is found. */ template struct IsDrawable : std::false_type {}; template struct IsDrawable> : std::true_type {}; template struct IsDrawable> : std::true_type {}; template struct IsDrawable> : std::true_type {}; template struct IsDrawable> : std::true_type {}; template struct IsDrawable> : std::true_type {}; template struct IsDrawable>: std::true_type {}; /// Traits class, used in intersects() for Polyline template struct IsShape : std::false_type {}; template struct IsShape> : std::true_type {}; template struct IsShape> : std::true_type {}; template struct IsShape> : std::true_type {}; template struct IsShape> : std::true_type {}; template struct IsShape>: std::true_type {}; //template struct IsShape>: std::true_type {}; /// Traits class, used to determine if we can use some "isInside()" function template struct HasArea : std::false_type {}; template struct HasArea> : std::true_type {}; template struct HasArea> : std::true_type {}; template struct HasArea> : std::true_type {}; template struct HasArea>: std::true_type {}; /// This one is used in base;;PolylineBase::isInside() template struct PolIsClosed : std::false_type {}; template struct PolIsClosed>: std::true_type {}; /// Traits class used in operator * ( const Hmatrix_& h, const Cont& vin ), /// used to detect if container is valid template struct IsContainer : std::false_type { }; template struct IsContainer> : std::true_type { }; template struct IsContainer> : std::true_type { }; template struct IsContainer> : std::true_type { }; template struct IsArray : std::false_type { }; template struct IsArray> : std::true_type { }; /// Traits class used to detect if container \c T is a \c std::array /** (because allocation is different, see \ref alloc() ) */ template struct Is_std_array : std::false_type {}; template struct Is_std_array> : std::true_type {}; /// Traits class, used for getBB() set of functions template struct IsSegment : std::false_type {}; template struct IsSegment> : std::true_type {}; template struct IsPoint : std::false_type {}; template struct IsPoint> : std::true_type {}; /// Traits class, used to check if type has a Bounding Box /** Difference with \c trait::HasArea is that this one is true for OPolyline, whereas the other is not */ template struct HasBB : std::false_type {}; template struct HasBB> : std::true_type {}; template struct HasBB> : std::true_type {}; template struct HasBB> : std::true_type {}; template struct HasBB>: std::true_type {}; #ifdef HOMOG2D_ENABLE_VRTP template struct IsVariant : std::false_type {}; template struct IsVariant>: std::true_type {}; #endif } // namespace trait //------------------------------------------------------------------ namespace detail { /// Holds 9 parameters of Ellipse_ template struct EllParams { template friend struct EllParams; EllParams() = default; /// Copy-constructor is needed when converting from one type to another type template EllParams( const EllParams& p ) { x0 = p.x0; y0 = p.y0; theta= p.theta; sint = p.sint; cost = p.cost; a = p.a; b = p.b; a2 = p.a2; b2 = p.b2; } T x0, y0; ///< center T theta = 0.; ///< angle T sint, cost; ///< \f$ \sin( \theta), \cos( \theta) \f$ T a, b; T a2, b2; ///< squared values of a and b template friend std::ostream& operator << ( std::ostream& f, const EllParams& par ); }; template std::ostream& operator << ( std::ostream& f, const EllParams& par ) { f << "EllParams: origin=" << par.x0 << "," << par.y0 << " angle=" << par.theta *180./M_PI << " a=" << par.a << " b=" << par.b << ' '; return f; } } // namespace detail //------------------------------------------------------------------ /// Ellipse as a conic in matrix form. /** This enables its projection using homography See: - https://en.wikipedia.org/wiki/Ellipse#General_ellipse - https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections General equation of an ellipse: \f[ A x^2 + B x y + C y^2 + D x + E y + F = 0 \f] It can be written as a 3 x 3 matrix: \f[ \begin{bmatrix} A & B/2 & D/2 \\ B/2 & C & E/2 \\ D/2 & E/2 & F \end{bmatrix} \f] Matrix coefficients computed from center x0,y0, major and minor distances (a,b) and angle theta: \f[ \begin{aligned} A &= a^2 \sin^2\theta + b^2 \cos^2\theta \\ B &= 2\left(b^2 - a^2\right) \sin\theta \cos\theta \\ C &= a^2 \cos^2\theta + b^2 \sin^2\theta \\ D &= -2A x_\circ - B y_\circ \\ E &= - B x_\circ - 2C y_\circ \\ F &= A x_\circ^2 + B x_\circ y_\circ + C y_\circ^2 - a^2 b^2 \end{aligned} \f] Homography projection: https://math.stackexchange.com/a/2320082/133647 \f[ Q' = H^{-T} \cdot Q \cdot H^{-1} \f] */ template class Ellipse_: public detail::Matrix_ #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; using detail::Common::isInside; Type type() const { return Type::Ellipse; } template friend class Ellipse_; template friend Ellipse_ operator * ( const Homogr_&, const Circle_& ); template friend Ellipse_ operator * ( const Homogr_&, const Ellipse_& ); public: /// \name Constructors ///@{ /// Default constructor: centered at (0,0), major=2, minor=1 Ellipse_(): Ellipse_( 0., 0., 2., 1., 0. ) {} /// Constructor 1 template Ellipse_( const Point2d_& pt, T2 major=2., T2 minor=1., T3 angle=0. ) : Ellipse_( pt.getX(), pt.getY(), major, minor, angle ) {} /// Constructor 2 template explicit Ellipse_( T1 x, T1 y, T2 major=2., T2 minor=1., T3 angle=0. ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); HOMOG2D_CHECK_IS_NUMBER(T3); if( major& cir ) { p_init( cir.center().getX(), cir.center().getY(), cir.radius(), cir.radius(), 0. ); } /// Copy-Constructor template Ellipse_( const Ellipse_& other ) : detail::Matrix_( other ) {} ///@} /// Translate Ellipse template void translate( TX dx, TY dy ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); auto par = p_getParams(); p_init( par.x0+dx, par.y0+dy, par.a, par.b, par.theta ); } /// Translate Ellipse template void translate( const std::pair& ppt ) { this->translate( ppt.first, ppt.second ); } /// Move Ellipse to new location template void moveTo( TX x, TY y ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); auto par = p_getParams(); p_init( x, y, par.a, par.b, par.theta ); } /// Move Ellipse to new location, given by \c new_org template void moveTo( const Point2d_& new_org ) { this->moveTo( new_org.getX(), new_org.getY() ); } /// \name attributes of ellipse ///@{ bool isCircle( HOMOG2D_INUMTYPE thres=1.E-10 ) const; Point2d_ getCenter() const; CPolyline_ getOBB() const; FRect_ getBB() const; HOMOG2D_INUMTYPE angle() const; std::pair getMajMin() const; ///@} /// Area of ellipse HOMOG2D_INUMTYPE area() const { auto par = p_getParams(); return M_PI * par.a * par.b; } // Circumference HOMOG2D_INUMTYPE length() const; std::pair,Line2d_> getAxisLines() const; template bool pointIsInside( const Point2d_& ) const; ////////////////////////// // OPERATORS // ////////////////////////// template friend std::ostream& operator << ( std::ostream& f, const Ellipse_& ell ); /// Comparison operator. Does normalization if required bool operator == ( const Ellipse_& h ) const { auto& data = detail::Matrix_::_mdata; if( !detail::Matrix_::isNormalized() ) detail::Matrix_::p_normalize(2,2); if( !h.isNormalized() ) h.p_normalize(2,2); for( int i=0; i<3; i++ ) for( int j=0; j<3; j++ ) if( std::fabs( static_cast( data[i][j] ) - h.value(i,j) ) >= thr::nullDeter() ) return false; return true; } /// Comparison operator. Does normalization if required bool operator != ( const Ellipse_& e ) const { return !(*this == e); } ////////////////////////// // PRIVATE FUNCTIONS // ////////////////////////// private: /// private constructor, needed in friend function only /** This is not public, because: 1-useless, 2-no guarantee would be given that the input is indeed a valid ellipse */ explicit Ellipse_( const detail::Matrix_& mat ): detail::Matrix_( mat ) {} template detail::EllParams p_getParams() const; template detail::EllParams p_computeParams() const; /// Called by all the constructors, fills the matrix. void p_init( double x0, double y0, double a, double b, double theta=0. ) { auto& data = detail::Matrix_::_mdata; HOMOG2D_INUMTYPE sin1 = std::sin(theta); HOMOG2D_INUMTYPE cos1 = std::cos(theta); HOMOG2D_INUMTYPE sin2 = sin1 * sin1; HOMOG2D_INUMTYPE cos2 = cos1 * cos1; HOMOG2D_INUMTYPE a2 = a*a; HOMOG2D_INUMTYPE b2 = b*b; HOMOG2D_INUMTYPE A = a2*sin2 + b2*cos2; HOMOG2D_INUMTYPE B = (HOMOG2D_INUMTYPE)2.*(b2-a2) * std::sin(theta) * std::cos(theta); HOMOG2D_INUMTYPE C = a2 * cos2 + b2 * sin2; HOMOG2D_INUMTYPE D = (HOMOG2D_INUMTYPE)(-2.)*A * x0 - B * y0; HOMOG2D_INUMTYPE E = - B * x0 - (HOMOG2D_INUMTYPE)2.*C * y0; HOMOG2D_INUMTYPE F = A*x0*x0 + B*x0*y0 + C*y0*y0 - a2*b2; data[0][0] = A; data[1][1] = C; data[2][2] = F; data[0][1] = data[1][0] = B / 2.; data[0][2] = data[2][0] = D / 2.; data[1][2] = data[2][1] = E / 2.; #ifdef HOMOG2D_OPTIMIZE_SPEED _epHasChanged = false; _par.a = a; _par.b = b; _par.a2 = a2; _par.b2 = b2; _par.theta = theta; _par.sint = sin1; _par.cost = cos1; _par.x0 = x0; _par.y0 = y0; #endif } public: #ifdef HOMOG2D_USE_OPENCV void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; #endif void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; ////////////////////////// // DATA SECTION // ////////////////////////// private: // (matrix holding the data is inherited from base class) #ifdef HOMOG2D_OPTIMIZE_SPEED mutable bool _epHasChanged = true; ///< if true, means we need to recompute parameters mutable detail::EllParams _par; #endif }; // class Ellipse //------------------------------------------------------------------ namespace detail { // forward declaration of template instanciation template base::LPBase crossProduct( const base::LPBase&, const base::LPBase& ); class Inters_1 {}; class Inters_2 {}; /// Common stuff for intersection code struct IntersectCommon { protected: bool _doesIntersect = false; public: bool operator()() const { return _doesIntersect; } }; /// Base class for intersection, gets specialized template class Intersect {}; //------------------------------------------------------------------ /// One point intersection template class Intersect: public IntersectCommon { template friend class ::h2d::Segment_; public: Point2d_ get() const { if( !_doesIntersect ) HOMOG2D_THROW_ERROR_1( "No intersection points" ); return _ptIntersect; } void set( const Point2d_& pt ) { _ptIntersect = pt; _doesIntersect = true; } size_t size() const { return _doesIntersect?1:0; } Intersect() {} Intersect( const Point2d_& ptInter ) : _ptIntersect(ptInter) { _doesIntersect = true; } /// To enable conversions from different floating-point types template Intersect( const Intersect& other ) : IntersectCommon( other ) { _ptIntersect = other._ptIntersect; } private: Point2d_ _ptIntersect; }; //------------------------------------------------------------------ /// Two points intersection template class Intersect: public IntersectCommon { template friend class ::h2d::base::LPBase; public: Intersect() {} Intersect( const Point2d_& p1, const Point2d_& p2 ) : _ptIntersect_1(p1), _ptIntersect_2(p2) { _doesIntersect = true; } /// To enable conversions from different floating-point types template Intersect( const Intersect& other ) : IntersectCommon( other ) { auto ppts = other.get(); _ptIntersect_1 = ppts.first; _ptIntersect_2 = ppts.second; } size_t size() const { return _doesIntersect?2:0; } PointPair1_ get() const { if( !_doesIntersect ) HOMOG2D_THROW_ERROR_1( "No intersection points" ); return std::make_pair( _ptIntersect_1, _ptIntersect_2 ); } private: Point2d_ _ptIntersect_1, _ptIntersect_2; friend std::ostream& operator << ( std::ostream& f, const Intersect& inters ) { f << "bool=" << inters._doesIntersect << " p1:" << inters._ptIntersect_1 << " p2:" << inters._ptIntersect_2; return f; } }; //------------------------------------------------------------------ /// Multiple points intersections template class IntersectM { private: mutable std::vector> _vecInters; ///< mutable, because it can get sorted in const functions public: IntersectM() {} /// To enable conversions from different floating-point types template IntersectM( const IntersectM& other ) { _vecInters.resize( other.size() ); auto it = _vecInters.begin(); for( const auto& elem: other.get() ) *it++ = elem; // automatic type conversion } bool operator()() const { return !_vecInters.empty(); } size_t size() const { return _vecInters.size(); } void add( const Point2d_& pt ) { _vecInters.push_back(pt); } void add( const std::vector>& vpt ) { for( const auto& pt: vpt ) _vecInters.push_back(pt); } /// Returns the intersection points, sorted std::vector> get() const { std::sort( std::begin(_vecInters), std::end(_vecInters) ); return _vecInters; } friend std::ostream& operator << ( std::ostream& f, const IntersectM& i ) { f << "IntersectM: size=" << i.size() << '\n' << i._vecInters; return f; } }; //------------------------------------------------------------------ /// Helper class, holds result of intersections of two FRect_ /// \sa FRect_::intersectArea() /** FPT: Floating Point Type, is defined by the rectangle on which member function \c intersectArea() is called. */ template class RectArea { private: bool _success = false; FRect_ _area; public: RectArea() = default; RectArea( const FRect_& r ) : _success(true), _area(r) {} bool operator()() const { return _success; } FRect_ get() const { if( !_success ) HOMOG2D_THROW_ERROR_1( "unable, no intersection between the two rectangles" ); return _area; } }; } // namespace detail //------------------------------------------------------------------ /// A Flat Rectangle, modeled by its two opposite points template class FRect_: public detail::Common #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; using detail::Common::isInside; Type type() const { return Type::FRect; } template friend class FRect_; private: Point2d_ _ptR1,_ptR2; public: /** \name Constructors */ ///@{ /// Default constructor, initialize rectangle to (0,0)-(1,1) FRect_() { _ptR2.set( 1., 1. ); } /// Constructor from 2 points template FRect_( const Point2d_& pa, const Point2d_& pb ) { set( pa, pb ); } /// Constructor from pair of points template FRect_( const PointPair1_& ppts ) { set( ppts.first, ppts.second ); } /// Constructor from center point, width and height template FRect_( const Point2d_& p0, T1 w, T2 h ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); set( Point2d_( p0.getX()-0.5L*w, p0.getY()-0.5L*h ), Point2d_( p0.getX()+0.5L*w, p0.getY()+0.5L*h ) ); } /// Constructor from x1, y1, x2, y2 (need to be all the same type) template FRect_( T x1, T y1, T x2, T y2 ) { HOMOG2D_CHECK_IS_NUMBER(T); set( Point2d_(x1,y1), Point2d_(x2,y2) ); } /// Copy-Constructor template FRect_( const FRect_& other ) : _ptR1(other._ptR1),_ptR2(other._ptR2) {} ///@} private: /// Private constructor from 4 points, used in intersectArea( const FRect_& ) template FRect_( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& pt3, const Point2d_& pt4 ) { // HOMOG2D_LOG( "pt1" << pt1 << " pt2="<< pt2 << " pt3=" << pt3 << " pt4=" << pt4 ); auto x0 = std::min( pt1.getX(), pt2.getX() ); x0 = std::min( x0, pt3.getX() ); x0 = std::min( x0, pt4.getX() ); auto y0 = std::min( pt1.getY(), pt2.getY() ); y0 = std::min( y0, pt3.getY() ); y0 = std::min( y0, pt4.getY() ); auto x1 = std::max( pt1.getX(), pt2.getX() ); x1 = std::max( x1, pt3.getX() ); x1 = std::max( x1, pt4.getX() ); auto y1 = std::max( pt1.getY(), pt2.getY() ); y1 = std::max( y1, pt3.getY() ); y1 = std::max( y1, pt4.getY() ); _ptR1 = Point2d_(x0,y0); _ptR2 = Point2d_(x1,y1); // HOMOG2D_LOG( "ptR1=" <<_ptR1 << " _ptR2=" << _ptR2 ); } public: /// Assigns points \c pa and \c pb to rectangle template void set( const Point2d_& pa, const Point2d_& pb ) { auto ppts = detail::getCorrectPoints( pa, pb ); _ptR1 = ppts.first; _ptR2 = ppts.second; } /// Assigns points (x1,y1) and (x2,y2) to rectangle template void set( T x1, T y1, T x2, T y2 ) { set( Point2d_(x1,y1), Point2d_(x2,y2) ); } /// \name Attributes access ///@{ HOMOG2D_INUMTYPE height() const { return _ptR2.getY() - _ptR1.getY(); } HOMOG2D_INUMTYPE width() const { return _ptR2.getX() - _ptR1.getX(); } HOMOG2D_INUMTYPE area() const { return height() * width(); } HOMOG2D_INUMTYPE length() const { return 2.*height() + 2.*width(); } /// Return size of rectangle in a std::pair std::pair size() const { return std::make_pair( width(), height() ); } /// Needed for getBB( pair of objects ) FRect_ getBB() const { return *this; } /// Returns the 2 major points of the rectangle /// \sa getPts( const FRect_& ) PointPair1_ getPts() const { return std::make_pair( _ptR1, _ptR2 ); } /// Returns center of rectangle Point2d_ getCenter() const { return Point2d_( (static_cast(_ptR1.getX() ) + _ptR2.getX() ) * 0.5, (static_cast(_ptR1.getY() ) + _ptR2.getY() ) * 0.5 ); } /// Return circle passing through 4 points of flat rectangle /// \sa h2d::getBoundingCircle() Circle_ getBoundingCircle() const { auto middle_pt = getCenter(); return Circle_( middle_pt, middle_pt.distTo( _ptR1 ) ); } /// Return circle inscribed in rectangle /// \sa h2d::getInscribedCircle() Circle_ getInscribedCircle() const { auto segs = getSegs(); auto center = getCenter(); Circle_ cir( center ); cir.radius() = std::min( center.distTo( segs[0] ), center.distTo( segs[1] ) ); return cir; } ///@} /// \name Modifying functions ///@{ /// Translate FRect template void translate( TX dx, TY dy ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); _ptR1.set( _ptR1.getX() + dx, _ptR1.getY() + dy ); _ptR2.set( _ptR2.getX() + dx, _ptR2.getY() + dy ); } /// Translate FRect template void translate( const std::pair& pa ) { this->translate( pa.first, pa.second ); } /// Move FRect to other location template void moveTo( TX x, TY y ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); auto s = size(); _ptR1.set(x,y); _ptR2.set( _ptR1.getX() + s.first, _ptR1.getY() + s.second ); } /// Move FRect to other location, given by \c pt template void moveTo( const Point2d_& pt ) { auto s = size(); _ptR1 = pt; _ptR2.set( _ptR1.getX() + s.first, _ptR1.getY() + s.second ); } template void rotate( Rotate, const Point2d_& ); void rotate( Rotate ); ///@} FRect_ getExtended() const { auto x1 = _ptR1.getX() * 2. - _ptR2.getX(); auto y1 = _ptR1.getY() * 2. - _ptR2.getY(); auto x2 = _ptR2.getX() * 2. - _ptR1.getX(); auto y2 = _ptR2.getY() * 2. - _ptR1.getY(); return FRect_( x1, y1, x2, y2 ); } std::pair,Segment_> getDiagonals() const { auto x1 = _ptR1.getX(); auto y1 = _ptR1.getY(); auto x2 = _ptR2.getX(); auto y2 = _ptR2.getY(); return std::make_pair( Segment_(x1,y1,x2,y2), Segment_(x1,y2,x2,y1) ); } /// \name Union/intersection area ///@{ template CPolyline_ unionArea( const FRect_& other ) const; template detail::RectArea intersectArea( const FRect_& other ) const; template CPolyline_ operator | ( const FRect_& other ) const { return this->unionArea( other ); } template detail::RectArea operator & ( const FRect_& other ) const { return this->intersectArea( other ); } ///@} /// Returns the 4 points of the rectangle, starting from "smallest" one, and /// in clockwise order /** \verbatim p1 +------+ p2 | | | | | | p0 +------+ p3 \endverbatim \sa get4Pts( const FRect_& ) */ std::array,4> get4Pts() const { std::array,4> arr; arr[0] = _ptR1; arr[1] = Point2d_( _ptR1.getX(), _ptR2.getY() ); arr[2] = _ptR2; arr[3] = Point2d_( _ptR2.getX(), _ptR1.getY() ); return arr; } /// Returns the 4 segments of the rectangle, starting with the first vertical one /** \verbatim s1 +------+p2 | | s0 | | s2 | | p1+------+ s3 \endverbatim \sa \ref h2d::getSegs( const FRect_& ) */ std::array,4> getSegs() const { auto pts = get4Pts(); std::array,4> out; out[0] = Segment_( pts[0], pts[1] ); out[1] = Segment_( pts[1], pts[2] ); out[2] = Segment_( pts[2], pts[3] ); out[3] = Segment_( pts[3], pts[0] ); return out; } private: /// Returns as a pair the two main segments of the rectangle, with the largest one in "first" /** - If one of the 3 considered points is at infinity, order is not checked */ std::pair,Segment_> p_getPairSegs() const { auto pts = get4Pts(); std::pair,Segment_> out; out.first = Segment_( pts[0], pts[1] ); out.second = Segment_( pts[1], pts[2] ); if( !pts[0].isInf() && !pts[1].isInf() && !pts[2].isInf() ) if( out.first.length() < out.second.length() ) std::swap( out.first, out.second ); return out; } public: /// \name Enclosing functions ///@{ /// Returns true if rectangle is inside \c shape (Circle_ or FRect_ or base::Polyline) /// \todo add some SFINAE to enable only for allowed types? template bool isInside( const T& shape ) const { HOMOG2D_START; for( const auto& pt: get4Pts() ) if( !pt.isInside( shape ) ) return false; return true; } /// A FRect is never inside an open polyline template constexpr bool isInside( const OPolyline_& ) const { return false; } /// For a rectangle to be inside a closed Polyline, two conditions are necessary: /** - all the points must be inside - no intersections */ template bool isInside( const CPolyline_& poly ) const { for( const auto& seg: getSegs() ) if( seg.intersects(poly)() ) return false; for( const auto& pt: get4Pts() ) if( !pt.isInside( poly ) ) return false; return true; } ///@} /// \name Intersection functions ///@{ /// FRect/Line intersection template detail::IntersectM intersects( const Line2d_& line ) const { return line.intersects( *this ); } /// FRect/Segment intersection template detail::IntersectM intersects( const Segment_& seg ) const { detail::IntersectM out; for( const auto& rseg: getSegs() ) { auto inters = rseg.intersects( seg ); // call of Segment/Segment if( inters() ) { auto pt = inters.get(); bool addPoint = true; if( out.size() == 1 ) // if we have already one if( out.get()[0] == pt ) addPoint = false; if( addPoint ) out.add( pt ); if( out.size() == 2 ) break; } } return out; } /// FRect/Circle intersection template detail::IntersectM intersects( const Circle_& circle ) const { // HOMOG2D_START; return p_intersects_R_C( circle ); } /// FRect/Polyline intersection template detail::IntersectM intersects( const base::PolylineBase& pl ) const { // HOMOG2D_START; return pl.intersects( *this ); } /// FRect/FRect intersection template detail::IntersectM intersects( const FRect_& rect ) const { // HOMOG2D_START; if( *this == rect ) // if rectangles are the same, return detail::IntersectM(); // no intersection return p_intersects_R_C( rect ); } ///@} /// \name Operators ///@{ template bool operator == ( const FRect_& other ) const { if( _ptR1 != other._ptR1 ) return false; if( _ptR2 != other._ptR2 ) return false; return true; } template bool operator != ( const FRect_& other ) const { return !( *this == other ); } ///@} template friend std::ostream& operator << ( std::ostream& f, const FRect_& r ); private: template std::vector> p_pointsInside( const FRect_& other ) const { std::vector> out; for( const auto& pt: get4Pts() ) if( pt.isInside( other ) ) out.push_back( pt ); return out; } /// Intersection of FRect vs FRect or Circle /** - We use a \c std::set to avoid having multiple times the same point. - second arg is used to fetch the indexes of intersecting segments, when needed */ template detail::IntersectM p_intersects_R_C( const T& other ) const { std::set> pts; for( const auto& rseg: getSegs() ) { auto inters = rseg.intersects( other ); // call of Segment/FRect => FRect/Segment, or Segment/Circle if( inters() ) { auto vpts = inters.get(); assert( vpts.size() < 3 ); if( vpts.size() > 0 ) pts.insert( vpts[0] ); if( vpts.size() > 1 ) pts.insert( vpts[1] ); } } detail::IntersectM out; // convert std::set to output container for( const auto& elem: pts ) out.add( elem ); return out; } public: #ifdef HOMOG2D_USE_OPENCV void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; #endif void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; }; // class FRect_ //------------------------------------------------------------------ /// A circle template class Circle_: public detail::Common #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; using detail::Common::isInside; Type type() const { return Type::Circle; } template friend class Circle_; private: FPT _radius; Point2d_ _center; public: /// \name Constructors ///@{ /// Default constructor, unit-radius circle at (0,0) Circle_() : _radius(1.) {} /// 1-arg constructor 1, given radius circle at (0,0) template< typename T, typename std::enable_if< std::is_arithmetic::value ,T >::type* = nullptr > explicit Circle_( T rad ) : Circle_( Point2d_(), rad ) {} /// 1-arg constructor 2, given center point, radius = 1.0 template explicit Circle_( Point2d_ center ) : Circle_( center, 1. ) {} /// 1-arg constructor 3, build circle from a set of 2, 3, or more points /// (Minimum Enclosing Circle, aka MEC) template< typename T, typename std::enable_if< trait::IsContainer::value ,T >::type* = nullptr > explicit Circle_( const T& pts ) { set( pts ); } /// 2-arg constructor 1: point and radius template Circle_( const Point2d_& center, T2 rad=1.0 ) : _radius(rad), _center(center) { #ifndef HOMOG2D_NOCHECKS if( homog2d_abs(rad) < thr::nullDistance() && !h2d::thr::doNotCheckRadius() ) HOMOG2D_THROW_ERROR_1( "radius value too small: " << std::scientific << homog2d_abs(rad) ); if( rad < 0. ) HOMOG2D_THROW_ERROR_1( "radius must not be <0" ); #endif } /// 2-arg constructor 2: circle from 2 points (may be of different types) template Circle_( const Point2d_& pt1, const Point2d_& pt2 ) { set( pt1, pt2 ); } /// 3-arg constructor 1: build circle from 3 floating-point values: x, y, radius /** We need Sfinae because there is another 3-args constructor (circle from 3 points) */ template< typename T1, typename T2, typename std::enable_if< (std::is_arithmetic::value && !std::is_same::value) ,T1 >::type* = nullptr > Circle_( T1 x, T1 y, T2 rad ) : Circle_( Point2d_(x,y), rad ) { // HOMOG2D_CHECK_IS_NUMBER(T1); // not needed, done by sfinae above HOMOG2D_CHECK_IS_NUMBER(T2); } /// 3-arg constructor 2: builds a circle from 3 points /** We need Sfinae because there is another 3-args constructor (x, y, radius as floating point values) */ template< typename PT, typename std::enable_if< !std::is_arithmetic::value ,PT >::type* = nullptr > Circle_( const PT& pt1, const PT& pt2, const PT& pt3 ) { set( pt1, pt2, pt3 ); } /// Copy-Constructor template Circle_( const Circle_& other ) : _radius(other._radius), _center(other._center) {} ///@} /// \name Attributes access ///@{ FPT& radius() { return _radius; } const FPT& radius() const { return _radius; } Point2d_& center() { return _center; } const Point2d_& center() const { return _center; } const Point2d_& getCenter() const { return _center; } HOMOG2D_INUMTYPE area() const { return static_cast(_radius) * _radius * M_PI; } HOMOG2D_INUMTYPE length() const { return static_cast(_radius) * M_PI * 2.0; } /// Returns Bounding Box FRect_ getBB() const { return FRect_( static_cast( _center.getX() ) - _radius, static_cast( _center.getY() ) - _radius, static_cast( _center.getX() ) + _radius, static_cast( _center.getY() ) + _radius ); } ///@} /// \name Edit values ///@{ /// Set circle center point, radius unchanged template void set( const Point2d_& center ) { _center = center; } /// Set circle radius, center point unchanged /** Use of Sfinae so it can be selected only for arithmetic types */ template< typename T, typename std::enable_if< (std::is_arithmetic::value && !std::is_same::value) ,T >::type* = nullptr > void set( T rad ) { _radius = rad; } /// Set circle from center point and radius template void set( const Point2d_& center, FPT3 rad ) { Circle_ c( center, rad ); std::swap( c, *this ); /// \todo 20211216: replace with move } /// Set circle from 3 values (x0,y0,radius) template< typename FPT2, typename std::enable_if< std::is_arithmetic::value ,FPT2 >::type* = nullptr > void set( FPT2 x, FPT2 y, FPT2 rad ) { set( Point2d_(x,y), rad ); } /// Set circle from 2 points template void set( const Point2d_& pt1, const Point2d_& pt2 ); // Set circle from 3 points template void set( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& pt3 ); // set Minimum Enclosing Circle (MEC) from a set of points template< typename T, typename std::enable_if< trait::IsContainer::value ,T >::type* = nullptr > void set( const T& ); private: template // helper function Circle_ p_WelzlHelper( std::vector&, std::vector, size_t ) const; template // helper function Circle_ p_min_circle_trivial( const std::vector& P ) const; public: /// Translate Circle template void translate( TX dx, TY dy ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); _center.translate( dx, dy ); } /// Translate Circle template void translate( const std::pair& pa ) { this->translate( pa.first, pa.second ); } /// Move Circle to other location template void moveTo( TX x, TY y ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); set( Point2d_(x, y) ); } /// Move Circle to other location, geiven by \c pt template void moveTo( const Point2d_& pt ) { set( pt ); } ///@} /// \name Enclosing functions ///@{ /// Returns true if circle is inside \c other circle template bool isInside( const Circle_& other ) const { return( _radius + _center.distTo( other.center() ) < other.radius() ); } /// Returns true if circle is inside rectangle defined by \c p1 and \c p2 template bool isInside( const Point2d_& p1, const Point2d_& p2 ) const { return implC_isInside( detail::getCorrectPoints( p1, p2 ) ); } /// Returns true if circle is inside flat rectangle \c rect template bool isInside( const FRect_& rect ) const { return implC_isInside( rect.getPts() ); } // Returns true if circle is inside close polyline \c poly template bool isInside( const base::PolylineBase& poly ) const; ///@} /// \name Intersection functions ///@{ /// Circle/Line intersection template detail::Intersect intersects( const Line2d_& li ) const { return li.intersects( *this ); } /// Circle/Segment intersection template detail::IntersectM intersects( const Segment_& seg ) const { return seg.intersects( *this ); } // Circle/Circle intersection template detail::IntersectM intersects( const Circle_& ) const; /// Circle/FRect intersection template detail::IntersectM intersects( const FRect_& rect ) const { return rect.intersects( * this ); } /// Circle/Polyline intersection template detail::IntersectM intersects( const base::PolylineBase& pl ) const { return pl.intersects( * this ); } ///@} private: template bool implC_isInside( const PointPair1_& ppts ) const { const auto& p1 = ppts.first; const auto& p2 = ppts.second; HOMOG2D_INUMTYPE rad = _radius; // convert to highest precision if( _center.getX() + rad < p2.getX() ) if( _center.getX() - rad > p1.getX() ) if( _center.getY() + rad < p2.getY() ) if( _center.getY() - rad > p1.getY() ) return true; return false; } public: /// \name Operators ///@{ template bool operator == ( const Circle_& other ) const { if( _radius != other._radius ) return false; if( _center != other._center ) return false; return true; } template bool operator != ( const Circle_& other ) const { return !( *this == other ); } ///@} template friend std::ostream& operator << ( std::ostream& f, const Circle_& r ); #ifdef HOMOG2D_USE_OPENCV void draw( img::Image&, img::DrawParams=img::DrawParams() ) const; #endif void draw( img::Image&, img::DrawParams=img::DrawParams() ) const; }; // class Circle_ //------------------------------------------------------------------ /// Holds private stuff namespace priv { /// Private free function, swap the points so that \c ptA.x <= \c ptB.x, and if equal, sorts on y template void fix_order( Point2d_& ptA, Point2d_& ptB ) { if( !(ptA < ptB) ) std::swap( ptA, ptB ); } //------------------------------------------------------------------ /// Free function, squared distance between points (sqrt not needed for comparisons, and can save some time) /// \sa Point2d_::distTo() /// \sa dist( const Point2d_&, const Point2d_& ) template HOMOG2D_INUMTYPE sqDist( const Point2d_& pt1, const Point2d_& pt2 ) { auto dx = (HOMOG2D_INUMTYPE)pt1.getX() - pt2.getX(); auto dy = (HOMOG2D_INUMTYPE)pt1.getY() - pt2.getY(); return dx*dx + dy*dy; } //------------------------------------------------------------------ /// Helper function, used to check for colinearity of three points /** This will return the same points as given in input but ordered as: - the pair that has the largest distance in [0] and [1] - the third point in [2] - the point closest to the third point in [1], the farthest in [0] \sa bool areCollinear() \todo 20220520: needs some optimization, once it has been extensively tested We have theses 6 situations described on below diagrams A through F, with the desired output order: \verbatim 1 + 1 + | A => 3,1,2 | B => 2,1,3 | | 2 +---------------+ 3 3 +---------------+ 2 2 + 2 + | C => 3,2,1 | D => 1,2,3 | | 1 +---------------+ 3 3 +---------------+ 1 3 + 3 + | E => 2,3,1 | F => 1,3,2 | | 1 +---------------+ 2 2 +---------------+ 1 \endverbatim */ template std::array getLargestDistancePoints( PT pt1, PT pt2, PT pt3 ) { auto d12 = sqDist( pt1, pt2 ); auto d13 = sqDist( pt1, pt3 ); auto d23 = sqDist( pt2, pt3 ); PT* pA = 0; PT* pB = 0; PT* pM = 0; if( d12 > d13 ) // case B, D, E { pA = &pt2; if( d12 > d23 ) // case B, D { pB = &pt1; pM = &pt3; if( d13 > d23 ) std::swap( *pA, *pB ); } else // case E { pB = &pt3; pM = &pt1; } } else // case A, C, F { pA = &pt3; if( d13 > d23 ) // A, F { pB = &pt1; pM = &pt2; if( d12 > d23 ) std::swap( *pA, *pB ); } else // case C { pB = &pt2; pM = &pt1; } } return std::array{ *pA, *pB, *pM }; } } // namespace priv /// Free function template, product of two points, returns a line template Line2d_ operator * ( const Point2d_& lhs, const Point2d_& rhs ) { #ifndef HOMOG2D_NOCHECKS if( lhs == rhs ) HOMOG2D_THROW_ERROR_1( "points are identical, unable to compute product:" << lhs ); #endif Line2d_ line = detail::crossProduct(lhs, rhs); line.p_normalizePL(); return line; } //------------------------------------------------------------------ /// Returns true if the 3 points are on the same line /** \todo at present, defined by the distance between third point and line. Need to change that, and replace by computation of the angle between the two lines */ template bool areCollinear( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& pt3 ) { if( pt1 == pt2 || pt2 == pt3 || pt1 == pt3 ) return true; auto pt_arr = priv::getLargestDistancePoints( pt1, pt2, pt3 ); auto li = pt_arr[0] * pt_arr[1]; if( li.distTo( pt_arr[2] ) < thr::nullDistance() ) return true; return false; } /* /// Constructor: build Ellipse from Circle /// \todo finish this template Ellipse_::Ellipse_( const Circle_& cir ) { p_init( cir.center().getX(), cir.center().getY(), radius(), radius(), 0. ); } */ //------------------------------------------------------------------ /// Holds private stuff namespace priv { /// Helper function, factorized here for the two impl_getPoints_A() implementations template auto getPoints_B2( const Point2d_& pt, FPT2 dist, const Line2d_& li ) { auto arr = li.get(); const HOMOG2D_INUMTYPE a = static_cast(arr[0]); const HOMOG2D_INUMTYPE b = static_cast(arr[1]); auto coeff = static_cast(dist) / homog2d_sqrt( a*a + b*b ); Point2d_ pt1( pt.getX() - b * coeff, pt.getY() + a * coeff ); Point2d_ pt2( pt.getX() + b * coeff, pt.getY() - a * coeff ); fix_order( pt1, pt2 ); return std::make_pair( pt1, pt2 ); } /// Helper function for impl_getOrthogonalLine_A() and impl_getOrthogonalLine_B(), /// Compute orthogonal line to \c li at point \c pt (that must lie on the line) template Line2d_ getOrthogonalLine_B2( const Point2d_& pt, const Line2d_& li ) { auto arr = li.get(); // get array of 3 values Line2d_ out( -arr[1], arr[0], arr[1] * pt.getX() - arr[0] * pt.getY() ); out.p_normalizePL(); return out; } #ifdef HOMOG2D_DEBUGMODE template void printVector( const std::vector& v, std::string msg=std::string(), bool linefeed=false ) { std::cout << "vector: "; if( !msg.empty() ) std::cout << msg; std::cout << " #=" << v.size() << '\n'; size_t c=0; for( const auto& elem: v ) { if( linefeed ) std::cout << c++ << ": "; std::cout << elem << (linefeed?'\n':'-'); } std::cout << '\n'; } template void printArray( const std::array& v, std::string msg=std::string() ) { std::cout << "array: " << msg << " #=" << N<< '\n'; for( const auto& elem: v ) std::cout << elem << "-"; std::cout << '\n'; } template void printVectorPairs( const std::vector>& v ) { std::cout << "vector of pairs: #=" << v.size() << '\n'; for( const auto& elem: v ) std::cout << " [" << (int)elem.first << "-" << (int)elem.second << "] "; std::cout << '\n'; } #endif } // namespace priv // forward declaration template Line2d_ operator * ( const Point2d_&, const Point2d_& ); // forward declaration template Point2d_ operator * ( const Line2d_&, const Line2d_& ); namespace base { //------------------------------------------------------------------ /// Base class, will be instanciated as \ref Point2d_ or \ref Line2d_ /** Type parameters: - LP: typ::IsPoint or typ::IsLine - FPT: Floating Point Type (float, double or long double) */ template class LPBase: public detail::Common #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; using detail::Common::isInside; private: template friend class Hmatrix_; // This is needed so we can convert from, say, Point2d_ to Point2d_ template friend class LPBase; template friend Point2d_ h2d::operator * ( const h2d::Line2d_&, const h2d::Line2d_& ); template friend auto h2d::operator * ( const h2d::Point2d_&, const h2d::Point2d_& ) -> h2d::Line2d_; template friend auto h2d::operator * ( const h2d::Homogr_&, const h2d::Line2d_& ) -> h2d::Line2d_; template friend base::LPBase detail::crossProduct( const base::LPBase&, const base::LPBase& ); template friend auto operator << ( std::ostream& f, const h2d::base::LPBase& r ) -> std::ostream&; template friend void detail::product( base::LPBase&, const detail::Matrix_&, const base::LPBase& ); template friend Line2d_ priv::getOrthogonalLine_B2( const Point2d_&, const Line2d_& ); public: /// Constructor: build a point from two lines template LPBase( const Line2d_& v1, const Line2d_& v2 ) { #ifndef HOMOG2D_NOCHECKS if( v1.isParallelTo(v2) ) HOMOG2D_THROW_ERROR_1( "unable to build point from these two lines, are parallel" ); #endif *this = detail::crossProduct( v1, v2 ); p_normalizePL(); } /// Constructor: build a line from two points template LPBase( const Point2d_& v1, const Point2d_& v2 ) { #ifndef HOMOG2D_NOCHECKS if( v1 == v2 ) HOMOG2D_THROW_ERROR_1( "unable to build line from these two points, are the same: " << v1 ); #endif *this = detail::crossProduct( v1, v2 ); p_normalizePL(); } /// Constructor: copy-constructor for lines /** \todo We should be able to declare this "explicit". This fails at present when attempting to convert a line (or point) from double to float, but I don't get why... */ template // explicit LPBase( const Line2d_& li ) { impl_init_1_Line( li, detail::BaseHelper() ); } /// Constructor with single arg of type "Point" /** This will call one of the two overloads of \c impl_init_1_Point(), depending on type of object: - if type is a point, then it can be seen as a copy-constructor - if type is a line, this will build a line from (0,0] to \c pt */ template LPBase( const Point2d_& pt ) { impl_init_1_Point( pt, detail::BaseHelper() ); } /// Constructor: build from two numerical values, depends on the type template LPBase( const T1& v1, const T2& v2 ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); impl_init_2( v1, v2, detail::BaseHelper() ); } /// Constructor of line/point from 3 values template LPBase( T0 v0, T1 v1, T2 v2 ) { set( v0, v1, v2 ); } /// Assign homogeneous values template void set( T0 v0, T1 v1, T2 v2 ) { HOMOG2D_CHECK_IS_NUMBER(T0); HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); _v[0] = v0; _v[1] = v1; _v[2] = v2; p_normalizePL(); } /// Constructor of line from 4 values x1,y1,x2,y2 template LPBase( T x1, T y1, T x2, T y2 ) { HOMOG2D_CHECK_IS_NUMBER(T); impl_init_4( x1, y1, x2, y2, detail::BaseHelper() ); } #if 0 /// Constructor of Point/Line from random type holding x,y values, see manual, section template< typename T, typename std::enable_if< ! trait::IsContainer::value, T >::type* = nullptr > LPBase( T val ) { impl_init_2( val.HOMOG2D_BIND_X, val.HOMOG2D_BIND_Y, detail::BaseHelper() ); } #endif /// Constructor from an array holding 3 values of same type (a direct copy can be done) template< typename T, typename std::enable_if< std::is_same>::value, T >::type* = nullptr > LPBase( const T& arr ) { _v = arr; p_normalizePL(); } /// Constructor from an array/vector holding 3 values of different type template< typename T, typename std::enable_if< ( !std::is_same>::value && ( std::is_same>::value || std::is_same>::value ) ), T >::type* = nullptr > LPBase( const T& arr ) { _v[0] = static_cast(arr[0]); _v[1] = static_cast(arr[1]); _v[2] = static_cast(arr[2]); p_normalizePL(); } /// Default constructor, depends on the type LPBase() { impl_init( detail::BaseHelper() ); } /// Constructor of horizontal/vertical line template LPBase( LineDir orient, T value ) { // HOMOG2D_CHECK_IS_NUMBER(T); impl_init_or( orient, value, detail::BaseHelper() ); } #ifdef HOMOG2D_USE_BOOSTGEOM public: /// Constructor from boost::geometry point type /** \note Although this one should work also for the other point type (`bg::model::d2::point_xy`), as that latter one is inherited from the first one), it does not, because there is another truely generic single-arg constructor, and the compiler will select that one first, leading to a build error. Thus, we need the second one. */ template // Boost Floating Point Type LPBase( const boost::geometry::model::point& pt ) { impl_init_BoostGeomPoint( pt, detail::BaseHelper() ); } /// Constructor from boost::geometry second point type template // Boost Floating Point Type LPBase( const boost::geometry::model::d2::point_xy& pt ) { impl_init_BoostGeomPoint( pt, detail::BaseHelper() ); } /// Set from boost::geometry point type template // Boost Floating Point Type void set( const boost::geometry::model::point& pt ) { set( boost::geometry::get<0>(pt), boost::geometry::get<1>(pt), 1.0 ); } #endif private: template void p_copyFrom( const LPBase& other ) { _v[0] = static_cast(other._v[0]); _v[1] = static_cast(other._v[1]); _v[2] = static_cast(other._v[2]); } /// Arg is a point, object is a point => copy-constructor template void impl_init_1_Point( const Point2d_& pt, const detail::BaseHelper& ) { p_copyFrom( pt ); } /// Arg is a point, object is a line: we build the line passing though (0,0) ant the given point template void impl_init_1_Point( const Point2d_& pt, const detail::BaseHelper& ) { *this = detail::crossProduct( pt, Point2d_() ); p_normalizePL(); } /// Arg is a line, object is a point: ILLEGAL INSTANCIATION template constexpr void impl_init_1_Line( const Line2d_&, const detail::BaseHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot build a point from a line" ); } /// Arg is a line, object is a line => copy-constructor template void impl_init_1_Line( const Line2d_& li, const detail::BaseHelper& ) { p_copyFrom( li ); } template constexpr void impl_init_or( LineDir, T, const detail::BaseHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot build a horiz/vertical point" ); } template< typename T, typename std::enable_if< std::is_arithmetic::value ,T >::type* = nullptr > void impl_init_or( LineDir dir, T value, const detail::BaseHelper& ) { _v[2] = -value; if( dir == LineDir::V ) { _v[0] = 1.; _v[1] = 0.; } else // = LineDir::H { _v[0] = 0.; _v[1] = 1.; } } template void impl_init_or( LineDir dir, const Point2d_& pt, const detail::BaseHelper& ) { if( dir == LineDir::V ) { _v[2] = -pt.getX(); _v[0] = 1.; _v[1] = 0.; } else // = LineDir::H { _v[2] = -pt.getY(); _v[0] = 0.; _v[1] = 1.; } } template constexpr void impl_init_4( T, T, T, T, const detail::BaseHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot build a point from 4 values" ); } template void impl_init_4( T x1, T y1, T x2, T y2, const detail::BaseHelper& ) { *this = Point2d_(x1, y1) * Point2d_(x2, y2); } public: Type type() const { return impl_type( detail::BaseHelper() ); } private: Type impl_type( const detail::BaseHelper& ) const { return Type::Point2d; } Type impl_type( const detail::BaseHelper& ) const { return Type::Line2d; } public: FPT getCoord( GivenCoord gc, FPT other ) const { return impl_getCoord( gc, other, detail::BaseHelper() ); } Point2d_ getPoint( GivenCoord gc, FPT other ) const { return impl_getPoint( gc, other, detail::BaseHelper() ); } /// Returns a pair of points that are lying on line at distance \c dist from a point defined by one of its coordinates. template PointPair1_ getPoints( GivenCoord gc, FPT coord, FPT2 dist ) const { return impl_getPoints_A( gc, coord, dist, detail::BaseHelper() ); } /// Returns a pair of points that are lying on line at distance \c dist from point \c pt, assuming that one is lying on the line. template PointPair1_ getPoints( const Point2d_& pt, FPT2 dist ) const { return impl_getPoints_B( pt, dist, detail::BaseHelper() ); } /// Returns an orthogonal line to the one it is called on, at a point defined by one of its coordinates. Line2d_ getOrthogonalLine( GivenCoord gc, FPT other ) const { return impl_getOrthogonalLine_A( gc, other, detail::BaseHelper() ); } /// Returns an orthogonal line to the one it is called on, at point \c pt, assuming that one is lying on the line. Line2d_ getOrthogonalLine( const Point2d_& pt ) const { return impl_getOrthogonalLine_B( pt, detail::BaseHelper() ); } /// Returns a line rotated at point \c pt with angle \c theta template Line2d_ getRotatedLine( const Point2d_& pt, T theta ) const { HOMOG2D_CHECK_IS_NUMBER(T); return impl_getRotatedLine( pt, theta, detail::BaseHelper() ); } /// Returns the segment from the point (not on line) to the line, shortest path Segment_ getOrthogSegment( const Point2d_& pt ) const { return impl_getOrthogSegment( pt, detail::BaseHelper() ); } /// Returns an parallel line to the one it is called on, with \c pt lying on it. /// \todo clarify orientation: on wich side will that line appear? Line2d_ getParallelLine( const Point2d_& pt ) const { return impl_getParallelLine( pt, detail::BaseHelper() ); } /// Returns the pair of parallel lines at a distance \c dist from line. template std::pair,Line2d_> getParallelLines( T dist ) const { return impl_getParallelLines( dist, detail::BaseHelper() ); } HOMOG2D_INUMTYPE getX() const { return impl_getX( detail::BaseHelper() ); } HOMOG2D_INUMTYPE getY() const { return impl_getY( detail::BaseHelper() ); } /// Translate Point2d, does nothing for Line2d ( \sa impl_translate() ) template void translate( T1 dx, T2 dy ) { HOMOG2D_CHECK_IS_NUMBER( T1 ); HOMOG2D_CHECK_IS_NUMBER( T2 ); impl_translate( dx, dy, detail::BaseHelper() ); } /// Translate Point2d, does nothing for Line2d ( \sa impl_translate() ) template void translate( const std::pair& pa ) { this->translate( pa.first, pa.second ); } /// Move point to other location (same as set(), but this one will be virtual). /// Does nothing for lines template void moveTo( const Point2d_& pt ) { impl_moveTo( pt, detail::BaseHelper() ); } #ifdef HOMOG2D_ENABLE_VRTP /// Needed because of variant (\sa CommonType) FRect_ getBB() const { HOMOG2D_THROW_ERROR_1( "invalid call, Point/Line has no Bounding Box" ); } #endif private: template ANY impl_getPt( const detail::BaseHelper& ) const { return ANY( getX(), getY() ); } template ANY impl_getPt( const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid call for lines" ); } public: /// Generic transformation into any other point type, as long as it provides a 2-args constructor /// (is the case for Opencv and Boost Geometry). /// \sa h2d::getPt() template ANY getPt() const { return impl_getPt( detail::BaseHelper() ); } std::array get() const { return std::array { _v[0], _v[1], _v[2] }; } template void set( T1 x, T2 y ) { impl_set( x, y, detail::BaseHelper() ); } template HOMOG2D_INUMTYPE distTo( const Point2d_& pt ) const { return impl_distToPoint( pt, detail::BaseHelper() ); } template HOMOG2D_INUMTYPE distTo( const Line2d_& li ) const { return impl_distToLine( li, detail::BaseHelper() ); } template HOMOG2D_INUMTYPE distTo( const Segment_& seg ) const { return impl_distToSegment( seg, detail::BaseHelper() ); } template bool isParallelTo( const LPBase& li ) const { return impl_isParallelTo( li, detail::BaseHelper() ); } template bool isParallelTo( const Segment_& seg ) const { return impl_isParallelTo( seg.getLine(), detail::BaseHelper() ); } /// Returns angle in rad. between the lines. \sa h2d::getAngle() /** Please check out warning described in impl_getAngle() */ template HOMOG2D_INUMTYPE getAngle( const LPBase& other ) const { return impl_getAngle( other, detail::BaseHelper() ); } /// Returns angle in rad. between line and segment \c seg. \sa h2d::getAngle() template HOMOG2D_INUMTYPE getAngle( const Segment_& seg ) const { return impl_getAngle( seg.getLine(), detail::BaseHelper() ); } /// Returns true if point is at infinity (third value less than thr::nullDenom() ) bool isInf() const { return impl_isInf( detail::BaseHelper() ); } HOMOG2D_INUMTYPE length() const { return impl_length( detail::BaseHelper() ); } /// Neither lines nor points have an area HOMOG2D_INUMTYPE area() const { return 0.; } private: /// A point has a null length HOMOG2D_INUMTYPE impl_length( const detail::BaseHelper& ) const { return 0.; } /// A line has an infinite length HOMOG2D_INUMTYPE impl_length( const detail::BaseHelper& ) const { HOMOG2D_THROW_ERROR_1( "unable, a line has an infinite length" ); } HOMOG2D_INUMTYPE impl_getX( const detail::BaseHelper& ) const { return static_cast(_v[0])/_v[2]; } HOMOG2D_INUMTYPE impl_getY( const detail::BaseHelper& ) const { return static_cast(_v[1])/_v[2]; } template void impl_set( T1 x, T2 y, const detail::BaseHelper& ) { HOMOG2D_CHECK_IS_NUMBER( T1 ); HOMOG2D_CHECK_IS_NUMBER( T2 ); _v[0] = x; _v[1] = y; _v[2] = 1.; p_normalizePL(); } template constexpr void impl_set( T1, T2, const detail::BaseHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid call for lines" ); } template void impl_translate( T1 dx, T2 dy, const detail::BaseHelper& ) { _v[0] = static_cast(_v[0]) / _v[2] + dx; _v[1] = static_cast(_v[1]) / _v[2] + dy; _v[2] = 1.; p_normalizePL(); } template constexpr void impl_translate( T1, T2, const detail::BaseHelper& ) { // changed on 20231216: this does now nothing. That way, it can become a virtual function, if needed // static_assert( detail::AlwaysFalse::value, "Invalid call for lines" ); } template void impl_moveTo( const Point2d_& pt, const detail::BaseHelper& ) { *this = pt; } template void impl_moveTo( const Point2d_& pt, const detail::BaseHelper& ) { } template HOMOG2D_INUMTYPE impl_distToPoint( const Point2d_&, const detail::BaseHelper& ) const; template HOMOG2D_INUMTYPE impl_distToPoint( const Point2d_&, const detail::BaseHelper& ) const; template HOMOG2D_INUMTYPE impl_distToLine( const Line2d_&, const detail::BaseHelper& ) const; template constexpr HOMOG2D_INUMTYPE impl_distToLine( const Line2d_&, const detail::BaseHelper& ) const; template HOMOG2D_INUMTYPE impl_distToSegment( const Segment_&, const detail::BaseHelper& ) const; template constexpr HOMOG2D_INUMTYPE impl_distToSegment( const Segment_&, const detail::BaseHelper& ) const; HOMOG2D_INUMTYPE impl_getAngle( const LPBase&, const detail::BaseHelper& ) const; constexpr HOMOG2D_INUMTYPE impl_getAngle( const LPBase&, const detail::BaseHelper& ) const; constexpr bool impl_isInf( const detail::BaseHelper& ) const { return false; } bool impl_isInf( const detail::BaseHelper& ) const { return homog2d_abs( _v[2] ) < thr::nullDenom(); } template bool impl_isParallelTo( const LPBase&, const detail::BaseHelper& ) const; template constexpr bool impl_isParallelTo( const LPBase&, const detail::BaseHelper& ) const; FPT impl_getCoord( GivenCoord gc, FPT other, const detail::BaseHelper& ) const; constexpr FPT impl_getCoord( GivenCoord gc, FPT other, const detail::BaseHelper& ) const; Point2d_ impl_getPoint( GivenCoord gc, FPT other, const detail::BaseHelper& ) const; constexpr Point2d_ impl_getPoint( GivenCoord gc, FPT other, const detail::BaseHelper& ) const; template PointPair1_ impl_getPoints_A( GivenCoord, FPT, FPT2, const detail::BaseHelper& ) const; template constexpr PointPair1_ impl_getPoints_A( GivenCoord, FPT, FPT2, const detail::BaseHelper& ) const; template PointPair1_ impl_getPoints_B( const Point2d_&, FPT2, const detail::BaseHelper& ) const; template constexpr PointPair1_ impl_getPoints_B( const Point2d_&, FPT2, const detail::BaseHelper& ) const; void impl_op_stream( std::ostream&, const Point2d_& ) const; void impl_op_stream( std::ostream&, const Line2d_& ) const; public: /// Line/Line intersection template detail::Intersect intersects( const Line2d_& other ) const { detail::Intersect out; if( this->isParallelTo( other ) ) return out; out.set( *this * other ); return out; } /// Line/FRect intersection (rectangle defined by \c pt1 and \c pt2) template detail::IntersectM intersects( const Point2d_& pt1, const Point2d_& pt2 ) const { return intersects( FRect_( pt1, pt2 ) ) ; } /// Line/FRect intersection template detail::IntersectM intersects( const FRect_& rect ) const { return impl_intersectsFRect( rect, detail::BaseHelper() ); } /// Line/Segment intersection /** \warning no implementation for points */ template detail::Intersect intersects( const Segment_& seg ) const { return seg.intersects( *this ); } /// Line/Circle intersection /**
The Sfinae below is needed to avoid ambiguity with the other 2 args "intersects()" function (with 2 points defining a FRect, see above) */ template< typename T, typename std::enable_if< (std::is_arithmetic::value && !std::is_same::value) ,T >::type* = nullptr > detail::Intersect intersects( const Point2d_& pt0, T radius ) const { return impl_intersectsCircle( pt0, radius, detail::BaseHelper() ); } /// Line/Circle intersection template detail::Intersect intersects( const Circle_& cir ) const { return impl_intersectsCircle( cir.center(), cir.radius(), detail::BaseHelper() ); } /// Line/Polyline intersection template detail::IntersectM intersects( const base::PolylineBase& pl ) const { return pl.intersects( *this ); } /// Point is inside flat rectangle bool isInside( const Point2d_& pt1, const Point2d_& pt2 ) const { HOMOG2D_START; return impl_isInsideRect( FRect_(pt1, pt2), detail::BaseHelper() ); } /// Point is inside FRect template bool isInside( const FRect_& rect ) const { HOMOG2D_START; return impl_isInsideRect( rect, detail::BaseHelper() ); } /// Point is inside circle defined by center and radius template bool isInside( const Point2d_& center, T radius ) const { HOMOG2D_CHECK_IS_NUMBER(T); HOMOG2D_START; return impl_isInsideCircle( center, radius, detail::BaseHelper() ); } /// Point is inside Circle template bool isInside( const Circle_& cir ) const { HOMOG2D_START; return impl_isInsideCircle( cir.center(), cir.radius(), detail::BaseHelper() ); } /// Point is inside Ellipse template bool isInside( const Ellipse_& ell ) const { HOMOG2D_START; return impl_isInsideEllipse( ell, detail::BaseHelper() ); } /// Point or line is inside Polyline template bool isInside( const base::PolylineBase& poly ) const { HOMOG2D_START; return impl_isInsidePoly( poly, detail::BaseHelper() ); } ////////////////////////// // OPERATORS // ////////////////////////// bool operator == ( const base::LPBase& other ) const { return impl_op_equal( other, detail::BaseHelper() ); } bool operator != ( const base::LPBase& other ) const { return !(*this == other); } bool operator < ( const base::LPBase& other ) const { return impl_op_sort( other, detail::BaseHelper() ); } /// SVG draw function void draw( img::Image& im, img::DrawParams dp=img::DrawParams() ) const { impl_draw_LP( im, dp, detail::BaseHelper() ); } #ifdef HOMOG2D_USE_OPENCV private: template OPENCVT impl_getCvPt( const detail::BaseHelper&, const OPENCVT& ) const { return OPENCVT( getX(),getY() ); } /// Build point from Opencv point template void impl_init_opencv( cv::Point_ pt, const detail::BaseHelper& ) { impl_init_2( pt.x, pt.y, detail::BaseHelper() ); } /// Build line from Opencv point template void impl_init_opencv( cv::Point_ pt, const detail::BaseHelper& ) { Point2d_ p(pt); impl_init_1_Point( p, detail::BaseHelper() ); } public: /// Opencv draw function void draw( img::Image& im, img::DrawParams dp=img::DrawParams() ) const { impl_draw_LP( im, dp, detail::BaseHelper() ); } cv::Point2i getCvPti() const { return impl_getPt( detail::BaseHelper() ); } cv::Point2i getCvPtd() const { return impl_getPt( detail::BaseHelper() ); } cv::Point2i getCvPtf() const { return impl_getPt( detail::BaseHelper() ); } /// Constructor: build from a single OpenCv point. template LPBase( cv::Point_ pt ) { impl_init_opencv( pt, detail::BaseHelper() ); } #endif // HOMOG2D_USE_OPENCV ////////////////////////// // DATA SECTION // ////////////////////////// private: std::array _v; ///< data, uses the template parameter FPT (for "Floating Point Type") ////////////////////////// // PRIVATE FUNCTIONS // ////////////////////////// private: void p_normalizePL() const { impl_normalize( detail::BaseHelper() ); } void impl_normalize( const detail::BaseHelper& ) const; void impl_normalize( const detail::BaseHelper& ) const; template detail::IntersectM impl_intersectsFRect( const FRect_&, const detail::BaseHelper& ) const; template constexpr detail::IntersectM impl_intersectsFRect( const FRect_&, const detail::BaseHelper& ) const; template detail::Intersect impl_intersectsCircle( const Point2d_&, T r, const detail::BaseHelper& ) const; template constexpr detail::Intersect impl_intersectsCircle( const Point2d_&, T r, const detail::BaseHelper& ) const; template bool impl_isInsideRect( const FRect_&, const detail::BaseHelper& ) const; template constexpr bool impl_isInsideRect( const FRect_&, const detail::BaseHelper& ) const; template bool impl_isInsideEllipse( const Ellipse_&, const detail::BaseHelper& ) const; template constexpr bool impl_isInsideEllipse( const Ellipse_&, const detail::BaseHelper& ) const; template bool impl_isInsideCircle( const Point2d_&, T r, const detail::BaseHelper& ) const; template constexpr bool impl_isInsideCircle( const Point2d_&, T r, const detail::BaseHelper& ) const; template bool impl_isInsidePoly( const base::PolylineBase&, const detail::BaseHelper& ) const; template constexpr bool impl_isInsidePoly( const base::PolylineBase&, const detail::BaseHelper& ) const; Line2d_ impl_getOrthogonalLine_A( GivenCoord, FPT, const detail::BaseHelper& ) const; constexpr Line2d_ impl_getOrthogonalLine_A( GivenCoord, FPT, const detail::BaseHelper& ) const; Line2d_ impl_getOrthogonalLine_B( const Point2d_&, const detail::BaseHelper& ) const; constexpr Line2d_ impl_getOrthogonalLine_B( const Point2d_&, const detail::BaseHelper& ) const; template Line2d_ impl_getRotatedLine( const Point2d_&, T, const detail::BaseHelper& ) const; template constexpr Line2d_ impl_getRotatedLine( const Point2d_&, T, const detail::BaseHelper& ) const; Segment_ impl_getOrthogSegment( const Point2d_&, const detail::BaseHelper& ) const; constexpr Segment_ impl_getOrthogSegment( const Point2d_&, const detail::BaseHelper& ) const; Line2d_ impl_getParallelLine( const Point2d_&, const detail::BaseHelper& ) const; constexpr Line2d_ impl_getParallelLine( const Point2d_&, const detail::BaseHelper& ) const; template std::pair,Line2d_> impl_getParallelLines( T, const detail::BaseHelper& ) const; template constexpr std::pair,Line2d_> impl_getParallelLines( T, const detail::BaseHelper& ) const; bool impl_op_equal( const LPBase&, const detail::BaseHelper& ) const; bool impl_op_equal( const LPBase&, const detail::BaseHelper& ) const; bool impl_op_sort( const LPBase&, const detail::BaseHelper& ) const; constexpr bool impl_op_sort( const LPBase&, const detail::BaseHelper& ) const; Point2d_ impl_op_product( const Line2d_& , const Line2d_& , const detail::BaseHelper& ) const; Line2d_ impl_op_product( const Point2d_&, const Point2d_&, const detail::BaseHelper& ) const; template void impl_draw_LP( img::Image&, img::DrawParams, const detail::BaseHelper& ) const; template void impl_draw_LP( img::Image&, img::DrawParams, const detail::BaseHelper& ) const; /// Called by default constructor, overload for lines void impl_init( const detail::BaseHelper& ) { _v[0] = 1.; _v[1] = 0.; _v[2] = 0.; } /// Called by default constructor, overload for points. Initialize to (0,0) void impl_init( const detail::BaseHelper& ) { _v[0] = 0.; _v[1] = 0.; _v[2] = 1.; } template void impl_init_2( const T1&, const T2&, const detail::BaseHelper& ); template void impl_init_2( const T1&, const T2&, const detail::BaseHelper& ); #ifdef HOMOG2D_USE_BOOSTGEOM template void impl_init_BoostGeomPoint( const boost::geometry::model::point&, const detail::BaseHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot build a Line2d using a boost::geometry point" ); } template // Boost Floating Point Type void impl_init_BoostGeomPoint( const boost::geometry::model::point& pt, const detail::BaseHelper& ) { set( pt ); } #endif }; // class LPBase } // namespace base ///////////////////////////////////////////////////////////////////////////// // SECTION - OPENCV API CODE ///////////////////////////////////////////////////////////////////////////// #ifdef HOMOG2D_USE_OPENCV template FRect_ getFRect( cv::Mat& mat ) { if( mat.cols == 0 || mat.rows == 0 ) HOMOG2D_THROW_ERROR_1( "Illegal values: cols=" << mat.cols << ", rows=" << mat.rows ); return FRect_( Point2d_(), // (0,0) Point2d_( mat.cols, mat.rows ) // (w,h) ); } /// Free function to return an OpenCv point (double) template cv::Point2d getCvPtd( const Point2d_& pt ) { return pt.getCvPtd(); } /// Free function to return an OpenCv point (float) template cv::Point2f getCvPtf( const Point2d_& pt ) { return pt.getCvPtf(); } /// Free function to return an OpenCv point (integer) template cv::Point2i getCvPti( const Point2d_& pt ) { return pt.getCvPti(); } #endif // HOMOG2D_USE_OPENCV ///////////////////////////////////////////////////////////////////////////// // SECTION - FREE FUNCTIONS ///////////////////////////////////////////////////////////////////////////// /// Generic free function to return a point of other type /** - RT: return type - FPT: Floating Point Type User code needs to provide the requested type as template argument, for example: \code auto p1 = getPt( pt ); // opencv type auto p2 = getPt>( pt ); // boost geometry type \endcode \sa LPBase::getPt() */ template RT getPt( const Point2d_& pt ) { return pt.template getPt(); } /// Free function, returns a vector of points of other type from a vector of h2d points /** - RT: return type - FPT: Floating Point Type User code needs to provide the requested type as template argument. \sa h2d::getPt<>() \todo 20230219: sfinae this to accept other containers using trait::IsContainer */ template std::vector getPts( const std::vector>& vpt ) { std::vector vout( vpt.size() ); auto it = vout.begin(); for( const auto& pt: vpt ) *it++ = getPt(pt); return vout; } //------------------------------------------------------------------ /// This namespace holds some private stuff, types here are not to be used directly by end-user code. namespace detail { //------------------------------------------------------------------ /// Private free function, returns true if point \c pt is inside the rectangle defined by (\c p00 , \c p11) template bool ptIsInside( const Point2d_& pt, const Point2d_& p00, const Point2d_& p11 ) { if( pt.getX() > p00.getX() && pt.getX() < p11.getX() ) if( pt.getY() > p00.getY() && pt.getY() < p11.getY() ) return true; return false; } #ifdef HOMOG2D_USE_EIGEN /// Build Homography from 2 sets of 4 points, using Eigen /** See - https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html - https://eigen.tuxfamily.org/dox/group__DenseMatrixManipulation__chapter.html */ template Homogr_ buildFrom4Points_Eigen( const std::vector>& vpt1, ///< source points const std::vector>& vpt2 ///< destination points ) { Eigen::MatrixXd A = Eigen::MatrixXd::Zero(8,8); Eigen::VectorXd b(8); for( int i=0; i<4; i++ ) { auto u1 = vpt1[i].getX(); auto v1 = vpt1[i].getY(); auto u2 = vpt2[i].getX(); auto v2 = vpt2[i].getY(); b(2*i) = u2; b(2*i+1) = v2; A(i*2,0) = A(i*2+1,3) = u1; A(i*2,1) = A(i*2+1,4) = v1; A(i*2,2) = A(i*2+1,5) = 1.; A(i*2,6) = - u1 * u2; A(i*2,7) = - v1 * u2; A(i*2+1,6) = - u1 * v2; A(i*2+1,7) = - v1 * v2; } #if 0 Eigen::VectorXd X = A.ldlt().solve(b); // for some reason this does not work... #else Eigen::MatrixXd Ai = A.inverse(); Eigen::VectorXd X = Ai * b; #endif Homogr_ H; for( int i=0; i<8; i++ ) H.set( i/3, i%3, X(i) ); H.set(2, 2, 1.); return H; } #endif // HOMOG2D_USE_EIGEN //------------------------------------------------------------------ #ifdef HOMOG2D_USE_OPENCV /// Build Homography from 2 sets of 4 points, using Opencv /** - see https://docs.opencv.org/master/d9/d0c/group__calib3d.html#ga4abc2ece9fab9398f2e560d53c8c9780 \note With current Opencv installed on current machine, it seems that \c cv::getPerspectiveTransform() requires that the points are "CV_32F" (\c float), and NOT double. */ template Homogr_ buildFrom4Points_Opencv ( const std::vector>& vpt1, ///< source points const std::vector>& vpt2 ///< destination points ) { const auto& src = getPts( vpt1 ); const auto& dst = getPts( vpt2 ); return cv::getPerspectiveTransform( src, dst ); // automatic type conversion to Hmatrix_ } #endif } // namespace detail //------------------------------------------------------------------ /// Build Homography from 2 sets of 4 points (free function) /// \sa Homogr_::buildFrom4Points() template Homogr_ buildFrom4Points( const std::vector>& vpt1, ///< source points const std::vector>& vpt2, ///< destination points int method=1 ///< 0: Eigen, 1: Opencv ) { Homogr_ H; H.buildFrom4Points( vpt1, vpt2, method ); return H; } //------------------------------------------------------------------ /// Build Homography from 2 sets of 4 points /** - Requires either Eigen or Opencv - we build a 8x8 matrix A and a 8x1 vector B, and get the solution from X = A^-1 B - see this for details: https://skramm.lautre.net/files/misc/Kramm_compute_H_from_4pts.pdf \sa free function: h2d::buildFrom4Points() \todo fix this so that user can provide a std::array of points */ template void Hmatrix_::buildFrom4Points( const std::vector>& vpt1, ///< source points const std::vector>& vpt2, ///< destination points int method ///< 0: Eigen, 1: Opencv (default) ) { if( vpt1.size() != 4 ) HOMOG2D_THROW_ERROR_1( "invalid vector size for source points, should be 4, value=" << vpt1.size() ); if( vpt2.size() != 4 ) HOMOG2D_THROW_ERROR_1( "invalid vector size for dest points, should be 4, value=" << vpt2.size() ); assert( method == 0 || method == 1 ); if( method == 0 ) { #ifdef HOMOG2D_USE_EIGEN *this = detail::buildFrom4Points_Eigen( vpt1, vpt2 ); #else throw std::runtime_error( "Unable, build without Eigen support" ); #endif } else { #ifdef HOMOG2D_USE_OPENCV *this = detail::buildFrom4Points_Opencv( vpt1, vpt2 ); #else throw std::runtime_error( "Unable, build without Opencv support" ); #endif } } //------------------------------------------------------------------ /// A line segment, defined by two points /** - Storage: "smallest" point is always stored as first element (see constructor) */ template class Segment_: public detail::Common #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; using detail::Common::isInside; template friend class Segment_; Type type() const { return Type::Segment; } private: Point2d_ _ptS1, _ptS2; public: /// \name Constructors ///@{ /// Default constructor: initializes segment to (0,0)--(1,1) Segment_(): _ptS2(1.,1.) {} /// Constructor 2: build segment from two points Segment_( Point2d_ p1, Point2d_ p2 ) : _ptS1(p1), _ptS2(p2) { #ifndef HOMOG2D_NOCHECKS if( p1 == p2 ) HOMOG2D_THROW_ERROR_1( "cannot build a segment with two identical points: " << p1 << " and " << p2 ); #endif priv::fix_order( _ptS1, _ptS2 ); } /// Constructor 3: build segment from two points coordinates, call constructor 2 template Segment_( T x1, T y1, T x2, T y2 ) : Segment_( Point2d_(x1,y1), Point2d_(x2,y2) ) { HOMOG2D_CHECK_IS_NUMBER(T); } /// Constructor 4: build segment from pair of points Segment_( const PointPair1_& ppts ) : Segment_(ppts.first, ppts.second) {} /// Copy-Constructor template Segment_( const Segment_& other ) : _ptS1(other._ptS1), _ptS2(other._ptS2) { // priv::fix_order( _ptS1, _ptS2 ); no need to call this, because the source segment was already OK } ///@} /// \name Modifying functions ///@{ /// Setter template void set( const Point2d_& p1, const Point2d_& p2 ) { #ifndef HOMOG2D_NOCHECKS if( p1 == p2 ) HOMOG2D_THROW_ERROR_1( "cannot define a segment with two identical points" << p1 << " and " << p2 ); #endif _ptS1 = p1; _ptS2 = p2; priv::fix_order( _ptS1, _ptS2 ); } /// Setter from a std::pair (points need to be of same underlying type) template void set( const PointPair1_& ppts ) { set( ppts.first, ppts.second ); } /// Setter from 4 raw point coordinates template void set( FPT2 x1, FPT2 y1, FPT2 x2, FPT2 y2 ) { set( Point2d_(x1,y1), Point2d_(x2,y2) ); } /// Translate Segment template void translate( T1 dx, T2 dy ) { HOMOG2D_CHECK_IS_NUMBER( T1 ); HOMOG2D_CHECK_IS_NUMBER( T2 ); _ptS1.translate( dx, dy ); _ptS2.translate( dx, dy ); } /// Translate Segment, using pair of values template void translate( const std::pair& pa ) { this->translate( pa.first, pa.second ); } /// Move Segment to other location template void moveTo( TX x, TY y ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); moveTo( Point2d_(x,y) ); } /// Move Segment to other location, given by \c pt template void moveTo( const Point2d_& pt ) { auto w = _ptS2.getX() - _ptS1.getX(); auto h = _ptS2.getY() - _ptS1.getY(); _ptS1 = pt; _ptS2.set( _ptS1.getX() + w, _ptS1.getY() + h ); } ///@} /// \name Attributes access ///@{ /// Get segment length HOMOG2D_INUMTYPE length() const { return _ptS1.distTo( _ptS2 ); } /// A segment always has a null area // /*constexpr*/ HOMOG2D_INUMTYPE area() const constexpr HOMOG2D_INUMTYPE area() const { return 0.; } /// Get angle between segment and other segment/line /** This will call the line angle function, thus the returned value will be in the range \f$ [0,\pi/2] \f$ */ template HOMOG2D_INUMTYPE getAngle( const U& other ) const { return other.getAngle( this->getLine() ); } #ifdef HOMOG2D_ENABLE_VRTP /// Needed because of variant (\sa CommonType) FRect_ getBB() const { HOMOG2D_THROW_ERROR_1( "invalid call, segment has no Bounding Box" ); } #endif ///@} /// \name Operators ///@{ bool operator == ( const Segment_& s2 ) const { if( _ptS1 != s2._ptS1 ) return false; if( _ptS2 != s2._ptS2 ) return false; return true; } bool operator != ( const Segment_& s2 ) const { return !(*this == s2); } bool operator < ( const Segment_& other ) const { return _ptS1 < other._ptS1; } ///@} private: template bool p_bothPtsAreInside( const S& shape ) const { if( !_ptS1.isInside( shape ) ) return false; if( !_ptS2.isInside( shape ) ) return false; return true; } public: /// \name Enclosing functions ///@{ /// Segment is inside Circle template bool isInside( const Circle_& shape ) const { return p_bothPtsAreInside( shape ); } /// Segment is inside FRect template bool isInside( const FRect_& shape ) const { return p_bothPtsAreInside( shape ); } /// Segment is inside Ellipse template bool isInside( const Ellipse_& shape ) const { return p_bothPtsAreInside( shape ); } /// Segment is inside OPolyline_ template constexpr bool isInside( const OPolyline_& ) const { return false; } /// Segment is inside CPolyline_ /** Requires both points inside AND no intersections */ template bool isInside( const CPolyline_& cpoly ) const { if( !p_bothPtsAreInside( cpoly ) ) return false; for( auto poly_seg: cpoly.getSegs() ) if( poly_seg.intersects( *this )() ) return false; return true; } ///@} /// \name Other const functions ///@{ /// Returns the points as a std::pair /** The one with smallest x coordinate will be returned as "first". If x-coordinate are equal, then the one with smallest y-coordinate will be returned first */ PointPair1_ getPts() const { return std::make_pair( _ptS1, _ptS2 ); } template std::pair getParallelSegs( T dist ) const { HOMOG2D_CHECK_IS_NUMBER( T ); auto plines = getLine().getParallelLines( dist ); auto lo1 = getLine().getOrthogonalLine( _ptS1 ); auto lo2 = getLine().getOrthogonalLine( _ptS2 ); auto pA1 = lo1 * plines.first; auto pA2 = lo2 * plines.first; auto pB1 = lo1 * plines.second; auto pB2 = lo2 * plines.second; return std::make_pair( Segment_( pA1, pA2 ), Segment_( pB1, pB2 ) ); } private: /// Computes the 4 points orthogonal to the segment std::array,4> p_getOrthog() const { auto x1 = static_cast(_ptS1.getX()); auto x2 = static_cast(_ptS2.getX()); auto y1 = static_cast(_ptS1.getY()); auto y2 = static_cast(_ptS2.getY()); auto dx = x1 - x2; auto dy = y1 - y2; std::array,4> out; out[0] = Point2d_( x1-dy, y1+dx ); out[1] = Point2d_( x1+dy, y1-dx ); out[3] = Point2d_( x2-dy, y2+dx ); out[2] = Point2d_( x2+dy, y2-dx ); return out; } public: /// Returns the 4 points orthogonal to the segment std::array,4> getOrthogPts() const { return p_getOrthog(); } /// Returns the 4 segments orthogonal to the segment std::array,4> getOrthogSegs() const { std::array,4> out; auto pts = p_getOrthog(); out[0] = Segment_( _ptS1, pts[0] ); out[1] = Segment_( _ptS1, pts[1] ); out[2] = Segment_( _ptS2, pts[2] ); out[3] = Segment_( _ptS2, pts[3] ); return out; } /// Returns supporting line Line2d_ getLine() const { Point2d_ pt1( _ptS1 ); Point2d_ pt2( _ptS2 ); auto li = pt1 * pt2; return Line2d_(li); } /// Returns a pair of segments split by the middle std::pair,Segment_> split() const { auto pt_mid = getCenter(); return std::make_pair( Segment_( _ptS1, pt_mid ), Segment_( _ptS2, pt_mid ) ); } template HOMOG2D_INUMTYPE distTo( const Point2d_&, int* segDistCase=0 ) const; template constexpr bool isParallelTo( const T& other ) const { static_assert( std::is_same>::value || std::is_same>::value, "type needs to be a segment or a line" ); return getLine().isParallelTo( other ); } /// Returns point that at middle distance between \c p1 and \c p2 Point2d_ getCenter() const { return Point2d_( ( static_cast(_ptS1.getX()) + _ptS2.getX() ) / 2., ( static_cast(_ptS1.getY()) + _ptS2.getY() ) / 2. ); } Segment_ getExtended() const; /// Returns the bisector line of the segment /// \sa free function h2d::getBisector() Line2d_ getBisector() const { Segment_ seg2 = *this; // convert to (possibly) enhance precision return seg2.getLine().getOrthogonalLine( seg2.getCenter() ); } ///@} /// \name Intersection functions ///@{ template detail::Intersect intersects( const Segment_& ) const; template detail::Intersect intersects( const Line2d_& ) const; template detail::IntersectM intersects( const Circle_& ) const; /// Segment/FRect intersection template detail::IntersectM intersects( const FRect_& r ) const { // HOMOG2D_START; return r.intersects( *this ); } /// Segment/Polyline intersection template detail::IntersectM intersects( const base::PolylineBase& other ) const { // HOMOG2D_START; return other.intersects( *this ); } ///@} template friend std::ostream& operator << ( std::ostream& f, const Segment_& seg ); #ifdef HOMOG2D_USE_OPENCV void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; #endif void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; }; // class Segment_ //------------------------------------------------------------------ /// Returns a segment with same support line but tripled length. /** With (1,0)-(2,0) as input, will return the segment (0,0)-(3,0) */ template Segment_ Segment_::getExtended() const { Segment_ seg(*this); // to get highest precision auto li = seg.getLine(); auto rad1 = seg.length(); auto c1 = Circle_( _ptS1, rad1 ); auto c2 = Circle_( _ptS2, rad1 ); auto int1 = li.intersects( c1 ); auto int2 = li.intersects( c2 ); assert( int1() ); assert( int2() ); auto ppt1 = int1.get(); auto ppt2 = int2.get(); return Segment_( ppt1.first, ppt2.second ); } //------------------------------------------------------------------ /// Distance from point to segment /** source: https://stackoverflow.com/a/6853926/193789 Temp implementation, until we get into this a bit more deeper */ template template HOMOG2D_INUMTYPE Segment_::distTo( const Point2d_& pt, int* segDistCase ) const { auto ppts = getPts(); auto x1 = static_cast( ppts.first.getX() ); auto y1 = static_cast( ppts.first.getY() ); auto x2 = static_cast( ppts.second.getX() ); auto y2 = static_cast( ppts.second.getY() ); auto A = pt.getX() - x1; auto B = pt.getY() - y1; auto C = x2 - x1; auto D = y2 - y1; auto dot = A * C + B * D; auto len_sq = C * C + D * D; auto param = dot / len_sq; HOMOG2D_INUMTYPE xx, yy; if( param < 0. ) { if( segDistCase ) *segDistCase = -1; xx = x1; yy = y1; } else { if( param > 1. ) { if( segDistCase ) *segDistCase = +1; xx = x2; yy = y2; } else { if( segDistCase ) *segDistCase = 0; xx = x1 + param * C; yy = y1 + param * D; } } auto dx = pt.getX() - xx; auto dy = pt.getY() - yy; return homog2d_sqrt( dx * dx + dy * dy ); } //------------------------------------------------------------------ /// Set circle from 2 points template template void Circle_::set( const Point2d_& pt1, const Point2d_& pt2 ) { #ifndef HOMOG2D_NOCHECKS if( pt1 == pt2 ) HOMOG2D_THROW_ERROR_1( "Unable, some points are identical" ); #endif Segment_ seg( pt1, pt2 ); _center = seg.getCenter(); _radius = seg.length() / 2.0; } //------------------------------------------------------------------ /// Set circle from 3 points /** Algorithm: we get the two largest segments and compute their bisector line. Their intersection point will be the center of circle, and the radius is the distance between center and any of the three points. We consider the the largest segments to improve numerical stability. Will throw if unable (numerical issue) One could think that the first checking would be enough, but experience shows that some times, the point are not colinear, but the two bisector lines are still parallel. Thus the second checking. \todo Check this other technique: https://www.johndcook.com/blog/2023/06/18/circle-through-three-points/ */ template template void Circle_::set( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& pt3 ) { #ifndef HOMOG2D_NOCHECKS if( areCollinear( pt1, pt2, pt3 ) ) HOMOG2D_THROW_ERROR_1( "Unable, points are colinear" ); #endif auto pts = priv::getLargestDistancePoints( pt1, pt2, pt3 ); auto seg1 = Segment_( pts[0], pts[1] ); auto seg2 = Segment_( pts[0], pts[2] ); auto li1 = seg1.getBisector(); auto li2 = seg2.getBisector(); #ifndef HOMOG2D_NOCHECKS if( li1.isParallelTo(li2) ) HOMOG2D_THROW_ERROR_1( "unable, bisector lines are parallel" ) #endif center() = li1 * li2; radius() = center().distTo(pt1); } //------------------------------------------------------------------ namespace priv { /// Free Function to check whether a circle encloses the given points template bool is_valid_circle( const Circle_& circ, const std::vector& pts ) { for( const auto& pt: pts ) if( !pt.isInside( circ ) ) return false; return true; } } // namespace priv //------------------------------------------------------------------ /// Helper function for Circle_::set( const T& ) /** Returns the minimum enclosing circle for N <= 3 */ template template Circle_ Circle_::p_min_circle_trivial( const std::vector& P ) const { if( P.empty() ) return Circle_(0.,0.,0.); // dummy circle of null radius if( P.size() == 1 ) return Circle_( P[0], 0. ); if( P.size() == 2 ) return Circle_( P[0], P[1] ); // To check if MEC can be determined by 2 points only for( int i=0; i<3; i++ ) for( int j=i+1; j<3; j++ ) { auto c = Circle_( P[i], P[j] ); if( priv::is_valid_circle( c, P ) ) return c; } return Circle_( P[0], P[1], P[2] ); // circle from 3 points } //------------------------------------------------------------------ /// Helper function for Circle_::set( const T& ), recursive /** Returns the MEC using Welzl's algorithm. Takes a set of input points P and a set R points on the circle boundary. \c T is the point type */ template template Circle_ Circle_::p_WelzlHelper( std::vector& P, ///< input set of points std::vector R, size_t n ///< Number of points in P that are not yet processed ) const { // Base case when all points processed or |R| = 3 if( n == 0 || R.size() == 3 ) return p_min_circle_trivial( R ); // Pick a random point randomly int idx = std::rand() % n; auto p = P[idx]; // Put the picked point at the end of P // since it's more efficient than // deleting from the middle of the vector std::swap( P[idx], P[n - 1] ); // Get the MEC circle d from the set of points P - {p} auto d = p_WelzlHelper( P, R, n - 1 ); if( p.isInside( d ) ) return d; // Otherwise, must be on the boundary of the MEC R.push_back( p ); // Return the MEC for P - {p} and R U {p} auto c = p_WelzlHelper( P, R, n - 1 ); return c; } //------------------------------------------------------------------ /// Compute circle from a set of points (Minimum Enclosing Circle, aka MEC) /// using the Welzl algorithm /** \c T may be std::vector, std::array or std::list holding points References: - https://en.wikipedia.org/wiki/Smallest-circle_problem - https://www.geeksforgeeks.org/minimum-enclosing-circle-using-welzls-algorithm/ */ template template< typename T, typename std::enable_if< trait::IsContainer::value ,T >::type* > void Circle_::set( const T& pts ) { if( pts.size() < 2 ) HOMOG2D_THROW_ERROR_1( "unable to build a circle from a single point" ); if( pts.size() == 2 ) { set( pts[0], pts[1] ); return; } if( pts.size() == 3 ) // todo: convert pts to HOMOG2D_INUMTYPE { set( pts[0], pts[1], pts[2] ); return; } std::vector> P_copy( std::begin(pts), std::end(pts) ); // std::random_shuffle( P_copy.begin(), P_copy.end() ); // ? check: what happens if removed? std::vector> R; h2d::thr::doNotCheckRadius() = true; auto cir = p_WelzlHelper( P_copy, R, P_copy.size() ); set( cir.center(), cir.radius() ); h2d::thr::doNotCheckRadius() = false; } //------------------------------------------------------------------ /// Returns true if circle is inside polyline /** Will be true if all these four conditions are met: - the polyline object must be of type "closed" AND a polygon (no intersection points) - center point is inside the polygon - no intersection points - any point of the polygon must be outside of the circle */ template template bool Circle_::isInside( const base::PolylineBase& poly ) const { HOMOG2D_START; if( !poly.isPolygon() ) return false; if( poly.getPts()[0].isInside(*this) ) // if a point of the polygon is inside the circle, return false; // then the circle cannot be inside the polygon if( !_center.isInside( poly ) ) return false; auto inters = intersects( poly ); return( inters() == 0 ); } //------------------------------------------------------------------ /// Circle/Circle intersection /** Ref: - https://stackoverflow.com/questions/3349125/ Can return 0, 1, or 2 intersection points \todo 20230219: in some situation, the difference below (x2 - x1) can be numericaly instable. Check if things would get improved by multiplying first (by a/d and h/d), before proceeding the difference. */ template template detail::IntersectM Circle_::intersects( const Circle_& other ) const { if( *this == other ) return detail::IntersectM(); HOMOG2D_INUMTYPE r1 = _radius; HOMOG2D_INUMTYPE r2 = other._radius; Point2d_ pt1 = _center; Point2d_ pt2 = other._center; HOMOG2D_INUMTYPE x1 = pt1.getX(); HOMOG2D_INUMTYPE y1 = pt1.getY(); HOMOG2D_INUMTYPE x2 = pt2.getX(); HOMOG2D_INUMTYPE y2 = pt2.getY(); HOMOG2D_INUMTYPE d_squared = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2); if( d_squared > r1*r1 + r2*r2 + (HOMOG2D_INUMTYPE)2.*r1*r2 ) // no intersection return detail::IntersectM(); if( d_squared < ( r1*r1 + r2*r2 - (HOMOG2D_INUMTYPE)2.*r1*r2 ) ) // no intersection: one circle inside the other return detail::IntersectM(); auto d = homog2d_sqrt( d_squared ); auto a = (r1*r1 - r2*r2 + d_squared) / (HOMOG2D_INUMTYPE)2. / d; auto h = homog2d_sqrt( r1*r1 - a*a ); Point2d_ P0( ( x2 - x1 ) * a / d + x1, ( y2 - y1 ) * a / d + y1 ); Point2d_ pt3( P0.getX() + ( y1 - y2 ) * h / d, P0.getY() - ( x1 - x2 ) * h / d ); Point2d_ pt4( P0.getX() - ( y1 - y2 ) * h / d, P0.getY() + ( x1 - x2 ) * h / d ); detail::IntersectM out; out.add( pt3 ); if( pt3 != pt4 ) out.add( pt4 ); return out; } namespace priv { //------------------------------------------------------------------ /// A POD value that needs some computing, associated with its flag. /// Used to be able to retain a value needing complex calculation. template class ValueFlag { private: T _value; bool _valIsCorrect = false; public: ValueFlag() = default; void set( T v ) { _value = v; _valIsCorrect = true; } T value() const { return _value; } void setBad() { _valIsCorrect = false; } bool isBad() const { return !_valIsCorrect; } }; //------------------------------------------------------------------ /// Holds attribute of a Polyline, allows storage of last computed value through the use of ValueFlag struct PolylineAttribs { priv::ValueFlag _length; priv::ValueFlag _area; priv::ValueFlag _isPolygon; priv::ValueFlag> _centroid; void setBad() { _length.setBad(); _area.setBad(); _isPolygon.setBad(); _centroid.setBad(); } }; //------------------------------------------------------------------ /// Returns the bounding box of points in vector/list/array of points \c vpts (free function) /** \todo This loops twice on the points. Maybe some improvement here. */ template< typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > PointPair1_ getBB_Points( const T& vpts ) { HOMOG2D_START; using FPT = typename T::value_type::FType; HOMOG2D_DEBUG_ASSERT( vpts.size(), "cannot run with no points" ); auto mm_x = std::minmax_element( std::begin( vpts ), std::end( vpts ), [] // lambda ( const Point2d_& pt1, const Point2d_& pt2 ) { return pt1.getX() < pt2.getX(); } ); auto mm_y = std::minmax_element( std::begin( vpts ), std::end( vpts ), [] // lambda ( const Point2d_& pt1, const Point2d_& pt2 ) { return pt1.getY() < pt2.getY(); } ); auto p1 = Point2d_( mm_x.first->getX(), mm_y.first->getY() ); auto p2 = Point2d_( mm_x.second->getX(), mm_y.second->getY() ); #ifndef HOMOG2D_NOCHECKS if( p1.distTo( p2 ) < thr::nullDistance() ) HOMOG2D_THROW_ERROR_1( "unable to compute bounding box of set, identical points:\n -p1:"<< p1 << "\n -p2:" << p2 ); if( shareCommonCoord( p1, p2 ) ) HOMOG2D_THROW_ERROR_1( "unable to compute bounding box of set, points share common coordinate:\n -p1:" << p1 << "\n -p2:" << p2 ); #endif // HOMOG2D_NOCHECKS return std::make_pair( p1, p2 ); } //------------------------------------------------------------------ /// Returns the bounding box of segments in vector/list/array of points \c vsegs template< typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > FRect_ getBB_Segments( const T& vsegs ) { HOMOG2D_START; using FPT = typename T::value_type::FType; HOMOG2D_DEBUG_ASSERT( vsegs.size(), "cannot compute bounding box of empty set of segments" ); std::vector> vpts( vsegs.size()*2 ); auto it = vpts.begin(); for( const auto& seg: vsegs ) { auto ppts = seg.getPts(); *it++ = ppts.first; *it++ = ppts.second; } return FRect_( getBB_Points( vpts ) ); } /// get BB for a set of FRect_ objects /// \todo same as getBB_Segments() ??? template auto getBB_FRect( const std::vector>& v_rects ) { HOMOG2D_START; HOMOG2D_DEBUG_ASSERT( v_rects.size(), "cannot compute bounding box of empty set of rectangles" ); std::vector> vpts( v_rects.size()*2 ); auto it = vpts.begin(); for( const auto& seg: v_rects ) { auto ppts = seg.getPts(); *it++ = ppts.first; *it++ = ppts.second; } return FRect_( getBB_Points( vpts ) ); } } // namespace priv // Forward declaration template CPolyline_ operator * ( const Homogr_&, const FRect_& ); // Forward declaration /// \todo 20230215: check if this has to be moved in the 'base' namespace. Can it be called **without** /// `using namespace std;`. This was the case for the streaming operator... template auto operator * ( const Homogr_&, const base::PolylineBase& ) -> base::PolylineBase; namespace base { //------------------------------------------------------------------ /// Polyline, will be instanciated either as \ref OPolyline_ (open polyline) or \ref CPolyline_ /** \warning When closed, in order to be able to compare two objects describing the same structure but potentially in different order, the comparison operator will proceed a sorting.
The consequence is that when adding points, if you have done a comparison before, you might not add point after the one you thought! template args: - PLT: PolyLine Type: type::IsClosed or type::IsOpen - FPT: Floating Point Type */ template class PolylineBase: public detail::Common #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; using PType = PLT; using detail::Common::isInside; template friend class PolylineBase; template friend class h2d::Ellipse_; template friend CPolyline_ h2d::operator * ( const h2d::Homogr_&, const h2d::FRect_& ); template friend auto h2d::operator * ( const Homogr_&, const base::PolylineBase& ) -> base::PolylineBase; Type type() const { return impl_Poly_type( detail::PlHelper() ); } private: Type impl_Poly_type( const detail::PlHelper& ) const { return Type::CPolyline; } Type impl_Poly_type( const detail::PlHelper& ) const { return Type::OPolyline; } private: mutable std::vector> _plinevec; mutable bool _plIsNormalized = false; mutable priv::PolylineAttribs _attribs; ///< Attributes. Will get stored upon computing. public: /// \name Constructors ///@{ /// Default constructor PolylineBase() = default; /// Constructor from FRect. Enabled for "closed" type, disabled for "open" type template PolylineBase( const FRect_& rect ) { imp_constrFRect( rect, detail::PlHelper() ); } template PolylineBase( const Segment_& seg ) { const auto& ppts = seg.getPts(); p_addPoint( ppts.first ); p_addPoint( ppts.second ); } /// Constructor from a vector of points. template< typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > PolylineBase( const T& vec ) { set( vec ); } /// Constructor that builds a regular convex polygon of \c n points at a distance \c rad, centered at (0,0). template PolylineBase( FPT2 rad, T n ) { impl_set_RCP( rad, n, detail::PlHelper() ); } /// Copy-Constructor from Closed Polyline template constexpr PolylineBase( const CPolyline_& other ) { static_assert( !(std::is_same::value), "Error, cannot build an Open Polyline from a closed one" ); set( other._plinevec ); } /// Copy-Constructor from Open Polyline template PolylineBase( const OPolyline_& other ) { set( other._plinevec ); } #ifdef HOMOG2D_USE_BOOSTGEOM /// Constructor from a boost geometry polygon, see misc/test_files/bg_test_1.cpp /** \note Only imports the "outer" envelope, homog2d does not handle polygons with holes Requirements: the Boost polygon must have 2-cartesian coordinates points, either `bg::model::point` or `bg::model::d2::point_xy` but the underlying numerical type is free. \note At present, the 3th template parameter (bool) (Closed or Open) is ignored, because it is unclear how this relates to the actual fact that last point is equal to first point. \todo 20230216: add some checking that the type BPT needs to fit certain requirements (must have 2-dimensions, and use cartesian coordinates). Maybe we should add some Sfinae to check this. */ template // PolylineBase( const boost::geometry::model::polygon< BPT, ///< Boost Point Type (either `bg::model::point` or `bg::model::d2::point_xy`) CLKW, ///< this one is ignored here CLOSED ///< true: closed, false: open >& bgpol ///< input boost geometry polygon ) { const auto& outer = bgpol.outer(); bool isClosed = false; auto ptf = outer.front(); auto ptb = outer.back(); Point2d_ pt_front( boost::geometry::get<0>(ptf), boost::geometry::get<1>(ptf) ); Point2d_ pt_back( boost::geometry::get<0>(ptb), boost::geometry::get<1>(ptb) ); /* this should work !!! (but doesn't...) Point2d_ pt_front( outer.front() ); Point2d_ pt_back( outer.back() ); */ if( pt_front == pt_back ) // means it's closed isClosed = true; if( isClosed && std::is_same::value ) // cannot build an open polyline from a closed one HOMOG2D_THROW_ERROR_1( "unable to convert a closed boost::polygon into an OPolyline" ); _plinevec.reserve( outer.size() - isClosed ); for( auto it=outer.begin(); it!=outer.end()-isClosed; it++ ) { const auto& bgpt = *it; _plinevec.emplace_back( Point2d_( boost::geometry::get<0>(bgpt), boost::geometry::get<1>(bgpt) ) ); } } #endif ///@} private: template std::pair impl_set_RCP( FPT2 rad, T n, const typename detail::PlHelper& ); template constexpr std::pair impl_set_RCP( FPT2, T, const typename detail::PlHelper& ) { static_assert( detail::AlwaysFalse::value, "cannot build an regular convex polygon for a OPolyline object" ); return std::make_pair(0.,0.); // to avoid a compiler warning } template void imp_constrFRect( const FRect_& rect, const detail::PlHelper& ) { for( const auto& pt: rect.get4Pts() ) _plinevec.push_back( pt ); } template constexpr void imp_constrFRect( const FRect_&, const detail::PlHelper& ) { static_assert( detail::AlwaysFalse::value, "cannot build an OPolyline object from a FRect"); } public: /// \name Attributes access ///@{ /// Returns the number of points size_t size() const { return _plinevec.size(); } constexpr bool isClosed() const { return impl_isClosed( detail::PlHelper()); } HOMOG2D_INUMTYPE length() const; HOMOG2D_INUMTYPE area() const; bool isPolygon() const; /// Returns Bounding Box of Polyline FRect_ getBB() const { HOMOG2D_START; #ifndef HOMOG2D_NOCHECKS if( size() < 2 ) HOMOG2D_THROW_ERROR_1( "cannot compute bounding box of empty Polyline" ); #endif auto ppts = priv::getBB_Points( getPts() ); #ifndef HOMOG2D_NOCHECKS if( shareCommonCoord( ppts.first, ppts.second ) ) HOMOG2D_THROW_ERROR_1( "unable, points share common coordinate" ); #endif return FRect_( ppts ); } LPBase centroid() const; /// Returns the number of segments. If "closed", /** the last segment (going from last to first point) is counted */ size_t nbSegs() const { if( size() < 2 ) return 0; return impl_nbSegs( detail::PlHelper() ); } Point2d_ getExtremePoint( CardDir ) const; Point2d_ getTmPoint() const; Point2d_ getBmPoint() const; Point2d_ getLmPoint() const; Point2d_ getRmPoint() const; bool isConvex() const; ///@} private: HOMOG2D_INUMTYPE p_ComputeSignedArea() const; constexpr bool impl_isClosed( const detail::PlHelper& ) const { return false; } constexpr bool impl_isClosed( const detail::PlHelper& ) const { return true; } size_t impl_nbSegs( const detail::PlHelper& ) const { return size() - 1; } size_t impl_nbSegs( const detail::PlHelper& ) const { return size(); } constexpr bool impl_isPolygon( const detail::PlHelper& ) const; bool impl_isPolygon( const detail::PlHelper& ) const; /// Add single point. private, because only to be used from other member functions /** \warning This function was discarded from public API in dec. 2021, because this will add the new point after the previous one \b only if the object has \b not been normalized. This normalizing operation will happen if you do a comparison (== or !=) */ template void p_addPoint( const Point2d_& pt ) { #ifndef HOMOG2D_NOCHECKS if( size() ) if( pt == _plinevec.back() ) HOMOG2D_THROW_ERROR_1( "cannot add a point identical to previous one: pt=" << pt << " size=" << size() ); #endif _attribs.setBad(); _plIsNormalized=false; _plinevec.push_back( pt ); } public: /// \name Data access ///@{ /// Returns the points (reference) std::vector>& getPts() { _attribs.setBad(); // because we send the non-const reference _plIsNormalized=false; return _plinevec; } /// Returns the points (const reference) const std::vector>& getPts() const { return _plinevec; } /// Returns (as a copy) the segments of the polyline std::vector> getSegs() const { if( size() < 2 ) // nothing to return return std::vector>(); std::vector> out; out.reserve( size() ); for( size_t i=0; i(pt1,pt2) ); } impl_getSegs( out, detail::PlHelper() ); return out; } /// Returns one point of the polyline. Point2d_ getPoint( size_t idx ) const { #ifndef HOMOG2D_NOCHECKS if( idx >= size() ) HOMOG2D_THROW_ERROR_1( "requesting point " << idx << ", only has " << size() ); #endif return _plinevec[idx]; } /// Returns one segment of the polyline. /** Segment \c n is the one between point \c n and point \c n+1 */ Segment_ getSegment( size_t idx ) const { #ifndef HOMOG2D_NOCHECKS if( idx >= nbSegs() ) HOMOG2D_THROW_ERROR_1( "requesting segment " << idx << ", only has " << nbSegs() ); if( size() < 2 ) HOMOG2D_THROW_ERROR_1( "no segment " << idx ); #endif return impl_getSegment( idx, detail::PlHelper() ); } CPolyline_ convexHull() const; ///@} private: /// empty implementation void impl_getSegs( std::vector>&, const detail::PlHelper& ) const {} /// that one is for closed Polyline, adds the last segment void impl_getSegs( std::vector>& out, const detail::PlHelper& ) const { out.push_back( Segment_(_plinevec.front(),_plinevec.back() ) ); } Segment_ impl_getSegment( size_t idx, const detail::PlHelper& ) const { return Segment_( _plinevec[idx], _plinevec[idx+1==nbSegs()?0:idx+1] ); } Segment_ impl_getSegment( size_t idx, const detail::PlHelper& ) const { return Segment_( _plinevec[idx], _plinevec[idx+1] ); } public: /// \name Modifiers (non-const member functions) ///@{ /// Clear all (does not change the "open/close" status). void clear() { _plinevec.clear(); _plIsNormalized=false; _attribs.setBad(); } template void rotate( Rotate, const Point2d_& ); void rotate( Rotate ); /// Miminize the PolylineBase: remove all points that lie in the middle of two segments with same angle. /** For example, if we have the following points ("Open" polyline): (0,0)-(1,0)-(2,0)
This will be replaced by: (0,0)--(2,0) \note We cannot use the Segment_::getAngle() function because it returns a value in [0,PI/2], so we would'nt be able to detect a segment going at 180° of the previous one. \todo 20230217: implement these: - https://en.wikipedia.org/wiki/Visvalingam%E2%80%93Whyatt_algorithm - https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm Also use the areCollinear() function */ void minimize() { if( size()<3 ) return; impl_minimizePL( detail::PlHelper() ); } /// Translate Polyline using \c dx, \c dy template void translate( TX dx, TY dy ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); for( auto& pt: _plinevec ) pt.translate( dx, dy ); } /// Translate Polyline, using a pair of numerical values template void translate( const std::pair& ppt ) { translate( ppt.first, ppt.second ); } /// Move Polyline to new origin template void moveTo( TX x, TY y ) { HOMOG2D_CHECK_IS_NUMBER( TX ); HOMOG2D_CHECK_IS_NUMBER( TY ); moveTo( Point2d_(x,y) ); } /// Move Polyline to new origin, given by \c new_org /** \warning The polygon origin is NOT its center but the point No 0! */ template void moveTo( const Point2d_& new_org ) { if( size() == 0 ) HOMOG2D_THROW_ERROR_1( "Invalid call, Polyline is empty" ); auto dx = new_org.getX() - getPoint(0).getX(); auto dy = new_org.getY() - getPoint(0).getY(); for( auto& pt: _plinevec ) pt.translate( dx, dy ); } /// Set from vector/array/list of points (discards previous points) /** - nb of elements must be 0 or 2 or more */ template< typename CONT, typename std::enable_if< trait::IsContainer::value, CONT >::type* = nullptr > void set( const CONT& vec ) { #ifndef HOMOG2D_NOCHECKS if( vec.size() == 1 ) HOMOG2D_THROW_ERROR_1( "Invalid: number of points must be 0, 2 or more" ); if( vec.size() > 1 ) p_checkInputData( vec ); #endif _attribs.setBad(); _plIsNormalized=false; _plinevec.resize( vec.size() ); auto it = std::begin( _plinevec ); for( const auto& pt: vec ) // copying one by one will *it++ = pt; // allow type conversions (std::copy implies same type) } /// Build a parallelogram (4 points) from 3 points template void setParallelogram( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& pt3 ) { impl_setParallelogram( pt1, pt2, pt3, detail::PlHelper() ); } /// Set from FRect template void set( const FRect_& rec ) { impl_setFromFRect( rec, detail::PlHelper() ); } /// Build RCP (Regular Convex Polygon), and return distance between consecutive points template std::pair set( FPT2 rad, T n ) { return impl_set_RCP( rad, n, detail::PlHelper() ); } ///@} private: template void impl_setFromFRect( const FRect_&, const detail::PlHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid: cannot set a OPolyline from a FRect" ); } template void impl_setFromFRect( const FRect_& rect, const detail::PlHelper& ) { CPolyline_ tmp(rect); std::swap( *this, tmp ); } template void impl_setParallelogram( const Point2d_&, const Point2d_&, const Point2d_&, const detail::PlHelper& ) { static_assert( detail::AlwaysFalse::value, "Invalid: cannot set a OPolyline as a parallelogram" ); } template void impl_setParallelogram( const Point2d_& p1, const Point2d_& p2, const Point2d_& p3, const detail::PlHelper& ) { Point2d_ pt1(p1); Point2d_ pt2(p2); Point2d_ pt3(p3); std::vector> vpts(4); vpts[0] = pt1; vpts[1] = pt2; vpts[2] = pt3; auto li_21 = pt1 * pt2; auto li_23 = pt3 * pt2; auto li_34 = li_21.getParallelLine( pt3 ); auto li_14 = li_23.getParallelLine( pt1 ); vpts[3] = li_34 * li_14; set( vpts ); } /// Checks that no contiguous identical points are stored template void p_checkInputData( const CONT& pts ) { for( auto it=pts.begin(); it!=std::prev(pts.end()); it++ ) { const auto& elem = *it; const auto& next = *std::next(it); if( elem == next ) HOMOG2D_THROW_ERROR_1( "cannot add two consecutive identical points:\npt:" << elem << " and pt:" << next << " in set of size " << pts.size() ); } } void impl_minimizePL( const detail::PlHelper& ); void impl_minimizePL( const detail::PlHelper& ); void p_minimizePL( PolylineBase&, size_t istart, size_t iend ); public: /// \name Operators ///@{ template bool operator == ( const PolylineBase& other ) const { if( size() != other.size() ) // for quick exit return false; return impl_operatorComp( other, detail::PlHelper() ); } template bool operator != ( const PolylineBase& other ) const { return !( *this == other ); } ///@} private: /// This one is guaranteed to operate on same 'PLT' types, is called by the four others. /** Assumes the size is the same */ template bool impl_operatorComp_root( const PolylineBase& other ) const { p_normalizePoly(); other.p_normalizePoly(); auto it = other._plinevec.begin(); for( const auto& elem: _plinevec ) if( *it++ != elem ) return false; return true; } template bool impl_operatorComp( const OPolyline_& other, const detail::PlHelper& ) const { return impl_operatorComp_root( other ); } template bool impl_operatorComp( const CPolyline_& other, const detail::PlHelper& ) const { return impl_operatorComp_root( other ); } /// Comparing open and closed polyline objects returns always \c false template constexpr bool impl_operatorComp( const OPolyline_&, const detail::PlHelper& ) const { return false; } /// Comparing open and closed polyline objects returns always \c false template constexpr bool impl_operatorComp( const CPolyline_&, const detail::PlHelper& ) const { return false; } public: /// Polyline intersection with Line, Segment, FRect, Circle template< typename T, typename std::enable_if< trait::IsShape::value, T >::type* = nullptr > detail::IntersectM intersects( const T& other ) const { detail::IntersectM out; for( const auto& pseg: getSegs() ) { auto inters = pseg.intersects( other ); if( inters() ) out.add( inters.get() ); } return out; } /// Polyline isInside other primitive /** Sfinae should resolve for T=Circle,FRect,Ellipse but not for CPolyline */ template< typename T, typename std::enable_if< ( trait::HasArea::value && !trait::PolIsClosed::value ) ,T >::type* = nullptr > bool isInside( const T& prim ) const { HOMOG2D_START; if( size() == 0 ) return false; for( const auto& pt: getPts() ) if( !pt.isInside( prim ) ) return false; return true; } /// Polyline isInside other CPolyline /** Sfinae should resolve ONLY for CPolyline */ template< typename T, typename std::enable_if< trait::PolIsClosed::value,T >::type* = nullptr > bool isInside( const T& cpol ) const { HOMOG2D_START; if( size() == 0 ) return false; for( const auto& pt: getPts() ) if( !pt.isInside( cpol ) ) return false; if( intersects(cpol)() ) return false; return true; } template friend std::ostream& h2d::base::operator << ( std::ostream&, const h2d::base::PolylineBase& ); public: #ifdef HOMOG2D_TEST_MODE /// this is only needed for testing bool isNormalized() const { return _plIsNormalized; } #endif private: void impl_normalizePoly( const detail::PlHelper& ) const { auto minpos = std::min_element( _plinevec.begin(), _plinevec.end() ); std::rotate( _plinevec.begin(), minpos, _plinevec.end() ); const auto& p1 = _plinevec[1]; const auto& p2 = _plinevec.back(); if( p2 < p1 ) { std::reverse( _plinevec.begin(), _plinevec.end() ); minpos = std::min_element( _plinevec.begin(), _plinevec.end() ); std::rotate( _plinevec.begin(), minpos, _plinevec.end() ); } } void impl_normalizePoly( const detail::PlHelper& ) const { if( _plinevec.back() < _plinevec.front() ) std::reverse( _plinevec.begin(), _plinevec.end() ); } /// Normalization of CPolyline_ /** Two tasks: - rotating so that the smallest one is first - reverse if needed, so that the second point is smaller than the last one */ void p_normalizePoly() const { if( size() == 0 ) return; if( !_plIsNormalized ) { impl_normalizePoly( detail::PlHelper() ); _plIsNormalized=true; } } public: void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; #ifdef HOMOG2D_USE_OPENCV void draw( img::Image&, img::DrawParams dp=img::DrawParams() ) const; #endif }; // class PolylineBase /// Build a Regular Convex Polygon of radius \c rad with \c n points, centered at (0,0) /** \return segment distance, inscribed circle radius */ template template std::pair PolylineBase::impl_set_RCP( FPT2 rad, ///< Radius T ni, ///< Nb of points const typename detail::PlHelper& ) { static_assert( std::is_integral::value, "2nd argument type must be integral type" ); if( ni < 3 ) HOMOG2D_THROW_ERROR_1( "unable, nb of points must be > 2" ); if( rad <= 0 ) HOMOG2D_THROW_ERROR_1( "unable, radius must be >= 0" ); size_t n(ni); std::vector> v_pts(n); auto it = std::begin( v_pts ); auto angle0 = (HOMOG2D_INUMTYPE)2. * M_PI / n; Point2d_ pt0( rad, 0. ); // initial point HOMOG2D_INUMTYPE radius = 0.; for( size_t i=0; i( x * rad, y * rad ); Segment_ seg( pt0, pt1 ); Line2d_ li( seg.getCenter() ); // line passing through (0,0) and middle point of segment auto inters_pt = seg.intersects( li ) ; // intersection point assert( inters_pt() ); radius = Point2d_(0,0).distTo( inters_pt.get() ); } *it = Point2d_( x * rad, y * rad ); it++; } *this = PolylineBase( v_pts ); return std::make_pair( getPoint(0).distTo( getPoint(1) ), radius ); } } // namespace base // FREE FUNCTIONS // (need to be here because the 5 below are being called by some member functions // of class PolylineBase, that is implemented further below namespace priv { //------------------------------------------------------------------ /// Return iterator on Bottom-most point of container holding points /** Used by - getBmPoint() - priv::chull::getPivotPoint() */ template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > auto getBmPoint_helper( const T& t ) { using FPT = typename T::value_type::FType; #ifndef HOMOG2D_NOCHECKS if( t.size() == 0 ) HOMOG2D_THROW_ERROR_1( "invalid call, container is empty" ); #endif return std::min_element( std::begin(t), std::end(t), [] // lambda ( const Point2d_& pt1, const Point2d_& pt2 ) { if( pt1.getY() < pt2.getY() ) return true; if( pt1.getY() > pt2.getY() ) return false; return( pt1.getX() < pt2.getX() ); } ); } } // namespace priv //------------------------------------------------------------------ /// Return Bottom-most point of container holding points template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > Point2d_ getBmPoint( const T& t ) { #ifndef HOMOG2D_NOCHECKS if( t.size() == 0 ) HOMOG2D_THROW_ERROR_1( "invalid call, container is empty" ); #endif return *priv::getBmPoint_helper( t ); } //------------------------------------------------------------------ /// Return Top-most point of container template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > Point2d_ getTmPoint( const T& t ) { using FPT = typename T::value_type::FType; #ifndef HOMOG2D_NOCHECKS if( t.size() == 0 ) HOMOG2D_THROW_ERROR_1( "invalid call, container is empty" ); #endif return *std::min_element( std::begin(t), std::end(t), [] // lambda ( const Point2d_& pt1, const Point2d_& pt2 ) { if( pt1.getY() > pt2.getY() ) return true; if( pt1.getY() < pt2.getY() ) return false; return( pt1.getX() < pt2.getX() ); } ); } //------------------------------------------------------------------ /// Return Left-most point of container template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > Point2d_ getLmPoint( const T& t ) { using FPT = typename T::value_type::FType; #ifndef HOMOG2D_NOCHECKS if( t.size() == 0 ) HOMOG2D_THROW_ERROR_1( "invalid call, container is empty" ); #endif return *std::min_element( std::begin(t), std::end(t), [] // lambda ( const Point2d_& pt1, const Point2d_& pt2 ) { if( pt1.getX() < pt2.getX() ) return true; if( pt1.getX() > pt2.getX() ) return false; return( pt1.getY() < pt2.getY() ); } ); } //------------------------------------------------------------------ /// Return Right-most point of container holding points template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > Point2d_ getRmPoint( const T& t ) { using FPT = typename T::value_type::FType; #ifndef HOMOG2D_NOCHECKS if( t.size() == 0 ) HOMOG2D_THROW_ERROR_1( "invalid call, container is empty" ); #endif return *std::min_element( std::begin(t), std::end(t), [] // lambda ( const Point2d_& pt1, const Point2d_& pt2 ) { if( pt1.getX() > pt2.getX() ) return true; if( pt1.getX() < pt2.getX() ) return false; return( pt1.getY() < pt2.getY() ); } ); } //------------------------------------------------------------------ /// Return Bottom-most point of Polyline (free function) template Point2d_ getBmPoint( const base::PolylineBase& poly ) { return getBmPoint( poly.getPts() ); } /// Return Top-most point of Polyline (free function) template Point2d_ getTmPoint( const base::PolylineBase& poly ) { return getTmPoint( poly.getPts() ); } /// Return Left-most point of Polyline (free function) template Point2d_ getLmPoint( const base::PolylineBase& poly ) { return getLmPoint( poly.getPts() ); } /// Return Right-most point of Polyline (free function) template Point2d_ getRmPoint( const base::PolylineBase& poly ) { return getRmPoint( poly.getPts() ); } //------------------------------------------------------------------ /// Get Top-most / Bottom-most / Left-most / Right-most point, depending on \c dir /// (free function) /** Type \c T can be either a Polyline (open or closed), or a container holding points, or a FRect_ */ template Point2d_ getExtremePoint( CardDir dir, const T& t ) { switch( dir ) { case CardDir::Top: return getTmPoint(t); break; case CardDir::Bottom: return getBmPoint(t); break; case CardDir::Left: return getLmPoint(t); break; case CardDir::Right: return getRmPoint(t); break; default: assert(0); } } namespace base { //------------------------------------------------------------------ /// Get Top-most / Bottom-most / Left-most / Right-most point template Point2d_ PolylineBase::getExtremePoint( CardDir dir ) const { return getExtremePoint( dir, *this ); } /// Return Bottom-most point of Polyline template Point2d_ PolylineBase::getBmPoint() const { return h2d::getBmPoint( *this ); } /// Return Top-most point of Polyline template Point2d_ PolylineBase::getTmPoint() const { return h2d::getTmPoint( *this ); } /// Return Left-most point of Polyline template Point2d_ PolylineBase::getLmPoint() const { return h2d::getLmPoint( *this ); } /// Return Right-most point of Polyline template Point2d_ PolylineBase::getRmPoint() const { return h2d::getRmPoint( *this ); } //------------------------------------------------------------------ /// Rotates the polyline by either 90°, 180°, 270° (-90°) at point \c refpt /** For an arbitrary angle \c alpha (rad.), you can write: ``` auto poly2 = Homogr(alpha) * poly; ``` \sa PolylineBase::rotate( Rotate ) */ template template void PolylineBase::rotate( Rotate rot, const Point2d_& refpt ) { translate( -refpt.getX(), -refpt.getY() ); this->rotate( rot ); translate( refpt.getX(), refpt.getY() ); } //------------------------------------------------------------------ /// Rotates the polyline by either 90°, 180°, 270° (-90°) at point (0,0) /** For an arbitrary angle \c alpha (rad.), you can write: ``` auto poly2 = Homogr(alpha) * poly; ``` \sa PolylineBase::rotate( Rotate, const Point2d_& ) */ template void PolylineBase::rotate( Rotate rot ) { switch( rot ) { case Rotate::CCW: for( auto& pt: getPts() ) pt.set( -pt.getY(), +pt.getX() ); break; case Rotate::CW: for( auto& pt: getPts() ) pt.set( +pt.getY(), -pt.getX() ); break; case Rotate::Full: for( auto& pt: getPts() ) pt.set( -pt.getX(), -pt.getY() ); break; case Rotate::VMirror: for( auto& pt: getPts() ) pt.set( -pt.getX(), pt.getY() ); break; case Rotate::HMirror: for( auto& pt: getPts() ) pt.set( pt.getX(), -pt.getY() ); break; default: assert(0); } } //------------------------------------------------------------------ /// Private member function, called by PolylineBase::impl_minimizePL() template void PolylineBase::p_minimizePL( PolylineBase& pl, size_t istart, size_t iend ) { auto nbpts = pl.size(); // step 1: check each point to see if it is the middle point of two segments with same angle std::vector ptset; for( size_t i=istart; i out; size_t vec_idx = 0; for( size_t i=0; i void PolylineBase::impl_minimizePL( const detail::PlHelper& ) { assert( size() > 2 ); p_minimizePL( *this, 1, size()-1 ); } template void PolylineBase::impl_minimizePL( const detail::PlHelper& ) { p_minimizePL( *this, 0, size() ); } //------------------------------------------------------------------ /// Returns true if object is a polygon (closed, and no segment crossing) template bool PolylineBase::isPolygon() const { if( size()<3 ) // needs at least 3 points to be a polygon return false; return impl_isPolygon( detail::PlHelper() ); } /// If open, then not a polygon template constexpr bool PolylineBase::impl_isPolygon( const detail::PlHelper& ) const { return false; } /// If closed, we need to check for crossings template bool PolylineBase::impl_isPolygon( const detail::PlHelper& ) const { if( _attribs._isPolygon.isBad() ) { auto nbs = nbSegs(); size_t i=0; bool notDone = true; bool hasIntersections = false; do { auto seg1 = getSegment(i); auto lastone = i==0?nbs-1:nbs; for( auto j=i+2; j bool PolylineBase::isConvex() const { if( !isPolygon() ) return false; int8_t sign = 0; const auto& vpts = getPts(); if( vpts.size() == 3 ) // 3 pts => always convex return true; for( size_t i=0; i0 ? +1 : -1); else if (sign != (crossproduct>0 ? +1 : -1) ) return false; } return true; } //------------------------------------------------------------------ /// Returns length of Polyline template HOMOG2D_INUMTYPE PolylineBase::length() const { if( _attribs._length.isBad() ) { HOMOG2D_INUMTYPE sum = 0.; for( const auto& seg: getSegs() ) sum += static_cast( seg.length() ); _attribs._length.set( sum ); } return _attribs._length.value(); } //------------------------------------------------------------------ /// Returns area of polygon (computed only if necessary) template HOMOG2D_INUMTYPE PolylineBase::area() const { if( !isPolygon() ) // implies that is both closed and has no intersections return 0.; if( _attribs._area.isBad() ) _attribs._area.set( homog2d_abs( p_ComputeSignedArea() ) ); return _attribs._area.value(); } //------------------------------------------------------------------ /// Compute and returns signed area (used in area() and in centroid() ) template HOMOG2D_INUMTYPE PolylineBase::p_ComputeSignedArea() const { HOMOG2D_INUMTYPE area = 0.; for( size_t i=0; i(pt1.getX()) * pt2.getY(); area -= static_cast(pt1.getY()) * pt2.getX(); } return area / 2.; } //------------------------------------------------------------------ /// Compute centroid of polygon /** ref: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon */ template Point2d_ base::PolylineBase::centroid() const { if( !isPolygon() ) HOMOG2D_THROW_ERROR_2( "PolylineBase::centroid", "unable, is not a polygon" ); if( _attribs._centroid.isBad() ) { HOMOG2D_INUMTYPE cx = 0.; HOMOG2D_INUMTYPE cy = 0.; for( size_t i=0; i( cx, cy ); _attribs._centroid.set( c ); } return _attribs._centroid.value(); } } // namespace base //------------------------------------------------------------------ /// Rotates the rectangle by either 90°, 180°, 270° (-90°) at point \c refpt /** For an arbitrary angle \c alpha (rad.), you can write: ``` auto r2 = Homogr(alpha) * rect; ``` \sa FRect_::rotate( Rotate ) \sa h2d::rotate() */ template template void FRect_::rotate( Rotate rot, const Point2d_& refpt ) { translate( -refpt.getX(), -refpt.getY() ); this->rotate( rot ); translate( refpt.getX(), refpt.getY() ); } //------------------------------------------------------------------ /// Rotates the rectangle by either 90°, 180°, 270° (-90°) at point (0,0) /** For an arbitrary angle \c alpha (rad.), you can write: ``` auto r2 = Homogr(alpha) * rect; ``` \sa FRect_::rotate( Rotate, const Point2d_& ) \sa h2d::rotate() \note This function converts the FRect to a CPolyline_, rotates it, and builds the output rectangle by getting the bounding box. This is because we cannot directly rotate the points, because the two points of a FRect_ are switched so that the smallest point is always in _ptR1. */ template void FRect_::rotate( Rotate rot ) { CPolyline_ pol( *this ); pol.rotate( rot ); *this = pol.getBB(); } //------------------------------------------------------------------ /// Returns Rectangle of the intersection area of two rectangles /** \return An object of type detail::RectArea, that can be checked for success using the () operator, and if success, will hold the resulting FRect_ Algorithm: 4 situations need to be considered, depending on the number of intersection points: - A: 2 intersection points on same segment => we have 2 points inside. \verbatim +------+ +------+ | | | | | +--+-----+ +--+----+ | | | | | | | | | | | | | or | | | | or ... | +--+-----+ +--+----+ | | | | | +------+ +------+ \endverbatim - B: 2 intersection points on different segments => for each rectangle, 1 point is inside the other. \verbatim +------+ | | | | | +--+-----+ | | | | or ... +---+--+ | | | +--------+ \endverbatim - C: 3 points \verbatim +------+-----+ | | | | | | +------+-----+ | | +------+ \endverbatim - D: 4 intersection points: the intersection rectangle is made of these 4 points \verbatim +------+ | | +--+------+-----+ | | | | | | | | +--+------+-----+ | | +------+ \endverbatim */ template template detail::RectArea FRect_::intersectArea( const FRect_& other ) const { if( *this == other ) // if same rectangles, then return detail::RectArea(other); // the intersection area is the rectangle auto inter = this->intersects( other ); if( !inter() ) // rectangles do not intersect { if( this->isInside( other ) ) return detail::RectArea(*this); if( other.isInside( *this ) ) return detail::RectArea(other); return detail::RectArea(); } if( inter.size() < 2 ) // only one intersection point return detail::RectArea(); // => no intersection area! if( inter.size() == 4 ) // 4 intersection points => case "D" return detail::RectArea(FRect_( inter.get().at(0), inter.get().at(3) ) ); if( inter.size() == 3 ) // 3 intersection points => case "C" { auto v = inter.get(); auto xmin = std::min( v[0].getX(), std::min( v[1].getX(),v[2].getX() ) ); auto ymin = std::min( v[0].getY(), std::min( v[1].getY(),v[2].getY() ) ); auto xmax = std::max( v[0].getX(), std::max( v[1].getX(),v[2].getX() ) ); auto ymax = std::max( v[0].getY(), std::max( v[1].getY(),v[2].getY() ) ); #ifndef HOMOG2D_DEBUGMODE assert( xmax-xmin > thr::nullOrthogDistance() ); assert( ymax-ymin > thr::nullOrthogDistance() ); #else HOMOG2D_DEBUG_ASSERT( ( ( xmax-xmin > thr::nullOrthogDistance() ) && ( ymax-ymin > thr::nullOrthogDistance() )) , std::scientific << "this=" << *this << " other=" << other << "\nxmax=" << xmax << " xmin=" << xmin << "\nymax=" << ymax << " ymin=" << ymin << "\nnod=" << thr::nullOrthogDistance() ); #endif return detail::RectArea( FRect_( xmin, ymin, xmax,ymax ) ); } /*------------------------ If we get to this point, we have 2 intersection points ------------------------*/ #ifndef HOMOG2D_DEBUGMODE assert( inter.size() == 2 ); #else /* HOMOG2D_DEBUG_ASSERT( inter.size() == 2, "inter.size()=" << inter.size() << "\n-this=" << *this << "\n-other=" << other <<"\n-inter:" << inter );*/ #endif const auto& r1 = *this; const auto& r2 = other; // HOMOG2D_LOG( "r1="<::unionArea() */ struct Cell { bool _isCorner = false; Cell() = default; #ifdef HOMOG2D_DEBUGMODE Cell( bool b ) : _isCorner(b) {} #endif template Cell( const Index& ix, const Index& iy ) { if( ix._rect_idx == iy._rect_idx ) _isCorner = true; } }; using Table = std::array,4>; using PCoord = std::pair; enum class Direction: uint8_t { N, E, S, W }; /* NEEDED TO DEBUG const char* getString( Direction dir ) { switch( dir ) { case Direction::W: return "W"; break; case Direction::S: return "S"; break; case Direction::E: return "E"; break; case Direction::N: return "N"; break; } return 0; // to avoid a warning }*/ enum class Turn: uint8_t { Left, Right }; inline Direction turn( Direction dir, Turn turn ) { switch( dir ) { case Direction::W: return turn == Turn::Left ? Direction::S : Direction::N; case Direction::S: return turn == Turn::Left ? Direction::E : Direction::W; case Direction::E: return turn == Turn::Left ? Direction::N : Direction::S; case Direction::N: return turn == Turn::Left ? Direction::W : Direction::E; } return Direction::N; // to avoid a warning } inline void moveToNextCell( uint8_t& row, uint8_t& col, const Direction& dir ) { switch( dir ) { case Direction::N: row--; break; case Direction::S: row++; break; case Direction::E: col++; break; case Direction::W: col--; break; } } /// Helper function for FRect_::unionArea() /** - Start from 0,0, direction East */ inline std::vector parseTable( Table& table ) { bool firstTime = true; bool done = false; uint8_t row = 0; uint8_t col = 0; Direction dir = Direction::E; std::vector out; do { if( table[row][col]._isCorner ) { auto new_pair = std::make_pair(row,col); if( out.size() > 0 ) if( new_pair == out.front() && out.size() > 2 ) done = true; if( !done ) out.push_back( new_pair ); if( firstTime ) firstTime = false; else dir = turn( dir, Turn::Right ); } else { if( ( row==2 || row==1 ) && ( col==2 || col==1 ) ) { out.push_back( std::make_pair(row,col) ); dir = turn( dir, Turn::Left ); } } moveToNextCell( row, col, dir ); } while( !done ); return out; } #ifdef HOMOG2D_DEBUGMODE inline void printTable( const Table& t, std::string msg ) { std::cout << "Table: " << msg << "\n | "; for( uint8_t r=0;r<4; r++ ) std::cout << (int)r << " "; std::cout << "\n--|---------|\n"; for( uint8_t r=0;r<4; r++ ) { std::cout << (int)r << " | "; for( uint8_t c=0;c<4; c++ ) std::cout << (t[r][c]._isCorner?'F':'.') << " "; std::cout << "|\n"; } } #endif //------------------------------------------------------------------ /// Helper function, used in FRect_::unionArea() template inline CPolyline_ convertToCoord( const std::vector& v_coord, ///< vector of coordinate indexes const std::array,4>& v_x, ///< holds x-coordinates const std::array,4>& v_y ///< holds y-coordinates ) { std::vector> v_pts; for( const auto& elem: v_coord ) { auto id_x = elem.first; auto id_y = elem.second; assert( id_x<4 && id_y<4 ); auto pt = Point2d_( v_x[id_x]._value, v_y[id_y]._value ); if( v_pts.empty() ) v_pts.push_back( pt ); else // add to vector only if not same as previous { // and not same as first if( v_pts.back() != pt && v_pts.front() != pt ) v_pts.push_back( pt ); } } // printVector( v_pts, "polyline" ); return CPolyline_( v_pts );//, IsClosed::Yes ); } } // namespace priv } // namespace runion //------------------------------------------------------------------ /// Computes the CPolyline of the union of two rectangles /** Algorithm: -# build vectors of x and y coordinates (4 elements) -# build table x-y (4x4), with corners tagged -# parse the table by turning right at each corner, and left if position is not one the outside row/col -# convert back indexes to real coordinates, to build the final CPolyline object Two examples: \verbatim 9 +----+ +-------+ | | | | 8 +---+----+---+ +---+---+ | | | | | | | | | 7 +---+----+---+ +---+---+ | | | | | 6 +----+ +-------+ 1 2 3 4 1 2 3 4 \endverbatim Step 1 will build (for both situations above): - vx: 1,2,3,4 - vy: 6,7,8,9 Step 2 will build a table showing where the corners are: \verbatim | 0 1 2 3 | 0 1 2 3 --|---------| --|---------| 0 | . F F . | 0| . F . F 1 | F F F F | 1| F F . . 2 | F F F F | 2| F F . . 3 | . F F . | 3| . F . F \endverbatim Step 3: parse that table and turn on each corner: - left if "inner" corner (row and col = 1 or = 2) - right if "outer" corner (row and col = 0 or = 3) Final step: convert indexes to real coordinates Special note: if the rectangles have an identical coordinate, as in this example: \verbatim 9 +----+ | | 8 +---+----+ | | | 7 +---+----+ | | 6 +----+ 1 2 3 \endverbatim Then the vectors are: - vx: 1,2,3,3 - vy: 6,7,8,9 (notice the duped coordinate) This will produce a Polyline_ with 2 extra points:
(1,7)-(1,8)-(2,8)-(2,9)-(3,9)-(3,8)-(3,7)-(3,6)-(2,6)-(2,7)
instead of:
(1,7)-(1,8)-(2,8)-(2,9)-(3,9)-(3,6)-(2,6)-(2,7)
We solve this by proceeding an extra Polyline minimization, see PolylineBase::minimize() */ template template CPolyline_ FRect_::unionArea( const FRect_& other ) const { using namespace priv::runion; if( *this == other ) // if same rectangles, then return CPolyline_( other ); // returns one of them if( !this->intersects(other)() ) // if no intersection, { if( this->isInside( other ) ) return CPolyline_(other); if( other.isInside( *this ) ) return CPolyline_(*this); return CPolyline_(); // return empty polygon } /* step 0: make sure the rect with highest x is first. This is needed to avoid this kind of situation in table: F F . . . . F F . . F F F F . . That would happen if we have an identical x in the two rect (thus, they DO intersect), because when sorting, the rect with smallest index (1 or 2) is placed first in the vector of coordinates */ const auto* pr1 = this; const auto* pr2 = &other; if( pr1->getPts().first.getX() < pr2->getPts().first.getX() ) std::swap( pr1, pr2 ); const auto& r1 = *pr1; const auto& r2 = *pr2; // step 1: build vectors of coordinates and sort them std::array,4> vx, vy; int i=0; vx[i++] = Index( r1.getPts().first.getX(), 1 ); vx[i++] = Index( r1.getPts().second.getX(), 1 ); vx[i++] = Index( r2.getPts().first.getX(), 2 ); vx[i++] = Index( r2.getPts().second.getX(), 2 ); i=0; vy[i++] = Index( r1.getPts().first.getY(), 1 ); vy[i++] = Index( r1.getPts().second.getY(), 1 ); vy[i++] = Index( r2.getPts().first.getY(), 2 ); vy[i++] = Index( r2.getPts().second.getY(), 2 ); std::sort( vx.begin(), vx.end() ); std::sort( vy.begin(), vy.end() ); // priv::printArray( vx, "vx" ); priv::printArray( vy, "vy"); // step 2: fill table\n"; Table table; for( uint8_t r=0;r<4; r++ ) for( uint8_t c=0;c<4; c++ ) table[r][c] = Cell( vx[r], vy[c] ); // printTable( table, "after step 2" ); // step 3: parse table auto vpts = parseTable( table ); // priv::printVectorPairs( vpts ); // step 4: convert back vector of coordinates indexes into vector of coordinates auto res1 = convertToCoord( vpts, vx, vy ); res1.minimize(); // remove unecessary points return res1; } //------------------------------------------------------------------ namespace detail { //------------------------------------------------------------------ /// See getPtLabel( const Point2d_& pt, const Circle_& circle ) enum class PtTag: uint8_t { Inside, Outside, OnEdge }; /// Returns a label characterizing point \c pt, related to \c circle: inside, outside, or on edge of circle template PtTag getPtLabel( const Point2d_& pt, const Circle_& circle ) { if( pt.isInside( circle ) ) return PtTag::Inside; if( homog2d_abs( pt.distTo( circle.center() ) - circle.radius() ) < thr::nullDistance() ) return PtTag::OnEdge; return PtTag::Outside; } #if 0 /// Debug, can be removed after void printTag( std::string txt, PtTag tag ) { std::cout << "point " << txt << ": "; switch( tag ) { case PtTag::Inside: std::cout << "Inside\n"; break; case PtTag::Outside: std::cout << "Outside\n"; break; case PtTag::OnEdge: std::cout << "OnEdge\n"; break; } } #endif } // namespace detail //------------------------------------------------------------------ /// Segment/Segment intersection /** Algorithm:
- first compute the intersection point - then check if the intersection point lies in between the range of both segments, both on x and on y */ template template detail::Intersect Segment_::intersects( const Segment_& s2 ) const { // HOMOG2D_START; if( *this == s2 ) // same segment => no intersection return detail::Intersect(); Line2d_ l1 = getLine(); Line2d_ l2 = s2.getLine(); if( l1.isParallelTo( l2 ) ) // if parallel, return detail::Intersect(); // then, no intersection const auto& ptA1 = getPts().first; const auto& ptA2 = getPts().second; const auto& ptB1 = s2.getPts().first; const auto& ptB2 = s2.getPts().second; auto ptInter = l1 * l2; // intersection point auto dA1 = ptA1.distTo( ptInter ); auto dA2 = ptA2.distTo( ptInter ); if( homog2d_abs(dA1+dA2 - length()) < thr::nullDistance() ) { auto dB1 = ptB1.distTo( ptInter ); auto dB2 = ptB2.distTo( ptInter ); if( homog2d_abs(dB1+dB2 - s2.length()) < thr::nullDistance() ) return detail::Intersect( ptInter ); } return detail::Intersect(); // no intersection } //------------------------------------------------------------------ /// Segment/Line intersection /** Algorithm:
We check if the intersection point lies in between the range of the segment, both on x and on y */ template template detail::Intersect Segment_::intersects( const Line2d_& li1 ) const { // HOMOG2D_START; // HOMOG2D_LOG( "seg=" << *this << " line=" << li1 ); detail::Intersect out; auto li2 = getLine(); if( li1.isParallelTo( li2 ) ) // if parallel, no intersection return out; out._ptIntersect = li1 * li2; // intersection point const auto& pi = out._ptIntersect; const auto& pt1 = getPts().first; const auto& pt2 = getPts().second; auto d1 = pt1.distTo( pi ); auto d2 = pt2.distTo( pi ); if( homog2d_abs(d1+d2-length()) < thr::nullDistance() ) out._doesIntersect = true; return out; } //------------------------------------------------------------------ /// Segment/Circle intersection /** For each point of the segment, we need to consider 3 different situations - point is inside (PI) - point is outside (PO) - point is on the edge (PE) That makes 6 different situations to handle: - S1: PI-PI => no intersection - S2: PI-PO => 1 intersection - S3: PI-PE => 1 intersection - S4: PO-PO => depends on the support line: - S4A: if line does NOT intersects circle, no intersection pts - S4B: if line does intersects circle, and intersections point in segment area => 2 intersection pts - S4C: if line does intersects circle, and intersections point NOT in segment area => no intersection pts - S5: PO-PE => 1 intersection - S6: PE-PE => 2 intersections */ template template detail::IntersectM Segment_::intersects( const Circle_& circle ) const { // HOMOG2D_START; using detail::PtTag; auto tag_ptS1 = detail::getPtLabel( _ptS1, circle ); // get status of segment points related to circle (inside/outside/on-edge) auto tag_ptS2 = detail::getPtLabel( _ptS2, circle ); if( tag_ptS1 == PtTag::Inside ) if( tag_ptS2 == PtTag::Inside ) return detail::IntersectM(); auto int_lc = getLine().intersects( circle ); if( !int_lc() ) return detail::IntersectM(); auto p_pts = int_lc.get(); // get the line intersection points const auto& p1 = p_pts.first; const auto& p2 = p_pts.second; if( // one inside, the other outside ( tag_ptS1 == PtTag::Inside && tag_ptS2 == PtTag::Outside ) || ( tag_ptS1 == PtTag::Outside && tag_ptS2 == PtTag::Inside ) ) { detail::IntersectM out; auto d1 = _ptS1.distTo( p1 ); auto d2 = _ptS2.distTo( p1 ); if( homog2d_abs( d1+d2-length() ) < thr::nullDistance() ) out.add( p1 ); // points is inside else out.add( p2 ); return out; } detail::IntersectM out; if( tag_ptS1 == PtTag::Outside && tag_ptS2 == PtTag::Outside ) // both outside { auto d1 = _ptS1.distTo( p1 ); auto d2 = _ptS2.distTo( p1 ); if( homog2d_abs( d1+d2-length() ) >= thr::nullDistance() ) // if sum of the two distances is higher than length return detail::IntersectM(); out.add( p1 ); out.add( p2 ); return out; } if( tag_ptS1 == PtTag::OnEdge ) out.add( _ptS1 ); if( tag_ptS2 == PtTag::OnEdge ) out.add( _ptS2 ); return out; } ///////////////////////////////////////////////////////////////////////////// // SECTION - STREAMING OPERATORS ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Overload for points /// \todo 20221217: add a global switch (static function) to select a printing mode: either [x, y], either [a,b,c] template void base::LPBase::impl_op_stream( std::ostream& f, const Point2d_& pt ) const { if( pt.isInf() ) f << '[' << pt._v[0] << ',' << pt._v[1] << ',' << pt._v[2] << "]"; else f // << std::scientific << std::setprecision(25) << '[' << pt.getX() << ',' << pt.getY() << "]"; } /// Overload for lines template void base::LPBase::impl_op_stream( std::ostream& f, const Line2d_& r ) const { f << '[' << r._v[0] << ',' << r._v[1] << ',' << r._v[2] << "]"; } namespace base { /// Stream operator, free function, call member function pseudo operator impl_op_stream() template std::ostream& operator << ( std::ostream& f, const h2d::base::LPBase& r ) { r.impl_op_stream( f, r ); return f; } } /// Stream operator for a container of points/lines, free function template< typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > std::ostream& operator << ( std::ostream& f, const T& vec ) { f << "#="<< vec.size() << '\n'; for( const auto& elem: vec ) f << elem << '\n'; return f; } /// Stream operator for a pair of points/lines, free function template std::ostream& operator << ( std::ostream& f, const std::pair,base::LPBase>& pr ) { f << "std::pair (" << getString(pr.first.type()) << "-" << getString(pr.second.type()) << "):\n -first=" << pr.first << "\n -second=" << pr.second << ' '; return f; } template std::ostream& operator << ( std::ostream& f, const FRect_& r ) { f << "pt1: " << r._ptR1 << " pt2: " << r._ptR2; return f; } template std::ostream& operator << ( std::ostream& f, const Circle_& r ) { f << "center: " << r._center << ", radius=" << r._radius; return f; } template std::ostream& operator << ( std::ostream& f, const Segment_& seg ) { f << seg._ptS1 << "-" << seg._ptS2; return f; } namespace base { template std::ostream& operator << ( std::ostream& f, const h2d::base::PolylineBase& pl ) { f << (pl.isClosed() ? 'C' : 'O' ) << "Polyline: "; if( !pl.size() ) f << "empty"; else { for( const auto& pt: pl._plinevec ) f << pt << "-"; } return f; } } // namespace base template std::ostream& operator << ( std::ostream& f, const Ellipse_& ell ) { auto par = ell.template p_getParams(); f << par; return f; } //------------------------------------------------------------------ /// Returns the 2 parallel lines at distance \c dist from \c li (free function) template std::pair,Line2d_> getParallelLines( const Line2d_& li, T dist ) { return li.getParallelLines( dist ); } ///////////////////////////////////////////////////////////////////////////// // SECTION - MEMBER FUNCTION IMPLEMENTATION: CLASS Ellipse_ ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Returns standard parameters from matrix coeffs template template detail::EllParams Ellipse_::p_getParams() const { #ifdef HOMOG2D_OPTIMIZE_SPEED if( _epHasChanged ) { _par = p_computeParams(); _epHasChanged = false; } return _par; #else return p_computeParams(); #endif // HOMOG2D_OPTIMIZE_SPEED } //------------------------------------------------------------------ /// Compute and returns standard parameters from matrix coeffs /** \note In the checking below, we **cannot** print the matrix using the << operator, because it call this function, thus it would enter an infinite loop (an eventually SO). */ template template detail::EllParams Ellipse_::p_computeParams() const { const auto& m = detail::Matrix_::_mdata; HOMOG2D_INUMTYPE A = m[0][0]; HOMOG2D_INUMTYPE C = m[1][1]; HOMOG2D_INUMTYPE F = m[2][2]; HOMOG2D_INUMTYPE B = 2. * m[0][1]; HOMOG2D_INUMTYPE D = 2. * m[0][2]; HOMOG2D_INUMTYPE E = 2. * m[1][2]; detail::EllParams par; // theta already set to zero auto denom = B*B - (HOMOG2D_INUMTYPE)4. * A * C; #ifndef HOMOG2D_NOCHECKS if( homog2d_abs(denom) < thr::nullDenom() ) HOMOG2D_THROW_ERROR_1( "unable to compute parameters, denom=" << std::scientific << std::setprecision(15) << denom ); #endif par.x0 = ( (HOMOG2D_INUMTYPE)2.*C*D - B*E ) / denom; par.y0 = ( (HOMOG2D_INUMTYPE)2.*A*E - B*D ) / denom; auto common_ab = (HOMOG2D_INUMTYPE)2. * ( A*E*E + C*D*D - B*D*E + denom*F ); auto AmC = A-C; auto AmC2 = AmC*AmC; auto sqr = homog2d_sqrt(AmC2+B*B); par.a = -homog2d_sqrt( common_ab * ( A+C+sqr ) )/ denom; par.b = -homog2d_sqrt( common_ab * ( A+C-sqr ) )/ denom; par.a2 = par.a * par.a; par.b2 = par.b * par.b; if( homog2d_abs(B) < thr::nullDenom() ) { if( A > C ) par.theta = 90.; } else { auto t = (C - A - sqr) / B; par.theta = homog2d_atan( t ); } par.sint = homog2d_sin( par.theta ); par.cost = homog2d_cos( par.theta ); return par; } //------------------------------------------------------------------ /// Returns true if ellipse is a circle /** Using the matrix representation, if A = C and B = 0, then the ellipse is a circle You can provide the 0 threshold as and argument */ template bool Ellipse_::isCircle( HOMOG2D_INUMTYPE thres ) const { const auto& m = detail::Matrix_::_mdata; HOMOG2D_INUMTYPE A = m[0][0]; HOMOG2D_INUMTYPE C = m[1][1]; HOMOG2D_INUMTYPE B2 = m[0][1]; if( homog2d_abs(A-C) < thres ) if( homog2d_abs(B2)*2. < thres ) return true; return false; } //------------------------------------------------------------------ /// Returns center of ellipse /// \sa center( const T& ) template Point2d_ Ellipse_::getCenter() const { auto par = p_getParams(); return Point2d_( par.x0, par.y0 ); } //------------------------------------------------------------------ template std::pair Ellipse_::getMajMin() const { auto par = p_getParams(); return std::make_pair( par.a, par.b ); } //------------------------------------------------------------------ /// Returns angle of ellipse /// \sa angle( const Ellipse_& ) template HOMOG2D_INUMTYPE Ellipse_::angle() const { auto par = p_getParams(); return par.theta; } //------------------------------------------------------------------ /// Returns (approximate) perimeter of ellipse using the Ramanujan second formulae /** See https://en.wikipedia.org/wiki/Ellipse#Circumference */ template HOMOG2D_INUMTYPE Ellipse_::length() const { auto par = p_getParams(); auto ab_sum = par.a + par.b; auto ab_diff = par.a - par.b; auto h = ab_diff * ab_diff / (ab_sum * ab_sum); auto denom = (HOMOG2D_INUMTYPE)10. + homog2d_sqrt(4. - 3. * h); return (par.a + par.b) * M_PI * ( 1. + 3. * h / denom ); } //------------------------------------------------------------------ /// Returns pair of axis lines of ellipse template std::pair,Line2d_> Ellipse_::getAxisLines() const { auto par = p_getParams(); auto dy = static_cast(par.sint) * par.a; auto dx = static_cast(par.cost) * par.a; Point2d_ ptA( par.x0 + dx, par.y0 + dy ); auto pt0 = Point2d_( par.x0, par.y0 ); auto li_H = ptA * pt0; HOMOG2D_LOG( std::scientific << "ptA=" << ptA << " pt0=" << pt0 << " li_H=" << li_H ); auto li_V = li_H.getOrthogonalLine( pt0 ); return std::make_pair( li_H, li_V ); } //------------------------------------------------------------------ /// Returns axis-aligned bounding box of ellipse /** - see https://math.stackexchange.com/questions/91132/how-to-get-the-limits-of-rotated-ellipse */ template FRect_ Ellipse_::getBB() const { auto par = p_getParams(); auto vx = par.a2 * par.cost * par.cost + par.b2 * par.sint * par.sint; auto vy = par.a2 * par.sint * par.sint + par.b2 * par.cost * par.cost; auto vx_sq = homog2d_sqrt( vx ); auto vy_sq = homog2d_sqrt( vy ); return FRect_( Point2d_( par.x0 - vx_sq, par.y0 - vy_sq ), Point2d_( par.x0 + vx_sq, par.y0 + vy_sq ) ); } //------------------------------------------------------------------ /// Returns oriented bounding box of ellipse as a closed Polyline /** Algorithm: - build line \c liH going through major axis, by using center point and point on semi-major axis, intersecting ellipse - get opposite point \c ptB, lying on line and at distance \c a - get the two parallel lines to \c liH, at a distance \c b - get the two orthogonal lines at \c ptA and \c ptB \todo 20240330: unclear, the text above does not match what is done below (or does it?). Clarify that, and build a gif showing how this is done. */ template CPolyline_ Ellipse_::getOBB() const { // step 1: build ptA using angle auto par = p_getParams(); auto dy = par.sint * par.a; auto dx = par.cost * par.a; Point2d_ ptA( par.x0 + dx, par.y0 + dy ); auto pt0 = Point2d_( par.x0, par.y0 ); // step 2: build main-axis line, going through center and ptA auto li_H = ptA * pt0; // step 3: get ptB, using line and distance auto ppts = li_H.getPoints( pt0, par.a ); auto ptB = ppts.first; if( ptB == ptA ) ptB = ppts.second; auto para = getParallelLines( li_H, par.b ); auto li_V1 = li_H.getOrthogonalLine( ptA ); auto li_V2 = li_H.getOrthogonalLine( ptB ); CPolyline_ out; #ifndef HOMOG2D_DEBUGMODE out.p_addPoint( para.first * li_V1 ); out.p_addPoint( para.second * li_V1 ); out.p_addPoint( para.second * li_V2 ); out.p_addPoint( para.first * li_V2 ); #else auto p1 = para.first * li_V1; auto p2 = para.second * li_V1; auto p3 = para.second * li_V2; auto p4 = para.first * li_V2; HOMOG2D_DEBUG_ASSERT( ( p2!=p1 && p3!=p2 && p4!=p3 ), "p1=" << p1 << " p2=" << p2 << " p3=" << p3 << " p4=" << p4 << "\n para.1=" << para.first << "\n li_V1=" << li_V1 << "\n ptA=" << ptA << "\n ptB=" << ptB << "\n " << ppts ); out.p_addPoint( p1 ); out.p_addPoint( p2 ); out.p_addPoint( p3 ); out.p_addPoint( p4 ); #endif return out; } //------------------------------------------------------------------ /// Returns true if point is inside ellipse /** taken from https://stackoverflow.com/a/16814494/193789 */ template template bool Ellipse_::pointIsInside( const Point2d_& pt ) const { HOMOG2D_INUMTYPE x = pt.getX(); HOMOG2D_INUMTYPE y = pt.getY(); auto par = p_getParams(); const auto& x0 = par.x0; const auto& y0 = par.y0; auto v1 = par.cost * (x-x0) + par.sint * (y-y0); HOMOG2D_INUMTYPE sum = v1*v1 / par.a2; auto v2 = par.sint * (x-x0) - par.cost * (y-y0); sum += v2*v2 / par.b2; if( sum < 1. ) return true; return false; } ///////////////////////////////////////////////////////////////////////////// // SECTION - MEMBER FUNCTION IMPLEMENTATION: CLASS LPBase ///////////////////////////////////////////////////////////////////////////// namespace base { //------------------------------------------------------------------ /// Normalize to unit length, and make sure \c a is always >0 /** \todo Checkout: in what situation will we be unable to normalize? Is the test below relevant? Clarify. */ template void LPBase::impl_normalize( const detail::BaseHelper& ) const { // not sure which one is the best... #if 0 auto sq = homog2d_sqrt( _v[0]*_v[0] + _v[1]*_v[1] ); #else #ifdef HOMOG2D_USE_TTMATH auto sq = ttmath::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] ); #else auto sq = std::hypot( _v[0], _v[1] ); #endif #endif #ifndef HOMOG2D_NOCHECKS if( sq <= std::numeric_limits::epsilon() ) HOMOG2D_THROW_ERROR_1( "unable to normalize line, sq=" << sq << " values: a=" << _v[0] << " b=" << _v[1] << " c=" << _v[2] ); #endif for( int i=0; i<3; i++ ) const_cast*>(this)->_v[i] /= sq; // needed to remove constness if( _v[0] < 0. ) // a always >0 for( int i=0; i<3; i++ ) const_cast*>(this)->_v[i] = -_v[i]; if( homog2d_abs(_v[0]) < thr::nullDenom() ) // then, change sign so that b>0 if( _v[1] < 0 ) { const_cast*>(this)->_v[1] = - _v[1]; const_cast*>(this)->_v[2] = - _v[2]; } } //------------------------------------------------------------------ /// Point normalization. Just make sure a>0 template void LPBase::impl_normalize( const detail::BaseHelper& ) const { if( _v[0] < 0. ) { const_cast*>(this)->_v[0] = -_v[0]; const_cast*>(this)->_v[1] = -_v[1]; const_cast*>(this)->_v[2] = -_v[2]; } #ifndef HOMOG2D_NOCHECKS if( homog2d_abs(_v[2]) < thr::nullDenom() && ( _v[0] < thr::nullOrthogDistance() && homog2d_abs(_v[1]) < thr::nullOrthogDistance() ) ) HOMOG2D_THROW_ERROR_1( "invalid point values" ); #endif } //------------------------------------------------------------------ template constexpr FPT LPBase::impl_getCoord( GivenCoord, FPT, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getCoord() on a point" ); } template FPT LPBase::impl_getCoord( GivenCoord gc, FPT other, const detail::BaseHelper& ) const { const auto a = static_cast( _v[0] ); const auto b = static_cast( _v[1] ); auto denom = ( gc == GivenCoord::X ? b : a ); #ifndef HOMOG2D_NOCHECKS if( homog2d_abs(denom) < thr::nullDenom() ) HOMOG2D_THROW_ERROR_2( "getCoord", "null denominator encountered" ); #endif if( gc == GivenCoord::X ) return ( -a * other - _v[2] ) / b; else return ( -b * other - _v[2] ) / a; } template constexpr Point2d_ LPBase::impl_getPoint( GivenCoord, FPT, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getPoint() on a point" ); } template Point2d_ LPBase::impl_getPoint( GivenCoord gc, FPT other, const detail::BaseHelper& ) const { auto coord = impl_getCoord( gc, other, detail::BaseHelper() ); if( gc == GivenCoord::X ) return Point2d_( other, coord ); return Point2d_( coord, other ); } //------------------------------------------------------------------ /// ILLEGAL INSTANCIATION template template constexpr PointPair1_ LPBase::impl_getPoints_A( GivenCoord, FPT, FPT2, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getPoints() on a point" ); } template template constexpr PointPair1_ LPBase::impl_getPoints_B( const Point2d_&, FPT2, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getPoints() on a point" ); } /// Returns pair of points on line at distance \c dist from point on line at coord \c coord. Implementation for lines template template PointPair1_ LPBase::impl_getPoints_A( GivenCoord gc, FPT coord, FPT2 dist, const detail::BaseHelper& ) const { const auto pt = impl_getPoint( gc, coord, detail::BaseHelper() ); return priv::getPoints_B2( pt, dist, *this ); } /// Returns pair of points on line at distance \c dist from point on line at coord \c coord. Implementation for lines template template PointPair1_ LPBase::impl_getPoints_B( const Point2d_& pt, FPT2 dist, const detail::BaseHelper& ) const { #ifndef HOMOG2D_NOCHECKS if( this->distTo( pt ) > thr::nullDistance() ) { std::cerr << "homog2d: distance=" << std::scientific << this->distTo( pt ) << " > null distance (" << thr::nullDistance() << ")\n"; HOMOG2D_THROW_ERROR_2( "getPoints", "point is not on line" ); } #endif return priv::getPoints_B2( pt, dist, *this ); } //------------------------------------------------------------------ /// Illegal instanciation template constexpr Line2d_ LPBase::impl_getOrthogonalLine_A( GivenCoord, FPT, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getOrthogonalLine() on a point" ); } /// Illegal instanciation template constexpr Line2d_ LPBase::impl_getOrthogonalLine_B( const Point2d_&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getOrthogonalLine() on a point" ); } /// Illegal instanciation template template constexpr Line2d_ LPBase::impl_getRotatedLine( const Point2d_&, T, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getRotatedLine() on a point" ); } /// Returns an orthogonal line, implementation of getOrthogonalLine(). template Line2d_ LPBase::impl_getOrthogonalLine_A( GivenCoord gc, FPT val, const detail::BaseHelper& ) const { auto other_val = impl_getCoord( gc, val, detail::BaseHelper() ); Point2d_ pt( other_val, val ) ; if( gc == GivenCoord::X ) pt.set( val, other_val ); return priv::getOrthogonalLine_B2( pt, *this ); } /// Returns an orthogonal line, implementation of getOrthogonalLine(). template Line2d_ LPBase::impl_getOrthogonalLine_B( const Point2d_& pt, const detail::BaseHelper& ) const { #ifndef HOMOG2D_NOCHECKS if( this->distTo( pt ) > thr::nullDistance() ) { std::cerr << "homog2d: distance=" << std::scientific << this->distTo( pt ) << "> null distance (" << thr::nullDistance() << ")\n"; HOMOG2D_THROW_ERROR_2( "getOrthogonalLine", "point is not on line" ); } #endif return priv::getOrthogonalLine_B2( pt, *this ); } template template Line2d_ LPBase::impl_getRotatedLine( const Point2d_& pt, T theta, const detail::BaseHelper& ) const { #ifndef HOMOG2D_NOCHECKS if( this->distTo( pt ) > thr::nullDistance() ) { std::cerr << "homog2d: distance=" << std::scientific << this->distTo( pt ) << "> null distance (" << thr::nullDistance() << ")\n"; HOMOG2D_THROW_ERROR_2( "getLineAngle", "point is not on line" ); } #endif Homogr_ H; H.addTranslation( -pt.getX(), -pt.getY() ).addRotation(theta).addTranslation( pt.getX(), pt.getY() ); return H * *this; } /// Returns the shortest segment that joins a point and a line template Segment_ LPBase::impl_getOrthogSegment( const Point2d_& pt, const detail::BaseHelper& ) const { Line2d_ src = *this; // copy to highest precision auto dist = src.distTo(pt); #ifndef HOMOG2D_NOCHECKS if( dist < thr::nullDistance() ) // sanity check HOMOG2D_THROW_ERROR_1( "unable to compute segment" ); #endif auto pair_lines = getParallelLines( dist ); // determine on which of the two parallel lines does the point lie? Line2d_* pline = &pair_lines.first; if( pt.distTo(pair_lines.second) < thr::nullDistance() ) pline = &pair_lines.second; auto oline = pline->getOrthogonalLine( pt ); auto p2 = *this * oline; return Segment_( pt, p2 ); } template constexpr Segment_ LPBase::impl_getOrthogSegment( const Point2d_&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid call" ); } //------------------------------------------------------------------ /// Illegal instanciation template constexpr Line2d_ LPBase::impl_getParallelLine( const Point2d_&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getParallelLine() on a point" ); } /// Returns an parallel line, implementation of getParallelLine(). template Line2d_ LPBase::impl_getParallelLine( const Point2d_& pt, const detail::BaseHelper& ) const { Line2d_ out = *this; out._v[2] = static_cast(-_v[0]) * pt.getX() - _v[1] * pt.getY(); out.p_normalizePL(); return out; } //------------------------------------------------------------------ /// Illegal instanciation template template constexpr std::pair,Line2d_> LPBase::impl_getParallelLines( T, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call getParallelLines() on a point" ); } /// Implementation for lines template template std::pair,Line2d_> LPBase::impl_getParallelLines( T dist, const detail::BaseHelper& ) const { Line2d_ l1 = *this; Line2d_ l2 = *this; l1._v[2] = static_cast(this->_v[2]) + dist; l2._v[2] = static_cast(this->_v[2]) - dist; return std::make_pair( l1, l2 ); } //------------------------------------------------------------------ /// Comparison operator, used for lines /** Definition: two lines will be equal: - if they are parallel AND - if their offset (3 third value) is less than nullDistance() */ template bool LPBase::impl_op_equal( const LPBase& other, const detail::BaseHelper& ) const { if( !this->isParallelTo( other ) ) return false; if( std::fabs( _v[2] - other._v[2] ) > thr::nullDistance() ) return false; return true; } /// Comparison operator, used for points template bool LPBase::impl_op_equal( const LPBase& other, const detail::BaseHelper& ) const { return ( this->distTo( other ) < thr::nullDistance() ); } //------------------------------------------------------------------ /// Sorting operator, for points template bool LPBase::impl_op_sort( const LPBase& other, const detail::BaseHelper& ) const { if( getX() < other.getX() ) return true; if( getX() > other.getX() ) return false; if( getY() < other.getY() ) return true; return false; } /// Sorting operator, for lines template constexpr bool LPBase::impl_op_sort( const LPBase&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid < operator: you cannot sort lines" ); return false; // to avoid a warning } } // namespace base //------------------------------------------------------------------ /// Inner implementation details namespace detail { /// Cross product, see https://en.wikipedia.org/wiki/Cross_product#Coordinate_notation template base::LPBase crossProduct( const base::LPBase& r1, const base::LPBase& r2 ) { auto r1_a = static_cast(r1._v[0]); auto r1_b = static_cast(r1._v[1]); auto r1_c = static_cast(r1._v[2]); auto r2_a = static_cast(r2._v[0]); auto r2_b = static_cast(r2._v[1]); auto r2_c = static_cast(r2._v[2]); base::LPBase res; res._v[0] = static_cast( r1_b * r2_c - r1_c * r2_b ); res._v[1] = static_cast( r1_c * r2_a - r1_a * r2_c ); res._v[2] = static_cast( r1_a * r2_b - r1_b * r2_a ); return res; } } // namespace detail //------------------------------------------------------------------ /////////////////////////////////////////// // CONSTRUCTORS /////////////////////////////////////////// namespace base { /// Points overload: generic init from two numeric args template template void LPBase::impl_init_2( const T1& v1, const T2& v2, const detail::BaseHelper& ) { _v[0] = v1; _v[1] = v2; _v[2] = 1.; } /// Lines overload: generic init from two numeric args template template void LPBase::impl_init_2( const T1& v1, const T2& v2, const detail::BaseHelper& ) { Point2d_ pt1; // 0,0 Point2d_ pt2(v1,v2); pt2.p_normalizePL(); *this = detail::crossProduct( pt1, pt2 ); p_normalizePL(); } /// Overload for point to point distance template template HOMOG2D_INUMTYPE LPBase::impl_distToPoint( const Point2d_& pt, const detail::BaseHelper& ) const { #ifndef HOMOG2D_USE_TTMATH return std::hypot( static_cast( getX() ) - static_cast( pt.getX() ), static_cast( getY() ) - static_cast( pt.getY() ) ); #else auto v1 = static_cast( getX() ) - static_cast( pt.getX() ); auto v2 = static_cast( getY() ) - static_cast( pt.getY() ); return ttmath::Sqrt( v1*v1 + v2*v2 ); #endif } //------------------------------------------------------------------ /// Returns distance between the line and point \b pt. overload for line to point distance. /** http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
        | a.x0 + b.y0 + c |
  d = -----------------------
         sqrt( a*a + b*b )
\todo Do we really require computation of hypot ? (because the line is supposed to be normalized, i.e. h=1 ?) */ template template HOMOG2D_INUMTYPE LPBase::impl_distToPoint( const Point2d_& pt, const detail::BaseHelper& ) const { return std::fabs( static_cast( _v[0] ) * pt.getX() + static_cast( _v[1] ) * pt.getY() + static_cast( _v[2] ) ) #ifndef HOMOG2D_USE_TTMATH / std::hypot( _v[0], _v[1] ); #else ; #endif } /// overload for line to point distance template template HOMOG2D_INUMTYPE LPBase::impl_distToLine( const Line2d_& li, const detail::BaseHelper& ) const { return li.distTo( *this ); } /// overload for line to line distance. Aborts build if instanciated (distance between two lines makes no sense). template template constexpr HOMOG2D_INUMTYPE LPBase::impl_distToLine( const Line2d_&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot compute distance between two lines" ); return 0.; // to avoid warning message on build } template template constexpr HOMOG2D_INUMTYPE LPBase::impl_distToSegment( const Segment_&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot compute distance between a line and a segment" ); return 0.; // to avoid warning message on build } template template HOMOG2D_INUMTYPE LPBase::impl_distToSegment( const Segment_& seg, const detail::BaseHelper& ) const { return seg.distTo( *this ); } } // namespace base //------------------------------------------------------------------ /// Free function, returns the angle (in Rad) between two lines/ or segments /// \sa Segment_::getAngle() /// \sa Line2d_::getAngle() template HOMOG2D_INUMTYPE getAngle( const T1& li1, const T2& li2 ) { return li1.getAngle( li2 ); } namespace base { template template bool LPBase::impl_isParallelTo( const LPBase& li, const detail::BaseHelper& ) const { if( getAngle(li) < thr::nullAngleValue() ) return true; return false; } template template constexpr bool LPBase::impl_isParallelTo( const LPBase&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot use IsParallel() with a point" ); return false; // to avoid a warning message on build } //------------------------------------------------------------------ /// Returns the angle (in rad.) between the line and the other one. /** - The returned value will be in the range [0, M_PI/2] If the lines \f$ (a_0,a_1,a_2) \f$ and \f$ (b_0,b_1,b_2) \f$ are correctly normalized (should be, but...) then the angle between them is \f$ \alpha = acos( a_0*b_0 + a_1*b_1) \f$.
However, in "some situations", even if the lines have been previously normalized (which is the case here) we can encounter numerical issues, so here we "reinforce the normalization and compute: \f[ \alpha = acos \left( \frac{a_0*b_0 + a_1*b_1} { \sqrt{ a_0^2 + a_1^2 } * \sqrt{ b_0^2 + b_1^2 } } \right) \f] In some situations, the value inside the parenthesis "may" be equal to \f$ 1+\epsilon \f$ (typically, something like "1.0000000000123"). This is out of bounds for the \f$ acos() \f$ function, that will then return "nan" (thus induce some failure further on).
To avoid this, a checking is done, and any value higher than 1 will be truncated. This is logged on \c std::cerr so that the user may take that into consideration. \todo more investigation needed ! : what are the exact situation that will lead to this event? */ template HOMOG2D_INUMTYPE LPBase::impl_getAngle( const LPBase& li, const detail::BaseHelper& ) const { HOMOG2D_INUMTYPE l1_a = _v[0]; HOMOG2D_INUMTYPE l1_b = _v[1]; HOMOG2D_INUMTYPE l2_a = li._v[0]; HOMOG2D_INUMTYPE l2_b = li._v[1]; HOMOG2D_INUMTYPE res = l1_a * l2_a + l1_b * l2_b; res /= homog2d_sqrt( (l1_a*l1_a + l1_b*l1_b) * (l2_a*l2_a + l2_b*l2_b) ); auto fres = homog2d_abs(res); if( fres > 1.0 ) { HOMOG2D_LOG_WARNING( "homog2d Warning: angle computation overflow detected, value " << std::scientific << std::setprecision(20) << fres << ", truncated to 1.0\n input lines:\n l1: " << *this << "\n l2: " << li ); fres = 1.0; } return homog2d_acos( fres ); } template constexpr HOMOG2D_INUMTYPE LPBase::impl_getAngle( const LPBase&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "cannot get angle of a point" ); return 0.; // to avoid a warning } //------------------------------------------------------------------ /// Returns true if point is inside (or on the edge) of a flat rectangle defined by (p0,p1) template template bool LPBase::impl_isInsideRect( const FRect_& rect, const detail::BaseHelper& ) const { auto pair_pts = rect.getPts(); const auto& p00 = pair_pts.first; const auto& p11 = pair_pts.second; return detail::ptIsInside( *this, p00, p11 ); } template template constexpr bool LPBase::impl_isInsideRect( const FRect_&, const detail::BaseHelper& ) const { return false; } //------------------------------------------------------------------ template template bool LPBase::impl_isInsideCircle( const Point2d_& center, T radius, const detail::BaseHelper& ) const { if( distTo( center ) < radius ) return true; return false; } template template constexpr bool LPBase::impl_isInsideCircle( const Point2d_&, T, const detail::BaseHelper& ) const { return false; } //------------------------------------------------------------------ #if 0 // deprecated. kept here just in case it would be useful in the future (who knows?) namespace sub { /// Returns the index of the segment (among the ones in \c bbox) /// that is the farthest from point \c pt template int getFarthestSegment( const Point2d_& pt, const FRect_& bbox ) { auto segs = bbox.getSegs(); using localPair = std::pair; std::array dist_array; for( int i=0; i<4; i++ ) { dist_array[i].first = pt.distTo( segs[i].getLine() ); dist_array[i].second = i; } std::sort( dist_array.begin(), dist_array.end(), [] // lambda ( const localPair& i, const localPair& j ) { return i.first < j.first; } ); return dist_array[3].second; } } // namespace sub #endif //------------------------------------------------------------------ /// Returns true if point is inside closed Polyline /** See https://en.wikipedia.org/wiki/Point_in_polygon Implementing this algorithm in a naive way may encounter some issue: if the main segment used to compute the number of intersections is colinear with one of the points of the polyline, we might run into a numerical issue. See page \ref md_docs_homog2d_algorithms for details. */ template template bool LPBase::impl_isInsidePoly( const base::PolylineBase& poly, const detail::BaseHelper& ) const { if( !poly.isPolygon() ) return false; // step 1: if point is outside bounding box, return false auto bbox = poly.getBB(); if( !this->isInside(bbox) ) return false; // step 2: check if point is lying on one of the segments. If so, return false; for( auto seg: poly.getSegs() ) if( seg.getLine().distTo( *this ) < thr::nullDistance() ) return false; auto segarr = bbox.getExtended().getSegs(); std::vector> seg_bb( segarr.begin(), segarr.end() ); int iter = 0; do { bool tooClose = false; for( size_t i=0; i seg_ref( *this, seg_bb[i].getCenter() ); auto poly_pts = poly.getPts(); for( size_t j=0; j( c%2 ); } } iter++; if( tooClose ) // if unable to find a suitable reference segment, { // we double the number of segments, by splitting each of them if( iter < HOMOG2D_MAXITER_PIP ) // (only if not max iterations!) { std::vector> seg_bb2( seg_bb.size() * 2 ); for( size_t j=0; j template constexpr bool LPBase::impl_isInsidePoly( const base::PolylineBase&, const detail::BaseHelper& ) const { return false; } //------------------------------------------------------------------ template template bool LPBase::impl_isInsideEllipse( const Ellipse_& ell, const detail::BaseHelper& ) const { return ell.pointIsInside( *this ); } template template constexpr bool LPBase::impl_isInsideEllipse( const Ellipse_&, const detail::BaseHelper& ) const { return false; } //------------------------------------------------------------------ /// Intersection of line and circle: implementation for points template template constexpr detail::Intersect LPBase::impl_intersectsCircle( const Point2d_&, T, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "cannot use intersects(Circle) with a point" ); } /// Line/Circle intersection: implementation /// For computation details, checkout http://skramm.lautre.net/files/misc/intersect_circle_line.pdf /// \todo 20230124: change return type: intersection of a circle with a line can be a single point template template detail::Intersect LPBase::impl_intersectsCircle( const Point2d_& pt, ///< circle origin T radius, ///< radius const detail::BaseHelper& ///< dummy arg, needed so that this overload is only called for lines ) const { detail::Intersect out; HOMOG2D_CHECK_IS_NUMBER(T); const HOMOG2D_INUMTYPE a = static_cast(_v[0]); // just to lighten a bit... const HOMOG2D_INUMTYPE b = static_cast(_v[1]); const HOMOG2D_INUMTYPE c = static_cast(_v[2]); // step 1: translate to origin auto cp = pt.getX() * a + pt.getY() * b + c; // step 2: compute distance between center (origin) and middle point auto a2b2 = a * a + b * b; auto d0 = homog2d_abs(cp) / homog2d_sqrt( a2b2 ); if( radius < d0 ) // if less than radius, return out; // no intersection auto d2 = radius*radius - d0*d0; // step 3: compute coordinates of middle point B auto xb = - a * cp / a2b2; auto yb = - b * cp / a2b2; // step 4: compute coordinates of intersection points, with center at (0,0) auto m = homog2d_sqrt( d2 / a2b2 ); auto x1 = xb + m*b; auto y1 = yb - m*a; auto x2 = xb - m*b; auto y2 = yb + m*a; // last step: translate back out._ptIntersect_1.set( x1 + pt.getX(), y1 + pt.getY() ); out._ptIntersect_2.set( x2 + pt.getX(), y2 + pt.getY() ); out._doesIntersect = true; priv::fix_order( out._ptIntersect_1, out._ptIntersect_2 ); return out; } //------------------------------------------------------------------ /// Overload used when attempting to use that on a point template template constexpr detail::IntersectM LPBase::impl_intersectsFRect( const FRect_&, const detail::BaseHelper& ) const { static_assert( detail::AlwaysFalse::value, "Invalid: you cannot call intersects(FRect) on a point" ); } /// Line/FRect intersection /** Returns 0, 1, or 2 intersections (1 when the line touches the edge of rectangle) \todo 20230118: check if there is a need to fix the order of points in returned object (this was done previously) */ template template detail::IntersectM LPBase::impl_intersectsFRect( const FRect_& rect, const detail::BaseHelper& ) const { std::vector> pti; for( const auto& seg: rect.getSegs() ) // get segment of rectangle { auto ppts_seg = seg.getPts(); auto inters = seg.intersects( *this ); if( inters() ) { bool storePoint(true); auto pt = inters.get(); if( pt == ppts_seg.first || pt == ppts_seg.second ) // if intersection point is one of the segment pts if( pti.size() == 1 ) { // AND if there is already one if( pti[0] == pt ) // AND that one is already stored storePoint = false; // THEN don't store it } if( storePoint ) pti.push_back( pt ); if( pti.size() == 2 ) // already got 2, done break; } } if( pti.empty() ) return detail::IntersectM(); detail::IntersectM out; out.add( pti ); return out; } } // namespace base ///////////////////////////////////////////////////////////////////////////// // SECTION - PRODUCT OPERATORS DEFINITIONS (HELPER FUNCTIONS) ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Apply homography to a vector/array/list (type T) of points or lines. template template void Hmatrix_::applyTo( T& vin ) const { #if 0 std::transform( std::begin( vin ), std::end( vin ), std::begin( vin ), [](){} ); #else for( auto& elem: vin ) elem = *this * elem; #endif } namespace detail { /// Implementation of product 3x3 by 3x3 template void product( Matrix_& out, const Matrix_& h1, const Matrix_& h2 ) { out.p_fillZero(); for( int i=0; i<3; i++ ) for( int j=0; j<3; j++ ) for( int k=0; k<3; k++ ) out.value(i,j) += static_cast( h1.value(i,k) ) * h2.value(k,j); } /// Implementation of product 3x3 by 3x1 /** - T1 and T2: type::IsLine or type::IsPoint (same but also different) */ template void product( base::LPBase& out, const Matrix_& h, const base::LPBase& in ) { for( int i=0; i<3; i++ ) { auto sum = static_cast(h._mdata[i][0]) * in._v[0]; sum += h._mdata[i][1] * in._v[1]; sum += h._mdata[i][2] * in._v[2]; out._v[i] = sum; } } } // namespace detail ///////////////////////////////////////////////////////////////////////////// // SECTION - PRODUCT OPERATORS DEFINITIONS ///////////////////////////////////////////////////////////////////////////// /// Free function template, product of two lines, returns a point template Point2d_ operator * ( const Line2d_& lhs, const Line2d_& rhs ) { #ifndef HOMOG2D_NOCHECKS if( lhs.isParallelTo(rhs) ) HOMOG2D_THROW_ERROR_1( "lines are parallel, unable to compute product:\nlhs=" << lhs << " rhs=" << rhs ); #endif return detail::crossProduct(lhs, rhs); } /// Free function template, product of two segments, returns a point template Point2d_ operator * ( const Segment_& lhs, const Segment_& rhs ) { return lhs.getLine() * rhs.getLine(); } #ifdef HOMOG2D_FUTURE_STUFF /// Apply Epipolar matrix to a point or line, this will return the opposite type. /// Free function, templated by point or line template base::LPBase::OtherType,V> operator * ( const Hmatrix_& h, const base::LPBase& in ) { LPBase::OtherType,V> out; detail::product( out, h._data, in ); return out; } #endif /// Free function, apply homography to a point. template Point2d_ operator * ( const Homogr_& h, const Point2d_& in ) { Point2d_ out; detail::product( out, h, in ); return out; } /// Free function, apply homography to a line. template Line2d_ operator * ( const Homogr_& h, const Line2d_& in ) { if( h._hmt == nullptr ) // if H^-T not allocated yet, do it { h._hmt = std::unique_ptr>( new detail::Matrix_() ); h._hasChanged = true; } if( h._hasChanged ) // if homography has changed, recompute inverse and transposed { auto hi = h; auto mat_inv = hi.inverse().transpose(); *h._hmt = mat_inv; h._hasChanged = false; } Line2d_ out; detail::product( out, *h._hmt, in ); out.p_normalizePL(); return out; } /// Apply homography to a Segment template Segment_ operator * ( const Homogr_& h, const Segment_& seg ) { const auto& pts = seg.getPts(); Point2d_ pt1 = h * pts.first; Point2d_ pt2 = h * pts.second; return Segment_( pt1, pt2 ); } /// Apply homography to a Polyline template base::PolylineBase operator * ( const Homogr_& h, const base::PolylineBase& pl ) { base::PolylineBase out; const auto& pts = pl.getPts(); for( const auto pt: pts ) out.p_addPoint( h * pt ); return out; } /// Apply homography to a flat rectangle produces a closed polyline template CPolyline_ operator * ( const h2d::Homogr_& h, const h2d::FRect_& rin ) { h2d::CPolyline_ out; for( const auto& pt: rin.get4Pts() ) out.p_addPoint( h * pt ); return out; } /// Apply homography to a Ellipse, produces an Ellipse /** \f[ Q' = H^{-T} \cdot Q \cdot H^{-1} \f] */ template Ellipse_ operator * ( const Homogr_& h, const Ellipse_& ell_in ) { auto hm = static_cast>(h); hm.inverse(); auto hmt = hm; hmt.transpose(); const auto& ell_in2 = static_cast>(ell_in); auto prod = hmt * ell_in2 * hm; Ellipse_ out( prod ); return out; } /// Apply homography to a Circle, produces an Ellipse /** Converts the circle to an ellipse, then calls the corresponding code */ template Ellipse_ operator * ( const Homogr_& h, const Circle_& cir ) { Ellipse_ ell_in( cir ); return h * ell_in; } //------------------------------------------------------------------ namespace priv { /// Allocation for \c std::array container /** \sa operator * ( const Hmatrix_&, const Cont& ); */ template< typename Cont, typename std::enable_if< trait::Is_std_array::value, Cont >::type* = nullptr > Cont alloc( std::size_t /* unused here */ ) { return Cont(); } /// Allocation for \c std::vector and \c std::list container template< typename Cont, typename std::enable_if< !trait::Is_std_array::value, Cont >::type* = nullptr > Cont alloc( std::size_t nb ) { return Cont(nb); } } // namespace priv //------------------------------------------------------------------ /// Used to proceed multiple products, whatever the container (\c std::list, \c std::vector, or \c std::array). /// Returned container is of same type as given input template typename std::enable_if::value,Cont>::type operator * ( const Hmatrix_& h, ///< Matrix const Cont& vin ///< Input container ) { Cont vout = priv::alloc( vin.size() ); auto it = std::begin( vout ); for( const auto& elem: vin ) *it++ = h * elem; return vout; } ///////////////////////////////////////////////////////////////////////////// // SECTION - FREE FUNCTIONS ///////////////////////////////////////////////////////////////////////////// template FPT getX( const Point2d_& pt) { return pt.getX(); } template FPT getY( const Point2d_& pt) { return pt.getY(); } //------------------------------------------------------------------ /// Free function, returns -1 or +1 depending on the side of \c pt, related to \c li template int side( const Point2d_& pt, const Line2d_& li ) { const auto& arr = li.get(); HOMOG2D_INUMTYPE a = arr[0]; HOMOG2D_INUMTYPE b = arr[1]; HOMOG2D_INUMTYPE c = arr[2]; auto dist = a * pt.getX() + b * pt.getY() + c; if( homog2d_abs(dist) < thr::nullDistance() ) return 0; return (std::signbit( dist ) ? -1 : +1); } /// Free function, distance between points /// \sa Point2d_::distTo() template HOMOG2D_INUMTYPE dist( const Point2d_& pt1, const Point2d_& pt2 ) { return pt1.distTo( pt2 ); } template HOMOG2D_INUMTYPE IoU( const FRect_& r1, const FRect_& r2 ) { auto ia = intersectArea( r1, r2 ); if( ia() ) return ia.get().area() / unionArea( r1, r2 ).area(); return 0.; } /// Free function, see FRect_::unionArea() template CPolyline_ unionArea( const FRect_& r1, const FRect_& r2 ) { return r1.unionArea(r2); } /// Free function, see FRect_::intersectArea() template detail::RectArea intersectArea( const FRect_& r1, const FRect_& r2 ) { return r1.intersectArea(r2); } /// Returns circle passing through 4 points of flat rectangle (free function) /// \sa FRect_::getBoundingCircle() template Circle_ getBoundingCircle( const FRect_& rect ) { return rect.getBoundingCircle(); } /// Return circle inscribed in rectangle /// \sa hFRect_::getInscribedCircle() template Circle_ getInscribedCircle( const FRect_& rect ) { return rect.getInscribedCircle(); } /// Returns true if is a polygon (free function) ///  \sa PolylineBase::isPolygon() template bool isPolygon( const base::PolylineBase& pl ) { return pl.isPolygon(); } /// Returns true if polygon is convex (free function) /// \sa PolylineBase::isConvex() template bool isConvex( const base::PolylineBase& poly ) { return poly.isConvex(); } /// Returns the number of segments (free function) /// \sa PolylineBase::nbSegs() template size_t nbSegs( const base::PolylineBase& pl ) { return pl.nbSegs(); } /// Returns the segments of the polyline (free function) /// \sa PolylineBase::getSegs() template std::vector> getSegs( const base::PolylineBase& pl ) { return pl.getSegs(); } /// Returns the number of points (free function) /// \sa PolylineBase::size() template size_t size( const base::PolylineBase& pl ) { return pl.size(); } /// Rotates the primitive (only available for Polyline and FRect) around (0,0) template void rotate( T& prim, Rotate rot ) { prim.rotate( rot ); } /// Rotates the primitive (only available for Polyline and FRect) around \c refpt template void rotate( T& prim, Rotate rot, const Point2d_& refpt ) { prim.rotate( rot, refpt ); } /// Returns width, height of rectangle (free function) /// \sa FRect_::size() template std::pair size( const FRect_& rect ) { return rect.size(); } /// Returns Bounding Box of Ellipse_ (free function) /// \sa Ellipse_::getBB() template CPolyline_ getOBB( const Ellipse_& ell ) { return ell.getOBB(); } /// Holds function returning a pair of points namespace ppair { /// Returns pair of points defining a BB of a primitive /// Overload 1/4, private free function /** Arg is a neither a Point2d, a Segment, a Line2d or a polyline */ template< typename T, typename std::enable_if< ( // arg is !std::is_same>::value // not a line, && !std::is_same>::value // not a point, && !std::is_same>::value // not a segment. && !std::is_same>::value // not a CPolyline && !std::is_same>::value // not a COolyline ) ,T >::type* = nullptr > auto getPointPair( const T& elem ) { HOMOG2D_START; return elem.getBB().getPts(); } /// Returns pair of points defining a BB of a Polyline /// Overload 2/4, private free function template< typename T, typename std::enable_if< ( std::is_same>::value || std::is_same>::value ) ,T >::type* = nullptr > auto getPointPair( const T& poly ) { HOMOG2D_START; if( poly.size() == 0 ) HOMOG2D_THROW_ERROR_1( "cannot compute point pair of empty Polyline" ); if( poly.size() == 2 ) { const auto& pts = poly.getPts(); return std::make_pair( pts[0], pts[1] ); } return priv::getBB_Points( poly.getPts() ); } /// Returns pair of points defining a BB of a Point2d /// Overload 3/4, private free function /** This seems useless at first glance, but it is used to get the common bounding box of two objects when one of them (or both) is a point. */ template< typename T, typename std::enable_if< std::is_same>::value ,T >::type* = nullptr > auto getPointPair( const T& elem ) { HOMOG2D_START; return std::make_pair( Point2d_(elem), Point2d_(elem) ); } /// Returns pair of points defining a BB of Segment /// Overload 4/4, private free function template< typename T, typename std::enable_if< std::is_same>::value ,T >::type* = nullptr > auto getPointPair( const T& elem ) { HOMOG2D_START; return elem.getPts(); } /// Needed because of variant template PointPair1_ getPointPair( const Line2d_& ) { HOMOG2D_START; HOMOG2D_THROW_ERROR_1( "Unable to get pair of points for a Line2d" ); } } // namespace ppair ///////////////////////////////////////////////////////////////////////////// // SECTION - FUNCTORS USED IN FOLLOWING FREE FUNCTIONS ///////////////////////////////////////////////////////////////////////////// #ifdef HOMOG2D_ENABLE_VRTP namespace fct { //------------------------------------------------------------------ /// A functor to get pair of points englobing the element /** \sa CommonType_ \sa area() */ struct PtPairFunct { template PointPair1_ operator ()(const T& a) { return ppair::getPointPair(a); } }; //------------------------------------------------------------------ /// A functor to get the Bounding Box /// \sa getBB() struct BBFunct { template // FRect_ operator ()(const T& a) FRect_ operator ()(const T& a) { return a.getBB(); } }; } // namespace fct #endif // HOMOG2D_ENABLE_VRTP ///////////////////////////////////////////////////////////////////////////// // SECTION - FREE FUNCTIONS HANDLING VARIANT TYPE ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Returns the type of object or variant /** Can be printed with `getString()` \sa CommonType_ */ template Type type( const T& elem ) { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( trait::IsVariant::value ) return std::visit( fct::TypeFunct{}, elem ); else #endif return elem.type(); } //------------------------------------------------------------------ /// Returns the underlying data type of object or variant /** Can be printed with `getString()` */ template Dtype dtype( const T& elem ) { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( trait::IsVariant::value ) return std::visit( fct::DTypeFunct{}, elem ); else #endif return elem.dtype(); } //------------------------------------------------------------------ /// Returns length of element or variant (free function) template HOMOG2D_INUMTYPE length( const T& elem ) { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( trait::IsVariant::value ) return std::visit( fct::LengthFunct{}, elem ); else #endif return elem.length(); } //------------------------------------------------------------------ /// Returns area of element or variant (free function) template HOMOG2D_INUMTYPE area( const T& elem ) { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( trait::IsVariant::value ) return std::visit( fct::AreaFunct{}, elem ); else #endif return elem.area(); } #ifdef HOMOG2D_ENABLE_VRTP namespace priv { /// Get Bounding Box for a container holding variant objects /** \note Here, we DO NOT preallocate the vector, because if the container holds lines, then no point pair will be added for these elements. If we would have preallocate with 2 x size of input vector, then some default-constructed points could sneak in, and introduce an error. \todo 20240513: At present, this is only implemented for std::vector. Let it handle std::list and std::array. */ template auto getBB_CommonType( const std::vector>& v_var ) { HOMOG2D_START; HOMOG2D_DEBUG_ASSERT( v_var.size(), "cannot compute bounding box of empty set of variant" ); std::vector> vpts; vpts.reserve( v_var.size()*2 ); PointPair1_ ppair; for( const auto& elem: v_var ) { try { ppair = std::visit( fct::PtPairFunct{}, elem ); } catch( const std::exception& err ) { HOMOG2D_LOG_WARNING( "unable to compute point pair\n -msg=" << err.what() ); } vpts.push_back( ppair.first ); vpts.push_back( ppair.second ); } return FRect_( getBB_Points( vpts ) ); } } // namespace priv #endif // HOMOG2D_ENABLE_VRTP //------------------------------------------------------------------ /// Return Bounding Box of primitive or container holding primitives (free function) /** tests: [BB-cont] */ template FRect_ getBB( const T& t ) { HOMOG2D_START; if constexpr( !trait::IsContainer::value ) // if not a container, return t.getBB(); // then call the member function else { if( t.empty() ) HOMOG2D_THROW_ERROR_1( "unable, can't compute BB of empty container" ); if constexpr( trait::IsPoint::value ) { if( t.size() < 2 ) HOMOG2D_THROW_ERROR_1( "unable, need at least two points" ); return FRect_( priv::getBB_Points( t ) ); } else { if constexpr( trait::HasBB::value ) { using FPT = typename T::value_type::FType; std::vector> v_bb( t.size() ); auto it = v_bb.begin(); for( const auto& elem: t ) // compute bounding box of each element *it++ = elem.getBB(); return priv::getBB_FRect( v_bb ); // compute BB of all the BB } else { if constexpr( trait::IsSegment::value ) { if( t.empty() ) HOMOG2D_THROW_ERROR_1( "unable, need at least one segment" ); return priv::getBB_Segments( t ); } else { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( !trait::IsVariant::value ) { HOMOG2D_THROW_ERROR_1( "Unable, cannot compute BoundingBox of a container holding Line2d objects" ); } else return priv::getBB_CommonType( t ); #else HOMOG2D_THROW_ERROR_1( "Unable, cannot compute BoundingBox of a container holding Line2d objects" ); #endif } } } } } #ifdef HOMOG2D_ENABLE_VRTP /// Apply homography to primitive /** \warning The floating-point type of the returned object (variant) will be the one of the homography \c h, NOT the one of the input element. */ template CommonType_ transform( const Homogr_& h, const T& elem ) { if constexpr( trait::IsVariant::value ) return std::visit( fct::TransformFunct(h), elem ); else return h * elem; } #endif ///////////////////////////////////////////////////////////////////////////// // SECTION - FREE FUNCTIONS ///////////////////////////////////////////////////////////////////////////// /// Returns a pair of points holding min and max coordinates of the two input pair of points /** Tests: search for [minmax] */ template auto getMinMax( const PointPair2_& pp1, const PointPair2_& pp2 ) { HOMOG2D_START; std::array,4> arr; arr[0] = pp1.first; arr[1] = pp2.first; arr[2] = pp1.second; arr[3] = pp2.second; return priv::getBB_Points( arr ); } //------------------------------------------------------------------ /// Returns bounding box of two pairs of points /** \warning May throw in case of identical coordinates (for example with [0,0]-[1,0] and [10,0]-[20,0] */ template auto getBB( const PointPair2_& pp1, const PointPair2_& pp2 ) { HOMOG2D_START; return FRect_( getMinMax( pp1, pp2 ) ); } //------------------------------------------------------------------ /// Overload 1/3. This one is called if NONE of the args are a Line2d template< typename T1, typename T2, typename std::enable_if< (!std::is_same>::value && !std::is_same>::value) ,T1 >::type* = nullptr > auto getBB( const T1& elem1, const T2& elem2 ) { HOMOG2D_START; FRect_ out; try { out = getBB( ppair::getPointPair( elem1 ), ppair::getPointPair( elem2 ) ); } catch( const std::exception& err ) { HOMOG2D_THROW_ERROR_1( "unable to compute bounding box:\n -arg1=" << elem1 << "\n -arg2=" << elem2 << "\n -err=" << err.what() ); } return out; } /// Overload 2/3. Called if template< typename T1, typename T2, typename PLT1, typename PLT2, typename std::enable_if< ( std::is_same>::value || std::is_same>::value ) ,T1 >::type* = nullptr > auto getBB( const T1& p1, const T2& p2 ) { HOMOG2D_START; if( p1.size() == 0 && p2.size() == 0 ) HOMOG2D_THROW_ERROR_1( "unable to compute bounding box, both polylines are empty" ); if( p1.size() != 0 && p2.size() == 0 ) return FRect_( ppair::getPointPair( p1 ) ); if( p1.size() == 0 && p2.size() != 0 ) return FRect_( ppair::getPointPair( p2 ) ); return getBB( ppair::getPointPair( p1 ), ppair::getPointPair( p2 ) ); } //------------------------------------------------------------------ /// Overload 3/3. Called if one of the args is a Line2d (=> no build!) template< typename T1, typename T2, typename std::enable_if< (std::is_same>::value || std::is_same>::value) ,T1 >::type* = nullptr > auto getBB( const T1&, const T2& ) { HOMOG2D_START; static_assert( detail::AlwaysFalse::value, "fallback: undefined function" ); return FRect_(); // to avoid a compile warning } //------------------------------------------------------------------ // forward declaration namespace priv { template class ClosestPoints; } // forward declaration template priv::ClosestPoints getClosestPoints( const base::PolylineBase& poly1, const base::PolylineBase& poly2 ); namespace priv { /// Used in getClosestPoints() template class ClosestPoints { friend priv::ClosestPoints h2d::getClosestPoints<>( const base::PolylineBase&,const base::PolylineBase& ); private: size_t _pt1_min = 0; size_t _pt2_min = 0; HOMOG2D_INUMTYPE _minDist; Point2d_ _pt1; Point2d_ _pt2; const base::PolylineBase& _poly1; const base::PolylineBase& _poly2; ClosestPoints( const base::PolylineBase& poly1, const base::PolylineBase& poly2 ) : _poly1(poly1), _poly2(poly2) { _minDist = dist( poly1.getPoint(0), poly2.getPoint(0) ); } void store( HOMOG2D_INUMTYPE d, size_t idx1, size_t idx2 ) { _minDist = d; _pt1_min = idx1; _pt2_min = idx2; } public: PointPair1_ getPoints() const { return std::make_pair( _poly1.getPoint(_pt1_min), _poly2.getPoint(_pt2_min) ); } HOMOG2D_INUMTYPE getMinDist() const { return _minDist; } std::pair getIndexes() const { return std::make_pair( _pt1_min, _pt2_min ); } void print() const { std::cout << "ClosestPoints:\n _pt1_min=" << _pt1_min << " _pt2_min=" << _pt2_min << "\n_minDist=" << _minDist << "\npt1=" << _pt1 << " pt2=" << _pt2 << std::endl; } }; } // namespace priv //------------------------------------------------------------------ /// Computes the closest points between two polylines (types can be different) template priv::ClosestPoints getClosestPoints( const base::PolylineBase& poly1, const base::PolylineBase& poly2 ) { #ifndef HOMOG2D_NOCHECKS if( poly1.size() == 0 ) HOMOG2D_THROW_ERROR_1( "arg 1 is empty" ); if( poly2.size() == 0 ) HOMOG2D_THROW_ERROR_1( "arg 2 is empty" ); #endif priv::ClosestPoints out( poly1, poly2 ); for( size_t i=0; i auto getPts( const Segment_& seg ) { return seg.getPts(); } /// Returns Segment supporting line (free function) /// \sa Segment_::getLine() template Line2d_ getLine( const Segment_& seg ) { return seg.getLine(); } /// Returns Extended Segment (free function) /// \sa Segment_::geExtended() template Segment_ getExtended( const Segment_& seg ) { return seg.getExtended(); } /// Returns a pair of segments, corresponding to the input segment split by middle point (free function) /// \sa Segment_::split() template std::pair,Segment_> split( const Segment_& seg ) { return seg.split(); } /// Returns orthogonal segments /// \sa Segment_::getOrthogSegs() template std::array,4> getOrthogSegs( const Segment_& seg ) { return seg.getOrthogSegs(); } /// Returns orthogonal points /// \sa Segment_::getOrthogPts() template std::array,4> getOrthogPts( const Segment_& seg ) { return seg.getOrthogPts(); } //------------------------------------------------------------------ /// Free function, returns segment between two circle centers (or ellipse) template Segment_ getSegment( const T1& c1, const T2& c2 ) { return Segment_( c1.getCenter(), c2.getCenter() ); } /// Free function, returns line between two circle centers template Line2d_ getLine( const Circle_& c1, const Circle_& c2 ) { return Line2d_( c1.center(), c2.center() ); } /// Free function, returns middle point of segment /// \sa Segment_::getCenter() template Point2d_ getCenter( const Segment_& seg ) { return seg.getCenter(); } /// Free function, returns bisector line of segment /// \sa Segment_::getBisector() template Line2d_ getBisector( const Segment_& seg ) { return seg.getBisector(); } /// Free function, returns a pair of parallel segments at a distance \c dist /// \sa Segment_::getBisector() template std::pair,Segment_> getParallelSegs( const Segment_& seg, T dist ) { return seg.getParallelSegs(dist); } /// Free function, returns middle points of set of segments/circles /** \sa Segment_::getCenter() \sa Center_::getCenter() - input: set of segments - output: set of points (same container) */ template< typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > auto getCenters( const T& vsegs ) { std::vector> vout( vsegs.size() ); auto it = std::begin( vout ); for( const auto& seg: vsegs ) *it++ = getCenter( seg ); return vout; } /// Free function, returns a set of lines from a set of segments /** \sa Segment_::getLine() */ template< typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > auto getLines( const T& vsegs ) { std::vector> vout( vsegs.size() ); auto it = std::begin( vout ); for( const auto& seg: vsegs ) *it++ = seg.getLine(); return vout; } /// Free function, returns segments of the rectangle /// \sa FRect_::getSegs() template std::array,4> getSegs( const FRect_& seg ) { return seg.getSegs(); } /// Free function, returns the pair of segments tangential to the two circles /** Sources: - https://math.stackexchange.com/questions/719758/ - https://math.stackexchange.com/a/211854/133647 - https://en.wikipedia.org/wiki/Tangent_lines_to_circles - https://gieseanw.wordpress.com/2012/09/12/finding-external-tangent-points-for-two-circles/ */ template std::pair,Segment_> getTanSegs( const Circle_& c1, const Circle_& c2 ) { #ifndef HOMOG2D_NOCHECKS if( c1 == c2 ) HOMOG2D_THROW_ERROR_1( "c1 and c2 identical" ); #endif // if same radius, return the two segments parallel to the one joining the centers if( homog2d_abs( c1.radius() - c2.radius() ) < thr::nullDistance() ) { Segment_ seg_center( c1.center(), c2.center() ); return seg_center.getParallelSegs( c1.radius() ); } // check which one is the smallest Circle_ cA = c1; Circle_ cB = c2; if( c1.radius() < c2.radius() ) std::swap( cA, cB ); auto h = dist( cA.center(), cB.center() ); auto theta = homog2d_asin( ( cA.radius() - cB.radius() ) / h ) ; auto hcost = h * homog2d_cos(theta); // get rotated lines at center of CB auto l0 = cA.center() * cB.center(); auto l1 = l0.getRotatedLine( cB.center(), theta ); auto l2 = l0.getRotatedLine( cB.center(), -theta ); // build segments by getting the opposite point auto ppts1 = l1.getPoints( cB.center(), hcost ); auto p1 = ppts1.first; if( ppts1.second.distTo( cA.center()) < p1.distTo( cA.center()) ) p1 = ppts1.second; auto ppts2 = l2.getPoints( cB.center(), hcost ); auto p2 = ppts2.first; if( ppts2.second.distTo( cA.center()) < p2.distTo( cA.center()) ) p2 = ppts2.second; // std::cout <<"p1=" << p1 << " p2=" << p2 << '\n'; auto seg1 = Segment_( p1, cB.center() ); auto seg2 = Segment_( p2, cB.center() ); auto psegs1 = seg1.getParallelSegs( cB.radius() ); if( psegs1.first.distTo(cA.center()) < psegs1.second.distTo(cA.center()) ) std::swap( psegs1.first, psegs1.second ); auto psegs2 = seg2.getParallelSegs( cB.radius() ); if( psegs2.second.distTo(cA.center()) < psegs2.first.distTo(cA.center()) ) std::swap( psegs2.first, psegs2.second ); return std::make_pair( psegs1.first, psegs2.second ); } /// Returns the 4 points of the rectangle (free function) /// \sa FRect_::get4Pts() template std::array,4> get4Pts( const FRect_& rect ) { return rect.get4Pts(); } /// Returns the 2 major points of the rectangle (free function) /// \sa FRect_::getPts() template auto getPts( const FRect_& rect ) { return rect.getPts(); } /// Free function template HOMOG2D_INUMTYPE height( const FRect_& rect ) { return rect.height(); } /// Free function template HOMOG2D_INUMTYPE width( const FRect_& rect ) { return rect.width(); } /// Free function template FRect_ getExtended( const FRect_& rect ) { return rect.getExtended(); } /// Free function template std::pair,Segment_> getDiagonals( const FRect_& rect ) { return rect.getDiagonals(); } /// Returns centroid of Polyline (free function) /// \sa PolylineBase::centroid() template base::LPBase centroid( const base::PolylineBase& pl ) { return pl.centroid(); } /// Returns reference on radius of circle (free function), non-const version /// \sa Circle_::radius() template FPT& radius( Circle_& cir ) { return cir.radius(); } /// Returns reference on radius of circle (free function), const version /// \sa Circle_::radius() template const FPT& radius( const Circle_& cir ) { return cir.radius(); } /// Returns reference on center of circle (free function), non-const version /// \sa Circle_::center() template Point2d_& center( Circle_& cir ) { return cir.center(); } /// Returns reference on center of circle (free function), const version /// \sa Circle_::center() template const Point2d_& center( const Circle_& cir ) { return cir.center(); } /// Get data size expressed as number of bits for, respectively, mantissa and exponent /// \sa detail::Common::dsize() template std::pair dsize( const T& t ) { return t.dsize(); } //------------------------------------------------------------------ /// Returns the distance between 2 parallel lines (free function) /** - ref: https://en.wikipedia.org/wiki/Distance_between_two_parallel_lines We first check that the two lines are indeed parallel. This should ensure that a1=a2 and b1=b2. But due to numeric issues they could by slightly different. In order to gain accuracy, we use the geometric mean of these. */ template HOMOG2D_INUMTYPE getParallelDistance( const Line2d_& li1, const Line2d_& li2 ) { #ifndef HOMOG2D_NOCHECKS if( !li1.isParallelTo(li2) ) HOMOG2D_THROW_ERROR_1( "lines are not parallel" ); #endif const auto ar1 = li1.get(); const auto ar2 = li2.get(); const HOMOG2D_INUMTYPE a1 = ar1[0]; const HOMOG2D_INUMTYPE b1 = ar1[1]; const HOMOG2D_INUMTYPE c1 = ar1[2]; const HOMOG2D_INUMTYPE a2 = ar2[0]; const HOMOG2D_INUMTYPE b2 = ar2[1]; const HOMOG2D_INUMTYPE c2 = ar2[2]; HOMOG2D_INUMTYPE a = homog2d_sqrt( a1*a2 ); HOMOG2D_INUMTYPE b = homog2d_sqrt( b1*b2 ); return homog2d_abs( c1 - c2 ) / homog2d_sqrt( a*a + b*b ); } /// Return angle of ellipse (free function) /// \sa Ellipse_::angle() template HOMOG2D_INUMTYPE angle( const Ellipse_& ell ) { return ell.angle(); } /// Return center of Segment_, FRect_, or Ellipse_ (free function) /// \sa Ellipse_::getCenter() /// \sa FRect_::getCenter() /// \sa Segment_::getCenter() template Point2d_ getCenter(const T& other ) { return other.getCenter(); } /// Translate primitive \c prim (free function) /** Calls the member function. Type checking is done there. */ template void translate( T& prim, FP1 dx, FP2 dy ) { prim.translate( dx, dy ); } /// Translate primitive \c prim with values in a \c std::pair (free function) /** Calls the member function. Type checking is done there. */ template void translate( T& prim, const std::pair& ppt ) { prim.translate( ppt.first, ppt.second ); } /// Move primitive to other location (free function) /** Calls the member function. Type checking is done there. */ template void moveTo( T& prim, const Point2d_& pt ) { prim.moveTo( pt ); } /// Move primitive to other location (free function) /** Calls the member function. Type checking is done there. */ template void moveTo( T& prim, FP1 x, FP2 y ) { prim.moveTo( x, y ); } /// Returns true if ellipse is a circle /// \sa Ellipse_::isCircle() template bool isCircle( const Ellipse_& ell, HOMOG2D_INUMTYPE thres=1.E-10 ) { return ell.isCircle( thres ); } /// Returns ellipse axis lines /// \sa Ellipse_::getAxisLines() template std::pair,Line2d_> getAxisLines( const Ellipse_& ell ) { return ell.getAxisLines(); } //------------------------------------------------------------------ /// Returns index of point in container \c cont that is the nearest to \c pt /// \todo add some sfinae and/or checking on type T /// \todo perform some speed analysis, and check if usage of squared distance is really better /// \todo add other distance computations (Manhattan?) /// \todo Would it be better to return an iterator? template size_t findNearestPoint( const Point2d_& pt, const T& cont ) { if( cont.empty() ) HOMOG2D_THROW_ERROR_1( "container is empty" ); auto minDist = priv::sqDist( pt, cont[0] ); size_t resIdx = 0; for( size_t i=1; i size_t findFarthestPoint( const Point2d_& pt, const T& cont ) { if( cont.empty() ) HOMOG2D_THROW_ERROR_1( "container is empty" ); auto maxDist = priv::sqDist( pt, cont[0] ); size_t resIdx = 0; for( size_t i=1; i maxDist ) { resIdx = i; maxDist = currentDist; } } return resIdx; } //------------------------------------------------------------------ /// Returns index of point in container \c cont that is the farthest to \c pt template auto findNearestFarthestPoint( const Point2d_& pt, const T& cont ) { if( cont.empty() ) HOMOG2D_THROW_ERROR_1( "container is empty" ); auto maxDist = priv::sqDist( pt, cont[0] ); auto minDist = maxDist; size_t idxMin = 0; size_t idxMax = 0; for( size_t i=1; i maxDist ) { idxMax = i; maxDist = currentDist; } if( currentDist < minDist ) { idxMin = i; minDist = currentDist; } } return std::make_pair(idxMin, idxMax); } ///////////////////////////////////////////////////////////////////////////// // SECTION - GENERIC DRAWING FREE FUNCTIONS (BACK-END INDEPENDENT) ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Free function, draws any of the primitives by calling the relevant member function /** \todo 20240504: see note on the IsDrawable trait class */ template< typename U, typename Prim, typename std::enable_if< trait::IsDrawable::value, Prim >::type* = nullptr > void draw( img::Image& img, const Prim& prim, const img::DrawParams& dp=img::DrawParams() ) { prim.draw( img, dp ); } /// Free function, draws text \c str at position \c pt template void drawText( img::Image& im, std::string str, Point2d_ pt, img::DrawParams dp=img::DrawParams() ) { im.drawText( str, pt, dp ); } namespace priv { /// Draw indexes for points template void impl_drawIndexes( img::Image& img, size_t c, const img::DrawParams& dp, const Point2d_& pt ) { if( dp._dpValues._showIndex ) drawText( img, std::to_string(c), pt, dp ); } /// Draw indexes for segment template void impl_drawIndexes( img::Image& img, size_t c, const img::DrawParams& dp, const Segment_& seg ) { if( dp._dpValues._showIndex ) drawText( img, std::to_string(c), seg.getCenter(), dp ); } /// Default signature, will be instanciated if no other fits (and does nothing) template void impl_drawIndexes( img::Image&, size_t, const img::DrawParams&, const DUMMY& ) {} } // namespace priv /// Free function, draws a set of primitives /** Type \c T can be \c std::array or \c std::vector, with \c type being anything drawable */ template< typename U, typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > void draw( img::Image& img, const T& cont, const img::DrawParams& dp=img::DrawParams() ) { size_t c=0; for( const auto& elem: cont ) { elem.draw( img, dp ); priv::impl_drawIndexes( img, c++, dp, elem ); } } /// This version holds a \c std::function as 3th parameter. It can be used to pass a function /// that will return a different img::DrawParams for a given index of the container. template< typename U, typename T, typename std::enable_if< trait::IsContainer::value, T >::type* = nullptr > void draw( img::Image& img, const T& cont, std::function& func ) { size_t c=0; for( const auto& elem: cont ) { auto dp = func(c); elem.draw( img, dp ); priv::impl_drawIndexes( img, c++, dp, elem ); } } //------------------------------------------------------------------ /// Free function, draws a pair of objects /** Type \c T1 and \c T2 can be anything drawable */ template void draw( img::Image& img, const std::pair& pa, const img::DrawParams& dp=img::DrawParams() ) { pa.first.draw( img, dp ); pa.second.draw( img, dp ); } namespace detail { /// Private helper function, used by LPBase::draw(). Draw point on image. template void drawPt( img::Image& img, img::PtStyle ps, std::vector> vpt, const img::DrawParams& dp, bool drawDiag=false ) { auto delta = dp._dpValues._ptDelta; auto delta2 = std::round( 0.7 * delta); switch( ps ) { case img::PtStyle::Times: vpt[0].translate( -delta2, +delta2 ); vpt[1].translate( +delta2, -delta2 ); vpt[2].translate( +delta2, +delta2 ); vpt[3].translate( -delta2, -delta2 ); break; case img::PtStyle::Plus: case img::PtStyle::Diam: vpt[0].translate( -delta, 0. ); vpt[1].translate( +delta, 0. ); vpt[2].translate( 0., -delta ); vpt[3].translate( 0., +delta ); break; default: assert(0); } auto dp2(dp); // we need to do this, because this is called by the dp2.showPoints(false); // segment drawing function. If not, infinite recursion if( !drawDiag ) { Segment_ s1( vpt[0], vpt[1] ); Segment_ s2( vpt[2], vpt[3] ); s1.draw( img, dp2 ); s2.draw( img, dp2 ); } else // draw 4 diagonal lines { Segment_( vpt[0], vpt[2] ).draw( img, dp2 ); Segment_( vpt[2], vpt[1] ).draw( img, dp2 ); Segment_( vpt[1], vpt[3] ).draw( img, dp2 ); Segment_( vpt[0], vpt[3] ).draw( img, dp2 ); } } } // namespace detail //------------------------------------------------------------------ namespace priv { /// Holds convex hull code namespace chull { //------------------------------------------------------------------ /// Used int the convex hull algorithm template size_t getPivotPoint( const std::vector>& in ) { auto pmin = h2d::priv::getBmPoint_helper( in ); return static_cast( pmin - in.begin() ); } //------------------------------------------------------------------ /// Sorts points by angle between the lines with horizontal axis template std::vector sortPoints( const std::vector>& in, size_t piv_idx ) { assert( in.size()>3 ); // step 1: create new vector holding the indexes of the points, including the pivot point (will be in first position) std::vector out( in.size() ); std::iota( out.begin(), out.end(), 0 ); // fill vector: [0,1,2,...] std::swap( out[piv_idx], out[0] ); auto pt0 = in[piv_idx]; // step 2: sort points by angle of lines between the current point and pivot point std::sort( out.begin()+1, out.end(), [&] // lambda ( size_t i1, size_t i2 ) { auto pt1 = in[i1]; auto pt2 = in[i2]; auto dx1 = pt1.getX() - pt0.getX(); auto dy1 = pt1.getY() - pt0.getY(); auto dx2 = pt2.getX() - pt0.getX(); auto dy2 = pt2.getY() - pt0.getY(); return ((dx1 * dy2 - dx2 * dy1) > 0); } ); return out; } //------------------------------------------------------------------ /// To find orientation of ordered triplet of points (p, q, r). /** The function returns following values - 0 --> p, q and r are colinear - 1 --> Clockwise - 2 --> Counterclockwise \todo 20240326: this is subject to numerical instability, as it is based on differences. \todo 20230212: replace const value HOMOG2D_THR_ZERO_DETER with related static function */ template int orientation( Point2d_ p, Point2d_ q, Point2d_ r ) { HOMOG2D_INUMTYPE px = p.getX(); HOMOG2D_INUMTYPE py = p.getY(); HOMOG2D_INUMTYPE qx = q.getX(); HOMOG2D_INUMTYPE qy = q.getY(); HOMOG2D_INUMTYPE rx = r.getX(); HOMOG2D_INUMTYPE ry = r.getY(); auto val = (qy - py) * (rx - qx) - (qx - px) * (ry - qy); if( homog2d_abs(val) < HOMOG2D_THR_ZERO_DETER ) return 0; // collinear return (val > 0 ? 1 : -1 ); // clock or counterclock wise } //------------------------------------------------------------------ /// Inherits std::stack<> and adds a member function to fetch the underlying std::vector. /// Used in h2d::convexHull() struct Mystack : std::stack> { const std::vector& getVect() const { return this->c; } }; //------------------------------------------------------------------ } // namespace chull } // namespace priv //------------------------------------------------------------------ /// Compute Convex Hull of a Polyline (free function) /** - type \c T: can be either OPolyline, CPolyline, or std::vector - Graham scan algorithm: https://en.wikipedia.org/wiki/Graham_scan */ template CPolyline_ convexHull( const base::PolylineBase& input ) { return convexHull( input.getPts() ); } namespace base { //------------------------------------------------------------------ /// Computes and returns the convex hull of a set of points (free function) /** - Graham scan algorithm: https://en.wikipedia.org/wiki/Graham_scan \todo 20230728: make this function accept also std::array and std::list (using Sfinae alogn with trait::IsContainer) */ template CPolyline_ convexHull( const std::vector>& input ) { if( input.size() < 4 ) // if 3 pts or less, then the hull is equal to input set return CPolyline_( input ); // step 1: find pivot (point with smallest Y coord) auto pivot_idx = priv::chull::getPivotPoint( input ); // step 2: sort points by angle of lines between the current point and pivot point auto v2 = priv::chull::sortPoints( input, pivot_idx ); std::stack firstPoint; priv::chull::Mystack hull; hull.push( v2[0] ); hull.push( v2[1] ); hull.push( v2[2] ); firstPoint.push( 0 ); firstPoint.push( 1 ); // step 3: iterate size_t idx1 = 1; size_t idx2 = 2; size_t idx3 = 3; do { // HOMOG2D_LOG( "** loop start, idx1=" << idx1 << ", idx2=" << idx2 << ", idx3=" << idx3 << ", hull " << hull.getVect() ); auto p = input.at( v2[idx1] ); auto q = input.at( v2[idx2] ); auto r = input.at( v2[idx3] ); auto orient = priv::chull::orientation( p, q, r ); // HOMOG2D_LOG( "considering pts: " << v2[idx1] << "," << v2[idx2] << "," << v2[idx3] << ": or = " << orient ); if( orient != 1 ) { hull.push( v2[idx3] ); idx1 = idx2; idx2 = idx3; idx3++; firstPoint.push( idx1 ); } else { hull.pop(); idx2=idx1; // idx3 stays the same firstPoint.pop(); idx1 = firstPoint.top(); } } while( idx3 < v2.size() ); // final step: copy hull indexes to vector of points const auto v = hull.getVect(); std::vector> vout( v.size() ); std::transform( v.begin(), v.end(), vout.begin(), [&input] // lambda (size_t idx) { return input.at( idx ); } ); return CPolyline_( vout ); } } // namespace base /// Return convex hull (member function implementation) template CPolyline_ base::PolylineBase::convexHull() const { return h2d::convexHull( *this ); } ///////////////////////////////////////////////////////////////////////////// // SECTION - OPENCV BINDING - GENERAL ///////////////////////////////////////////////////////////////////////////// #ifdef HOMOG2D_USE_OPENCV //------------------------------------------------------------------ /// Copy matrix to Opencv \c cv::Mat /** The output matrix is passed by reference to avoid issues with Opencv copy operator, and is allocated here. User can pass a type as second argument: CV_32F for \c float, CV_64F for \c double (default) */ template void Hmatrix_::copyTo( cv::Mat& mat, int type ) const { auto& data = detail::Matrix_::_mdata; #ifndef HOMOG2D_NOCHECKS if( type != CV_64F && type != CV_32F ) throw std::runtime_error( "invalid OpenCv matrix type" ); #endif mat.create( 3, 3, type ); // default:CV_64F size_t i=0; switch( type ) { case CV_64F: for( auto it = mat.begin(); it != mat.end(); it++, i++ ) *it = data[i/3][i%3]; break; case CV_32F: for( auto it = mat.begin(); it != mat.end(); it++, i++ ) *it = data[i/3][i%3]; break; default: assert(0); } } //------------------------------------------------------------------ /// Get homography from Opencv \c cv::Mat template Hmatrix_& Hmatrix_::operator = ( const cv::Mat& mat ) { auto type = mat.type(); #ifndef HOMOG2D_NOCHECKS if( mat.rows != 3 || mat.cols != 3 ) throw std::runtime_error( "invalid matrix size, rows=" + std::to_string(mat.rows) + " cols=" + std::to_string(mat.cols) ); if( mat.channels() != 1 ) throw std::runtime_error( "invalid matrix nb channels: " + std::to_string(mat.channels() ) ); if( type != CV_64F && type != CV_32F ) throw std::runtime_error( "invalid matrix type" ); #endif size_t i=0; auto& data = detail::Matrix_::_mdata; switch( type ) { case CV_64F: for( auto it = mat.begin(); it != mat.end(); it++, i++ ) data[i/3][i%3] = *it; break; case CV_32F: for( auto it = mat.begin(); it != mat.end(); it++, i++ ) data[i/3][i%3] = *it; break; default: assert(0); } return *this; } #endif // HOMOG2D_USE_OPENCV ///////////////////////////////////////////////////////////////////////////// // SECTION .1 - CLASS DRAWING MEMBER FUNCTIONS (backend-agnostic) ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Draw Line2d on image, backend independent /** Steps: -# find the intersection points between the line and the image rectangle, should find 2. (but 1 is possible) -# draw a line between these 2 points */ template template void base::LPBase::impl_draw_LP( img::Image& im, img::DrawParams dp, const detail::BaseHelper& ) const { assert( im.rows() > 2 ); assert( im.cols() > 2 ); Point2d_ pt1; // 0,0 Point2d_ pt2( im.cols()-1, im.rows()-1 ); auto ri = this->intersects( pt1, pt2 ); if( ri() ) { if( ri.size() == 2 ) // if only one intersection point, do nothing { auto ppts = ri.get(); h2d::Segment_ seg( ppts[0], ppts[1] ); dp.showPoints(false); seg.draw( im, dp ); } } } //------------------------------------------------------------------ /// Draw points on image implementation, backend independent. /// Returns false if point not in image template template void base::LPBase::impl_draw_LP( img::Image& im, img::DrawParams dp, const detail::BaseHelper& ///< dummy arg ) const { if( getX()<0 || getX()>=im.cols() ) // check if point is in image return; if( getY()<0 || getY()>=im.rows() ) return; std::vector> vpt( 4, *this ); switch( dp._dpValues._ptStyle ) { case img::PtStyle::Dot: Circle_( *this, dp._dpValues._pointSize ).draw( im, dp ); break; case img::PtStyle::Plus: // "+" symbol detail::drawPt( im, img::PtStyle::Plus, vpt, dp ); break; case img::PtStyle::Star: detail::drawPt( im, img::PtStyle::Plus, vpt, dp ); detail::drawPt( im, img::PtStyle::Times, vpt, dp ); break; case img::PtStyle::Diam: detail::drawPt( im, img::PtStyle::Plus, vpt, dp, true ); break; case img::PtStyle::Times: // "times" symbol detail::drawPt( im, img::PtStyle::Times, vpt, dp ); break; default: assert(0); } } ///////////////////////////////////////////////////////////////////////////// // SECTION .2 - CLASS DRAWING MEMBER FUNCTIONS (OpenCv) ///////////////////////////////////////////////////////////////////////////// #ifdef HOMOG2D_USE_OPENCV #if 0 // DEPRECATED (???) //------------------------------------------------------------------ /// Extension to PolylineBase::draw(), to draw point indexes /** (Opencv-dependent) */ template template void base::PolylineBase::impl_draw_pl( img::Image& im ) const { for( size_t i=0; i void FRect_::draw( img::Image& im, img::DrawParams dp ) const { cv::rectangle( im.getReal(), _ptR1.getCvPti(), _ptR2.getCvPti(), dp.cvColor(), dp._dpValues._lineThickness, dp._dpValues._lineType==1?cv::LINE_AA:cv::LINE_8 ); if( dp._dpValues._showPoints ) { _ptR1.draw( im, dp ); _ptR2.draw( im, dp ); } } //------------------------------------------------------------------ /// Draw \c Segment (Opencv implementation) template void Segment_::draw( img::Image& im, img::DrawParams dp ) const { cv::line( im.getReal(), _ptS1.getCvPtd(), _ptS2.getCvPtd(), dp.cvColor(), dp._dpValues._lineThickness, dp._dpValues._lineType==1?cv::LINE_AA:cv::LINE_8 ); if( dp._dpValues._showPoints ) { _ptS1.draw( im, dp ); _ptS2.draw( im, dp ); } } //------------------------------------------------------------------ /// Draw \c Circle (Opencv implementation) template void Circle_::draw( img::Image& im, img::DrawParams dp ) const { cv::circle( im.getReal(), _center.getCvPti(), static_cast(_radius), dp.cvColor(), dp._dpValues._lineThickness, dp._dpValues._lineType==1?cv::LINE_AA:cv::LINE_8 ); if( dp._dpValues._showPoints ) // WARNING: can't use the "Dot" style here, because it call this function, so infinite recursion _center.draw( im, dp.setPointStyle( img::PtStyle::Plus) ); } //------------------------------------------------------------------ /// Draw \c Ellipse (Opencv implementation) /** - see https://docs.opencv.org/3.4/d6/d6e/group__imgproc__draw.html#ga28b2267d35786f5f890ca167236cbc69 */ template void Ellipse_::draw( img::Image& im, img::DrawParams dp ) const { auto par = p_getParams(); cv::ellipse( im.getReal(), cv::Point( par.x0,par.y0 ), cv::Size( par.a, par.b ), par.theta*180./M_PI, 0., 360., dp.cvColor(), dp._dpValues._lineThickness, dp._dpValues._lineType==1?cv::LINE_AA:cv::LINE_8 ); } //------------------------------------------------------------------ namespace base { /// Draw PolylineBase (Opencv implementation) template void PolylineBase::draw( img::Image& im, img::DrawParams dp ) const { if( size() < 2 ) // nothing to draw return; for( size_t i=0; i inline void img::Image::drawText( std::string str, Point2d_ pt, img::DrawParams dp ) { getReal()._svgString << "" << str << "\n"; } #ifdef HOMOG2D_USE_OPENCV /// Free function, draw text on Opencv image template <> inline void img::Image::drawText( std::string str, Point2d_ pt, img::DrawParams dp ) //=img::DrawParams() ) { auto col = dp.color(); cv::putText( getReal(), str, pt.getCvPtd(), cv::FONT_HERSHEY_SIMPLEX, // font id, see https://docs.opencv.org/4.7.0/ 0.03 * dp._dpValues._fontSize, // scale factor (approximate...) cv::Scalar( col.b, col.g, col.r ) ); } #endif //------------------------------------------------------------------ /// Draw \c Circle (SVG implementation) template void Circle_::draw( img::Image& im, img::DrawParams dp ) const { im.getReal()._svgString << "\n"; } //------------------------------------------------------------------ /// Draw \c Ellipse (SVG implementation) template void Ellipse_::draw( img::Image& im, img::DrawParams dp ) const { im.getReal()._svgString << "\n"; } //------------------------------------------------------------------ /// Draw \c FRect (SVG implementation) template void FRect_::draw( img::Image& im, img::DrawParams dp ) const { im.getReal()._svgString << "\n"; } //------------------------------------------------------------------ /// Draw \c Segment (SVG implementation) template void Segment_::draw( img::Image& im, img::DrawParams dp ) const { if( dp._dpValues._showPoints ) im.getReal()._svgString << "\n"; auto pts = getPts(); im.getReal()._svgString << "\n"; if( dp._dpValues._showPoints ) { _ptS1.draw( im, dp ); _ptS2.draw( im, dp ); im.getReal()._svgString << "\n"; } } //------------------------------------------------------------------ #ifdef HOMOG2D_ENABLE_PRTP namespace rtp { /// Stream operator for \c Root type /** \todo replace this by a call to a virtual function `print()` (that needs to be defined in all the child classes as: \code void print( std::ostream& f ) { f << *this; } \endcode */ std::ostream& operator << ( std::ostream& f, const Root& p ) { f << "type="<< getString( p.type() ) << std::endl; switch( p.type() ) { case Type::Circle: { const Circle_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::Ellipse: { const Ellipse_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::FRect: { const FRect_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::Line2d: { const Line2d_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::Point2d: { const Point2d_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::Segment: { const Segment_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::OPolyline: { const OPolyline_* p2 = static_cast*>( &p ); f << *p2; } break; case Type::CPolyline: { const CPolyline_* p2 = static_cast*>( &p ); f << *p2; } break; default: assert(0); } return f; } } // namespace rtp #endif // HOMOG2D_ENABLE_PRTP //------------------------------------------------------------------ namespace base { /// Draw Polyline (SVG implementation) /** \note To show the points index, we don't use the svg "marker-start/marker-mid/marker-end" syntax so that the dots always have the same color as the segments */ template void PolylineBase::draw( img::Image& im, img::DrawParams dp ) const { if( size() < 2 ) // nothing to draw return; if( dp._dpValues._showIndex || dp._dpValues._showPoints ) im.getReal()._svgString << "\n"; im.getReal()._svgString << '<' << (isClosed() ? "polygon" : "polyline") << " stroke=\"" << dp.getSvgRgbColor() << "\" stroke-width=\"" << dp._dpValues._lineThickness << "\" "; if( !dp.holdsFill() ) im.getReal()._svgString << "fill=\"none\" "; im.getReal()._svgString << dp.getAttrString() << "points=\""; for( const auto& pt: getPts() ) im.getReal()._svgString << pt.getX() << ',' << pt.getY() << ' '; im.getReal()._svgString << "\"/>\n"; if( dp._dpValues._showIndex ) { im.getReal()._svgString << "\n"; for( size_t i=0; i" << i << "\n"; } im.getReal()._svgString << "\n"; } if( dp._dpValues._showPoints ) { im.getReal()._svgString << "\n"; for( size_t i=0; i\n"; } if( dp._dpValues._showIndex || dp._dpValues._showPoints ) im.getReal()._svgString << "\n"; } } // namespace base ///////////////////////////////////////////////////////////////////////////// // SECTION - TYPEDEFS ///////////////////////////////////////////////////////////////////////////// /// Default line type, uses \c double as numerical type using Line2d = Line2d_; /// Default point type, uses \c double as numerical type //using Point2d = base::LPBase; using Point2d = Point2d_; /// Default homography (3x3 matrix) type, uses \c double as numerical type using Homogr = Homogr_; /// Default homogeneous matrix, uses \c double as numerical type using Epipmat = Hmatrix_; /// Default segment type using Segment = Segment_; /// Default circle type using Circle = Circle_; /// Default rectangle type using FRect = FRect_; /// Default polyline type using CPolyline = base::PolylineBase; using OPolyline = base::PolylineBase; /// Default ellipse type using Ellipse = Ellipse_; // float types using Line2dF = Line2d_; using Point2dF = Point2d_; using HomogrF = Homogr_; using SegmentF = Segment_; using CircleF = Circle_; using FRectF = FRect_; using EllipseF = Ellipse_; // double types using Line2dD = Line2d_; using Point2dD = Point2d_; using HomogrD = Homogr_; using SegmentD = Segment_; using CircleD = Circle_; using FRectD = FRect_; using EllipseD = Ellipse_; // long double types using Line2dL = Line2d_; using Point2dL = Point2d_; using HomogrL = Homogr_; using SegmentL = Segment_; using CircleL = Circle_; using FRectL = FRect_; using EllipseL = Ellipse_; using CPolylineF = base::PolylineBase; using CPolylineD = base::PolylineBase; using CPolylineL = base::PolylineBase; using OPolylineF = base::PolylineBase; using OPolylineD = base::PolylineBase; using OPolylineL = base::PolylineBase; #ifdef HOMOG2D_ENABLE_VRTP // variant type using CommonType = CommonType_; using CommonTypeF = CommonType_; using CommonTypeD = CommonType_; using CommonTypeL = CommonType_; #endif #ifdef HOMOG2D_USE_SVG_IMPORT /// Holds stuff related to SVG import namespace svg { /// Private functions related to the SVG import code namespace svgp { //------------------------------------------------------------------- /// General string tokenizer, taken from http://stackoverflow.com/a/236803/193789 /** - see also this one: http://stackoverflow.com/a/53878/193789 */ inline std::vector tokenize( const std::string &s, char delim ) { std::vector velems; // std::stringstream ss( TrimSpaces(s) ); std::stringstream ss( s ); std::string item; while( std::getline( ss, item, delim ) ) velems.push_back(item); return velems; } //------------------------------------------------------------------ /// Importing rotated ellipse from SVG data inline std::pair,HOMOG2D_INUMTYPE> getEllipseRotateAttr( const char* rot_str ) { std::string s(rot_str); // std::cout << __FUNCTION__ << "(): " << s << "\n"; auto v1 = tokenize( s, '(' ); if( v1.size() == 2 ) if( v1[0] == "rotate" ) { auto v2 = v1[1].substr(0, v1[1].size() - 1 ); // std::cout << __FUNCTION__ << "(): v2=" << v2 << "\n"; auto v3 = tokenize( v2, ',' ); if( v3.size() == 3 ) { try // in case of an incorrect numerical string { auto angle = std::stod( v3[0] ); auto x0 = std::stod( v3[1] ); auto y0 = std::stod( v3[2] ); return std::make_pair( Point2d_(x0,y0),angle*M_PI/180. ); } catch( std::exception& err ) { HOMOG2D_THROW_ERROR_2( "invalid 'transform' attribute for svg ellipse import", err.what() ); } } } HOMOG2D_THROW_ERROR_1( "invalid 'transform' attribute for svg ellipse import" ); } //------------------------------------------------------------------ /// Svg import: Basic parsing of points that are in the format "10,20 30,40 50,60" /// \todo 20240326: this is used to import SVG polygon type. Maybe this can be replaced by /// the "path" import code? inline std::vector parsePoints( const char* pts ) { std::vector out; std::string s(pts); // std::cout << "processing " << s << '\n'; // trimString( s ); auto v1 = tokenize( s, ' ' ); for( const auto& pt: v1 ) { auto v2 = tokenize( pt, ',' ); if( v2.size() != 2 ) throw "h2d:img::svg: invalid point format in importing svg element: " + s; auto x = std::stod( v2[0] ); auto y = std::stod( v2[1] ); out.emplace_back( Point2d(x,y) ); } return out; } inline bool isDigit( char c ) { return ( (c >= '0' && c <= '9') || c == '.' || c == '-' ); } inline bool svgPathCommandIsAllowed(char c) { if( c == 'M' || c == 'L' || c == 'H'|| c == 'V' || c == 'Z' ) return true; return false; } /// Get next element in svg path string. inline std::string getNextElem( const std::string& str, std::string::const_iterator& it ) { std::string out; while( *it == ' ' || *it == ',' ) // skip spaces and commas it++; if( it >= str.cend() ) return out; // return empty string if( !isDigit(*it) ) // alpha character { out.push_back( *it ); it++; } else { while( isDigit(*it) && it < str.cend() ) { out.push_back( *it ); it++; } // std::cout << "getNextElem() DIGITS:-" << out << "- #=" << out.size() << '\n'; // return out; } return out; } /// Returns nb of expected values for a given SVG path command. /// Ref: https://www.w3.org/TR/SVG2/paths.html inline std::map& numberValues() { static std::map nbval; nbval['M'] = 2; // M-m nbval['L'] = 2; // L-l nbval['H'] = 1; // H-h nbval['V'] = 1; // V-v nbval['C'] = 6; // C-c !NOT HANDLED! nbval['S'] = 4; // S-s !NOT HANDLED! nbval['Q'] = 4; // Q-q !NOT HANDLED! nbval['T'] = 2; // T t !NOT HANDLED! nbval['A'] = 7; // A-a !NOT HANDLED! nbval['Z'] = 0; // Z-z return nbval; } enum class PathMode { Absolute, Relative }; /// Holds the current SVG "path" command, and the number of required numerical values /// \sa SvgValuesBuffer struct SvgPathCommand { PathMode _absRel = PathMode::Absolute; char _command = 'M'; uint8_t _nbValues = 2; bool isAbsolute() const { return _absRel == PathMode::Absolute ? true : false; } void setCommand( char c ) { _command = c; _nbValues = numberValues().at(c); } bool isNewPolyline() const { if( _command == 'm' || _command == 'M' || _command == 'z' || _command == 'Z' ) return true; return false; } }; /// Generate new point from current mode and previous point, handles absolute/relative coordinates inline Point2d_ generateNewPoint( SvgPathCommand mode, ///< SVG path command Point2d_ prevPt, ///< previous point, is needed if relative mode const std::vector& val ///< numerical values that have been stored ) { // std::cout << "generateNewPoint(): command=" << mode._command // << std::hex << " hex=" << (int)mode._command << std::dec << '\n'; auto nb = val.size(); // priv::printVector( val,"generateNewPoint()" ); HOMOG2D_LOG( "abs/rel=" << (mode.isAbsolute()?"ABS":"REL") ); // some checking to lighten up the following code, maybe can be removed afterwards switch( mode._command ) { case 'M': case 'L': assert( nb == 2 ); break; case 'H': case 'V': assert( nb == 1 ); break; case 'Z': assert( nb == 0 ); break; default: assert(0); } Point2d_ out; switch( mode._command ) { case 'M': case 'L': if( mode.isAbsolute() ) out.set( val[0], val[1] ); else out.set( prevPt.getX() + val[0], prevPt.getY() + val[1] ); break; case 'H': // Horizontal if( mode.isAbsolute() ) out.set( val[0], prevPt.getY() ); else out.set( val[0] + prevPt.getX() , prevPt.getY() ); break; case 'V': // Vertical if( mode.isAbsolute() ) out.set( prevPt.getX(), val[0] ); else out.set( prevPt.getX(), val[0] + prevPt.getY() ); break; default: assert(0); } return out; } inline SvgPathCommand getCommand( char c ) { HOMOG2D_LOG( "search command for " << c ); static std::string commands( "MLHVCSQTAZ" ); // allowed SVG path commands (and their counterparts relative) SvgPathCommand out; std::string str( 1, c ); HOMOG2D_LOG( " str=" << str << " #=" << str.size() ); auto pos = commands.find( str ); // HOMOG2D_LOG( "pos=" << pos << " str=" << str ); bool invalid = false; if( pos == std::string::npos ) { if( c>='a' && c <= 'z' ) // check if lowercase { c = (char)((int)c+'A'-'a'); // std::cout << "lowercase, new c=" << c << '\n'; str[0] = c; // HOMOG2D_LOG( " str2=" << str << "#=" << str.size() ); pos = commands.find( str ); if( pos != std::string::npos ) out._absRel = PathMode::Relative; else invalid = true; } else invalid = true; } if( invalid ) HOMOG2D_THROW_ERROR_1( "Illegal character in SVG path element:-" << str << "- ascii=" << std::hex << str[0] << std::dec ); out.setCommand( commands[pos] ); // std::cout << "pos=" << pos << " _command=" << out._command << " _nbValues=" << (int)out._nbValues << '\n'; return out; } /// Removes dupes in set of points. Needed when importing SVG files using a "path" command, /// because sometimes they hold duplicates points, and that can't be in polylines template std::vector> purgeSetDupes( const std::vector>& pts ) { std::vector> out; out.reserve( pts.size() ); for( auto it=pts.begin(); it!=std::prev(pts.end()); it++ ) { const auto& elem = *it; const auto& next = *std::next(it); if( elem != next ) out.push_back( elem ); } out.push_back( pts.back() ); // add last one return out; } /// This will hold the values read from the SVG Path parsing code, before they are /// converted to points /// \sa SvgPathCommand struct SvgValuesBuffer { std::vector _values; Point2d _previousPt; size_t size() const { return _values.size(); } void storeValues( std::vector& out, SvgPathCommand mode ) { if( _values.size() != (size_t)mode._nbValues ) HOMOG2D_THROW_ERROR_1( "SVG path command: inconsistency with stored values, expected " << (size_t)mode._nbValues << ", got " << _values.size() ); auto pt = generateNewPoint( mode, _previousPt, _values ); HOMOG2D_LOG( "new point added: " << pt ); out.push_back( pt ); _previousPt = pt; _values.clear(); } void addValue( std::string elem ) { _values.push_back( std::stod(elem) ); } }; /// Parse a SVG "path" string and convert it to a vector holding a set (vector) of points /** Input string example: \verbatim m 261.68497,138.79393 2.57,3.15 -0.72,1.27 2.18,1.94 -0.7,4.93 1.88,0.9 \endverbatim The return value holds as 'second' a bool value, will be true if closed polyline */ inline auto parsePath( const char* s ) { SvgPathCommand mode; std::vector> vout(1); SvgValuesBuffer values; std::string str(s); size_t idx = 0; HOMOG2D_LOG( "parsing string -" << str << "- #=" << str.size() ); if( str.size() == 0 ) HOMOG2D_THROW_ERROR_1( "SVG path string is empty" ); auto it = str.cbegin(); Point2d previousPt; do { auto e = getNextElem( str, it ); // HOMOG2D_LOG( "parsing element -" << e << "- #=" << e.size() ); if( e.size() == 1 && !isDigit(e[0]) ) // we have a command ! { if( values.size() != 0 ) // if we have some values stored, values.storeValues( vout[idx], mode ); // first process them and add new point mode = getCommand( e[0] ); HOMOG2D_LOG( "command=" << e[0] ); if( !svgPathCommandIsAllowed(mode._command) ) HOMOG2D_THROW_ERROR_1( "SVG path command -" << mode._command << "- not handled" ); if( mode.isNewPolyline() && !vout[idx].empty() ) { vout.push_back( std::vector() ); idx++; HOMOG2D_LOG( "NEW vector idx=" << idx ); } } else // not a command, but a value { // HOMOG2D_LOG( "process value, values size=" << values.size() ); if( values.size() == (size_t)mode._nbValues ) // already got enough values values.storeValues( vout[idx], mode ); values.addValue( e ); } } while( it < str.cend() ); if( values.size() ) // process remaining values that have been stored values.storeValues( vout[idx], mode ); HOMOG2D_LOG( "Nb vectors=" << vout.size() ); for( auto& v: vout ) if( v.size() ) v = purgeSetDupes( v ); if( vout.back().empty() ) vout.pop_back(); return std::make_pair( vout, mode._command == 'Z' ? true : false ); } } // namespace svgp //------------------------------------------------------------------ /// Visitor class, derived from the tinyxml2 visitor class. Used to import SVG data. /** Holds the imported data through std::variant */ class Visitor: public tinyxml2::XMLVisitor { /// This type is used to provide a type that can be used in a switch (see VisitExit() ), /// as this cannot be done with a string |-( enum SvgType { T_circle, T_rect, T_line, T_polygon, T_polyline, T_ellipse, T_path, T_other ///< for other elements (\c ) or illegal ones, that will just be ignored }; /// A map holding correspondences between type as a string and type as a SvgType. /// Populated in constructor std::map _svgTypesTable; std::vector> _vecVar; ///< all the data is stored here public: /// Constructor, populates the table giving type from svg string Visitor() { // svg name local type id _svgTypesTable["circle"] = T_circle; _svgTypesTable["rect"] = T_rect; _svgTypesTable["line"] = T_line; _svgTypesTable["polyline"] = T_polyline; _svgTypesTable["polygon"] = T_polygon; _svgTypesTable["ellipse"] = T_ellipse; _svgTypesTable["path"] = T_path; } /// Returns the type as a member of enum SvgType, so the type can be used in a switch SvgType getSvgType( std::string s ) const { auto it = _svgTypesTable.find( s ); if( it == std::cend(_svgTypesTable) ) return T_other; return it->second; } const std::vector>& get() const { return _vecVar; } bool VisitExit( const tinyxml2::XMLElement& ) override; }; namespace svgp { //------------------------------------------------------------------ /// Fetch attribute from XML element. Tag \c e_name is there just in case of trouble. inline double getAttribValue( const tinyxml2::XMLElement& e, const char* str, std::string e_name ) { double value=0.; if( tinyxml2::XML_SUCCESS != e.QueryDoubleAttribute( str, &value ) ) HOMOG2D_THROW_ERROR_1( "h2d::svg::import error, failed to read attribute '" << std::string{str} << "' while reading element '" << e_name << "'" ); return value; } /// Helper function for SVG import /** \todo Who owns the data? Should we return a string and/or release the memory? */ inline const char* getAttribString( const char* attribName, const tinyxml2::XMLElement& e ) { const char *pts = e.Attribute( attribName ); if( !pts ) throw std::string("h2d::img::svg Error: unable to find attribute '") + attribName + "' in tag " + e.Name(); return pts; } } // namespace svgp /// This is the place where actual SVG data is converted and stored into vector /** (see manual, section "SVG import") Overload of the root class `VisitExit()` member function */ inline bool Visitor::VisitExit( const tinyxml2::XMLElement& e ) { std::string n = e.Name(); // std::cout << "element name:" << n << '\n'; try { switch( getSvgType( n ) ) { case T_circle: _vecVar.emplace_back( CircleD( svgp::getAttribValue( e, "cx", n ), svgp::getAttribValue( e, "cy", n ), svgp::getAttribValue( e, "r", n ) ) ); break; case T_rect: { auto x1 = svgp::getAttribValue( e, "x", n ); auto y1 = svgp::getAttribValue( e, "y", n ); auto w = svgp::getAttribValue( e, "width", n ); auto h = svgp::getAttribValue( e, "height", n ); _vecVar.emplace_back( FRectD( x1, y1, x1+w, y1+h ) ); } break; case T_line: _vecVar.emplace_back( SegmentD( svgp::getAttribValue( e, "x1", n ), svgp::getAttribValue( e, "y1", n ), svgp::getAttribValue( e, "x2", n ), svgp::getAttribValue( e, "y2", n ) ) ); break; case T_polygon: { auto pts_str = svgp::getAttribString( "points", e ); auto vec_pts = svgp::parsePoints( pts_str ); _vecVar.emplace_back( CPolylineD(vec_pts) ); } break; case T_polyline: { auto pts_str = svgp::getAttribString( "points", e ); auto vec_pts = svgp::parsePoints( pts_str ); _vecVar.emplace_back( OPolylineD(vec_pts) ); } break; case T_path: // a path can hold multiple polygons (because of the 'L' command) { auto pts_str = svgp::getAttribString( "d", e ); try { auto parse_res = svgp::parsePath( pts_str ); const auto& vec_vec_pts = parse_res.first; // for( const auto& vec_pts: vec_vec_pts ) if( parse_res.second == true ) _vecVar.emplace_back( CPolylineD(vec_pts) ); else _vecVar.emplace_back( OPolylineD(vec_pts) ); } catch( std::exception& err ) // an unhandled path command will just get the whole path command ignored { HOMOG2D_LOG_WARNING( "Unable to import SVG path command\n -msg=" << err.what() << "\n -input string=" << pts_str ); } } break; case T_ellipse: { auto x = svgp::getAttribValue( e, "cx", n ); auto y = svgp::getAttribValue( e, "cy", n ); auto rx = svgp::getAttribValue( e, "rx", n ); auto ry = svgp::getAttribValue( e, "ry", n ); auto rot = svgp::getEllipseRotateAttr( svgp::getAttribString( "transform", e ) ); auto ell = EllipseD( x, y, rx, ry ); auto H = Homogr().addTranslation(-x,-y).addRotation(rot.second).addTranslation(x,y); _vecVar.push_back( H * ell ); } break; default: // for T_other elements if( n != "svg" ) // because that one will be always there, so no need to show a warning HOMOG2D_LOG_WARNING( "found SVG element '" << n << "' in SVG file, left unprocessed" ); break; } } catch( std::string& msg ) { HOMOG2D_THROW_ERROR_1( "h2d: Tinyxml read error: " << msg ); return false; // to avoid a compile warning } return true; } inline void printFileAttrib( const tinyxml2::XMLDocument& doc ) { const tinyxml2::XMLElement* root = doc.RootElement(); const tinyxml2::XMLAttribute* pAttrib = root->FirstAttribute(); size_t i=0; while( pAttrib ) { std::cout << "Attrib " << i++ << ": Name=" << pAttrib->Name() << "; Value=" << pAttrib->Value() << '\n'; pAttrib=pAttrib->Next(); } } /// Fetch size of image in SVG file inline auto getImgSize( const tinyxml2::XMLDocument& doc ) { const tinyxml2::XMLElement* root = doc.RootElement(); const tinyxml2::XMLAttribute* pAttrib = root->FirstAttribute(); double w = -1., h = -1.; while( pAttrib ) { auto attr = std::string( pAttrib->Name() ); if( attr == "width" ) w = std::stod( pAttrib->Value() ); if( attr == "height" ) h = std::stod( pAttrib->Value() ); pAttrib=pAttrib->Next(); } if( w == -1. || h == -1. ) HOMOG2D_THROW_ERROR_1( "unable to find size in SVG file" ); return std::make_pair( w, h ); } } // namespace svg #endif // HOMOG2D_USE_SVG_IMPORT } // namespace h2d #endif // HG_HOMOG2D_HPP