diff --git a/bindings/python/src/image/core/raster_image.cpp b/bindings/python/src/image/core/raster_image.cpp index d614cf665..3fe3d7e51 100644 --- a/bindings/python/src/image/core/raster_image.cpp +++ b/bindings/python/src/image/core/raster_image.cpp @@ -29,8 +29,8 @@ #define PYTHON_RASTER_IMAGE( dimension ) \ const auto name##dimension = \ "RasterImage" + std::to_string( dimension ) + "D"; \ - pybind11::class_< RasterImage##dimension##D, CellArray##dimension##D >( \ - module, name##dimension.c_str() ) \ + pybind11::class_< RasterImage##dimension##D, CellArray##dimension##D, \ + Identifier >( module, name##dimension.c_str() ) \ .def( pybind11::init< std::array< index_t, dimension > >() ) \ .def( \ "native_extension", &RasterImage##dimension##D::native_extension ) \ diff --git a/bindings/python/src/mesh/core/light_regular_grid.cpp b/bindings/python/src/mesh/core/light_regular_grid.cpp index 34bee5daf..9b07354cd 100644 --- a/bindings/python/src/mesh/core/light_regular_grid.cpp +++ b/bindings/python/src/mesh/core/light_regular_grid.cpp @@ -33,8 +33,8 @@ #define PYTHON_LIGHT_REGULAR_GRID( dimension ) \ const auto name##dimension = \ "LightRegularGrid" + std::to_string( dimension ) + "D"; \ - pybind11::class_< LightRegularGrid##dimension##D, Grid##dimension##D >( \ - module, name##dimension.c_str() ) \ + pybind11::class_< LightRegularGrid##dimension##D, Grid##dimension##D, \ + Identifier >( module, name##dimension.c_str() ) \ .def( pybind11::init< Point< dimension >, \ std::array< index_t, dimension >, \ std::array< double, dimension > >() ) \ diff --git a/cmake/ConfigureEarcut.cmake b/cmake/ConfigureEarcut.cmake index 965776c1e..f5d43b5fc 100644 --- a/cmake/ConfigureEarcut.cmake +++ b/cmake/ConfigureEarcut.cmake @@ -26,7 +26,7 @@ ExternalProject_Add(earcut BINARY_DIR ${EARCUT_PATH}/build STAMP_DIR ${EARCUT_PATH}/stamp GIT_REPOSITORY https://github.com/mapbox/earcut.hpp - GIT_TAG a30c14b5676adabe4714ff4173dae8a5d568ab59 + GIT_TAG 7fa7aa30183849e988ae79ab2eef19f9ae97acf4 GIT_SHALLOW ON GIT_PROGRESS ON CMAKE_GENERATOR ${CMAKE_GENERATOR} diff --git a/include/geode/basic/console_progress_logger_client.hpp b/include/geode/basic/console_progress_logger_client.hpp index ea9711c8d..7379d3f4c 100644 --- a/include/geode/basic/console_progress_logger_client.hpp +++ b/include/geode/basic/console_progress_logger_client.hpp @@ -40,16 +40,33 @@ namespace geode private: void start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t nb_steps ) override; void update( const uuid& progress_logger_id, + Logger::LEVEL level, index_t current_step, index_t nb_steps ) override; - void completed( const uuid& progress_logger_id ) override; + void completed( + const uuid& progress_logger_id, Logger::LEVEL level ) override; - void failed( const uuid& progress_logger_id ) override; + void failed( + const uuid& progress_logger_id, Logger::LEVEL level ) override; + + [[deprecated]] void start( const uuid& progress_logger_id, + const std::string& message, + index_t nb_steps ) override; + + [[deprecated]] void update( const uuid& progress_logger_id, + index_t current_step, + index_t nb_steps ) override; + + [[deprecated]] void completed( + const uuid& progress_logger_id ) override; + + [[deprecated]] void failed( const uuid& progress_logger_id ) override; private: IMPLEMENTATION_MEMBER( impl_ ); diff --git a/include/geode/basic/logger.hpp b/include/geode/basic/logger.hpp index d5f3408bb..aacf0e49c 100644 --- a/include/geode/basic/logger.hpp +++ b/include/geode/basic/logger.hpp @@ -54,6 +54,12 @@ namespace geode static void set_level( LEVEL level ); + template < typename... Args > + static void log( LEVEL level, const Args &...args ) + { + log( level, absl::StrCat( args... ) ); + } + template < typename... Args > static void trace( const Args &...args ) { @@ -96,6 +102,7 @@ namespace geode [[nodiscard]] static Logger &instance(); + static void log( LEVEL level, const std::string &message ); static void log_trace( const std::string &message ); static void log_debug( const std::string &message ); static void log_info( const std::string &message ); diff --git a/include/geode/basic/progress_logger.hpp b/include/geode/basic/progress_logger.hpp index 080970437..aa2b4d0ff 100644 --- a/include/geode/basic/progress_logger.hpp +++ b/include/geode/basic/progress_logger.hpp @@ -26,6 +26,7 @@ #include #include +#include #include namespace absl @@ -38,7 +39,10 @@ namespace geode class opengeode_basic_api ProgressLogger { public: - ProgressLogger( const std::string& message, index_t nb_steps ); + [[deprecated]] ProgressLogger( + const std::string& message, index_t nb_steps ); + ProgressLogger( + Logger::LEVEL level, const std::string& message, index_t nb_steps ); ~ProgressLogger(); index_t increment(); diff --git a/include/geode/basic/progress_logger_client.hpp b/include/geode/basic/progress_logger_client.hpp index a3af419c2..4d44d9853 100644 --- a/include/geode/basic/progress_logger_client.hpp +++ b/include/geode/basic/progress_logger_client.hpp @@ -26,6 +26,7 @@ #include #include +#include namespace geode { @@ -40,16 +41,34 @@ namespace geode virtual ~ProgressLoggerClient() = default; virtual void start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t nb_steps ) = 0; virtual void update( const uuid& progress_logger_id, + Logger::LEVEL level, index_t current_step, index_t nb_steps ) = 0; - virtual void completed( const uuid& progress_logger_id ) = 0; + virtual void completed( + const uuid& progress_logger_id, Logger::LEVEL level ) = 0; - virtual void failed( const uuid& progress_logger_id ) = 0; + virtual void failed( + const uuid& progress_logger_id, Logger::LEVEL level ) = 0; + + [[deprecated]] virtual void start( const uuid& progress_logger_id, + const std::string& message, + index_t nb_steps ) = 0; + + [[deprecated]] virtual void update( const uuid& progress_logger_id, + index_t current_step, + index_t nb_steps ) = 0; + + [[deprecated]] virtual void completed( + const uuid& progress_logger_id ) = 0; + + [[deprecated]] virtual void failed( + const uuid& progress_logger_id ) = 0; protected: ProgressLoggerClient() = default; diff --git a/include/geode/basic/progress_logger_manager.hpp b/include/geode/basic/progress_logger_manager.hpp index 88fb935aa..ebe03121d 100644 --- a/include/geode/basic/progress_logger_manager.hpp +++ b/include/geode/basic/progress_logger_manager.hpp @@ -43,16 +43,32 @@ namespace geode std::unique_ptr< ProgressLoggerClient >&& client ); static void start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t nb_steps ); static void update( const uuid& progress_logger_id, + Logger::LEVEL level, index_t current_step, index_t nb_steps ); - static void completed( const uuid& progress_logger_id ); + static void completed( + const uuid& progress_logger_id, Logger::LEVEL level ); - static void failed( const uuid& progress_logger_id ); + static void failed( + const uuid& progress_logger_id, Logger::LEVEL level ); + + [[deprecated]] static void start( const uuid& progress_logger_id, + const std::string& message, + index_t nb_steps ); + + [[deprecated]] static void update( const uuid& progress_logger_id, + index_t current_step, + index_t nb_steps ); + + [[deprecated]] static void completed( const uuid& progress_logger_id ); + + [[deprecated]] static void failed( const uuid& progress_logger_id ); private: ProgressLoggerManager(); diff --git a/include/geode/basic/types.hpp b/include/geode/basic/types.hpp index e1a3fd032..46a5f06dc 100644 --- a/include/geode/basic/types.hpp +++ b/include/geode/basic/types.hpp @@ -28,6 +28,7 @@ #endif #include +#include #include namespace geode diff --git a/include/geode/geometry/basic_objects/ellipse.hpp b/include/geode/geometry/basic_objects/ellipse.hpp new file mode 100644 index 000000000..ce6379c1e --- /dev/null +++ b/include/geode/geometry/basic_objects/ellipse.hpp @@ -0,0 +1,110 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#pragma once + +#include +#include + +namespace geode +{ + FORWARD_DECLARATION_DIMENSION_CLASS( OwnerEllipse ); + FORWARD_DECLARATION_DIMENSION_CLASS( Point ); + FORWARD_DECLARATION_DIMENSION_CLASS( Frame ); + + template < index_t dimension > + using RefPoint = std::reference_wrapper< const Point< dimension > >; + ALIAS_2D_AND_3D( RefPoint ); + template < index_t dimension > + using RefFrame = std::reference_wrapper< const Frame< dimension > >; + ALIAS_2D_AND_3D( RefFrame ); +} // namespace geode + +namespace geode +{ + template < typename PointType, typename FrameType, index_t dimension > + class GenericEllipse + { + public: + GenericEllipse( PointType center, FrameType axis ); + + GenericEllipse( + const GenericEllipse< PointType, FrameType, dimension >& other ); + GenericEllipse< PointType, FrameType, dimension >& operator=( + const GenericEllipse< PointType, FrameType, dimension >& other ); + GenericEllipse( GenericEllipse< PointType, FrameType, dimension >&& + other ) noexcept; + GenericEllipse< PointType, FrameType, dimension >& operator=( + GenericEllipse< PointType, FrameType, dimension >&& + other ) noexcept; + + [[nodiscard]] const Point< dimension >& center() const; + [[nodiscard]] const Frame< dimension >& axes() const; + + private: + PointType center_; + FrameType axes_; + }; + + template < index_t dimension > + class OwnerEllipse : public GenericEllipse< Point< dimension >, + Frame< dimension >, + dimension > + { + using Base = + GenericEllipse< Point< dimension >, Frame< dimension >, dimension >; + + public: + explicit OwnerEllipse( + Point< dimension > center, Frame< dimension > axes ); + + OwnerEllipse( const OwnerEllipse< dimension >& other ); + OwnerEllipse< dimension >& operator=( + const OwnerEllipse< dimension >& other ); + OwnerEllipse( OwnerEllipse< dimension >&& other ) noexcept; + OwnerEllipse< dimension >& operator=( + OwnerEllipse< dimension >&& other ) noexcept; + }; + ALIAS_2D_AND_3D( OwnerEllipse ); + + template < index_t dimension > + class Ellipse : public GenericEllipse< RefPoint< dimension >, + RefFrame< dimension >, + dimension > + { + using Base = GenericEllipse< RefPoint< dimension >, + RefFrame< dimension >, + dimension >; + + public: + Ellipse( + const Point< dimension >& center, const Frame< dimension >& axes ); + Ellipse( const Ellipse< dimension >& other ); + Ellipse( const OwnerEllipse< dimension >& other ); + Ellipse< dimension >& operator=( const Ellipse< dimension >& other ); + Ellipse( Ellipse< dimension >&& other ) noexcept; + Ellipse< dimension >& operator=( + Ellipse< dimension >&& other ) noexcept; + }; + ALIAS_2D_AND_3D( Ellipse ); +} // namespace geode diff --git a/include/geode/geometry/detail/aabb_impl.hpp b/include/geode/geometry/detail/aabb_impl.hpp index 88810c1ff..4fc4592a2 100644 --- a/include/geode/geometry/detail/aabb_impl.hpp +++ b/include/geode/geometry/detail/aabb_impl.hpp @@ -95,9 +95,15 @@ namespace geode bboxes, ROOT_INDEX, 0, bboxes.size() ); } const auto grain = async::detail::auto_grain_size( bboxes.size() ); - const auto nb_async_depth = std::log2( grain ); - async_depth_ = - depth_ > nb_async_depth ? depth_ - nb_async_depth : depth_; + if( grain < 8 ) + { + async_depth_ = 0; + } + else + { + const auto nb_async_depth = std::log2( grain ); + async_depth_ = depth_ - nb_async_depth; + } } [[nodiscard]] index_t nb_bboxes() const diff --git a/include/geode/geometry/distance.hpp b/include/geode/geometry/distance.hpp index 18f3b4e11..742e80908 100644 --- a/include/geode/geometry/distance.hpp +++ b/include/geode/geometry/distance.hpp @@ -32,10 +32,12 @@ namespace geode FORWARD_DECLARATION_DIMENSION_CLASS( Segment ); FORWARD_DECLARATION_DIMENSION_CLASS( Triangle ); FORWARD_DECLARATION_DIMENSION_CLASS( Sphere ); + FORWARD_DECLARATION_DIMENSION_CLASS( Ellipse ); ALIAS_2D_AND_3D( Point ); ALIAS_2D_AND_3D( Triangle ); ALIAS_2D_AND_3D( InfiniteLine ); ALIAS_3D( Segment ); + ALIAS_2D_AND_3D( Ellipse ); class Plane; class Tetrahedron; class Circle; @@ -136,6 +138,34 @@ namespace geode opengeode_geometry_api segment_triangle_distance( const Segment3D& segment, const Triangle3D& triangle ); + /*! + * Compute the smallest distance between two triangles + * @return a tuple containing: + * - the smallest distance. + * - the closest point on the first triangle. + * - the closest point on the second triangle. + */ + [[nodiscard]] std::tuple< double, Point3D, Point3D > + opengeode_geometry_api triangle_triangle_distance( + const Triangle3D& triangle0, const Triangle3D& triangle1 ); + + /*! + * Compute the smallest distance between two triangles + * @details if the two triangles are the same, return nullopt. Only non + * conformal part of triangles are considered in computation of distance, + * i.e. if the triangle have a common point, it iterates on opposite + * segment, if the triangle have a common edge, it computes distance with + * the opposite point + * @return a tuple containing: + * - the smallest distance. + * - the closest point on the first triangle. + * - the closest point on the second triangle. + */ + [[nodiscard]] std::optional< std::tuple< double, Point3D, Point3D > > + opengeode_geometry_api + triangle_triangle_distance_between_non_conformal_parts( + const Triangle3D& triangle0, const Triangle3D& triangle1 ); + /*! * Compute the distance between a point and a tetrahedron * @return a tuple containing: @@ -229,4 +259,15 @@ namespace geode */ [[nodiscard]] std::tuple< double, Point3D > opengeode_geometry_api point_disk_distance( const Point3D& point, const Disk& disk ); + + /*! + * Compute the smallest distance between a point and an ellipse + * @return a tuple containing: + * - the smallest distance. + * - the closest point on the ellipse. + */ + template < index_t dimension > + [[nodiscard]] std::tuple< double, Point< dimension > > + point_ellipse_distance( const Point< dimension >& point, + const Ellipse< dimension >& ellipse ); } // namespace geode diff --git a/include/geode/geometry/intersection.hpp b/include/geode/geometry/intersection.hpp index 3b3aaef23..29231cb53 100644 --- a/include/geode/geometry/intersection.hpp +++ b/include/geode/geometry/intersection.hpp @@ -37,11 +37,15 @@ namespace geode FORWARD_DECLARATION_DIMENSION_CLASS( Segment ); FORWARD_DECLARATION_DIMENSION_CLASS( Sphere ); FORWARD_DECLARATION_DIMENSION_CLASS( Triangle ); + FORWARD_DECLARATION_DIMENSION_CLASS( Ellipse ); + FORWARD_DECLARATION_DIMENSION_CLASS( OwnerEllipse ); ALIAS_2D_AND_3D( InfiniteLine ); ALIAS_2D_AND_3D( OwnerInfiniteLine ); ALIAS_2D_AND_3D( Point ); ALIAS_2D_AND_3D( Segment ); ALIAS_3D( Triangle ); + ALIAS_2D_AND_3D( Ellipse ); + ALIAS_1D_AND_2D_AND_3D( OwnerEllipse ); class Circle; class Cylinder; class Plane; @@ -228,4 +232,27 @@ namespace geode [[nodiscard]] IntersectionResult< OwnerInfiniteLine3D > opengeode_geometry_api plane_plane_intersection( const Plane& plane0, const Plane& plane1 ); + + /*! + * Compute the intersection between a segment and an ellipse + * @return an optional of the intersection points. + */ + template < index_t dimension > + [[nodiscard]] IntersectionResult< + absl::InlinedVector< Point< dimension >, 2 > > + opengeode_geometry_api segment_ellipse_intersection( + const Segment< dimension >& segment, + const Ellipse< dimension >& ellipse ); + + /*! + * Compute the intersection between a line and an ellipse + * @return an optional of the intersection points. + */ + template < index_t dimension > + [[nodiscard]] IntersectionResult< + absl::InlinedVector< Point< dimension >, 2 > > + opengeode_geometry_api line_ellipse_intersection( + const InfiniteLine< dimension >& line, + const Ellipse< dimension >& ellipse ); + } // namespace geode \ No newline at end of file diff --git a/src/geode/basic/console_progress_logger_client.cpp b/src/geode/basic/console_progress_logger_client.cpp index 366338b4e..899d6aee0 100644 --- a/src/geode/basic/console_progress_logger_client.cpp +++ b/src/geode/basic/console_progress_logger_client.cpp @@ -48,34 +48,38 @@ namespace geode public: void start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t /*nb_steps */ ) { info_.emplace( progress_logger_id, message ); - Logger::info( message, " started" ); + Logger::log( level, message, " started" ); } - void update( - const uuid& progress_logger_id, index_t current, index_t nb_steps ) + void update( const uuid& progress_logger_id, + Logger::LEVEL level, + index_t current, + index_t nb_steps ) { const auto percent = std::floor( static_cast< double >( current ) / nb_steps * 100 ); - Logger::info( info_.at( progress_logger_id ).message, " ", current, - "/", nb_steps, " (", percent, "%)" ); + Logger::log( level, info_.at( progress_logger_id ).message, " ", + current, "/", nb_steps, " (", percent, "%)" ); } - void completed( const uuid& progress_logger_id ) + void completed( const uuid& progress_logger_id, Logger::LEVEL level ) { const auto& info = info_.at( progress_logger_id ); - Logger::info( - info.message, " completed in ", info.timer.duration() ); + Logger::log( + level, info.message, " completed in ", info.timer.duration() ); info_.erase( progress_logger_id ); } - void failed( const uuid& progress_logger_id ) + void failed( const uuid& progress_logger_id, Logger::LEVEL level ) { const auto& info = info_.at( progress_logger_id ); - Logger::info( info.message, " failed in ", info.timer.duration() ); + Logger::log( + level, info.message, " failed in ", info.timer.duration() ); info_.erase( progress_logger_id ); } @@ -88,26 +92,56 @@ namespace geode ConsoleProgressLoggerClient::~ConsoleProgressLoggerClient() = default; void ConsoleProgressLoggerClient::start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t nb_steps ) { - impl_->start( progress_logger_id, message, nb_steps ); + impl_->start( progress_logger_id, level, message, nb_steps ); + } + + void ConsoleProgressLoggerClient::update( const uuid& progress_logger_id, + Logger::LEVEL level, + index_t current, + index_t nb_steps ) + { + impl_->update( progress_logger_id, level, current, nb_steps ); + } + + void ConsoleProgressLoggerClient::completed( + const uuid& progress_logger_id, Logger::LEVEL level ) + { + impl_->completed( progress_logger_id, level ); + } + + void ConsoleProgressLoggerClient::failed( + const uuid& progress_logger_id, Logger::LEVEL level ) + { + impl_->failed( progress_logger_id, level ); + } + + void ConsoleProgressLoggerClient::start( const uuid& progress_logger_id, + const std::string& message, + index_t nb_steps ) + { + impl_->start( + progress_logger_id, Logger::LEVEL::info, message, nb_steps ); } void ConsoleProgressLoggerClient::update( const uuid& progress_logger_id, index_t current, index_t nb_steps ) { - impl_->update( progress_logger_id, current, nb_steps ); + impl_->update( + progress_logger_id, Logger::LEVEL::info, current, nb_steps ); } void ConsoleProgressLoggerClient::completed( const uuid& progress_logger_id ) { - impl_->completed( progress_logger_id ); + impl_->completed( progress_logger_id, Logger::LEVEL::info ); } void ConsoleProgressLoggerClient::failed( const uuid& progress_logger_id ) { - impl_->failed( progress_logger_id ); + impl_->failed( progress_logger_id, Logger::LEVEL::info ); } } // namespace geode diff --git a/src/geode/basic/logger.cpp b/src/geode/basic/logger.cpp index f194e762a..350b205a4 100644 --- a/src/geode/basic/logger.cpp +++ b/src/geode/basic/logger.cpp @@ -23,6 +23,8 @@ #include #include +#include + #include #include @@ -41,6 +43,11 @@ namespace geode level_ = level; } + void log( LEVEL level, const std::string &message ) + { + leveled_log.at( level )( message ); + } + void log_trace( const std::string &message ) { if( level_ <= LEVEL::trace ) @@ -90,6 +97,17 @@ namespace geode } private: + const absl::flat_hash_map< geode::Logger::LEVEL, + std::function< void( const std::string & ) > > + leveled_log{ + { geode::Logger::LEVEL::trace, geode::Logger::log_trace }, + { geode::Logger::LEVEL::debug, geode::Logger::log_debug }, + { geode::Logger::LEVEL::info, geode::Logger::log_info }, + { geode::Logger::LEVEL::warn, geode::Logger::log_warn }, + { geode::Logger::LEVEL::err, geode::Logger::log_error }, + { geode::Logger::LEVEL::critical, geode::Logger::log_critical } + }; + LEVEL level_{ LEVEL::info }; }; @@ -113,6 +131,11 @@ namespace geode instance().impl_->set_level( level ); } + void Logger::log( LEVEL level, const std::string &message ) + { + instance().impl_->log( level, message ); + } + void Logger::log_trace( const std::string &message ) { instance().impl_->log_trace( message ); diff --git a/src/geode/basic/progress_logger.cpp b/src/geode/basic/progress_logger.cpp index 1f12c57a4..3c941fbd2 100644 --- a/src/geode/basic/progress_logger.cpp +++ b/src/geode/basic/progress_logger.cpp @@ -37,21 +37,24 @@ namespace geode class ProgressLogger::Impl { public: - Impl( const std::string& message, index_t nb_steps ) - : nb_steps_( nb_steps ), current_time_{ absl::Now() } + Impl( + Logger::LEVEL level, const std::string& message, index_t nb_steps ) + : nb_steps_( nb_steps ), + current_time_{ absl::Now() }, + level_( level ) { - ProgressLoggerManager::start( id_, message, nb_steps_ ); + ProgressLoggerManager::start( id_, level, message, nb_steps_ ); } ~Impl() { if( current_ == nb_steps_ ) { - ProgressLoggerManager::completed( id_ ); + ProgressLoggerManager::completed( id_, level_ ); } else { - ProgressLoggerManager::failed( id_ ); + ProgressLoggerManager::failed( id_, level_ ); } } @@ -63,7 +66,8 @@ namespace geode if( now - current_time_ > refresh_interval_ ) { current_time_ = now; - ProgressLoggerManager::update( id_, current_, nb_steps_ ); + ProgressLoggerManager::update( + id_, level_, current_, nb_steps_ ); } return current_; } @@ -87,11 +91,18 @@ namespace geode absl::Time current_time_; std::mutex lock_; absl::Duration refresh_interval_{ absl::Seconds( 1 ) }; + Logger::LEVEL level_; }; ProgressLogger::ProgressLogger( const std::string& message, index_t nb_steps ) - : impl_( message, nb_steps ) + : ProgressLogger( Logger::LEVEL::info, message, nb_steps ) + { + } + + ProgressLogger::ProgressLogger( + Logger::LEVEL level, const std::string& message, index_t nb_steps ) + : impl_( level, message, nb_steps ) { } diff --git a/src/geode/basic/progress_logger_manager.cpp b/src/geode/basic/progress_logger_manager.cpp index ce6097469..274fa90f6 100644 --- a/src/geode/basic/progress_logger_manager.cpp +++ b/src/geode/basic/progress_logger_manager.cpp @@ -39,42 +39,45 @@ namespace geode } void start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t nb_steps ) { const std::lock_guard< std::mutex > locking{ lock_ }; for( auto& logger : loggers_ ) { - logger->start( progress_logger_id, message, nb_steps ); + logger->start( progress_logger_id, level, message, nb_steps ); } } void update( const uuid& progress_logger_id, + Logger::LEVEL level, index_t current_step, index_t nb_steps ) { const std::lock_guard< std::mutex > locking{ lock_ }; for( auto& logger : loggers_ ) { - logger->update( progress_logger_id, current_step, nb_steps ); + logger->update( + progress_logger_id, level, current_step, nb_steps ); } } - void completed( const uuid& progress_logger_id ) + void completed( const uuid& progress_logger_id, Logger::LEVEL level ) { const std::lock_guard< std::mutex > locking{ lock_ }; for( auto& logger : loggers_ ) { - logger->completed( progress_logger_id ); + logger->completed( progress_logger_id, level ); } } - void failed( const uuid& progress_logger_id ) + void failed( const uuid& progress_logger_id, Logger::LEVEL level ) { const std::lock_guard< std::mutex > locking{ lock_ }; for( auto& logger : loggers_ ) { - logger->failed( progress_logger_id ); + logger->failed( progress_logger_id, level ); } } @@ -94,26 +97,32 @@ namespace geode } void ProgressLoggerManager::start( const uuid& progress_logger_id, + Logger::LEVEL level, const std::string& message, index_t nb_steps ) { - instance().impl_->start( progress_logger_id, message, nb_steps ); + instance().impl_->start( progress_logger_id, level, message, nb_steps ); } - void ProgressLoggerManager::update( - const uuid& progress_logger_id, index_t current_step, index_t nb_steps ) + void ProgressLoggerManager::update( const uuid& progress_logger_id, + Logger::LEVEL level, + index_t current_step, + index_t nb_steps ) { - instance().impl_->update( progress_logger_id, current_step, nb_steps ); + instance().impl_->update( + progress_logger_id, level, current_step, nb_steps ); } - void ProgressLoggerManager::completed( const uuid& progress_logger_id ) + void ProgressLoggerManager::completed( + const uuid& progress_logger_id, Logger::LEVEL level ) { - instance().impl_->completed( progress_logger_id ); + instance().impl_->completed( progress_logger_id, level ); } - void ProgressLoggerManager::failed( const uuid& progress_logger_id ) + void ProgressLoggerManager::failed( + const uuid& progress_logger_id, Logger::LEVEL level ) { - instance().impl_->failed( progress_logger_id ); + instance().impl_->failed( progress_logger_id, level ); } ProgressLoggerManager& ProgressLoggerManager::instance() diff --git a/src/geode/geometry/CMakeLists.txt b/src/geode/geometry/CMakeLists.txt index 4a2939cf3..954bdb19e 100644 --- a/src/geode/geometry/CMakeLists.txt +++ b/src/geode/geometry/CMakeLists.txt @@ -26,6 +26,7 @@ add_geode_library( "barycentric_coordinates.cpp" "basic_objects/circle.cpp" "basic_objects/cylinder.cpp" + "basic_objects/ellipse.cpp" "basic_objects/infinite_line.cpp" "basic_objects/plane.cpp" "basic_objects/polygon.cpp" @@ -61,6 +62,7 @@ add_geode_library( "barycentric_coordinates.hpp" "basic_objects/circle.hpp" "basic_objects/cylinder.hpp" + "basic_objects/ellipse.hpp" "basic_objects/infinite_line.hpp" "basic_objects/plane.hpp" "basic_objects/polygon.hpp" diff --git a/src/geode/geometry/basic_objects/ellipse.cpp b/src/geode/geometry/basic_objects/ellipse.cpp new file mode 100644 index 000000000..49f63af86 --- /dev/null +++ b/src/geode/geometry/basic_objects/ellipse.cpp @@ -0,0 +1,119 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#include + +namespace geode +{ + template < typename PointType, typename FrameType, index_t dimension > + GenericEllipse< PointType, FrameType, dimension >::GenericEllipse( + PointType center, FrameType axes ) + : center_( std::move( center ) ), axes_( std::move( axes ) ) + { + } + template < typename PointType, typename FrameType, index_t dimension > + GenericEllipse< PointType, FrameType, dimension >::GenericEllipse( + const GenericEllipse< PointType, FrameType, dimension >& ) = default; + template < typename PointType, typename FrameType, index_t dimension > + GenericEllipse< PointType, FrameType, dimension >& + GenericEllipse< PointType, FrameType, dimension >::operator=( + const GenericEllipse< PointType, FrameType, dimension >& ) = + default; + template < typename PointType, typename FrameType, index_t dimension > + GenericEllipse< PointType, FrameType, dimension >::GenericEllipse( + GenericEllipse< PointType, FrameType, dimension >&& ) noexcept = + default; + + template < typename PointType, typename FrameType, index_t dimension > + GenericEllipse< PointType, FrameType, dimension >& + GenericEllipse< PointType, FrameType, dimension >::operator=( + GenericEllipse< PointType, FrameType, dimension >&& ) noexcept = + default; + + template < typename PointType, typename FrameType, index_t dimension > + const Point< dimension >& + GenericEllipse< PointType, FrameType, dimension >::center() const + { + return center_; + } + template < typename PointType, typename FrameType, index_t dimension > + const Frame< dimension >& + GenericEllipse< PointType, FrameType, dimension >::axes() const + { + return axes_; + } + + template < index_t dimension > + OwnerEllipse< dimension >::OwnerEllipse( + Point< dimension > center, Frame< dimension > axes ) + : Base( std::move( center ), std::move( axes ) ) + { + } + template < index_t dimension > + OwnerEllipse< dimension >::OwnerEllipse( + const OwnerEllipse< dimension >& ) = default; + template < index_t dimension > + OwnerEllipse< dimension >& OwnerEllipse< dimension >::operator=( + const OwnerEllipse< dimension >& ) = default; + template < index_t dimension > + OwnerEllipse< dimension >::OwnerEllipse( + OwnerEllipse< dimension >&& ) noexcept = default; + template < index_t dimension > + OwnerEllipse< dimension >& OwnerEllipse< dimension >::operator=( + OwnerEllipse< dimension >&& ) noexcept = default; + + template < index_t dimension > + Ellipse< dimension >::Ellipse( + const Point< dimension >& center, const Frame< dimension >& axes ) + : Base( center, axes ) + { + } + template < index_t dimension > + Ellipse< dimension >::Ellipse( const Ellipse< dimension >& ) = default; + template < index_t dimension > + Ellipse< dimension >::Ellipse( const OwnerEllipse< dimension >& other ) + : Base( other.center(), other.axes() ) + { + } + template < index_t dimension > + Ellipse< dimension >& Ellipse< dimension >::operator=( + const Ellipse< dimension >& ) = default; + template < index_t dimension > + Ellipse< dimension >::Ellipse( Ellipse< dimension >&& ) noexcept = default; + template < index_t dimension > + Ellipse< dimension >& Ellipse< dimension >::operator=( + Ellipse< dimension >&& ) noexcept = default; + + template class opengeode_geometry_api + GenericEllipse< Point< 2 >, Frame< 2 >, 2 >; + template class opengeode_geometry_api + GenericEllipse< Point< 3 >, Frame< 3 >, 3 >; + template class opengeode_geometry_api + GenericEllipse< RefPoint< 2 >, RefFrame< 2 >, 2 >; + template class opengeode_geometry_api + GenericEllipse< RefPoint< 3 >, RefFrame< 3 >, 3 >; + template class opengeode_geometry_api OwnerEllipse< 2 >; + template class opengeode_geometry_api OwnerEllipse< 3 >; + template class opengeode_geometry_api Ellipse< 2 >; + template class opengeode_geometry_api Ellipse< 3 >; +} // namespace geode diff --git a/src/geode/geometry/basic_objects/sphere.cpp b/src/geode/geometry/basic_objects/sphere.cpp index 4752826c7..464a3a108 100644 --- a/src/geode/geometry/basic_objects/sphere.cpp +++ b/src/geode/geometry/basic_objects/sphere.cpp @@ -63,15 +63,9 @@ namespace geode BoundingBox< dimension > GenericSphere< PointType, dimension >::bounding_box() const { - Point< dimension > translation; - for( const auto i : LRange{ dimension } ) - { - translation.set_value( i, radius_ ); - } BoundingBox< dimension > bbox; - const Point< dimension >& origin = origin_; - bbox.add_point( origin + translation ); - bbox.add_point( origin - translation ); + bbox.add_point( origin_ ); + bbox.extends( radius_ ); return bbox; } diff --git a/src/geode/geometry/distance.cpp b/src/geode/geometry/distance.cpp index 7ed73baaa..6860d91cf 100644 --- a/src/geode/geometry/distance.cpp +++ b/src/geode/geometry/distance.cpp @@ -23,6 +23,7 @@ #include +#include #include #include @@ -339,6 +340,92 @@ namespace geode::point_point_distance( point, closest_point ); return std::make_tuple( distance, std::move( closest_point ) ); } + + std::pair< std::vector< geode::local_index_t >, + std::vector< geode::local_index_t > > + find_non_colocated_triangles_points( const geode::Triangle3D& triangle0, + const geode::Triangle3D& triangle1 ) + { + std::pair< std::array< bool, 3 >, std::array< bool, 3 > > colocated{ + { false, false, false }, { false, false, false } + }; + const auto& vertices0 = triangle0.vertices(); + const auto& vertices1 = triangle1.vertices(); + for( const auto vertex0 : geode::LRange{ 3 } ) + { + for( const auto vertex1 : geode::LRange{ 3 } ) + { + if( geode::point_point_distance( + vertices0[vertex0].get(), vertices1[vertex1].get() ) + <= geode::GLOBAL_EPSILON ) + { + colocated.first[vertex0] = true; + colocated.second[vertex1] = true; + } + } + } + + std::pair< std::vector< geode::local_index_t >, + std::vector< geode::local_index_t > > + result; + for( const auto v : geode::LRange{ 3 } ) + { + if( !colocated.first[v] ) + { + result.first.push_back( v ); + } + if( !colocated.second[v] ) + { + result.second.push_back( v ); + } + } + return result; + } + + std::tuple< double, geode::Point3D, geode::Point3D > test_close_triangles( + const std::vector< geode::local_index_t >& non_colocated_points, + const geode::Triangle3D& base_triangle, + const geode::Triangle3D& other_triangle ) + { + const auto& base_vertices = base_triangle.vertices(); + double min_distance{ std::numeric_limits< double >::max() }; + geode::Point3D point0; + geode::Point3D point1; + for( const auto vertex0 : non_colocated_points ) + { + for( const auto vertex1 : non_colocated_points ) + { + if( vertex0 == vertex1 ) + { + continue; + } + if( geode::point_point_distance( base_vertices[vertex0].get(), + base_vertices[vertex1].get() ) + <= geode::GLOBAL_EPSILON ) + { + DEBUG( vertex0 ); + DEBUG( vertex1 ); + SDEBUG( base_vertices[0].get() ); + SDEBUG( base_vertices[1].get() ); + SDEBUG( base_vertices[2].get() ); + SDEBUG( other_triangle.vertices()[0].get() ); + SDEBUG( other_triangle.vertices()[1].get() ); + SDEBUG( other_triangle.vertices()[2].get() ); + } + const geode::Segment3D edge{ base_vertices[vertex0], + base_vertices[vertex1] }; + auto [cur_distance, cur_pt0, cur_pt1] = + geode::segment_triangle_distance( edge, other_triangle ); + if( cur_distance < min_distance ) + { + min_distance = cur_distance; + point0 = cur_pt0; + point1 = cur_pt1; + } + } + } + return std::make_tuple( min_distance, point0, point1 ); + } } // namespace namespace geode @@ -923,6 +1010,112 @@ namespace geode closest_on_segment, closest_on_triangle ); } + std::tuple< double, Point3D, Point3D > + opengeode_geometry_api triangle_triangle_distance( + const Triangle3D& triangle0, const Triangle3D& triangle1 ) + { + const auto non_colocated_points = + find_non_colocated_triangles_points( triangle0, triangle1 ); + if( non_colocated_points.first.size() < 3 ) + { + for( const auto v : LRange{ 3 } ) + { + if( !absl::c_contains( non_colocated_points.first, v ) ) + { + auto& pt0 = triangle0.vertices()[v].get(); + auto [distance, pt1] = + point_triangle_distance( pt0, triangle1 ); + return std::tuple{ distance, pt0, pt1 }; + } + } + } + if( non_colocated_points.second.size() < 3 ) + { + for( const auto v : LRange{ 3 } ) + { + if( !absl::c_contains( non_colocated_points.second, v ) ) + { + auto& pt1 = triangle1.vertices()[v].get(); + auto [distance, pt0] = + point_triangle_distance( pt1, triangle0 ); + return std::tuple{ distance, pt0, pt1 }; + } + } + } + auto [cur_distance0, cur_pt00, cur_pt01] = test_close_triangles( + non_colocated_points.first, triangle0, triangle1 ); + auto [cur_distance1, cur_pt11, cur_pt10] = test_close_triangles( + non_colocated_points.second, triangle1, triangle0 ); + if( cur_distance0 < cur_distance1 ) + { + return std::tuple{ cur_distance0, cur_pt00, cur_pt01 }; + } + return std::tuple{ cur_distance1, cur_pt10, cur_pt11 }; + } + + std::optional< std::tuple< double, Point3D, Point3D > > + triangle_triangle_distance_between_non_conformal_parts( + const Triangle3D& triangle0, const Triangle3D& triangle1 ) + { + const auto [non_colocated_points0, non_colocated_points1] = + find_non_colocated_triangles_points( triangle0, triangle1 ); + if( non_colocated_points0.size() == 0 + || non_colocated_points1.size() == 0 ) + { + return std::nullopt; + } + const auto& vertices0 = triangle0.vertices(); + const auto& vertices1 = triangle1.vertices(); + double min_distance{ std::numeric_limits< double >::max() }; + Point3D point0; + Point3D point1; + if( non_colocated_points0.size() == 1 ) + { + auto [cur_distance, cur_pt] = point_triangle_distance( + vertices0[non_colocated_points0[0]].get(), triangle1 ); + if( cur_distance < min_distance ) + { + min_distance = cur_distance; + point0 = vertices0[non_colocated_points0[0]].get(); + point1 = cur_pt; + } + } + if( non_colocated_points1.size() == 1 ) + { + auto [cur_distance, cur_pt] = point_triangle_distance( + vertices1[non_colocated_points1[0]].get(), triangle0 ); + if( cur_distance < min_distance ) + { + min_distance = cur_distance; + point0 = cur_pt; + point1 = vertices1[non_colocated_points1[0]].get(); + } + } + if( non_colocated_points0.size() > 1 ) + { + auto [cur_distance, cur_pt0, cur_pt1] = test_close_triangles( + non_colocated_points0, triangle0, triangle1 ); + if( cur_distance < min_distance ) + { + min_distance = cur_distance; + point0 = cur_pt0; + point1 = cur_pt1; + } + } + if( non_colocated_points1.size() > 1 ) + { + auto [cur_distance, cur_pt1, cur_pt0] = test_close_triangles( + non_colocated_points1, triangle1, triangle0 ); + if( cur_distance < min_distance ) + { + min_distance = cur_distance; + point0 = cur_pt0; + point1 = cur_pt1; + } + } + return std::tuple{ min_distance, point0, point1 }; + } + std::tuple< double, Point3D > point_tetrahedron_distance( const Point3D& point, const Tetrahedron& tetra ) { @@ -1048,7 +1241,8 @@ namespace geode } OPENGEODE_ASSERT( !circle.plane().normal().inexact_equal( other_direction ), - "[point_circle_distance] Problem while getting circle nearest " + "[point_circle_distance] Problem while getting circle " + "nearest " "point" ); const Vector3D other_projected_on_plane = other_direction @@ -1098,6 +1292,262 @@ namespace geode return point_circle_distance( point, disk ); } + template < index_t dimension > + double Bisector( const geode::index_t number_of_components, + const std::array< double, dimension >& locE, + const std::array< double, dimension >& locY, + std::array< double, dimension >& locX ) + { + double sum_z_squared{ 0 }; + std::array< double, dimension > zPos; + for( const auto i : geode::LRange{ number_of_components } ) + { + zPos[i] = locY[i] / locE[i]; + sum_z_squared += zPos[i] * zPos[i]; + } + if( std::fabs( sum_z_squared - 1 ) < geode::GLOBAL_EPSILON ) + { + for( const auto i : geode::LRange{ number_of_components } ) + { + locX[i] = locY[i]; + } + return 0; + } + const auto emin = locE[number_of_components - 1]; + std::array< double, dimension > numerator; + std::array< double, dimension > pSqr; + pSqr.fill( 0 ); + numerator.fill( 0 ); + for( const auto i : geode::LRange{ number_of_components } ) + { + const auto p = locE[i] / emin; + pSqr[i] = p * p; + numerator[i] = p * zPos[i]; + } + double s{ 0 }; + auto smin = zPos[number_of_components - 1] - 1; + double smax{ 0 }; + if( sum_z_squared >= 1 ) + { + geode::Vector< dimension > v{ numerator }; + smax = v.length() - 1; + } + const geode::index_t jmax = 2048; + geode::index_t j{ 0 }; + while( j < jmax ) + { + s = ( smin + smax ) * 0.5; + if( s == smin || s == smax ) + { + break; + } + double g = -1; + for( const auto i : geode::LRange{ number_of_components } ) + { + const auto ratio = numerator[i] / ( s + pSqr[i] ); + g += ratio * ratio; + } + if( g > 0 ) + { + smin = s; + j++; + } + else if( g < 0 ) + { + smax = s; + j++; + } + else + { + break; + } + } + auto squared_distance = 0.; + for( const auto i : geode::LRange{ number_of_components } ) + { + locX[i] = pSqr[i] * locY[i] / ( s + pSqr[i] ); + const auto diff = locX[i] - locY[i]; + squared_distance += diff * diff; + } + return squared_distance; + } + + template < index_t dimension > + std::tuple< double, std::array< double, dimension > > SqrDistanceSpecial( + const std::array< double, dimension >& extents, + const std::array< double, dimension >& query_point_coordinates ) + { + std::tuple< double, std::array< double, dimension > > result; + auto& [squared_distance, closest_point_coordinates] = result; + std::array< double, dimension > ePos; + std::array< double, dimension > yPos; + std::array< double, dimension > xPos; + index_t numpos{ 0 }; + for( const auto i : geode::LRange{ dimension } ) + { + if( query_point_coordinates[i] > 0 ) + { + ePos.at( numpos ) = extents[i]; + yPos.at( numpos ) = query_point_coordinates[i]; + ++numpos; + continue; + } + closest_point_coordinates.at( numpos ) = 0.; + } + if( query_point_coordinates[dimension - 1] > 0 ) + { + squared_distance = + Bisector< dimension >( numpos, ePos, yPos, xPos ); + } + else + { + std::array< double, dimension - 1 > numer; + std::array< double, dimension - 1 > denom; + const auto extent_squared = + extents[dimension - 1] * extents[dimension - 1]; + for( const auto i : geode::LRange{ numpos } ) + { + numer[i] = ePos[i] * yPos[i]; + denom[i] = ePos[i] * ePos[i] - extent_squared; + } + bool inSubbox{ true }; + for( const auto i : geode::LRange{ numpos } ) + { + if( numer[i] >= denom[i] ) + { + inSubbox = false; + break; + } + } + bool inSubellipsoid{ false }; + double discr{ 1 }; + if( inSubbox ) + { + std::array< double, dimension - 1 > xde; + for( const auto i : geode::LRange{ numpos } ) + { + xde[i] = numer[i] / denom[i]; + discr -= xde[i] * xde[i]; + } + if( discr > 0 ) + { + squared_distance = 0; + for( const auto i : geode::LRange{ numpos } ) + { + xPos[i] = ePos[i] * xde[i]; + const auto diff = xPos[i] - yPos[i]; + squared_distance += diff * diff; + } + closest_point_coordinates[dimension - 1] = + extents[dimension - 1] * std::sqrt( discr ); + squared_distance += + closest_point_coordinates[dimension - 1] + * closest_point_coordinates[dimension - 1]; + inSubellipsoid = true; + } + } + if( !inSubellipsoid ) + { + closest_point_coordinates[dimension - 1] = 0; + squared_distance = + Bisector< dimension >( numpos, ePos, yPos, xPos ); + } + } + numpos = 0; + for( const auto i : geode::LRange{ dimension } ) + { + if( query_point_coordinates[i] > 0 ) + { + closest_point_coordinates[i] = xPos[numpos]; + ++numpos; + } + } + return result; + } + + template < index_t dimension > + std::tuple< double, Point< dimension > > SquaredDistance( + const Ellipse< dimension >& ellipse, + const std::array< double, dimension >& query_point_coordinates ) + { + std::array< bool, dimension > is_query_points_coordinates_negative; + std::tuple< double, Point< dimension > > result; + auto& [squared_distance, closest_point] = result; + for( const auto i : geode::LRange{ dimension } ) + { + is_query_points_coordinates_negative[i] = + ( query_point_coordinates[i] < 0. ); + } + std::array< std::pair< double, index_t >, dimension > + axis_sorted_by_decreasing_extent; + for( const auto i : geode::LRange{ dimension } ) + { + auto& [extent, axis] = axis_sorted_by_decreasing_extent.at( i ); + extent = ellipse.axes().direction( i ).length(); + axis = i; + }; + absl::c_sort( axis_sorted_by_decreasing_extent, + []( const std::pair< double, index_t >& a, + const std::pair< double, index_t >& b ) { + return a.second < b.second; + } ); + std::array< index_t, dimension > reverse_permutation; + for( const auto i : geode::LRange{ dimension } ) + { + reverse_permutation.at( + axis_sorted_by_decreasing_extent[i].second ) = i; + } + std::array< double, dimension > extents; + std::array< double, dimension > locY; + for( const auto i : geode::LRange{ dimension } ) + { + const auto j = axis_sorted_by_decreasing_extent[i].second; + extents[i] = ellipse.axes().direction( j ).length(); + locY[i] = std::fabs( query_point_coordinates[j] ); + } + auto [squared_distance_special, closest_point_coordinates] = + SqrDistanceSpecial< dimension >( extents, locY ); + squared_distance = squared_distance_special; + for( const auto i : geode::LRange{ dimension } ) + { + const auto j = reverse_permutation[i]; + if( is_query_points_coordinates_negative[i] ) + { + closest_point_coordinates[j] = -closest_point_coordinates[j]; + } + } + for( const auto i : geode::LRange{ dimension } ) + { + closest_point.set_value( i, closest_point_coordinates[i] ); + } + return result; + }; + + template < index_t dimension > + std::tuple< double, Point< dimension > > point_ellipse_distance( + const Point< dimension >& point, const Ellipse< dimension >& ellipse ) + { + std::tuple< double, Point< dimension > > result; + std::array< double, dimension > point_coordinates; + const Vector< dimension > center_to_point{ ellipse.center(), point }; + for( const auto i : geode::LRange{ dimension } ) + { + point_coordinates[i] = center_to_point.dot( + ellipse.axes().direction( i ).normalize() ); + } + auto& [distance, closest_point] = result; + const auto [squared_distance, closest_point_result] = + SquaredDistance< dimension >( ellipse, point_coordinates ); + distance = std::sqrt( squared_distance ); + closest_point = ellipse.center(); + for( const auto i : geode::LRange{ dimension } ) + { + closest_point += ellipse.axes().direction( i ).normalize() + * closest_point_result.value( i ); + } + return result; + } + template double opengeode_geometry_api point_point_distance( const Point1D&, const Point1D& ); template double opengeode_geometry_api point_segment_distance( @@ -1123,6 +1573,8 @@ namespace geode point_sphere_signed_distance( const Point2D&, const Sphere2D& ); template std::tuple< double, Point2D > opengeode_geometry_api point_ball_distance( const Point2D&, const Ball2D& ); + template std::tuple< double, Point2D > opengeode_geometry_api + point_ellipse_distance( const Point2D&, const Ellipse2D& ellipse ); template double opengeode_geometry_api point_point_distance( const Point3D&, const Point3D& ); @@ -1142,4 +1594,6 @@ namespace geode point_sphere_signed_distance( const Point3D&, const Sphere3D& ); template std::tuple< double, Point3D > opengeode_geometry_api point_ball_distance( const Point3D&, const Ball3D& ); + template std::tuple< double, Point3D > opengeode_geometry_api + point_ellipse_distance( const Point3D&, const Ellipse3D& ellipse ); } // namespace geode diff --git a/src/geode/geometry/intersection.cpp b/src/geode/geometry/intersection.cpp index fbf939153..3e0d8aafd 100644 --- a/src/geode/geometry/intersection.cpp +++ b/src/geode/geometry/intersection.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -34,6 +35,7 @@ #include #include #include +#include #include namespace @@ -1010,6 +1012,124 @@ namespace geode return { std::move( line ), { first_correctness, second_correctness } }; } + template < index_t dimension > + IntersectionResult< absl::InlinedVector< Point< dimension >, 2 > > + line_ellipse_intersection( const InfiniteLine< dimension >& line, + const Ellipse< dimension >& ellipse ) + { + geode::SquareMatrix< dimension > M( 0 ); + for( const auto i : geode::LRange{ dimension } ) + { + const auto ratio = ellipse.axes().direction( i ).normalize() + * ( 1 / ellipse.axes().direction( i ).length() ); + geode::SquareMatrix< dimension > M_i; + for( const auto j : geode::LRange{ dimension } ) + { + for( const auto k : geode::LRange{ dimension } ) + { + M_i.set_value( j, k, ratio.value( j ) * ratio.value( k ) ); + } + } + M += M_i; + } + const Vector< dimension > diff{ ellipse.center(), line.origin() }; + const auto matDir = M * line.direction(); + const auto matDiff = M * diff; + const auto a0 = diff.dot( matDiff ) - 1; + const auto a1 = line.direction().dot( matDiff ); + const auto a2 = line.direction().dot( matDir ); + const auto discr = a1 * a1 - a0 * a2; + if( discr > GLOBAL_EPSILON ) + { + absl::InlinedVector< Point< dimension >, 2 > results; + const auto root = std::sqrt( discr ); + results.reserve( 2 ); + results.emplace_back( + line.origin() + line.direction() * ( -a1 - root ) / a2 ); + results.emplace_back( + line.origin() + line.direction() * ( -a1 + root ) / a2 ); + typename CorrectnessInfo< absl::InlinedVector< Point< dimension >, + 2 > >::Correctness first_correctness; + first_correctness.first = + point_line_distance( results.front(), line ) <= GLOBAL_EPSILON; + first_correctness.second.push_back( + point_line_projection( results.front(), line ) ); + first_correctness.first = + first_correctness.first + && point_line_distance( results.back(), line ) + <= GLOBAL_EPSILON; + first_correctness.second.push_back( + point_line_projection( results.back(), line ) ); + typename CorrectnessInfo< absl::InlinedVector< Point< dimension >, + 2 > >::Correctness second_correctness; + const auto front_output = + point_ellipse_distance( results.front(), ellipse ); + second_correctness.first = + std::get< 0 >( front_output ) <= GLOBAL_EPSILON; + second_correctness.second.push_back( + std::get< 1 >( front_output ) ); + const auto back_output = + point_ellipse_distance( results.back(), ellipse ); + second_correctness.first = + second_correctness.first + && std::get< 0 >( back_output ) <= GLOBAL_EPSILON; + second_correctness.second.push_back( std::get< 1 >( back_output ) ); + return { std::move( results ), + { first_correctness, second_correctness } }; + } + if( discr > -GLOBAL_EPSILON ) + { + absl::InlinedVector< Point< dimension >, 2 > results; + results.reserve( 1 ); + results.emplace_back( line.origin() - line.direction() * a1 / a2 ); + typename CorrectnessInfo< absl::InlinedVector< Point< dimension >, + 2 > >::Correctness first_correctness; + first_correctness.first = + point_line_distance( results.front(), line ) <= GLOBAL_EPSILON; + first_correctness.second.push_back( + point_line_projection( results.front(), line ) ); + const auto output = + point_ellipse_distance( results.front(), ellipse ); + typename CorrectnessInfo< absl::InlinedVector< Point< dimension >, + 2 > >::Correctness second_correctness; + second_correctness.first = + std::get< 0 >( output ) <= GLOBAL_EPSILON; + second_correctness.second.push_back( std::get< 1 >( output ) ); + return { std::move( results ), + { first_correctness, second_correctness } }; + } + return { INTERSECTION_TYPE::none }; + } + + template < index_t dimension > + IntersectionResult< absl::InlinedVector< Point< dimension >, 2 > > + segment_ellipse_intersection( const Segment< dimension >& segment, + const Ellipse< dimension >& ellipse ) + { + auto line_intersections = line_ellipse_intersection< dimension >( + InfiniteLine< dimension >{ segment }, ellipse ); + if( line_intersections ) + { + absl::InlinedVector< Point< dimension >, 2 > segment_intersections; + segment_intersections.reserve( + line_intersections.result.value().size() ); + for( auto&& point : line_intersections.result.value() ) + { + if( point_segment_distance( point, segment ) <= GLOBAL_EPSILON ) + { + segment_intersections.emplace_back( point ); + } + } + if( segment_intersections.empty() ) + { + return { INTERSECTION_TYPE::none }; + } + return { std::move( segment_intersections ), + std::move( line_intersections.correctness.value() ) }; + } + return line_intersections.type; + } + template IntersectionResult< absl::InlinedVector< Point2D, 2 > > opengeode_geometry_api line_sphere_intersection( const InfiniteLine2D& segment, const Sphere2D& sphere ); @@ -1018,6 +1138,14 @@ namespace geode opengeode_geometry_api segment_sphere_intersection( const Segment2D& segment, const Sphere2D& sphere ); + template IntersectionResult< absl::InlinedVector< Point2D, 2 > > + opengeode_geometry_api line_ellipse_intersection( + const InfiniteLine2D& line, const Ellipse2D& sphere ); + + template IntersectionResult< absl::InlinedVector< Point2D, 2 > > + opengeode_geometry_api segment_ellipse_intersection( + const Segment2D& segment, const Ellipse2D& sphere ); + template IntersectionResult< absl::InlinedVector< Point3D, 2 > > opengeode_geometry_api line_sphere_intersection( const InfiniteLine3D& segment, const Sphere3D& sphere ); @@ -1026,4 +1154,12 @@ namespace geode opengeode_geometry_api segment_sphere_intersection( const Segment3D& segment, const Sphere3D& sphere ); + template IntersectionResult< absl::InlinedVector< Point3D, 2 > > + opengeode_geometry_api line_ellipse_intersection( + const InfiniteLine3D& line, const Ellipse3D& sphere ); + + template IntersectionResult< absl::InlinedVector< Point3D, 2 > > + opengeode_geometry_api segment_ellipse_intersection( + const Segment3D& segment, const Ellipse3D& sphere ); + } // namespace geode \ No newline at end of file diff --git a/src/geode/mesh/helpers/euclidean_distance_transform.cpp b/src/geode/mesh/helpers/euclidean_distance_transform.cpp index c7f96f93d..f5dfa586d 100644 --- a/src/geode/mesh/helpers/euclidean_distance_transform.cpp +++ b/src/geode/mesh/helpers/euclidean_distance_transform.cpp @@ -364,7 +364,8 @@ namespace geode template <> void EuclideanDistanceTransform< 2 >::compute_squared_distance_map() { - ProgressLogger logger{ "Compute 2D euclidian distance", 2 }; + ProgressLogger logger{ Logger::LEVEL::info, + "Compute 2D euclidian distance", 2 }; propagate_directional_squared_distance( 0 ); logger.increment(); combine_squared_distance_components( 1 ); @@ -373,7 +374,8 @@ namespace geode template <> void EuclideanDistanceTransform< 3 >::compute_squared_distance_map() { - ProgressLogger logger{ "Compute 3D euclidian distance", 3 }; + ProgressLogger logger{ Logger::LEVEL::info, + "Compute 3D euclidian distance", 3 }; propagate_directional_squared_distance( 0 ); logger.increment(); combine_squared_distance_components( 1 ); diff --git a/src/geode/model/helpers/component_mesh_edges.cpp b/src/geode/model/helpers/component_mesh_edges.cpp index f11d3a060..d78dae3d5 100644 --- a/src/geode/model/helpers/component_mesh_edges.cpp +++ b/src/geode/model/helpers/component_mesh_edges.cpp @@ -109,6 +109,11 @@ namespace geode const Model& model, const std::array< geode::index_t, 2 >& edge_unique_vertices ) { + if( edge_unique_vertices[0] == NO_ID + || edge_unique_vertices[1] == NO_ID ) + { + return {}; + } const auto line_pairs = model_edge_pairs( model, edge_unique_vertices, geode::Line< Model::dim >::component_type_static() ); @@ -146,6 +151,11 @@ namespace geode surface_component_mesh_edges( const Model& model, const std::array< geode::index_t, 2 >& edge_unique_vertices ) { + if( edge_unique_vertices[0] == NO_ID + || edge_unique_vertices[1] == NO_ID ) + { + return {}; + } const auto surface_pairs = model_edge_pairs( model, edge_unique_vertices, geode::Surface< Model::dim >::component_type_static() ); @@ -184,6 +194,11 @@ namespace geode const geode::BRep& model, const std::array< geode::index_t, 2 >& edge_unique_vertices ) { + if( edge_unique_vertices[0] == NO_ID + || edge_unique_vertices[1] == NO_ID ) + { + return {}; + } const auto block_pairs = model_edge_pairs( model, edge_unique_vertices, geode::Block3D::component_type_static() ); if( block_pairs.empty() ) diff --git a/tests/basic/test-progress-logger.cpp b/tests/basic/test-progress-logger.cpp index fc7e398e9..60055f76c 100644 --- a/tests/basic/test-progress-logger.cpp +++ b/tests/basic/test-progress-logger.cpp @@ -30,7 +30,8 @@ void test() { geode::index_t nb{ 30000 }; - geode::ProgressLogger logger{ "Cool message", nb }; + geode::ProgressLogger logger{ geode::Logger::LEVEL::info, "Cool message", + nb }; for( const auto i : geode::Range{ nb } ) { auto test = 0; diff --git a/tests/geometry/test-distance.cpp b/tests/geometry/test-distance.cpp index 5d70fcf53..609699f90 100644 --- a/tests/geometry/test-distance.cpp +++ b/tests/geometry/test-distance.cpp @@ -24,9 +24,11 @@ #include #include +#include #include #include +#include #include #include #include @@ -1087,19 +1089,131 @@ void test_segment_triangle_distance() "[Test] Wrong result for segment_triangle_distance with seg_hx" ); } +void test_point_ellipse_distance_2d() +{ + const geode::Point2D center{ { 0.0, 0.0 } }; + const geode::Vector2D first_axis{ { 3.0, 0.0 } }; + const geode::Vector2D second_axis{ { 0.0, 2.0 } }; + const geode::Frame2D frame{ { first_axis, second_axis } }; + const geode::Ellipse2D ellipse{ center, frame }; + const geode::Point2D on_boundary{ { 3.0, 0.0 } }; + const geode::Point2D outside1{ { 4.0, 0 } }; + const geode::Point2D outside2{ { 0, 5.0 } }; + const auto [on_boundary_distance, on_boundary_closest_point] = + geode::point_ellipse_distance( on_boundary, ellipse ); + OPENGEODE_EXCEPTION( + on_boundary_distance == 0. && on_boundary_closest_point == on_boundary, + "[Test] Wrong result for point_ellipse_distance_2d with on_boundary " + "point" ); + const auto [outside1_distance, outside1_closest_point] = + geode::point_ellipse_distance( outside1, ellipse ); + const geode::Point2D result_outside1{ { 3, 0 } }; + OPENGEODE_EXCEPTION( + outside1_distance == 1 + && outside1_closest_point.inexact_equal( result_outside1 ), + "[Test] Wrong result for point_ellipse_distance_2d with outside1 " + "point" ); + const auto [outside2_distance, outside2_closest_point] = + geode::point_ellipse_distance( outside2, ellipse ); + const geode::Point2D result_outside2{ { 0, 2 } }; + OPENGEODE_EXCEPTION( + outside2_distance == 3 + && outside2_closest_point.inexact_equal( result_outside2 ), + "[Test] Wrong result for point_ellipse_distance_2d with outside2 " + "point" ); +} + +void test_point_ellipse_distance_2d_not_aligned() +{ + const geode::Point2D center{ { 0.0, 0.0 } }; + const geode::Vector2D first_axis{ { 1.0, 1.0 } }; + const geode::Vector2D second_axis{ { -2.0, 2.0 } }; + const geode::Frame2D frame{ { first_axis, second_axis } }; + const geode::Ellipse2D ellipse{ center, frame }; + const geode::Point2D point{ { 2.0, 2.0 } }; + const auto [distance, closest_point] = + geode::point_ellipse_distance( point, ellipse ); + const geode::Point2D result{ { 1.0, 1.0 } }; + OPENGEODE_EXCEPTION( + std::fabs( std::sqrt( 2 ) - distance ) < geode::GLOBAL_EPSILON + && closest_point.inexact_equal( result ), + "[Test] Wrong result for point_ellipse_distance_2d_not_aligned" ); +} + +void test_point_ellipse_distance_3d() +{ + const geode::Point3D center{ { 0.0, 0.0, 0.0 } }; + const geode::Vector3D first_axis{ { 3.0, 0.0, 0.0 } }; + const geode::Vector3D second_axis{ { 0.0, 2.0, 0.0 } }; + const geode::Vector3D third_axis{ { 0.0, 0.0, 1.0 } }; + const geode::Frame3D frame{ { first_axis, second_axis, third_axis } }; + const geode::Ellipse3D ellipse{ center, frame }; + const geode::Point3D on_boundary{ { 3.0, 0.0, 0.0 } }; + const geode::Point3D outside1{ { 4.0, 0, 0 } }; + const geode::Point3D outside2{ { 0, 5.0, 0 } }; + const geode::Point3D outside3{ { 0, 0, 3 } }; + const auto [on_boundary_distance, on_boundary_closest_point] = + geode::point_ellipse_distance( on_boundary, ellipse ); + OPENGEODE_EXCEPTION( + on_boundary_distance == 0. && on_boundary_closest_point == on_boundary, + "[Test] Wrong result for point_ellipse_distance_3d with " + "on_boundary " + "point" ); + const auto [outside1_distance, outside1_closest_point] = + geode::point_ellipse_distance( outside1, ellipse ); + const geode::Point3D result_outside1{ { 3, 0, 0 } }; + OPENGEODE_EXCEPTION( + outside1_distance == 1 + && outside1_closest_point.inexact_equal( result_outside1 ), + "[Test] Wrong result for point_ellipse_distance_3d with outside1 " + "point" ); + const auto [outside2_distance, outside2_closest_point] = + geode::point_ellipse_distance( outside2, ellipse ); + const geode::Point3D result_outside2{ { 0, 2, 0 } }; + OPENGEODE_EXCEPTION( + outside2_distance == 3 + && outside2_closest_point.inexact_equal( result_outside2 ), + "[Test] Wrong result for point_ellipse_distance_3d with outside2 " + "point" ); + const auto [outside3_distance, outside3_closest_point] = + geode::point_ellipse_distance( outside3, ellipse ); + const geode::Point3D result_outside3{ { 0, 0, 1 } }; + OPENGEODE_EXCEPTION( + outside3_distance == 2 + && outside3_closest_point.inexact_equal( result_outside3 ), + "[Test] Wrong result for point_ellipse_distance_3d with outside3 " + "point" ); + const auto [inside_distance, inside_closest_point] = + geode::point_ellipse_distance( center, ellipse ); + const geode::Point3D result_inside{ { 0, 0, 1 } }; + OPENGEODE_EXCEPTION( + inside_distance == 1 + && inside_closest_point.inexact_equal( result_inside ), + "[Test] Wrong result for point_ellipse_distance_3d with inside " + "point" ); +} + +void test_point_ellipse_distance() +{ + test_point_ellipse_distance_2d(); + test_point_ellipse_distance_3d(); + test_point_ellipse_distance_2d_not_aligned(); +} + void test() { - test_point_segment_distance(); - test_segment_segment_distance(); - test_segment_line_distance(); - test_point_line_distance(); - test_point_triangle_distance(); - test_point_tetrahedron_distance(); - test_point_plane_distance(); - test_point_sphere_distance(); - test_point_circle_distance(); - test_line_triangle_distance(); - test_segment_triangle_distance(); + // test_point_segment_distance(); + // test_segment_segment_distance(); + // test_segment_line_distance(); + // test_point_line_distance(); + // test_point_triangle_distance(); + // test_point_tetrahedron_distance(); + // test_point_plane_distance(); + // test_point_sphere_distance(); + // test_point_circle_distance(); + // test_line_triangle_distance(); + // test_segment_triangle_distance(); + test_point_ellipse_distance(); } OPENGEODE_TEST( "distance" ) diff --git a/tests/geometry/test-intersection.cpp b/tests/geometry/test-intersection.cpp index 213a03bb7..b82e48987 100644 --- a/tests/geometry/test-intersection.cpp +++ b/tests/geometry/test-intersection.cpp @@ -26,12 +26,14 @@ #include #include +#include #include #include #include #include #include #include +#include #include #include @@ -750,6 +752,216 @@ void test_triangle_circle_intersection() !XY0_result, "[Test] Wrong intersection between planeZO, circleYZ0" ); } +void test_line_ellipse_intersection_2D() +{ + const geode::Point2D center{ { 0.0, 0.0 } }; + const geode::Vector2D first_axis{ { 3.0, 0.0 } }; + const geode::Vector2D second_axis{ { 0.0, 2.0 } }; + const geode::Frame2D frame{ { first_axis, second_axis } }; + const geode::Ellipse2D ellipse{ center, frame }; + const geode::Point2D parallel_origin{ { 0., 3 } }; + const geode::Point2D parallel_dest{ { 1, 3 } }; + const geode::Segment2D parallel_segment{ parallel_origin, parallel_dest }; + const geode::InfiniteLine2D parallel{ parallel_segment }; + const auto parallel_result = + geode::line_ellipse_intersection( parallel, ellipse ); + OPENGEODE_EXCEPTION( !parallel_result, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Point2D intersecting_once_origin{ { 3.0, 2.0 } }; + const geode::Point2D intersecting_once_dest{ { 3.0, 0.0 } }; + const geode::Segment2D intersecting_once_segment{ intersecting_once_origin, + intersecting_once_dest }; + const geode::InfiniteLine2D intersecting_once{ intersecting_once_segment }; + const auto intersecting_once_result = + geode::line_ellipse_intersection( intersecting_once, ellipse ); + OPENGEODE_EXCEPTION( intersecting_once_result.result->size() == 1, + "[Test] Wrong intersection between line and ellipse" ); + OPENGEODE_EXCEPTION( + intersecting_once_result.result->at( 0 ) == intersecting_once_dest, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Segment2D intersecting_twice_segment( + center, geode::Point2D{ { 3, 0 } } ); + const geode::InfiniteLine2D intersecting_twice_line{ + intersecting_twice_segment + }; + const auto intersecting_twice_result = + geode::line_ellipse_intersection( intersecting_twice_line, ellipse ); + OPENGEODE_EXCEPTION( intersecting_twice_result.result->size() == 2, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Point2D result1{ { -3, 0 } }; + OPENGEODE_EXCEPTION( intersecting_twice_result.result->at( 0 ) == result1, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Point2D result2{ { 3, 0 } }; + OPENGEODE_EXCEPTION( intersecting_twice_result.result->at( 1 ) == result2, + "[Test] Wrong intersection between line and ellipse" ); +} + +void test_line_ellipse_intersection_3D() +{ + const geode::Point3D center{ { 0.0, 0.0 } }; + const geode::Vector3D first_axis{ { 3.0, 0.0, 0.0 } }; + const geode::Vector3D second_axis{ { 0.0, 2.0, 0.0 } }; + const geode::Vector3D third_axis{ { 0.0, 0.0, 1.0 } }; + const geode::Frame3D frame{ { first_axis, second_axis, third_axis } }; + const geode::Ellipse3D ellipse{ center, frame }; + + const geode::Point3D parallel_origin{ { 0., 3, 1 } }; + const geode::Point3D parallel_dest{ { 1, 3, 1 } }; + const geode::Segment3D parallel_segment{ parallel_origin, parallel_dest }; + const geode::InfiniteLine3D parallel{ parallel_segment }; + const auto parallel_result = + geode::line_ellipse_intersection( parallel, ellipse ); + OPENGEODE_EXCEPTION( !parallel_result, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Point3D intersecting_once_origin{ { 3.0, 2.0, 0.0 } }; + const geode::Point3D intersecting_once_dest{ { 3.0, 0.0, 0.0 } }; + const geode::Segment3D intersecting_once_segment{ intersecting_once_origin, + intersecting_once_dest }; + const geode::InfiniteLine3D intersecting_once{ intersecting_once_segment }; + const auto intersecting_once_result = + geode::line_ellipse_intersection( intersecting_once, ellipse ); + OPENGEODE_EXCEPTION( intersecting_once_result.result->size() == 1, + "[Test] Wrong intersection between line and ellipse" ); + OPENGEODE_EXCEPTION( + intersecting_once_result.result->at( 0 ) == intersecting_once_dest, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Segment3D intersecting_twice_segment( + geode::Point3D{ { 0, 0, 1 } }, geode::Point3D{ { 3, 0, 0 } } ); + const geode::InfiniteLine3D intersecting_twice_line{ + intersecting_twice_segment + }; + const auto intersecting_twice_result = + geode::line_ellipse_intersection( intersecting_twice_line, ellipse ); + OPENGEODE_EXCEPTION( intersecting_twice_result.result->size() == 2, + "[Test] Wrong intersection between line and ellipse" ); + const geode::Point3D result1{ { 0, 0, 1 } }; + OPENGEODE_EXCEPTION( + intersecting_twice_result.result->at( 0 ).inexact_equal( result1 ), + "[Test] Wrong intersection between line and ellipse" ); + const geode::Point3D result2{ { 3, 0, 0 } }; + OPENGEODE_EXCEPTION( + intersecting_twice_result.result->at( 1 ).inexact_equal( result2 ), + "[Test] Wrong intersection between line and ellipse" ); +} + +void test_line_ellipse_intersection() +{ + test_line_ellipse_intersection_2D(); + test_line_ellipse_intersection_3D(); +} + +void test_segment_ellipse_intersection_2D() +{ + const geode::Point2D center{ { 0.0, 0.0 } }; + const geode::Vector2D first_axis{ { 3.0, 0.0 } }; + const geode::Vector2D second_axis{ { 0.0, 2.0 } }; + const geode::Frame2D frame{ { first_axis, second_axis } }; + const geode::Ellipse2D ellipse{ center, frame }; + const geode::Point2D parallel_origin{ { 0., 3 } }; + const geode::Point2D parallel_dest{ { 1, 3 } }; + const geode::Segment2D parallel_segment{ parallel_origin, parallel_dest }; + const auto parallel_result = + geode::segment_ellipse_intersection( parallel_segment, ellipse ); + OPENGEODE_EXCEPTION( !parallel_result, + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point2D intersecting_once_origin{ { 3.0, 2.0 } }; + const geode::Point2D intersecting_once_dest{ { 3.0, -2.0 } }; + const geode::Segment2D intersecting_once_segment{ intersecting_once_origin, + intersecting_once_dest }; + const auto intersecting_once_result = geode::segment_ellipse_intersection( + intersecting_once_segment, ellipse ); + OPENGEODE_EXCEPTION( intersecting_once_result.result->size() == 1, + "[Test] Wrong intersection between segment and ellipse" ); + geode::Point2D result{ { 3.0, 0.0 } }; + OPENGEODE_EXCEPTION( intersecting_once_result.result->at( 0 ) == result, + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point2D intersecting_twice_origin{ { -4.0, 0 } }; + const geode::Point2D intersecting_twice_dest{ { 4.0, 0 } }; + const geode::Segment2D intersecting_twice_segment{ + intersecting_twice_origin, intersecting_twice_dest + }; + const auto intersecting_twice_result = geode::segment_ellipse_intersection( + intersecting_twice_segment, ellipse ); + OPENGEODE_EXCEPTION( intersecting_twice_result.result->size() == 2, + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point2D result1{ { -3, 0 } }; + OPENGEODE_EXCEPTION( + intersecting_twice_result.result->at( 0 ).inexact_equal( result1 ), + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point2D result2{ { 3, 0 } }; + OPENGEODE_EXCEPTION( + intersecting_twice_result.result->at( 1 ).inexact_equal( result2 ), + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point2D not_intersecting_origin{ { 1.0, 0.0 } }; + const geode::Point2D not_intersecting_dest{ { 2.0, 0.0 } }; + const geode::Segment2D not_intersecting_segment{ not_intersecting_origin, + not_intersecting_dest }; + const auto not_intersecting_result = geode::segment_ellipse_intersection( + not_intersecting_segment, ellipse ); + OPENGEODE_EXCEPTION( !not_intersecting_result, + "[Test] Wrong intersection between segment and ellipse" ); +} + +void test_segment_ellipse_intersection_3D() +{ + const geode::Point3D center{ { 0.0, 0.0, 0.0 } }; + const geode::Vector3D first_axis{ { 3.0, 0.0, 0.0 } }; + const geode::Vector3D second_axis{ { 0.0, 2.0, 0.0 } }; + const geode::Vector3D third_axis{ { 0.0, 0.0, 1.0 } }; + const geode::Frame3D frame{ { first_axis, second_axis, third_axis } }; + const geode::Ellipse3D ellipse{ center, frame }; + + const geode::Point3D parallel_origin{ { 0., 3, 1 } }; + const geode::Point3D parallel_dest{ { 1, 3, 1 } }; + const geode::Segment3D parallel_segment{ parallel_origin, parallel_dest }; + const auto parallel_result = + geode::segment_ellipse_intersection( parallel_segment, ellipse ); + OPENGEODE_EXCEPTION( !parallel_result, + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point3D intersecting_once_origin{ { 3.0, 2.0, 0 } }; + const geode::Point3D intersecting_once_dest{ { 3.0, -2.0, 0 } }; + const geode::Segment3D intersecting_once_segment{ intersecting_once_origin, + intersecting_once_dest }; + const auto intersecting_once_result = geode::segment_ellipse_intersection( + intersecting_once_segment, ellipse ); + OPENGEODE_EXCEPTION( intersecting_once_result.result->size() == 1, + "[Test] Wrong intersection between segment and ellipse" ); + geode::Point3D result{ { 3.0, 0.0, 0 } }; + OPENGEODE_EXCEPTION( intersecting_once_result.result->at( 0 ) == result, + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point3D intersecting_twice_origin{ { -4.0, 0, 0 } }; + const geode::Point3D intersecting_twice_dest{ { 4.0, 0, 0 } }; + const geode::Segment3D intersecting_twice_segment{ + intersecting_twice_origin, intersecting_twice_dest + }; + const auto intersecting_twice_result = geode::segment_ellipse_intersection( + intersecting_twice_segment, ellipse ); + OPENGEODE_EXCEPTION( intersecting_twice_result.result->size() == 2, + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point3D result1{ { -3, 0 } }; + OPENGEODE_EXCEPTION( + intersecting_twice_result.result->at( 0 ).inexact_equal( result1 ), + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point3D result2{ { 3, 0 } }; + OPENGEODE_EXCEPTION( + intersecting_twice_result.result->at( 1 ).inexact_equal( result2 ), + "[Test] Wrong intersection between segment and ellipse" ); + const geode::Point3D not_intersecting_origin{ { 1.0, 0.0, 0.0 } }; + const geode::Point3D not_intersecting_dest{ { 2.0, 0.0, 0.5 } }; + const geode::Segment3D not_intersecting_segment{ not_intersecting_origin, + not_intersecting_dest }; + const auto not_intersecting_result = geode::segment_ellipse_intersection( + not_intersecting_segment, ellipse ); + OPENGEODE_EXCEPTION( !not_intersecting_result, + "[Test] Wrong intersection between segment and ellipse" ); +} + +void test_segment_ellipse_intersection() +{ + test_segment_ellipse_intersection_2D(); + test_segment_ellipse_intersection_3D(); +} + void test() { test_line_sphere_intersection(); @@ -764,6 +976,8 @@ void test() test_plane_plane_intersection(); test_plane_circle_intersection(); test_triangle_circle_intersection(); + test_line_ellipse_intersection(); + test_segment_ellipse_intersection(); } OPENGEODE_TEST( "intersection" )