/************************************************************************** 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-2025 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 // HOMOG2D_USE_SVG_IMPORT #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 // << std::scientific << std::setprecision(10) #ifdef HOMOG2D_DEBUGMODE #define HOMOG2D_LOG(a) \ std::cout << '-' << __FUNCTION__ << "(), line " << __LINE__ << ": " \ << 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 ) \ if( err::printWarnings() == true ) \ std::cerr << "homog2d warning (" << ++err::warningCount() << "), l. " << __LINE__ << ": " << 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.13.1" // 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 \todo 20250123: for some reason, the error count is always 42! Need to investigate this. */ inline size_t& errorCount() { static size_t c_err; return c_err; } /// Used in macro HOMOG2D_LOG_WARNING inline size_t& warningCount() { static size_t c_war; return c_war; } /// User can use this to silence warnings at runtime. inline bool& printWarnings() { static bool b_war = true; return b_war; } } //namespace err /// Holds the types needed for policy based design namespace typ { /// Used to determine if line or point, see base::LPBase struct IsLine {}; struct IsPoint {}; struct IsHomogr {}; struct IsEpipmat {}; /// Used to determine the type of "point pair (segment or vector), see base::SegVec struct IsSegment {}; /// \sa IsSegment struct IsOSeg {}; /// Used to determine the type of polyline (\ref CPolyline_ or \ref OPolyline_), see base::PolylineBase struct IsClosed {}; /// \sa IsClosed struct IsOpen {}; /// \name These are used as static member type SType ///@{ struct T_Point {}; struct T_Line {}; struct T_Circle {}; struct T_FRect {}; struct T_Segment {}; struct T_OSeg {}; struct T_OPol {}; struct T_CPol {}; struct T_Ellipse {}; ///@} } // 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 SegVec; } template class Hmatrix_; template using Homogr_ = Hmatrix_; #ifdef HOMOG2D_FUTURE_STUFF template using Epipmat_ = Hmatrix_; #endif 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 Segment_ = base::SegVec; template using OSegment_ = base::SegVec; template using PointPair2_ = std::pair,Point2d_>; template using PointPair_ = 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_, OSegment_, 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::FRect_; template friend class h2d::Ellipse_; template friend class h2d::base::PolylineBase; template friend class h2d::base::LPBase; template friend class h2d::base::SegVec; /// 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 bool _showAngles = false; 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 ) { if( ps%2 == 0 ) HOMOG2D_THROW_ERROR_1( "odd number required" ); _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; } /// Set or unset the drawing of angles of polylines DrawParams& showAngles( bool b=true ) { _dpValues._showAngles = 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 { friend std::ostream& operator << ( std::ostream&, const 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 getOrthogLine() 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() /// \sa type() enum class Type: uint8_t { Line2d, Point2d, Segment, OSegment, 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::OSegment: s="OSegment"; 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. /// See https://github.com/skramm/homog2d/blob/master/docs/homog2d_manual.md#section_rtp 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 to get the size (nb of points) of an object in a std::variant, call with std::visit() /** \sa CommonType_ \sa size() */ struct SizeFunct { template HOMOG2D_INUMTYPE operator ()(const T& a) { return a.size(); } }; //------------------------------------------------------------------ /// 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() /** The constructor has a third optional parameter that can be used to pass drawing parameters */ 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 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'; } } // 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 /// \todo 20250201: why do we need this and why isn't the same required for \c OSegment_ ? namespace base { template auto operator << ( std::ostream&, const h2d::base::LPBase& ) -> std::ostream&; template auto operator << ( std::ostream&, const h2d::base::PolylineBase& ) -> std::ostream&; } // namespace base // forward declaration, related to https://github.com/skramm/homog2d/issues/2 template Line2d_ operator * ( const Homogr_&, const Line2d_& ); // forward declaration template std::vector> getBisectorLines( const base::PolylineBase& ); template HOMOG2D_INUMTYPE getAngle( const T1&, const T2& ); namespace detail { // forward declaration template void product( Matrix_&, const Matrix_&, const Matrix_& ); //------------------------------------------------------------------ /// Returns true if one of the points share a common coordinate /// (thus making them unable to define a bounding box) 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; } //------------------------------------------------------------------ /// Returns true if one of the points share a common coordinate /// (thus making them unable to define a bounding box) template bool shareCommonCoord( const std::pair,Point2d_>& ppts ) { return shareCommonCoord( ppts.first , ppts.second ); } //------------------------------------------------------------------ /// 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 PointPair_ 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 ); } /// An alias used to hold data of a 3x3 matrix, see detail::Matrix_ 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; virtual size_t size() 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_normalizeMat( 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 /// (const because flag \c _hasChanged declared as \c mutable) void normalize() const { detail::Matrix_::p_normalizeMat(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 {}; 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 {}; //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() /// \todo 20250131: probably useless, check if this can be replaced by some "if constexpr" 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 {}; /// Traits class used to detect if container \c T is a \c std::array /** (because allocation is different, see \ref alloc() ) */ template struct IsArray : std::false_type {}; template struct IsArray> : 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; using SType = typ::T_Ellipse; 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 ///@{ /// Returns 1 constexpr size_t size() const { return 1; } bool isCircle( HOMOG2D_INUMTYPE thres=1.E-10 ) const; Point2d_ getCenter() const; CPolyline_ getOBB() const; auto 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_normalizeMat(2,2); if( !h.isNormalized() ) h.p_normalizeMat(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& ); template HOMOG2D_INUMTYPE crossProductV( const base::SegVec&, const base::SegVec& ); 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::base::SegVec; 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; } PointPair_ 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'; priv::printVector( 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 SType = typ::T_FRect; 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 PointPair_& 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 template FRect_( T1 x1, T2 y1, T3 x2, T4 y2 ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); HOMOG2D_CHECK_IS_NUMBER(T3); HOMOG2D_CHECK_IS_NUMBER(T4); 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(); } constexpr size_t size() const { return 4; } /// get BB of rectangle. Needed for getBB( pair of objects ) /// \todo 20250205: fix this so that it return a "full res" rectangle FRect_ getBB() const { //FRect_ out; // out = *this; return *this; // return out; } /// Returns the 2 major points of the rectangle /// \sa getPts( const FRect_& ) PointPair_ 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 SType = typ::T_Circle; 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 template Circle_( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& 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; } constexpr size_t size() const { return 1; } /// Area of circle HOMOG2D_INUMTYPE area() const { return static_cast(_radius) * _radius * M_PI; } /// Perimeter of circle HOMOG2D_INUMTYPE length() const { return static_cast(_radius) * M_PI * 2.0; } /// Returns Bounding Box of circle auto 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 ); } template void set( const Point2d_&, const Point2d_& ); template void set( const Point2d_&, const Point2d_&, const Point2d_& ); // 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 PointPair_& 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 ); } /// Get sign of value /** source: https://stackoverflow.com/a/4609795/ */ template int sign(T val) { return (T(0) < val) - (val < T(0)); } //------------------------------------------------------------------ /// 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; } } // 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; // The following line enables defining the SType based on how the class is instanciated: Line or Point using SType = LP; // using SType = std::conditional,typ::T_Point,typ::T_Line>; 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&, const h2d::base::LPBase& ) -> 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 ) { static_assert( std::is_same_v, "Invalid: you cannot build a point from a line" ); p_copyFrom( li ); } /// 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 ); } /// 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); p_init_2( v1, v2 ); } /// 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( T1 x1, T2 y1, T3 x2, T4 y2 ) { HOMOG2D_CHECK_IS_NUMBER(T1); HOMOG2D_CHECK_IS_NUMBER(T2); HOMOG2D_CHECK_IS_NUMBER(T3); HOMOG2D_CHECK_IS_NUMBER(T4); static_assert( std::is_same_v, "Invalid: you cannot build a point from 4 values" ); *this = Point2d_(x1, y1) * Point2d_(x2, y2); } #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 ) { p_init_2( val.HOMOG2D_BIND_X, val.HOMOG2D_BIND_Y ); } #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() { if constexpr( std::is_same_v ) { _v[0] = 1.; _v[1] = 0.; _v[2] = 0.; } else { _v[0] = 0.; _v[1] = 0.; _v[2] = 1.; } } /// Constructor of horizontal/vertical line template LPBase( LineDir orient, T value ) { impl_init_or( orient, value ); } #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 ) { init_BoostGeomPoint( pt ); } /// Constructor from boost::geometry second point type template // Boost Floating Point Type LPBase( const boost::geometry::model::d2::point_xy& pt ) { init_BoostGeomPoint( pt ); } /// 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]); } /// init function. if object is a point => copy-constructor, /// else we build the line passing though (0,0) ant the given point template void impl_init_1_Point( const Point2d_& pt ) { if constexpr( std::is_same_v ) p_copyFrom( pt ); else { *this = detail::crossProduct( pt, Point2d_() ); p_normalizePL(); } } template void impl_init_or( LineDir, T ); template void impl_init_or( LineDir, const Point2d_& ); public: Type type() const { if constexpr( std::is_same_v ) return Type::Point2d; else return Type::Line2d; } size_t size() const { if constexpr( std::is_same_v ) return 1; else return 0; } template HOMOG2D_INUMTYPE getCoord( GivenCoord gc, FPT2 other ) const; template Point2d_ getPoint( GivenCoord gc, FPT2 other ) const { static_assert( std::is_same_v, "Invalid: you cannot call on a point" ); HOMOG2D_CHECK_IS_NUMBER( FPT2 ); auto coord = getCoord( gc, other ); if( gc == GivenCoord::X ) return Point2d_( other, coord ); return Point2d_( coord, other ); } /// Returns a pair of points that are lying on line at distance \c dist from a point defined by one of its coordinates. template PointPair_ getPoints( GivenCoord gc, FPT2 coord, FPT3 dist ) const { static_assert( std::is_same_v, "Invalid: you cannot call on a point" ); HOMOG2D_CHECK_IS_NUMBER(FPT2); HOMOG2D_CHECK_IS_NUMBER(FPT3); const auto pt = getPoint( gc, coord ); return priv::getPoints_B2( pt, dist, *this ); } /// 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 PointPair_ getPoints( const Point2d_& pt, FPT3 dist ) const; template Line2d_ getOrthogLine( GivenCoord gc, FPT2 other ) const; template Line2d_ getOrthogLine( const Point2d_& pt ) const; /// Returns a line rotated at point \c pt with angle \c theta template Line2d_ getRotatedLine( const Point2d_& pt, T theta ) const; /// Returns the segment from the point (not on line) to the line, shortest path template OSegment_ getOrthogSegment( const Point2d_& pt ) const; /// Returns an parallel line to the one it is called on, with \c pt lying on it. /// \todo clarify orientation: on which side will that line appear? template Line2d_ getParallelLine( const Point2d_& ) const; /// Returns the pair of parallel lines at a distance \c dist from line. template std::pair,Line2d_> getParallelLines( T dist ) const; HOMOG2D_INUMTYPE getX() const { static_assert( std::is_same_v, "Invalid: cannot get x for a line" ); return static_cast(_v[0]) / _v[2]; } HOMOG2D_INUMTYPE getY() const { static_assert( std::is_same_v, "Invalid: cannot get y for a line" ); return static_cast(_v[1]) / _v[2]; } /// Translate Point2d, does nothing for Line2d template void translate( T1 dx, T2 dy ) { HOMOG2D_CHECK_IS_NUMBER( T1 ); HOMOG2D_CHECK_IS_NUMBER( T2 ); if constexpr( std::is_same_v ) { _v[0] = static_cast(_v[0]) / _v[2] + dx; _v[1] = static_cast(_v[1]) / _v[2] + dy; _v[2] = 1.; p_normalizePL(); } } /// Translate Point2d, does nothing for Line2d 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 => NO, FAILS TO BUILD ! template void moveTo( const Point2d_& pt ) { static_assert( std::is_same_v, "Invalid: cannot move a line" ); *this = pt; } #ifdef HOMOG2D_ENABLE_VRTP /// Needed so the function getBB(T1,T2) builds, whatever the types /// and 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] }; } /// Set point from 2 euclidean values template void set( T1 x, T2 y ) { HOMOG2D_CHECK_IS_NUMBER( T1 ); HOMOG2D_CHECK_IS_NUMBER( T2 ); static_assert( std::is_same_v, "Invalid call for lines" ); _v[0] = x; _v[1] = y; _v[2] = 1.; p_normalizePL(); } template HOMOG2D_INUMTYPE distTo( const Point2d_& pt ) const { return impl_distToPoint( pt, detail::BaseHelper() ); } template HOMOG2D_INUMTYPE distTo( const Line2d_& ) const; template HOMOG2D_INUMTYPE distTo( const Segment_& ) const; template bool isParallelTo( const Line2d_& ) const; template bool isParallelTo( const Segment_& seg ) const { return isParallelTo( seg.getLine() ); } template HOMOG2D_INUMTYPE getAngle( const LPBase& other ) const; /// Returns angle in rad. between line and segment \c seg. \sa h2d::getAngle() template HOMOG2D_INUMTYPE getAngle( const base::SegVec& seg ) const { return getAngle( seg.getLine() ); } /// Returns true if point is at infinity (third value less than thr::nullDenom() ) bool isInf() const { if constexpr( std::is_same_v ) return false; else return homog2d_abs( _v[2] ) < thr::nullDenom(); } /// A point has a null length, a line has an infinite length HOMOG2D_INUMTYPE length() const { if constexpr( std::is_same_v ) { HOMOG2D_THROW_ERROR_1( "unable, a line has an infinite length" ); } else return 0.; } /// Neither lines nor points have an area constexpr HOMOG2D_INUMTYPE area() const { return 0.; } private: template HOMOG2D_INUMTYPE impl_distToPoint( const Point2d_&, const detail::BaseHelper& ) const; template HOMOG2D_INUMTYPE impl_distToPoint( const Point2d_&, const detail::BaseHelper& ) 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 template bool isInside( const Point2d_& pt1, const Point2d_& pt2 ) const { HOMOG2D_START; return isInside( FRect_(pt1, pt2) ); } /// Point is inside FRect template bool isInside( const FRect_& rect ) const; /// Point is inside circle defined by center and radius template< typename T, typename std::enable_if< !std::is_same_v> ,T >::type* = nullptr > bool isInside( const Point2d_& center, T radius ) const { HOMOG2D_CHECK_IS_NUMBER(T); if constexpr( std::is_same_v ) return false; else { if( distTo( center ) < radius ) return true; return false; } } /// Point is inside Circle template bool isInside( const Circle_& cir ) const { return isInside( cir.center(), cir.radius() ); } /// Point is inside Ellipse template bool isInside( const Ellipse_& ell ) const; 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; /// 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() ); } 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 point or line from a single OpenCv point. template LPBase( cv::Point_ pt ) { if constexpr( std::is_same_v ) { p_init_2( pt.x, pt.y ); } else { Point2d_ p(pt); impl_init_1_Point( p ); } } #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_isInsidePoly( const base::PolylineBase&, const detail::BaseHelper& ) const; template constexpr bool impl_isInsidePoly( const base::PolylineBase&, const detail::BaseHelper& ) const; bool impl_op_equal( const LPBase&, const detail::BaseHelper& ) const; bool impl_op_equal( 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; template void p_init_2( const T1&, const T2& ); #ifdef HOMOG2D_USE_BOOSTGEOM template // Boost Floating Point Type void init_BoostGeomPoint( const boost::geometry::model::point& pt ) { static_assert( std::is_same_v, "Invalid: you cannot build a Line2d using a boost::geometry point" ); set( pt ); } #endif }; // class LPBase //------------------------------------------------------------------ /// Helper function for constructor of horizontal or vertical line /// (overload for argument of type Point2d) /// \sa LineDir template template void LPBase::impl_init_or( LineDir dir, const Point2d_& pt ) { static_assert( std::is_same_v, "Invalid: you cannot build a horiz/vertical point" ); 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.; } } //------------------------------------------------------------------ /// Helper function for constructor of horizontal or vertical line /// (overload for numerical argument) /// \sa LineDir template template void LPBase::impl_init_or( LineDir dir, T arg ) { static_assert( std::is_same_v, "Invalid: you cannot build a horiz/vertical point" ); HOMOG2D_CHECK_IS_NUMBER(T); _v[2] = -arg; if( dir == LineDir::V ) { _v[0] = 1.; _v[1] = 0.; } else // = LineDir::H { _v[0] = 0.; _v[1] = 1.; } } } // 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 ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ /// Returns true is lines (or segments) are parallel template bool areParallel( const T1& t1, const T2& t2 ) { return t1.isParallelTo(t2); } //------------------------------------------------------------------ /// 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 } } /// Point related to a OSegment ///\sa SegVec::getPointSide() enum class PointSide: uint8_t { Left,Right,Neither }; inline const char* getString( PointSide t ) { const char* s=0; switch( t ) { case PointSide::Left: s="Left"; break; case PointSide::Right: s="Right"; break; case PointSide::Neither: s="Neither"; break; assert(0); } return s; } namespace base { //------------------------------------------------------------------ /// A line segment, oriented (\ref OSegment_) or not (\ref Segment_). Holds the two points. /** This will get instanciated as \ref Segment_ or \ref OSegment_ The difference between theses two types is that with \c Segment_, the "smallest" point is always stored as first element (see constructor). */ template class SegVec: public detail::Common #ifdef HOMOG2D_ENABLE_PRTP , public rtp::Root #endif { public: using FType = FPT; // using SType = std::conditional,typ::T_Segment,typ::T_OSeg>; using SType = SV; //std::conditional,typ::T_Segment,typ::T_OSeg>; using detail::Common::isInside; template friend class SegVec; /// \todo 20250127: if this works, then generalize to all the other base type() member functions Type type() const { if constexpr( std::is_same_v ) return Type::Segment; else return Type::OSegment; } private: Point2d_ _ptS1, _ptS2; public: /// \name Constructors ///@{ /// Default constructor: initializes segment to (0,0)--(1,1) SegVec(): _ptS2(1.,1.) {} /// Constructor 2: build segment from two points, can hold different FP types /** Please note that the source (points) floating-point type is lost */ template SegVec( 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 if constexpr( std::is_same_v ) priv::fix_order( _ptS1, _ptS2 ); } /// Constructor 3: build segment from two points coordinates, calls constructor 2 template SegVec( FP1 x1, FP2 y1, FP3 x2, FP4 y2 ) : SegVec( Point2d_(x1,y1), Point2d_(x2,y2) ) { HOMOG2D_CHECK_IS_NUMBER(FP1); HOMOG2D_CHECK_IS_NUMBER(FP2); HOMOG2D_CHECK_IS_NUMBER(FP3); HOMOG2D_CHECK_IS_NUMBER(FP4); } /// Constructor 4: build segment from pair of points template SegVec( const PointPair2_& ppts ) : SegVec(ppts.first, ppts.second) {} /// Copy-Constructor, behavior depends on concrete types /** - OSegment(OSegment): OK - Segment(OSegment): OK, but loose orientation - Segment(Segment): OK - OSegment(Segment): No build, because orientation would be arbitrary */ template SegVec( const SegVec& other ) : _ptS1(other._ptS1), _ptS2(other._ptS2) { static_assert( ( std::is_same_v || ( std::is_same_v && std::is_same_v ) ), "Cannot build a OSegment from a Segment" ); if constexpr( std::is_same_v ) priv::fix_order( _ptS1, _ptS2 ); } ///@} /// \name Modifying functions ///@{ /// Reverse oriented segment SegVec operator - () { static_assert( !std::is_same_v, "cannot reverse non-oriented segment" ); std::swap( _ptS1, _ptS2 ); return *this; } /// 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; if constexpr( std::is_same_v ) priv::fix_order( _ptS1, _ptS2 ); } /// Setter from a std::pair template void set( const PointPair2_& ppts ) { set( ppts.first, ppts.second ); } /// Setter from 4 raw point coordinates template void set( FP1 x1, FP2 y1, FP3 x2, FP4 y2 ) { HOMOG2D_CHECK_IS_NUMBER(FP1); HOMOG2D_CHECK_IS_NUMBER(FP2); HOMOG2D_CHECK_IS_NUMBER(FP3); HOMOG2D_CHECK_IS_NUMBER(FP4); 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 ///@{ #ifdef HOMOG2D_ENABLE_VRTP /// Get Bounding Box of segment, always throws but needed so the function getBB(T1,T2) builds, whatever the types /// and because of variant (\sa CommonType) (HOMOG2D_ENABLE_VRTP) FRect_ getBB() const { HOMOG2D_THROW_ERROR_1( "invalid call, Segment has no Bounding Box" ); } #endif /// Returns 2 constexpr size_t size() const { return 2; } /// Get segment length HOMOG2D_INUMTYPE length() const { return _ptS1.distTo( _ptS2 ); } /// A segment always has a null area constexpr HOMOG2D_INUMTYPE area() const { return 0.; } /// Get vector of Oriented segment as a pair of values std::pair getVector() const { static_assert( std::is_same_v, "cannot return vector of unoriented segment" ); return std::make_pair( _ptS2.getX() - _ptS1.getX(), _ptS2.getY() - _ptS1.getY() ); } /// Get angle between segment and other segment /** - if either the object or the argument is not oriented, then this will return the line angles, in the range [0:+PI] - if both are oriented, will return a value in the range [-PI:+PI] */ template HOMOG2D_INUMTYPE getAngle( const SegVec& other ) const { // if one of the two is a (unoriented) segment, then we just return the lines angle if constexpr( std::is_same_v || std::is_same_v ) return other.getLine().getAngle( this->getLine() ); else { // both are oriented segments auto v1 = this->getVector(); auto v2 = other.getVector(); auto dx1 = v1.first; auto dx2 = v2.first; auto dy1 = v1.second; auto dy2 = v2.second; return std::atan2( dx1 * dy2 - dy1 * dx2, dx1 * dx2 + dy1 * dy2 ); } } /// Get angle between segment and line template HOMOG2D_INUMTYPE getAngle( const base::LPBase& other ) const { static_assert( std::is_same_v, "cannot compute angle between segment and point" ); return getLine().getAngle( other ); } ///@} /// \name Operators ///@{ bool operator == ( const SegVec& s2 ) const { if( _ptS1 != s2._ptS1 ) return false; if( _ptS2 != s2._ptS2 ) return false; return true; } bool operator != ( const SegVec& s2 ) const { return !(*this == s2); } bool operator < ( const SegVec& 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 of segment as a std::pair /** If "Segment", then 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 */ PointPair_ getPts() const { return std::make_pair( _ptS1, _ptS2 ); } template std::pair getParallelSegs( T dist ) const; template Point2d_ getPointAt( FPT2 dist ) const; 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/vector std::array,4> getOrthogPts() const { return p_getOrthog(); } /// Returns the 4 segments/vectors 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 auto getLine() const { Point2d_ pt1( _ptS1 ); Point2d_ pt2( _ptS2 ); auto li = pt1 * pt2; return Line2d_(li); } template PointSide getPointSide( const Point2d_& pt ) const; std::pair,SegVec> split() const; template std::pair,SegVec> split( T dist ) const; 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, a vector, or a line" ); return getLine().isParallelTo( other ); } /// Returns point at middle distance between \c p1 and \c p2 auto getCenter() const { return Point2d_( ( static_cast(_ptS1.getX()) + _ptS2.getX() ) / 2., ( static_cast(_ptS1.getY()) + _ptS2.getY() ) / 2. ); } SegVec getExtended() const; /// Returns the bisector line of the segment /// \sa free function h2d::getBisector() auto getBisector() const { SegVec seg2 = *this; // convert to (possibly) enhance precision return seg2.getLine().getOrthogLine( seg2.getCenter() ); } ///@} /// \name Intersection functions ///@{ template detail::Intersect intersects( const SegVec& ) 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 { return r.intersects( *this ); } /// Segment/Polyline intersection template detail::IntersectM intersects( const base::PolylineBase& other ) const { return other.intersects( *this ); } ///@} template friend std::ostream& operator << ( std::ostream&, const SegVec& ); #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 SegVec //------------------------------------------------------------------ /// Returns a pair of segments/vectors split by the middle template std::pair,SegVec> SegVec::split() const { auto pt_mid = getCenter(); return std::make_pair( SegVec( _ptS1, pt_mid ), SegVec( pt_mid, _ptS2 ) ); } //------------------------------------------------------------------ /// Returns a pair of segments/vectors split by point at arbitrary distance /// \note: distance value must be >0 and less than segment length template template std::pair,OSegment_> SegVec::split( T dist ) const { static_assert( std::is_same_v, "cannot use this on non-oriented segments" ); HOMOG2D_CHECK_IS_NUMBER(T); #ifndef HOMOG2D_NOCHECKS if( dist <= 0 ) HOMOG2D_THROW_ERROR_1( "distance value must be >=0, current value=" << dist ); if( dist >= length() ) HOMOG2D_THROW_ERROR_1( "distance value must less than length, current value=" << dist << " length=" << length() ); #endif auto pt = getPointAt(dist); return std::make_pair( SegVec( _ptS1, pt ), SegVec( pt, _ptS2 ) ); } //------------------------------------------------------------------ /// Return side of point \c pt relatively to the segment (Oriented only) template template PointSide SegVec::getPointSide( const Point2d_& pt ) const { static_assert( std::is_same_v, "unable to get side of point related to Segment" ); if( pt == _ptS1 || pt == _ptS2 ) return PointSide::Neither; OSegment_ other( _ptS1, pt ); HOMOG2D_INUMTYPE cp = detail::crossProductV( other, *this ); PointSide out; switch ( priv::sign(cp) ) { case 0: out = PointSide::Neither; break; case -1: out = PointSide::Left; break; case +1: out = PointSide::Right; break; default: assert(0); } return out; } //------------------------------------------------------------------ /// Returns the point at a distance \c dist from starting point of oriented vector template template Point2d_ SegVec::getPointAt( FPT2 dist ) const { static_assert( std::is_same_v, "Cannot apply on non-oriented segment" ); HOMOG2D_CHECK_IS_NUMBER(FPT2); #ifndef HOMOG2D_NOCHECKS if( dist < 0 ) HOMOG2D_THROW_ERROR_1( "distance value must be >0, current value=" << dist ); #endif Line2d_ li = getLine(); auto ppts = li.getPoints( _ptS1, dist ); if( _ptS2.distTo( ppts.first ) < _ptS2.distTo( ppts.second ) ) return ppts.first; else return ppts.second; } //------------------------------------------------------------------ /// Returns two parallel segments to the current one in a pair /** Behavior differs depending on type of segment: - if oriented (\ref OSegment_), the returned pair will hold as \c first the segment that is \b left of the original - if not oriented (\ref Segment_), the order will be unpredictable */ template template std::pair,SegVec> SegVec::getParallelSegs( T dist ) const { HOMOG2D_CHECK_IS_NUMBER( T ); #ifndef HOMOG2D_NOCHECKS if( dist <= 0 ) HOMOG2D_THROW_ERROR_1( "Invalid value for distance:" << dist ); #endif auto plines = getLine().getParallelLines( dist ); auto lo1 = getLine().getOrthogLine( _ptS1 ); auto lo2 = getLine().getOrthogLine( _ptS2 ); auto pA1 = lo1 * plines.first; auto pA2 = lo2 * plines.first; auto pB1 = lo1 * plines.second; auto pB2 = lo2 * plines.second; if constexpr( std::is_same_v ) if( getPointSide(pB1) == PointSide::Left ) { std::swap( pA1, pB1 ); std::swap( pA2, pB2 ); } return std::make_pair( SegVec( pA1, pA2 ), SegVec( pB1, pB2 ) ); } //------------------------------------------------------------------ /// 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 SegVec SegVec::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 \todo 20250219: pass an enum instead as an int, for \c segDistCase */ template template HOMOG2D_INUMTYPE SegVec::distTo( const Point2d_& pt, ///< the point int* segDistCase ///< optional, may be used to know if its the orthogonal distance or not ) 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 = static_cast(pt.getX()) - x1; auto B = static_cast(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 ); } } // namespace base //------------------------------------------------------------------ /// 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.isSimple() ) 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 _isSimplePolyg; priv::ValueFlag> _centroid; void setBad() { _length.setBad(); _area.setBad(); _isSimplePolyg.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 > PointPair_ getBB_Points( const T& vpts ) { HOMOG2D_START; using FPT = typename T::value_type::FType; HOMOG2D_DEBUG_ASSERT( static_cast(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 ) ); } /// Private helper function for base::PolylineBase::getSegs() and base::PolylineBase::getOSegs() template // Segment Type, PolyLine Type, Floating Point Type std::vector p_getSegs( const base::PolylineBase& pl, const ST& ) { auto siz = pl.size(); if( siz < 2 ) // nothing to return return std::vector(); std::vector out; out.reserve( siz ); for( size_t i=0; i ) out.push_back( ST(pl.getPoint(siz-1), pl.getPoint(0) ) ); return out; } } // 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; //------------------------------------------------------------------ /// Parameters for base::PolylineBase::getOffsetPoly() /** (Unused at present) */ struct OffsetPolyParams { bool _angleSplit = false; }; 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: typ::IsClosed or typ::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 SType = std::conditional,typ::T_CPol,typ::T_OPol>; 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 { if constexpr( std::is_same_v ) return Type::CPolyline; else 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 ) { static_assert( std::is_same_v, "cannot build an OPolyline object from a FRect"); for( const auto& pt: rect.get4Pts() ) _plinevec.push_back( pt ); } 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 ) { set( rad, n ); } /// 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_v ) // 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 ///@} public: /// \name Attributes access ///@{ /// Returns the number of points size_t size() const { return _plinevec.size(); } constexpr bool isClosed() const { if constexpr( std::is_same_v ) return true; else return false; } HOMOG2D_INUMTYPE length() const; HOMOG2D_INUMTYPE area() const; bool isSimple() const; /// Returns Bounding Box of Polyline auto 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; if constexpr( std::is_same_v ) return size(); else return size() - 1; } 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; /// 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 oriented segments of the polyline std::vector> getOSegs() const { return priv::p_getSegs( *this, OSegment_() ); } /// Returns (as a copy) the segments of the polyline std::vector> getSegs() const { return priv::p_getSegs( *this, Segment_() ); } /// 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]; } private: /// Private helper function template ST impl_getSegment( size_t idx, const ST& ) const { #ifndef HOMOG2D_NOCHECKS if( idx >= nbSegs() ) HOMOG2D_THROW_ERROR_1( "requesting segment " << idx << ", only has " << nbSegs() ); if( size() < 2 ) HOMOG2D_THROW_ERROR_1( "empty, no segment " << idx ); #endif if constexpr( std::is_same_v ) return ST( _plinevec[idx], _plinevec[idx+1==nbSegs()?0:idx+1] ); else return ST( _plinevec[idx], _plinevec[idx+1] ); } public: /// Returns one oriented segment of the polyline. OSegment_ getOSegment( size_t idx ) const { return impl_getSegment( idx, OSegment_() ); } /// 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 { return impl_getSegment( idx, Segment_() ); } CPolyline_ convexHull() const; ///@} 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 allow *it++ = pt; // type conversions (std::copy implies same FP type) } template void setParallelogram( const Point2d_& pt1, const Point2d_& pt2, const Point2d_& pt3 ); /// Set from FRect template void set( const FRect_& rect ) { static_assert( std::is_same_v, "Invalid: cannot set a OPolyline from a FRect" ); CPolyline_ tmp(rect); std::swap( *this, tmp ); } /// Build RCP (Regular Convex Polygon), and return distance between consecutive points template std::pair set( FPT2 rad, T n ); ///@} private: /// 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() ); } if( pts.front() == pts.back() ) HOMOG2D_THROW_ERROR_1( "cannot add first point equal to last point:\npt:" << pts.front() << " 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 constexpr( std::is_same_v ) { if( size() != other.size() ) // for quick exit return false; p_normalizePoly(); other.p_normalizePoly(); auto it = other._plinevec.begin(); for( const auto& elem: _plinevec ) if( *it++ != elem ) return false; return true; } else return false; } template bool operator != ( const PolylineBase& other ) const { return !( *this == other ); } ///@} /// 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& ); #ifdef HOMOG2D_TEST_MODE /// this is only needed for testing bool isNormalized() const { return _plIsNormalized; } #endif private: void p_normalizePoly() const; public: template PolylineBase getOffsetPoly( T value, OffsetPolyParams p=OffsetPolyParams{} ) const; std::vector> getBisectorLines() const { return h2d::getBisectorLines( *this ); } 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 //------------------------------------------------------------------ /// Return an "offsetted" closed polyline, requires simple polygon (CPolyline AND no crossings) as input /** On failure (for whatever reason), will return an empty CPolyline - If dist>0, returns the polyline "outside" the source one - If dist<0, returns the polyline "inside" the source one */ template template PolylineBase PolylineBase::getOffsetPoly( T dist, OffsetPolyParams params ) const { HOMOG2D_CHECK_IS_NUMBER(T); auto dist0 = (HOMOG2D_INUMTYPE)dist; // casting bool valid = true; if( homog2d_abs(dist) < thr::nullDistance() ) { HOMOG2D_LOG_WARNING( "Failure, distance value is null, returning empty CPolyline" ); valid = false; } if( size()<3 ) { HOMOG2D_LOG_WARNING( "Failure, computing offsetted Polyline requires at least 3 points, returning empty CPolyline" ); valid = false; } if( !isSimple() ) { HOMOG2D_LOG_WARNING( "Failure, Polyline is not a polygon, returning empty CPolyline" ); valid = false; } if( !valid ) return PolylineBase(); p_normalizePoly(); size_t current = 0; // bool paraLines = false; std::vector> v_out; auto osegs = getOSegs(); auto oseg1 = osegs.at(current); do { auto next = (current==size()-1 ? 0 : current+1); auto pt1 = getPoint(next); auto oseg2 = osegs.at(next); auto psegs1 = oseg1.getParallelSegs(homog2d_abs(dist0)); auto psegs2 = oseg2.getParallelSegs(homog2d_abs(dist0)); auto pseg1 = &psegs1.first; auto pseg2 = &psegs2.first; if( dist < 0 ) { pseg1 = &psegs1.second; pseg2 = &psegs2.second; } auto li1 = pseg1->getLine(); auto li2 = pseg2->getLine(); if( !areParallel(li1,li2) ) { auto pt_int = li1 * li2; if( !params._angleSplit || getAngle(oseg1,oseg2)>0 || dist < 0 ) { v_out.push_back( pt_int ); // add point intersection of the two lines } else { auto oseg = OSegment_( pt1, pt_int ); auto dist_cut = std::min( oseg.length(), dist0 ); auto pt_cut = oseg.getPointAt( dist_cut ); auto oli = oseg.getLine().getOrthogLine( pt_cut ); auto pt_cut1 = oli * li1; auto pt_cut2 = oli * li2; v_out.push_back( pt_cut1 ); if( pt_cut1 != pt_cut2 ) // so we don't add two identical pts (in case they were computed as such) v_out.push_back( pt_cut2 ); } } // else // std::cout << "parallel!\n"; current++; oseg1 = oseg2; } while( current( v_out ); } //------------------------------------------------------------------ /// Build a parallelogram (4 points) from 3 points template template void PolylineBase::setParallelogram( const Point2d_& p1, const Point2d_& p2, const Point2d_& p3 ) { static_assert( std::is_same_v, "Invalid: cannot set a OPolyline as a parallelogram" ); 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 ); } //------------------------------------------------------------------ /// 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::set( FPT2 rad, ///< Radius T ni ///< Nb of points ) { static_assert( std::is_integral::value, "2nd argument type must be integral type" ); static_assert( std::is_same_v, "Cannot build a RCP as open polyline" ); 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 > auto getBmPoint( const T& t ) { #ifndef HOMOG2D_NOCHECKS if( t.size() == 0 ) HOMOG2D_THROW_ERROR_1( "invalid call, container is empty" ); #endif auto it = priv::getBmPoint_helper( t ); return std::make_pair( *it, std::distance( std::begin(t), it ) ); } //------------------------------------------------------------------ /// Return Top-most point of container template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > auto 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 auto it = 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 std::make_pair( *it, std::distance( std::begin(t), it ) ); } //------------------------------------------------------------------ /// Return Left-most point of container as a pair holding: /** - first: the point - second: the distance (index) from begin position \todo 20250222: CHANGE RETURN TYPE! => should return an iterator, not a pair (so this allows usage of any non-random access container) */ template< typename T, typename std::enable_if< ( trait::IsContainer::value && std::is_same>::value ), T >::type* = nullptr > auto 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 auto it = 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 std::make_pair( *it, std::distance( std::begin(t), it ) ); } //------------------------------------------------------------------ /// 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 > auto 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 auto it = 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 std::make_pair( *it, std::distance( std::begin(t), it ) ); } //------------------------------------------------------------------ /// Return Bottom-most point of Polyline (free function) template Point2d_ getBmPoint( const base::PolylineBase& poly ) { return getBmPoint( poly.getPts() ).first; } /// Return Top-most point of Polyline (free function) template Point2d_ getTmPoint( const base::PolylineBase& poly ) { return getTmPoint( poly.getPts() ).first; } /// Return Left-most point of Polyline (free function) template Point2d_ getLmPoint( const base::PolylineBase& poly ) { return getLmPoint( poly.getPts() ).first; } /// Return Right-most point of Polyline (free function) template Point2d_ getRmPoint( const base::PolylineBase& poly ) { return getRmPoint( poly.getPts() ).first; } //------------------------------------------------------------------ /// 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 auto 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 ); } //------------------------------------------------------------------ #ifdef HOMOG2D_PRELIMINAR /// TMP struct TMP_DebugSquare { Line2d_ li0; Line2d_ line_ortho; Line2d_ li_L, li_R; Segment_ s1; Segment_ sL,sR; double dist; PointPair_ ppts; Point2d_ pt_resL, pt_resR; }; template auto buildSquare( const CONT& pts ) { TMP_DebugSquare dbg; if( pts.size() != 4 ) HOMOG2D_THROW_ERROR_1( "Error, size must be = 4, current size=" << pts.size() ); if( areCollinear( pts[0], pts[1], pts[2] ) ) HOMOG2D_THROW_ERROR_1( "Unable, points 012 are colinear" ); if( areCollinear( pts[1], pts[2], pts[3] ) ) HOMOG2D_THROW_ERROR_1( "Unable, points 123 are colinear" ); if( areCollinear( pts[1], pts[0], pts[3] ) ) HOMOG2D_THROW_ERROR_1( "Unable, points 013 are colinear" ); if( areCollinear( pts[0], pts[2], pts[3] ) ) HOMOG2D_THROW_ERROR_1( "Unable, points 023 are colinear" ); auto pt_L = getLmPoint( pts ); auto pt_R = getRmPoint( pts ); auto pt_T = getTmPoint( pts ); auto pt_B = getBmPoint( pts ); Segment_ s1( pt_L.first, pt_R.first ); dbg.s1 = s1; auto dist = s1.length(); dbg.dist = dist; auto line_ortho = s1.getLine().getOrthogLine( pt_T.first ); dbg.line_ortho = line_ortho; auto ppts = line_ortho.getPoints( pt_T.first, dist ); dbg.ppts = ppts; auto d1 = ppts.first.distTo( pt_L.first ); auto d2 = ppts.second.distTo( pt_L.first ); auto goodpt = ppts.first; if( d2 < d1 ) goodpt = ppts.second; auto li0 = goodpt * pt_B.first; dbg.li0 = li0; auto li_L = li0.getOrthogLine( pt_L.first ); auto li_R = li0.getOrthogLine( pt_R.first ); dbg.li_L = li_L; dbg.li_R = li_R; auto pt_resL = li_L * li0; auto pt_resR = li_R * li0; dbg.pt_resL = pt_resL; dbg.pt_resR = pt_resR; auto s_Top = li_L.getOrthogSegment( pt_T.first ); auto li_Top = s_Top.getLine(); auto pt_res3 = li_L * li_Top; auto pt_res4 = li_R * li_Top; std::vector> vout; vout.push_back( pt_resL ); vout.push_back( pt_res3 ); vout.push_back( pt_res4 ); vout.push_back( pt_resR ); return std::make_pair( PolylineBase( vout ), dbg ); } template PolylineBase buildSquare( const Point2d_& p1, const Point2d_& p2, const Point2d_& p3, const Point2d_& p4 ) { std::vector> v(4); v[0] = p1; v[1] = p2; v[2] = p3; v[3] = p4; return buildSquare( v ); } #endif // HOMOG2D_PRELIMINAR //------------------------------------------------------------------ /// 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::isSimple() const { if( size()<3 ) // needs at least 3 points to be a polygon return false; if constexpr( std::is_same_v ) // If open, then not a polygon return false; else // If closed, we need to check for crossings { if( _attribs._isSimplePolyg.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( !isSimple() ) 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( !isSimple() ) // 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 \warning: centroid CAN be outside of polygon! */ template Point2d_ base::PolylineBase::centroid() const { if( !isSimple() ) HOMOG2D_THROW_ERROR_1( "unable, Polyline object is not simple" ); 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(); } //------------------------------------------------------------------ /// Normalization of CPolyline_ /** Two tasks: - rotating so that the smallest one is first - for closed polylines: constant orientation \note: is \c const, although it moves the points, but all the attributes are the same (done by declaring the vector of points as \c mutable). Reference for the orientation: - http://www.faqs.org/faqs/graphics/algorithms-faq/ Subject 2.07: How do I find the orientation of a simple polygon? - https://en.wikipedia.org/wiki/Curve_orientation#Orientation_of_a_simple_polygon */ template void base::PolylineBase::p_normalizePoly() const { if( size() == 0 ) return; if( !_plIsNormalized ) { if constexpr ( std::is_same_v ) { auto minpos = std::min_element( _plinevec.begin(), _plinevec.end() ); std::rotate( _plinevec.begin(), minpos, _plinevec.end() ); auto ptB = _plinevec[0]; auto ptA = _plinevec[1]; auto ptC = _plinevec[size()-1]; HOMOG2D_INUMTYPE xA = ptA.getX(); HOMOG2D_INUMTYPE xB = ptB.getX(); HOMOG2D_INUMTYPE xC = ptC.getX(); HOMOG2D_INUMTYPE yA = ptA.getY(); HOMOG2D_INUMTYPE yB = ptB.getY(); HOMOG2D_INUMTYPE yC = ptC.getY(); auto det = (xB-xA)*(yC-yA) - (xC-xA)*(yB-yA); if( det < 0) { std::reverse( _plinevec.begin(), _plinevec.end() ); minpos = std::min_element( _plinevec.begin(), _plinevec.end() ); std::rotate( _plinevec.begin(), minpos, _plinevec.end() ); } } else { if( _plinevec.back() < _plinevec.front() ) std::reverse( _plinevec.begin(), _plinevec.end() ); } _plIsNormalized=true; } } } // 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 namespace base { //------------------------------------------------------------------ /// 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 SegVec::intersects( const SegVec& 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 SegVec::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 SegVec::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; } } // namespace base ///////////////////////////////////////////////////////////////////////////// // SECTION - STREAMING OPERATORS ///////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------ 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& pl ) { if constexpr( std::is_same_v ) f << '[' << pl._v[0] << ',' << pl._v[1] << ',' << pl._v[2] << "]"; else { if( pl.isInf() ) f << '[' << pl._v[0] << ',' << pl._v[1] << ',' << pl._v[2] << "]"; else f // << std::scientific << std::setprecision(25) << '[' << pl.getX() << ',' << pl.getY() << "]"; } return f; } } // namespace base /// 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; } namespace base { template std::ostream& operator << ( std::ostream& f, const h2d::base::SegVec& seg ) { f << seg._ptS1; if constexpr( std::is_same_v ) f << "=>"; else f << "-"; f << seg._ptS2; return f; } 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.getOrthogLine( 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 auto 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.getOrthogLine( ptA ); auto li_V2 = li_H.getOrthogLine( 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 template HOMOG2D_INUMTYPE LPBase::getCoord( GivenCoord gc, FPT2 other ) const { static_assert( std::is_same_v, "Invalid: you cannot call on a point" ); 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_1( "null denominator encountered" ); #endif if( gc == GivenCoord::X ) return ( -a * other - _v[2] ) / b; else return ( -b * other - _v[2] ) / a; } //------------------------------------------------------------------ /// Returns pair of points on line at distance \c dist from point on line /// at coord \c coord. template template PointPair_ LPBase::getPoints( const Point2d_& pt, FPT3 dist ) const { static_assert( std::is_same_v, "Invalid: you cannot call on a point" ); HOMOG2D_CHECK_IS_NUMBER(FPT2); HOMOG2D_CHECK_IS_NUMBER(FPT3); #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 ); } //------------------------------------------------------------------ /// Returns an orthogonal line to the one it is called on, at a point on the line defined by one of its coordinates. template template Line2d_ LPBase::getOrthogLine( GivenCoord gc, FPT2 val ) const { static_assert( std::is_same_v, "Invalid: you cannot call getOrthogLine() on a point" ); auto other_val = getCoord( gc, val ); Point2d_ pt( other_val, val ) ; if( gc == GivenCoord::X ) pt.set( val, other_val ); return priv::getOrthogonalLine_B2( pt, *this ); } /// Returns an orthogonal line to the one it is called on, at point \c pt template template Line2d_ LPBase::getOrthogLine( const Point2d_& pt ) const { static_assert( std::is_same_v, "Invalid: you cannot call getOrthogLine() on a point" ); /* #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_1( "point is not on line" ); } #endif */ return priv::getOrthogonalLine_B2( pt, *this ); } //------------------------------------------------------------------ template template Line2d_ LPBase::getRotatedLine( const Point2d_& pt, T theta ) const { HOMOG2D_CHECK_IS_NUMBER(T); static_assert( std::is_same_v, "Invalid: you cannot call getRotatedLine() on a point" ); #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 (oriented) segment that joins a point and a line /** Upon return, the first point will hold the intersection point (projection of point on line), and the second will hold the given point. */ template template OSegment_ LPBase::getOrthogSegment( const Point2d_& pt ) const { static_assert( !std::is_same_v, "Invalid call, cannot compute orthogonal segment between two points" ); 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, distance too small=" << dist ); #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->getOrthogLine( pt ); auto p2 = *this * oline; return OSegment_( p2, pt ); } //------------------------------------------------------------------ /// Returns an parallel line, implementation of getParallelLine(). template template Line2d_ LPBase::getParallelLine( const Point2d_& pt ) const { static_assert( std::is_same_v, "Invalid: you cannot call getParallelLine() on a point" ); Line2d_ out = *this; out._v[2] = static_cast(-_v[0]) * pt.getX() - _v[1] * pt.getY(); out.p_normalizePL(); return out; } //------------------------------------------------------------------ /// Implementation for lines template template std::pair,Line2d_> LPBase::getParallelLines( T dist ) const { static_assert( std::is_same_v, "Invalid: you cannot call getParallelLines() on a point" ); 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 (illegal for lines) template bool LPBase::operator < ( const LPBase& other ) const { static_assert( std::is_same_v, "Invalid < operator: you cannot sort lines" ); if( getX() < other.getX() ) return true; if( getX() > other.getX() ) return false; if( getY() < other.getY() ) return true; return false; } } // namespace base //------------------------------------------------------------------ /// Inner implementation details namespace detail { /// Cross product of points * points or line * line /** - 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; } /// Cross product of two oriented segments, return a scalar template HOMOG2D_INUMTYPE crossProductV( const base::SegVec& v1, const base::SegVec& v2 ) { auto ppts1 = v1.getPts(); auto ppts2 = v2.getPts(); auto dx1 = static_cast( ppts1.second.getX() ) - static_cast( ppts1.first.getX() ); auto dy1 = static_cast( ppts1.second.getY() ) - static_cast( ppts1.first.getY() ); auto dx2 = static_cast( ppts2.second.getX() ) - static_cast( ppts2.first.getX() ); auto dy2 = static_cast( ppts2.second.getY() ) - static_cast( ppts2.first.getY() ); return dx1 * dy2 - dy1 * dx2; } } // namespace detail //------------------------------------------------------------------ /////////////////////////////////////////// // CONSTRUCTORS /////////////////////////////////////////// namespace base { /// Generic init from two numeric args template template void LPBase::p_init_2( const T1& v1, const T2& v2 ) { if constexpr( std::is_same_v ) { _v[0] = v1; _v[1] = v2; _v[2] = 1.; } else { 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::distTo( const Line2d_& li ) const { static_assert( !std::is_same_v, "Invalid: you cannot compute distance between two lines" ); return li.distTo( *this ); } template template HOMOG2D_INUMTYPE LPBase::distTo( const Segment_& seg ) const { static_assert( !std::is_same_v, "Invalid: you cannot compute distance between a line and a segment" ); return seg.distTo( *this ); } } // namespace base //------------------------------------------------------------------ /// Free function, returns angle between two segments/lines /// \sa Segment_::getAngle() /// \sa Line2d_::getAngle() template HOMOG2D_INUMTYPE getAngle( const T1& t1, const T2& t2 ) { static_assert( ( std::is_same_v || std::is_same_v || std::is_same_v ) && ( std::is_same_v || std::is_same_v || std::is_same_v ), "Both types must be either a Line or a Segment" ); return t1.getAngle( t2 ); } //------------------------------------------------------------------ namespace base { template template bool LPBase::isParallelTo( const Line2d_& li ) const { static_assert( std::is_same_v, "Invalid: you cannot use IsParallel() with a point" ); if( this->getAngle(li) < thr::nullAngleValue() ) return true; return false; } //------------------------------------------------------------------ /// 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 template HOMOG2D_INUMTYPE LPBase::getAngle( const LPBase& li ) const { static_assert( !std::is_same_v, "cannot get angle of a point" ); 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 ); } //------------------------------------------------------------------ /// Returns true if point is inside (or on the edge) of a flat rectangle template template bool LPBase::isInside( const FRect_& rect ) const { if constexpr( std::is_same_v ) return false; else { auto pair_pts = rect.getPts(); const auto& p00 = pair_pts.first; const auto& p11 = pair_pts.second; return detail::ptIsInside( *this, p00, p11 ); } } //------------------------------------------------------------------ #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.isSimple() ) 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; } //------------------------------------------------------------------ /// Point is inside Polyline template template bool LPBase::isInside( const Ellipse_& ell ) const { if constexpr( std::is_same_v ) return ell.pointIsInside( *this ); else 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: typ::IsLine or typ::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/vector, returns a point template Point2d_ operator * ( const base::SegVec& lhs, const base::SegVec& 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 base::SegVec& seg ) { const auto& pts = seg.getPts(); Point2d_ pt1 = h * pts.first; Point2d_ pt2 = h * pts.second; return base::SegVec( 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::IsArray::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::IsArray::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 ); } //------------------------------------------------------------------ /// Intersection area over Union area (free function) 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::isSimple() template bool isSimple( const base::PolylineBase& pl ) { return pl.isSimple(); } //------------------------------------------------------------------ /// 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(); } //------------------------------------------------------------------ /// 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 Bounding Box of Ellipse_ (free function) /// \sa Ellipse_::getBB() template CPolyline_ getOBB( const Ellipse_& ell ) { return ell.getOBB(); } //------------------------------------------------------------------ /// Returns bisector lines of a Polyline /** Size of vector will be: - equal to size() of polyline if closed - equal to size()-2 of polyline if open \warning: - if open: the FIRST line (index 0) will correspond to SECOND point of the polyline */ template std::vector> getBisectorLines( const base::PolylineBase& pl ) { if( pl.size() < 3 ) { HOMOG2D_THROW_ERROR_1( "unable, minimum size is 3, currently=" << pl.size() ); } const auto& pts = pl.getPts(); std::vector> out; if constexpr ( std::is_same_v ) { out.reserve( pl.size()-2 ); OSegment_ seg1( pts[1], pts[0] ); for( size_t i=0; i seg2( pts[i+1], pts[i+2] ); auto angle = seg1.getAngle( seg2 ); out.emplace_back( seg1.getLine().getRotatedLine( pts[i+1], angle/2. ) ); seg1 = -seg2; } } else { out.reserve( pl.size() ); OSegment_ seg1( pts[0], pts[pl.size()-1] ); for( size_t i=0; i seg2( pts[i], pts[next] ); auto angle = seg1.getAngle( seg2 ); out.emplace_back( seg1.getLine().getRotatedLine( pts[i], angle/2. ) ); seg1 = -seg2; } } return out; } //------------------------------------------------------------------ /// Holds free functions 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< ( // type 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 vector && !std::is_same>::value // not a CPolyline && !std::is_same>::value // not a OPolyline ) ,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_v> || std::is_same_v> ) ,T >::type* = nullptr > auto getPointPair( const T& elem ) { HOMOG2D_START; return elem.getPts(); } /// Needed because of variant template PointPair_ 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 PointPair_ 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 ///////////////////////////////////////////////////////////////////////////// /** \addtogroup varff Free functions handling variant type defgroup These function can be used on a "regular" geometric object or on a `CommonType` object, holding a std::variant and allowing run-time polymorphism, See md_docs_homog2d_manual.html#section_rtp */ //------------------------------------------------------------------ /// Free function. Returns the type of object or variant /** Can be printed with `getString()` \sa CommonType_ \ingroup varff */ 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()` \ingroup varff */ 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) /// \ingroup varff 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) /// \ingroup varff 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(); } //------------------------------------------------------------------ /// Returns size of element or variant (free function) /// \ingroup varff template HOMOG2D_INUMTYPE size( const T& elem ) { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( trait::IsVariant::value ) return std::visit( fct::SizeFunct{}, elem ); else #endif return elem.size(); } #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 ); PointPair_ 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] \ingroup varff */ 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, \b NOT the one of the input element. \ingroup varff */ 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_v> && !std::is_same_v> ),T1 >::type* = nullptr > auto getBB( const T1& elem1, const T2& elem2 ) { HOMOG2D_START; bool isPoly_1 = false; if( elem1.type() == Type::OPolyline || elem1.type() == Type::CPolyline ) isPoly_1 = true; bool isPoly_2 = false; if( elem2.type() == Type::OPolyline || elem2.type() == Type::CPolyline ) isPoly_1 = true; auto siz_1 = elem1.size(); auto siz_2 = elem2.size(); if( isPoly_1 && isPoly_2 ) if( siz_1==0 && siz_2==0 ) HOMOG2D_THROW_ERROR_1( "unable to compute bounding box of two empty polyline objects" ); FRect_ out; // we need to do that in two steps: first, copy return value into `out`, then return // because the underlying floating point might not be the same for the two arguments. if( isPoly_2 && siz_1==0 ) { auto pp = ppair::getPointPair( elem2 ); if( shareCommonCoord(pp) ) { HOMOG2D_THROW_ERROR_1( "Unable to build common Bounding Box of empty polyline and elem2:" << elem2 ); } else return FRect_( pp ); } if( isPoly_1 && siz_2==0 ) { auto pp = ppair::getPointPair( elem1 ); if( shareCommonCoord(pp) ) { HOMOG2D_THROW_ERROR_1( "Unable to build common Bounding Box of empty polyline and elem1:" << elem1 ); } else return FRect_( pp ); } 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 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: PointPair_ 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,base::SegVec> split( const base::SegVec& seg ) { return seg.split(); } /// Returns a pair of segments, corresponding to the input segment split by point at distance \c dist (free function) /// \note: can only be used on oriented segments /// \sa Segment_::split( T ) template std::pair,OSegment_> split( const OSegment_& seg, T dist ) { HOMOG2D_CHECK_IS_NUMBER(T); return seg.split( dist ); } /// Returns the 4 orthogonal segments to the segment /// \sa Segment_::getOrthogSegs() template std::array,4> getOrthogSegs( const Segment_& seg ) { return seg.getOrthogSegs(); } /// Returns the 4 orthogonal points to the segment /// \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 Circle_::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(); } namespace priv { /// used in findPoint() struct F_MIN{}; struct F_MAX{}; /// Private. Common function for searching nearest of farthest point /* \sa findNearestPoint() \sa findFarthestPoint */ template size_t findPoint( const Point2d_& qpt, ///< query point const CONT& cont, ///< container const S_WHAT& ///< F_MIN or F_MAX // const CONT2* mask ///< optional vector or list or array of bool ) { if( cont.size() < 2 ) HOMOG2D_THROW_ERROR_1( "container holds " << cont.size() \ << " points, minimum is 2" ); decltype( priv::sqDist( qpt, qpt ) ) resDist; size_t startIdx = 1; size_t resIdx = 0; // if point is the first one of the container, then // initialize with the second one if( cont[0] == qpt ) { startIdx++; resIdx++; resDist = priv::sqDist( qpt, cont[1] ); } else // initialize with first point { resDist = priv::sqDist( qpt, cont[0] ); } for( size_t i=startIdx; i ) { if( currentDist < resDist ) { resIdx = i; resDist = currentDist; } } else { if( currentDist > resDist ) { resIdx = i; resDist = currentDist; } } } } return resIdx; } } //namespace priv //------------------------------------------------------------------ /// 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 ) { return priv::findPoint( pt, cont, priv::F_MIN() ); } //------------------------------------------------------------------ /// Returns index of point in container \c cont that is the farthest to \c pt template size_t findFarthestPoint( const Point2d_& pt, const T& cont ) { return priv::findPoint( pt, cont, priv::F_MAX() ); } //------------------------------------------------------------------ /// Returns indexes of points in container \c cont that are nearest/farthest /** - return pair: first is nearest, second is farthest */ template auto findNearestFarthestPoint( const Point2d_& pt, const T& cont ) { if( cont.size() < 2 ) HOMOG2D_THROW_ERROR_1( "container holds " << cont.size() \ << " points, minimum is 2" ); auto maxDist = priv::sqDist( pt, cont[0] ); auto minDist = maxDist; size_t idxMin = 0; size_t idxMax = 0; size_t startIdx = 1; if( pt == cont[0] ) { idxMin++; idxMax++; startIdx++; maxDist = priv::sqDist( pt, cont[1] ); minDist = maxDist; } for( size_t i=startIdx; i maxDist ) { idxMax = i; maxDist = currentDist; } if( currentDist < minDist ) { idxMin = i; minDist = currentDist; } } } return std::make_pair(idxMin, idxMax); } /// Returns set of points that are inside primitive \c prim /** Output container will be of same type as input container (vector, list or array \todo 20250123: maybe replace the bunch of static_asserts with something around trait::HasArea ? */ template CONT getPtsInside( const CONT& input_set, const PRIM& prim ///< geometrical primitive (FRect, Circle, ...?) ) { static_assert( !std::is_same_v, "Cannot find points inside a Open Polyline" ); static_assert( !std::is_same_v, "Cannot find points inside a Line" ); static_assert( !std::is_same_v, "Cannot find points inside a Point" ); static_assert( !std::is_same_v, "Cannot find points inside a Segment" ); static_assert( !trait::IsArray::value, "Cannot use std::array as container" ); CONT out; out.reserve( input_set.size() ); for( auto pt: input_set ) if( pt.isInside( prim ) ) out.push_back( pt ); return out; } ///////////////////////////////////////////////////////////////////////////// // 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. - The types inside the container can be either plain \c h2d types (\c FRect, \c Segment, ...) of variant types, using the \c CommonType class (requires HOMOG2D_ENABLE_VRTP symbol). \ingroup varff */ 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() ) { #ifdef HOMOG2D_ENABLE_VRTP if constexpr( trait::IsVariant::value ) { fct::DrawFunct vde( img, dp ); // functor for( auto& e: cont ) std::visit( vde, e ); } else #endif { 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. /** At present, cannot be used with container holding variants... */ 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: case img::PtStyle::Squ: 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 ) { if( ps == img::PtStyle::Squ ) { std::vector> vseg{ {vpt[0], vpt[1]}, {vpt[2], vpt[1]}, {vpt[2], vpt[3]}, {vpt[0], vpt[3]} }; for( const auto& s:vseg ) s.draw( img, dp2 ); } else { 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::Squ: detail::drawPt( im, img::PtStyle::Squ, vpt, dp, true ); break; case img::PtStyle::Times: // "times" symbol detail::drawPt( im, img::PtStyle::Times, vpt, dp ); break; default: assert(0); } } //------------------------------------------------------------------ /// Helper function used by the draw(OSegment) function, returns the 3 segment corresponding /// to the "arrows" as 3 pairs of points /** Used both in the SVG and the Opencv backends \verbatim /|\ | | | --- \endverbatim */ namespace priv { template std::array,3> getArrowSegments( const base::SegVec& vec ) { std::array,Point2d_>,3> out; auto ppts = vec.getPts(); auto pt1 = ppts.first; auto pt2 = ppts.second; const int arSize = 8; Line2d_ liA = vec.getLine().getOrthogLine( pt1 ); out[0] = liA.getPoints( pt1, arSize ); auto ppts_B = vec.getLine().getPoints( pt2, arSize ); Point2d_ p0 = ppts_B.first; if( dist( pt1, ppts_B.first) > dist(pt1, ppts_B.second) ) p0 = ppts_B.second; Line2d_ liB = vec.getLine().getOrthogLine( p0 ); auto ppts_C = liB.getPoints( p0, arSize ); out[1] = std::make_pair(ppts_C.first, pt2); out[2] = std::make_pair(ppts_C.second, pt2); return out; } } // namespace priv ///////////////////////////////////////////////////////////////////////////// // SECTION .2 - CLASS DRAWING MEMBER FUNCTIONS (OpenCv) ///////////////////////////////////////////////////////////////////////////// #ifdef HOMOG2D_USE_OPENCV //------------------------------------------------------------------ /// Draw \c FRect (Opencv implementation) template 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 / \c OSegment (Opencv implementation) namespace base { template void SegVec::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 constexpr( std::is_same_v ) { auto arrsegs = priv::getArrowSegments( *this ); for( auto ppts: arrsegs ) cv::line( im.getReal(), ppts.first.getCvPtd(), ppts.second.getCvPtd(), dp.cvColor(), dp._dpValues._lineThickness, dp._dpValues._lineType==1?cv::LINE_AA:cv::LINE_8 ); } else { if( dp._dpValues._showPoints ) { _ptS1.draw( im, dp ); _ptS2.draw( im, dp ); } } } } // namespace base //------------------------------------------------------------------ /// 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& im ) { f << "\n" << "\n"; f << im._realImg._svgString.str(); f << "\n"; return f; } } // namespace img /// Draw text on Svg image /// \todo 20230118: find a way to add a default parameter for dp (not allowed on explicit instanciation) template <> inline void img::Image::drawText( std::string str, Point2d_ pt, img::DrawParams dp ) { getReal()._svgString << "" << str << "\n"; } #ifdef HOMOG2D_USE_OPENCV /// Draw text on Opencv image template <> inline void img::Image::drawText( std::string str, Point2d_ pt, img::DrawParams dp ) { 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"; } //------------------------------------------------------------------ namespace priv { /// Helper function to draw SVG segment template void drawSvgSeg( img::Image& im, const Point2d_& pt1, const Point2d_& pt2, std::string color, int thickness, std::string attribs=std::string() ) { im.getReal()._svgString << "\n"; } template void drawSvgSeg( img::Image& im, const std::pair,Point2d_>& ppts, std::string color, int thickness, std::string attribs=std::string() ) { drawSvgSeg( im, ppts.first, ppts.second, color, thickness, attribs ); } } // namespace priv //------------------------------------------------------------------ /// Draw \c Segment / \c OSegment (SVG implementation) /// \todo 20250127: implement arrows for the Opencv/png version, and share the code between the two versions namespace base { template void SegVec::draw( img::Image& im, img::DrawParams dp ) const { auto group = false; if( dp._dpValues._showPoints || std::is_same_v ) group = true; if( group ) im.getReal()._svgString << ""; auto pts = getPts(); priv::drawSvgSeg( im, pts, dp.getSvgRgbColor(), dp._dpValues._lineThickness, dp.getAttrString() ); if constexpr( std::is_same_v ) { auto arrsegs = priv::getArrowSegments( *this ); for( auto ppts: arrsegs ) priv::drawSvgSeg( im, ppts, dp.getSvgRgbColor(), dp._dpValues._lineThickness ); } else { if( dp._dpValues._showPoints ) { _ptS1.draw( im, dp ); _ptS2.draw( im, dp ); } } if( group ) im.getReal()._svgString << "\n"; } } // namespace base //------------------------------------------------------------------ #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 */ inline 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::OSegment: { const OSegment_* 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"; /// \todo 20240215 Why don't we use the drawText() function ? 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._showAngles ) { auto osegs = getOSegs(); // std::cout << "osegs size=" << osegs.size() << "\n"; const auto& pts = getPts(); for( size_t i=0; i\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 = 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_; using OSegment = OSegment_; /// Default circle type using Circle = Circle_; /// Default rectangle type using FRect = FRect_; /// Default polyline type using CPolyline = CPolyline_; using OPolyline = OPolyline_; /// Default ellipse type using Ellipse = Ellipse_; // float types using Line2dF = Line2d_; using Point2dF = Point2d_; using HomogrF = Homogr_; using SegmentF = Segment_; using OSegmentF = OSegment_; using CircleF = Circle_; using FRectF = FRect_; using EllipseF = Ellipse_; // double types using Line2dD = Line2d_; using Point2dD = Point2d_; using HomogrD = Homogr_; using SegmentD = Segment_; using OSegmentD = OSegment_; 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 OSegmentL = OSegment_; using CircleL = Circle_; using FRectL = FRect_; using EllipseL = Ellipse_; using CPolylineF = CPolyline_; using CPolylineD = CPolyline_; using CPolylineL = CPolyline_; using OPolylineF = OPolyline_; using OPolylineD = OPolyline_; using OPolylineL = OPolyline_; using PointPairF = PointPair_; using PointPairD = PointPair_; using PointPairL = PointPair_; using PointPair = PointPair_; #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; } /// Used to access the data once the file has been read 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; } /// Helper function called by Visitor::VisitExit() to process Polyline/Polygons inline std::vector importSvgPoints( const tinyxml2::XMLElement& e ) { auto pts_str = svgp::getAttribString( "points", e ); auto vec_pts = svgp::parsePoints( pts_str ); if( vec_pts.front() == vec_pts.back() ) // if first point equal to last vec_pts.pop_back(); // point, remove last point return vec_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 vpts = svgp::importSvgPoints( e ); _vecVar.emplace_back( CPolylineD(vpts) ); } break; case T_polyline: { auto vpts = svgp::importSvgPoints( e ); _vecVar.emplace_back( OPolylineD(vpts) ); } 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( auto vec_pts: vec_vec_pts ) // we need a copy so we may remove last point if equal to first // for( const auto& vec_pts: vec_vec_pts ) // we need a copy so we may remove last point if equal to first { if( vec_pts.front() == vec_pts.back() ) // if first point equal to last vec_pts.pop_back(); // point, remove last point 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