diff --git a/CMakeLists.txt b/CMakeLists.txt index 498caee7e..f8b658be3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -99,6 +99,7 @@ option(LIBIGL_USE_STATIC_LIBRARY "Use libigl as static library" ${LIBIGL_TOPLEVE # Permissive modules. These modules are available under MPL2 license, and their dependencies are available # under a permissive or public domain license. +option(LIBIGL_CYCODEBASE "Build target igl::cycodebase" ${LIBIGL_TOPLEVEL_PROJECT}) option(LIBIGL_EMBREE "Build target igl::embree" ${LIBIGL_TOPLEVEL_PROJECT}) option(LIBIGL_GLFW "Build target igl::glfw" ${LIBIGL_TOPLEVEL_PROJECT}) option(LIBIGL_IMGUI "Build target igl::imgui" ${LIBIGL_TOPLEVEL_PROJECT}) diff --git a/cmake/igl/modules/cycodebase.cmake b/cmake/igl/modules/cycodebase.cmake new file mode 100644 index 000000000..db61abe38 --- /dev/null +++ b/cmake/igl/modules/cycodebase.cmake @@ -0,0 +1,26 @@ +# 1. Define module +igl_add_library(igl_cycodebase) + +# 2. Include headers +include(GNUInstallDirs) +target_include_directories(igl_cycodebase ${IGL_SCOPE} + $ + $ +) + +# 3. Target sources +file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/cycodebase/*.h") +file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/cycodebase/*.cpp") +igl_target_sources(igl_cycodebase ${INC_FILES} ${SRC_FILES}) + +# 4. Dependencies +include(cycodebase) +target_link_libraries(igl_cycodebase ${IGL_SCOPE} + igl::core + cyCodeBase::cyCodeBase +) + +# 5. Unit tests +file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/cycodebase/*.cpp") +igl_add_test(igl_cycodebase ${SRC_FILES}) + diff --git a/cmake/libigl.cmake b/cmake/libigl.cmake index 41b71e519..88569488d 100644 --- a/cmake/libigl.cmake +++ b/cmake/libigl.cmake @@ -13,6 +13,7 @@ include(igl_windows) # Libigl permissive modules igl_include(core) +igl_include_optional(cycodebase) igl_include_optional(embree) igl_include_optional(opengl) igl_include_optional(glfw) diff --git a/cmake/recipes/external/cycodebase.cmake b/cmake/recipes/external/cycodebase.cmake new file mode 100644 index 000000000..ebc951b17 --- /dev/null +++ b/cmake/recipes/external/cycodebase.cmake @@ -0,0 +1,18 @@ +if(TARGET cycodebase::cycodebase) + return() +endif() + +FetchContent_Declare( + cyCodeBase + GIT_REPOSITORY https://github.com/cemyuksel/cyCodeBase/ + GIT_TAG e36f3cffca65eb12a8a071f0443128b7de6ed75d +) +FetchContent_Populate(cyCodeBase) +add_library(cyCodeBase_interface INTERFACE) +target_include_directories(cyCodeBase_interface INTERFACE ${cycodebase_SOURCE_DIR}) + +if(NOT (CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|i[3-6]86")) + target_compile_definitions(cyCodeBase_interface INTERFACE CY_NO_INTRIN_H) +endif() + +add_library(cyCodeBase::cyCodeBase ALIAS cyCodeBase_interface) diff --git a/cmake/recipes/external/libigl_tutorial_data.cmake b/cmake/recipes/external/libigl_tutorial_data.cmake index c7492d69e..a3e1f17fb 100644 --- a/cmake/recipes/external/libigl_tutorial_data.cmake +++ b/cmake/recipes/external/libigl_tutorial_data.cmake @@ -8,7 +8,7 @@ include(FetchContent) FetchContent_Declare( libigl_tutorial_data GIT_REPOSITORY https://github.com/libigl/libigl-tutorial-data - GIT_TAG 6700bf49000cd64199835fb40323e5ca9c7796ab + GIT_TAG 644dd4104843b6d736745d9dafbd70bf8d175648 ) FetchContent_MakeAvailable(libigl_tutorial_data) diff --git a/include/igl/predicates/Orientation.h b/include/igl/Orientation.h similarity index 54% rename from include/igl/predicates/Orientation.h rename to include/igl/Orientation.h index 078cd9517..6c7c50c1f 100644 --- a/include/igl/predicates/Orientation.h +++ b/include/igl/Orientation.h @@ -7,23 +7,21 @@ // v. 2.0. If a copy of the MPL was not distributed with this file, You can // obtain one at http://mozilla.org/MPL/2.0/. #pragma once -#ifndef IGL_PREDICATES_ORIENTATION_H -#define IGL_PREDICATES_ORIENTATION_H +#ifndef IGL_ORIENTATION_H +#define IGL_ORIENTATION_H -#include "../igl_inline.h" +#include "igl_inline.h" #include namespace igl { - namespace predicates { - /// Types of orientations and other predicate results. - /// - /// \fileinfo - enum class Orientation { - POSITIVE=1, INSIDE=1, - NEGATIVE=-1, OUTSIDE=-1, - COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0 - }; - } + /// Types of orientations and other predicate results. + /// + /// \fileinfo + enum class Orientation { + POSITIVE=1, INSIDE=1, + NEGATIVE=-1, OUTSIDE=-1, + COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0 + }; } diff --git a/include/igl/box_simplices.cpp b/include/igl/box_simplices.cpp index 01a84336c..123737825 100644 --- a/include/igl/box_simplices.cpp +++ b/include/igl/box_simplices.cpp @@ -6,6 +6,7 @@ // v. 2.0. If a copy of the MPL was not distributed with this file, You can // obtain one at http://mozilla.org/MPL/2.0/. #include "box_simplices.h" +#include "parallel_for.h" template < typename DerivedV, @@ -20,7 +21,8 @@ IGL_INLINE void igl::box_simplices( { B1.setConstant(F.rows(),V.cols(),std::numeric_limits::infinity()); B2.setConstant(F.rows(),V.cols(),-std::numeric_limits::infinity()); - for(int f = 0; f < F.rows(); f++) + //for(int f = 0; f < F.rows(); f++) + igl::parallel_for(F.rows(),[&](const int f) { for(int c = 0; c < F.cols(); c++) { @@ -30,7 +32,7 @@ IGL_INLINE void igl::box_simplices( B2(f,d) = std::max(B2(f,d),V(F(f,c),d)); } } - } + },1000); } #ifdef IGL_STATIC_LIBRARY diff --git a/include/igl/colormap.cpp b/include/igl/colormap.cpp index ff9e433e4..e881fc867 100644 --- a/include/igl/colormap.cpp +++ b/include/igl/colormap.cpp @@ -1567,6 +1567,266 @@ static double parula_cm[256][3] = { { 0.9736, 0.9752, 0.0597 }, { 0.9763, 0.9831, 0.0538 } }; +static double zoe_cm[256][3] = { + {1.000000,1.000000,1.000000}, +{0.944653,0.947049,0.948510}, +{0.988743,0.993787,0.996863}, +{0.933958,0.941146,0.945529}, +{0.977485,0.987574,0.993725}, +{0.923264,0.935244,0.942549}, +{0.966228,0.981361,0.990588}, +{0.912569,0.929341,0.939568}, +{0.954970,0.975148,0.987451}, +{0.901874,0.923439,0.936588}, +{0.943713,0.968934,0.984313}, +{0.891180,0.917536,0.933608}, +{0.932455,0.962721,0.981176}, +{0.880485,0.911634,0.930627}, +{0.921198,0.956508,0.978039}, +{0.869791,0.905732,0.927647}, +{0.909941,0.950295,0.974901}, +{0.859096,0.899829,0.924666}, +{0.898683,0.944082,0.971764}, +{0.848402,0.893927,0.921686}, +{0.887426,0.937869,0.968627}, +{0.837707,0.888024,0.918705}, +{0.876168,0.931656,0.965490}, +{0.827013,0.882122,0.915725}, +{0.864911,0.925443,0.962352}, +{0.816318,0.876219,0.912744}, +{0.853653,0.919229,0.959215}, +{0.805623,0.870317,0.909764}, +{0.842396,0.913016,0.956078}, +{0.794929,0.864414,0.906783}, +{0.831138,0.906803,0.952940}, +{0.784234,0.858512,0.903803}, +{0.819881,0.900590,0.949803}, +{0.773540,0.852609,0.900823}, +{0.808624,0.894377,0.946666}, +{0.762845,0.846707,0.897842}, +{0.797366,0.888164,0.943528}, +{0.752151,0.840804,0.894862}, +{0.786109,0.881951,0.940391}, +{0.741456,0.834902,0.891881}, +{0.774851,0.875738,0.937254}, +{0.730761,0.829000,0.888901}, +{0.763594,0.869524,0.934116}, +{0.720067,0.823097,0.885920}, +{0.752336,0.863311,0.930979}, +{0.709372,0.817195,0.882940}, +{0.741079,0.857098,0.927842}, +{0.698678,0.811292,0.879959}, +{0.729822,0.850885,0.924704}, +{0.687983,0.805390,0.876979}, +{0.718564,0.844672,0.921567}, +{0.677289,0.799487,0.873998}, +{0.707307,0.838459,0.918430}, +{0.666594,0.793585,0.871018}, +{0.696049,0.832246,0.915292}, +{0.655900,0.787682,0.868038}, +{0.684792,0.826033,0.912155}, +{0.645205,0.781780,0.865057}, +{0.673534,0.819820,0.909018}, +{0.634510,0.775877,0.862077}, +{0.662277,0.813606,0.905880}, +{0.623816,0.769975,0.859096}, +{0.651020,0.807393,0.902743}, +{0.613121,0.764072,0.856116}, +{0.639762,0.801180,0.899606}, +{0.602427,0.758170,0.853135}, +{0.628505,0.794967,0.896469}, +{0.591732,0.752267,0.850155}, +{0.617247,0.788754,0.893331}, +{0.581038,0.746365,0.847174}, +{0.605990,0.782541,0.890194}, +{0.570343,0.740463,0.844194}, +{0.594732,0.776328,0.887057}, +{0.559648,0.734560,0.841214}, +{0.583475,0.770115,0.883919}, +{0.548954,0.728658,0.838233}, +{0.572218,0.763901,0.880782}, +{0.538259,0.722755,0.835253}, +{0.560960,0.757688,0.877645}, +{0.527565,0.716853,0.832272}, +{0.549703,0.751475,0.874507}, +{0.516870,0.710950,0.829292}, +{0.538445,0.745262,0.871370}, +{0.506176,0.705048,0.826311}, +{0.527188,0.739049,0.868233}, +{0.495481,0.699145,0.823331}, +{0.515930,0.732836,0.865095}, +{0.484787,0.693243,0.820350}, +{0.504673,0.726623,0.861958}, +{0.474092,0.687340,0.817370}, +{0.493415,0.720410,0.858821}, +{0.463397,0.681438,0.814389}, +{0.482158,0.714197,0.855683}, +{0.452703,0.675535,0.811409}, +{0.470901,0.707983,0.852546}, +{0.442008,0.669633,0.808429}, +{0.459643,0.701770,0.849409}, +{0.431314,0.663731,0.805448}, +{0.448386,0.695557,0.846271}, +{0.420619,0.657828,0.802468}, +{0.437128,0.689344,0.843134}, +{0.409925,0.651926,0.799487}, +{0.425871,0.683131,0.839997}, +{0.399230,0.646023,0.796507}, +{0.414613,0.676918,0.836859}, +{0.388535,0.640121,0.793526}, +{0.403356,0.670705,0.833722}, +{0.377841,0.634218,0.790546}, +{0.392099,0.664492,0.830585}, +{0.367146,0.628316,0.787565}, +{0.380841,0.658278,0.827448}, +{0.356452,0.622413,0.784585}, +{0.369584,0.652065,0.824310}, +{0.345757,0.616511,0.781604}, +{0.358326,0.645852,0.821173}, +{0.335063,0.610608,0.778624}, +{0.347069,0.639639,0.818036}, +{0.324368,0.604706,0.775644}, +{0.335811,0.633426,0.814898}, +{0.313674,0.598803,0.772663}, +{0.324554,0.627213,0.811761}, +{0.302979,0.592901,0.769683}, +{0.313297,0.621000,0.808624}, +{0.292284,0.586999,0.766702}, +{0.302039,0.614787,0.805486}, +{0.281590,0.581096,0.763722}, +{0.290782,0.608573,0.802349}, +{0.270895,0.575194,0.760741}, +{1.000000,0.978342,0.887274}, +{0.950000,0.925509,0.837329}, +{1.000000,0.970099,0.875524}, +{0.950000,0.917678,0.826167}, +{1.000000,0.961856,0.863775}, +{0.950000,0.909847,0.815005}, +{1.000000,0.953612,0.852025}, +{0.950000,0.902016,0.803843}, +{1.000000,0.945369,0.840275}, +{0.950000,0.894185,0.792681}, +{1.000000,0.937126,0.828526}, +{0.950000,0.886354,0.781519}, +{1.000000,0.928883,0.816776}, +{0.950000,0.878523,0.770356}, +{1.000000,0.920640,0.805027}, +{0.950000,0.870692,0.759194}, +{1.000000,0.912397,0.793277}, +{0.950000,0.862861,0.748032}, +{1.000000,0.904154,0.781528}, +{0.950000,0.855030,0.736870}, +{1.000000,0.895910,0.769778}, +{0.950000,0.847199,0.725708}, +{1.000000,0.887667,0.758028}, +{0.950000,0.839368,0.714546}, +{1.000000,0.879424,0.746279}, +{0.950000,0.831537,0.703384}, +{1.000000,0.871181,0.734529}, +{0.950000,0.823706,0.692222}, +{1.000000,0.862938,0.722780}, +{0.950000,0.815875,0.681060}, +{1.000000,0.854695,0.711030}, +{0.950000,0.808044,0.669898}, +{1.000000,0.846452,0.699281}, +{0.950000,0.800213,0.658736}, +{1.000000,0.838208,0.687531}, +{0.950000,0.792382,0.647574}, +{1.000000,0.829965,0.675782}, +{0.950000,0.784551,0.636411}, +{1.000000,0.821722,0.664032}, +{0.950000,0.776720,0.625249}, +{1.000000,0.813479,0.652282}, +{0.950000,0.768889,0.614087}, +{1.000000,0.805236,0.640533}, +{0.950000,0.761058,0.602925}, +{1.000000,0.796993,0.628783}, +{0.950000,0.753227,0.591763}, +{1.000000,0.788749,0.617034}, +{0.950000,0.745396,0.580601}, +{1.000000,0.780506,0.605284}, +{0.950000,0.737566,0.569439}, +{1.000000,0.772263,0.593535}, +{0.950000,0.729735,0.558277}, +{1.000000,0.764020,0.581785}, +{0.950000,0.721904,0.547115}, +{1.000000,0.755777,0.570035}, +{0.950000,0.714073,0.535953}, +{1.000000,0.747534,0.558286}, +{0.950000,0.706242,0.524791}, +{1.000000,0.739291,0.546536}, +{0.950000,0.698411,0.513629}, +{1.000000,0.731047,0.534787}, +{0.950000,0.690580,0.502466}, +{1.000000,0.722804,0.523037}, +{0.950000,0.682749,0.491304}, +{1.000000,0.714561,0.511288}, +{0.950000,0.674918,0.480142}, +{1.000000,0.706318,0.499538}, +{0.950000,0.667087,0.468980}, +{1.000000,0.698075,0.487789}, +{0.950000,0.659256,0.457818}, +{1.000000,0.689832,0.476039}, +{0.950000,0.651425,0.446656}, +{1.000000,0.681589,0.464289}, +{0.950000,0.643594,0.435494}, +{1.000000,0.673345,0.452540}, +{0.950000,0.635763,0.424332}, +{1.000000,0.665102,0.440790}, +{0.950000,0.627932,0.413170}, +{1.000000,0.656859,0.429041}, +{0.950000,0.620101,0.402008}, +{1.000000,0.648616,0.417291}, +{0.950000,0.612270,0.390846}, +{1.000000,0.640373,0.405542}, +{0.950000,0.604439,0.379683}, +{1.000000,0.632130,0.393792}, +{0.950000,0.596608,0.368521}, +{1.000000,0.623887,0.382042}, +{0.950000,0.588777,0.357359}, +{1.000000,0.615643,0.370293}, +{0.950000,0.580946,0.346197}, +{1.000000,0.607400,0.358543}, +{0.950000,0.573115,0.335035}, +{1.000000,0.599157,0.346794}, +{0.950000,0.565284,0.323873}, +{1.000000,0.590914,0.335044}, +{0.950000,0.557453,0.312711}, +{1.000000,0.582671,0.323295}, +{0.950000,0.549622,0.301549}, +{1.000000,0.574428,0.311545}, +{0.950000,0.541791,0.290387}, +{1.000000,0.566184,0.299796}, +{0.950000,0.533960,0.279225}, +{1.000000,0.557941,0.288046}, +{0.950000,0.526129,0.268063}, +{1.000000,0.549698,0.276296}, +{0.950000,0.518298,0.256901}, +{1.000000,0.541455,0.264547}, +{0.950000,0.510467,0.245738}, +{1.000000,0.533212,0.252797}, +{0.950000,0.502636,0.234576}, +{1.000000,0.524969,0.241048}, +{0.950000,0.494805,0.223414}, +{1.000000,0.516726,0.229298}, +{0.950000,0.486974,0.212252}, +{1.000000,0.508482,0.217549}, +{0.950000,0.479143,0.201090}, +{1.000000,0.500239,0.205799}, +{0.950000,0.471312,0.189928}, +{1.000000,0.491996,0.194050}, +{0.950000,0.463481,0.178766}, +{1.000000,0.483753,0.182300}, +{0.950000,0.455650,0.167604}, +{1.000000,0.475510,0.170550}, +{0.950000,0.447819,0.156442}, +{1.000000,0.467267,0.158801}, +{0.950000,0.439988,0.145280}, +{1.000000,0.459024,0.147051}, +{0.950000,0.432157,0.134118} +}; + + } template @@ -1602,6 +1862,9 @@ IGL_INLINE void igl::colormap( case COLOR_MAP_TYPE_VIRIDIS: colormap(viridis_cm, x_in, r, g, b); break; + case COLOR_MAP_TYPE_ZOE: + colormap(zoe_cm, x_in, r, g, b); + break; default: throw std::invalid_argument("igl::colormap(): Selected colormap is unsupported!"); break; diff --git a/include/igl/colormap.h b/include/igl/colormap.h index 68a304a97..fd8c97f2b 100644 --- a/include/igl/colormap.h +++ b/include/igl/colormap.h @@ -24,7 +24,8 @@ namespace igl { COLOR_MAP_TYPE_PLASMA = 4, COLOR_MAP_TYPE_VIRIDIS = 5, COLOR_MAP_TYPE_TURBO = 6, - NUM_COLOR_MAP_TYPES = 7 + COLOR_MAP_TYPE_ZOE = 7, + NUM_COLOR_MAP_TYPES = 8 }; /// Compute [r,g,b] values of the selected colormap for /// a given factor f between 0 and 1 diff --git a/include/igl/cubic.cpp b/include/igl/cubic.cpp new file mode 100644 index 000000000..59a0c8193 --- /dev/null +++ b/include/igl/cubic.cpp @@ -0,0 +1,39 @@ +#include "cubic.h" + template + < + typename DerivedC, + typename DerivedP> + IGL_INLINE void igl::cubic( + const Eigen::MatrixBase& C, + const typename DerivedP::Scalar & t, + Eigen::PlainObjectBase& P) +{ + // static assert that C has 4 rows or Dynamic + static_assert( + DerivedC::RowsAtCompileTime == 4 || + DerivedC::RowsAtCompileTime == Eigen::Dynamic, + "C must have 4 rows."); + // runtime assert that C has 4 rows + assert(C.rows() == 4 && "C must have 4 rows."); + using Scalar = typename DerivedP::Scalar; + // Evaluate cubic Bezier at parameter t + P = + (Scalar(1) - t) * (Scalar(1) - t) * (Scalar(1) - t) * C.row(0) + + Scalar(3) * (Scalar(1) - t) * (Scalar(1) - t) * t * C.row(1) + + Scalar(3) * (Scalar(1) - t) * t * t * C.row(2) + + t * t * t * C.row(3); +} + + + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template void igl::cubic, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic const, -1, 1, true>, Eigen::Matrix>(Eigen::MatrixBase const, -1, 1, true>> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic, Eigen::Block, 1, -1, false>, Eigen::internal::AllRange<-1>>, Eigen::Matrix>(Eigen::MatrixBase, Eigen::Block, 1, -1, false>, Eigen::internal::AllRange<-1>>> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cubic.h b/include/igl/cubic.h new file mode 100644 index 000000000..91a0993e5 --- /dev/null +++ b/include/igl/cubic.h @@ -0,0 +1,29 @@ +#ifndef IGL_CUBIC_H +#define IGL_CUBIC_H + +#include "igl_inline.h" +#include + +namespace igl +{ + /// Evaluate a cubic Bézier curve defined by control points C at parameter t. + /// + /// @param[in] C 4 by dimensions matrix of control points for a cubic + /// Bézier curve + /// @param[in] t Parameter at which to evaluate curve + /// @param[out] P 1 by dimensions point on curve C(t) + template + < + typename DerivedC, + typename DerivedP> + IGL_INLINE void cubic( + const Eigen::MatrixBase& C, + const typename DerivedP::Scalar & t, + Eigen::PlainObjectBase& P); +} + +#ifndef IGL_STATIC_LIBRARY +# include "cubic.cpp" +#endif + +#endif diff --git a/include/igl/cubic_is_flat.cpp b/include/igl/cubic_is_flat.cpp new file mode 100644 index 000000000..1a5cf3ce2 --- /dev/null +++ b/include/igl/cubic_is_flat.cpp @@ -0,0 +1,23 @@ +#include "cubic_is_flat.h" + +template +IGL_INLINE bool igl::cubic_is_flat( + const Eigen::MatrixBase& C, + const typename DerivedC::Scalar squared_distance_bound) +{ + using Scalar = typename DerivedC::Scalar; + const auto u = (Scalar(3) * C.row(1) - Scalar(2) * C.row(0) - C.row(3)).array().square().eval(); + const auto v = (Scalar(3) * C.row(2) - Scalar(2) * C.row(3) - C.row(0)).array().square().eval(); + const auto max_uv = u.cwiseMax(v).eval(); + const Scalar tao = max_uv.sum(); + const Scalar tolerance = Scalar(16) * squared_distance_bound; + return tao <= tolerance; +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template bool igl::cubic_is_flat>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar); +// generated by autoexplicit.sh +template bool igl::cubic_is_flat>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar); +#endif diff --git a/include/igl/cubic_is_flat.h b/include/igl/cubic_is_flat.h new file mode 100644 index 000000000..2a84221f6 --- /dev/null +++ b/include/igl/cubic_is_flat.h @@ -0,0 +1,29 @@ +#ifndef IGL_CUBIC_IS_FLAT_H +#define IGL_CUBIC_IS_FLAT_H +#include "igl_inline.h" +#include + +namespace igl +{ + /// "Piecewise Linear Approximation of Bézier Curves" [Fischer 2000]. The test + /// computes a scalar and checks if it is less than 16 times + /// `squared_distance_bound`. If so, the curve's maximum squared distance to + /// the chord from first to last point will be less than + /// `squared_distance_bound`. The test is well behaved for degenerate input + /// (e.g., all rows in C are the same point). + /// + /// @param[in] C 4 by dimensions matrix of control points for a cubic + /// Bézier curve + /// @param[in] squared_distance_bound Squared distance tolerance + /// @return True if the curve is flat within the given tolerance + template + IGL_INLINE bool cubic_is_flat( + const Eigen::MatrixBase& C, + const typename DerivedC::Scalar squared_distance_bound); +} + +#ifndef IGL_STATIC_LIBRARY +# include "cubic_is_flat.cpp" +#endif + +#endif diff --git a/include/igl/cubic_monomial_bases.cpp b/include/igl/cubic_monomial_bases.cpp new file mode 100644 index 000000000..86754eaf5 --- /dev/null +++ b/include/igl/cubic_monomial_bases.cpp @@ -0,0 +1,57 @@ +#include "cubic_monomial_bases.h" + +template < + typename DerivedC, + typename DerivedM, + typename DerivedD, + typename DerivedB> +IGL_INLINE void igl::cubic_monomial_bases( + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& M, + Eigen::PlainObjectBase& D, + Eigen::PlainObjectBase& B) +{ + // static assert that C has 4 rows or Dynamic + static_assert( + DerivedC::RowsAtCompileTime == 4 || + DerivedC::RowsAtCompileTime == Eigen::Dynamic, + "C must have 4 rows."); + // runtime assert that C has 4 rows + assert(C.rows() == 4 && "C must have 4 rows."); + const auto C0 = C.row(0); + const auto C1 = C.row(1); + const auto C2 = C.row(2); + const auto C3 = C.row(3); + // monomial coefficients for C(t) + M.resize(4,C.cols()); + M << + C.row(0).eval(), + 3 * (C1 - C0), + 3 * (C0 - 2 * C1 + C2), + C3 - C0 + 3 * (C1 - C2); + // Monomial coefficients for dC/dt + D.resize(3,C.cols()); + D << M.row(1), + 2 * M.row(2), + 3 * M.row(3); + // Should work for row and column vectors + B.resize(6); + B.setConstant(0); + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 4; ++j) + { + B(i + j) += D.row(i).dot(M.row(j)); + } + } +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template void igl::cubic_monomial_bases, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic_monomial_bases, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic_monomial_bases, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cubic_monomial_bases.h b/include/igl/cubic_monomial_bases.h new file mode 100644 index 000000000..e144cbe5a --- /dev/null +++ b/include/igl/cubic_monomial_bases.h @@ -0,0 +1,34 @@ +#ifndef IGL_CUBIC_MONOMIAL_BASES_H +#define IGL_CUBIC_MONOMIAL_BASES_H + +#include "igl_inline.h" +#include + +namespace igl +{ + /// Compute monomial basis representations for cubic Bézier curves. + /// + /// @param[in] C 4 by dimensions matrix of control points for a cubic + /// Bézier curve + /// @param[out] M 4 by dimensions matrix of monomial coefficients for C(t) + /// @param[out] D 3 by dimensions matrix of monomial coefficients for dC/dt + /// @param[out] B 6-vector of inner products of those basis functions for C(t) + /// + /// \see igl::cycodebase::point_cubic_squared_distance + template < + typename DerivedC, + typename DerivedM, + typename DerivedD, + typename DerivedB> + IGL_INLINE void cubic_monomial_bases( + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& M, + Eigen::PlainObjectBase& D, + Eigen::PlainObjectBase& B); +} + +#ifndef IGL_STATIC_LIBRARY +# include "cubic_monomial_bases.cpp" +#endif + +#endif diff --git a/include/igl/cubic_split.cpp b/include/igl/cubic_split.cpp new file mode 100644 index 000000000..4653f2021 --- /dev/null +++ b/include/igl/cubic_split.cpp @@ -0,0 +1,65 @@ +#include "cubic_split.h" + +template < + typename DerivedC, + typename DerivedK> +IGL_INLINE void igl::cubic_split( + const Eigen::MatrixBase& C, + const typename DerivedC::Scalar & t, + Eigen::PlainObjectBase& C01, + Eigen::PlainObjectBase& C012, + Eigen::PlainObjectBase& C0123, + Eigen::PlainObjectBase& C123, + Eigen::PlainObjectBase& C23) +{ + const auto C0 = C.row(0); + const auto C1 = C.row(1); + const auto C2 = C.row(2); + const auto C3 = C.row(3); + + C01 = (C1 - C0) * t + C0; + const auto C12 = ((C2 - C1) * t + C1).eval(); + C23 = (C3 - C2) * t + C2; + + C012 = (C12 - C01) * t + C01; + C123 = (C23 - C12) * t + C12; + + C0123 = (C123 - C012) * t + C012; +} + +template < + typename DerivedC, + typename DerivedK> +IGL_INLINE void igl::cubic_split( + const Eigen::MatrixBase& C, + const typename DerivedC::Scalar & t, + Eigen::PlainObjectBase& C1, + Eigen::PlainObjectBase& C2) +{ + using Scalar = typename DerivedC::Scalar; + typedef Eigen::Matrix RowVectorS; + RowVectorS C01,C012,C0123,C123,C23; + igl::cubic_split(C,t,C01,C012,C0123,C123,C23); + C1.resize(4,C.cols()); + C1 << C.row(0), + C01, + C012, + C0123; + C2.resize(4,C.cols()); + C2 << C0123, + C123, + C23, + C.row(3); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template void igl::cubic_split, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic_split, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cubic_split, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +template void igl::cubic_split, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +template void igl::cubic_split, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cubic_split.h b/include/igl/cubic_split.h new file mode 100644 index 000000000..205403267 --- /dev/null +++ b/include/igl/cubic_split.h @@ -0,0 +1,52 @@ +#ifndef IGL_CUBIC_SPLIT_H +#define IGL_CUBIC_SPLIT_H +#include "igl_inline.h" +#include + +namespace igl +{ + /// Splits a cubic Bézier curve defined by control points C at parameter t + /// into two cubic Bézier curves defined by control points C1 and C2. + /// C1 = [C0,C01,C012,C0123], C2 = [C0123,C123,C23,C3] + /// + /// @param[in] C 4 by dimensions matrix of control points for a cubic + /// Bézier curve + /// @param[in] t Parameter at which to split curve + /// @param[out] C01 1 by dim new control point + /// @param[out] C012 1 by dim new control point + /// @param[out] C0123 1 by dim new control point + /// @param[out] C123 1 by dim new control point + /// @param[out] C23 1 by dim new control point + /// + template < + typename DerivedC, + typename DerivedK> + IGL_INLINE void cubic_split( + const Eigen::MatrixBase& C, + const typename DerivedC::Scalar & t, + Eigen::PlainObjectBase& C01, + Eigen::PlainObjectBase& C012, + Eigen::PlainObjectBase& C0123, + Eigen::PlainObjectBase& C123, + Eigen::PlainObjectBase& C23); + /// + /// @param[out] C1 4 by dimensions matrix of control points for first cubic + /// Bézier curve from C(0) to C(t) + /// @param[out] C2 4 by dimensions matrix of control points for second cubic + /// Bézier curve from C(t) to C(1) + template < + typename DerivedC, + typename DerivedK> + IGL_INLINE void cubic_split( + const Eigen::MatrixBase& C, + const typename DerivedC::Scalar & t, + Eigen::PlainObjectBase& C1, + Eigen::PlainObjectBase& C2); +} + +#ifndef IGL_STATIC_LIBRARY +# include "cubic_split.cpp" +#endif + +#endif + diff --git a/include/igl/cycodebase/box_cubic.cpp b/include/igl/cycodebase/box_cubic.cpp new file mode 100644 index 000000000..752cfe888 --- /dev/null +++ b/include/igl/cycodebase/box_cubic.cpp @@ -0,0 +1,82 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2026 Alec Jacobson +// +// 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 http://mozilla.org/MPL/2.0/. +#include "box_cubic.h" +#include "../cubic.h" +#include "../parallel_for.h" +#include + +template < + typename DerivedC, + typename DerivedB> +IGL_INLINE void igl::cycodebase::box_cubic( + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& B1, + Eigen::PlainObjectBase& B2) +{ + using Scalar = typename DerivedC::Scalar; + typedef Eigen::Matrix VectorS3; + typedef Eigen::Matrix VectorSC; + + // Using the control points is a simple (but not tight) bound + B1 = C({0,3},Eigen::all).colwise().minCoeff(); + B2 = C({0,3},Eigen::all).colwise().maxCoeff(); + + // Better to find each of the t values where dC/dt = 0 and evaluate C there + for(int d = 0;d(r.data(), coef.data(),0.0,1.0); + for(int i = 0;i Cdt; + igl::cubic(C.col(d),t,Cdt); + B1(d) = std::min(B1(d),Cdt(0)); + B2(d) = std::max(B2(d),Cdt(0)); + } + } +} + +template < + typename DerivedP, + typename DerivedC, + typename DerivedB> +IGL_INLINE void igl::cycodebase::box_cubic( + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& B1, + Eigen::PlainObjectBase& B2) +{ + B1.resize(C.rows(),P.cols()); + B2.resize(C.rows(),P.cols()); + typedef Eigen::Matrix RowVectorP; + + igl::parallel_for(C.rows(),[&](const int c) + { + RowVectorP B1_c, B2_c; + // Eval copies, but is it really good to make a template for the non copy? + box_cubic(P(C.row(c),Eigen::all).eval(),B1_c,B2_c); + B1.row(c) = B1_c; + B2.row(c) = B2_c; + },1000); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template void igl::cycodebase::box_cubic, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +template void igl::cycodebase::box_cubic, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +template void igl::cycodebase::box_cubic, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cycodebase/box_cubic.h b/include/igl/cycodebase/box_cubic.h new file mode 100644 index 000000000..c96e07aee --- /dev/null +++ b/include/igl/cycodebase/box_cubic.h @@ -0,0 +1,46 @@ +#ifndef IGL_CYCODEBASE_BOX_CUBIC_H +#define IGL_CYCODEBASE_BOX_CUBIC_H + +#include "../igl_inline.h" +#include + +namespace igl { + namespace cycodebase { + /// Compute the min/max box corners tightly containing a given cubic bezier + /// curve. + /// + /// @param[in] C 4 by dim matrix of control points defining the cubic bezier + /// curve + /// @param[out] B1 1 by dim min corner of the bounding box + /// @param[out] B2 1 by dim max corner of the bounding box + /// + /// \see igl::box_simplices + template < + typename DerivedC, + typename DerivedB> + IGL_INLINE void box_cubic( + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& B1, + Eigen::PlainObjectBase& B2); + /// \brief overload + /// + /// @param[in] P #P by dim matrix of control point locations + /// @param[in] C #C by 4 matrix of indices into P defining the cubics + /// @param[out] B1 #C by dim matrix of min corners of the bounding boxes + /// @param[out] B2 #C by dim matrix of max corners of the bounding boxes + template < + typename DerivedP, + typename DerivedC, + typename DerivedB> + IGL_INLINE void box_cubic( + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& B1, + Eigen::PlainObjectBase& B2); + } +} + +#ifndef IGL_STATIC_LIBRARY + #include "box_cubic.cpp" +#endif +#endif diff --git a/include/igl/cycodebase/point_cubic_squared_distance.cpp b/include/igl/cycodebase/point_cubic_squared_distance.cpp new file mode 100644 index 000000000..4f6873802 --- /dev/null +++ b/include/igl/cycodebase/point_cubic_squared_distance.cpp @@ -0,0 +1,197 @@ +#include "point_cubic_squared_distance.h" +#include "../parallel_for.h" +#include "../cubic.h" +#include "../cubic_monomial_bases.h" +#include +#include + +template < + typename DerivedQ, + typename DerivedC, + typename DerivedsqrD, + typename DerivedS, + typename DerivedK> +void igl::cycodebase::point_cubic_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K) +{ + using Scalar = typename DerivedQ::Scalar; + constexpr int ColsAtCompileTime = DerivedQ::ColsAtCompileTime; + const int dim = C.cols(); + assert(C.rows() == 4 && "C should be 4 x dim."); + assert(Q.cols() == dim && "Q and C should have the same number of columns."); + // min_t ‖ Q - C(t) ‖^2 + // + // Necessary condition for a minimum: + // + // d/dt ‖ Q - C(t) ‖^2 = 0 + // => d/dt (Q - C(t)) · (Q - C(t)) = 0 + // => -2 (dC/dt) · (Q - C(t)) = 0 + // => (dC/dt) · (C(t) - Q) = f(t) = 0 + // + using RowVectorSD = Eigen::RowVector; + using MatrixS4D = Eigen::Matrix; + using MatrixS3D = Eigen::Matrix; + MatrixS4D M; + MatrixS3D D; + Eigen::RowVector B; + cubic_monomial_bases(C, M, D, B); + return point_cubic_squared_distance( Q, C, D, B, sqrD, S, K); +} + +template < + typename DerivedQ, + typename DerivedC, + typename DerivedD, + typename DerivedB, + typename DerivedsqrD, + typename DerivedS, + typename DerivedK> +void igl::cycodebase::point_cubic_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& D, + const Eigen::MatrixBase& B, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K) +{ + // static assert that C has 4 rows or Dynamic + static_assert( + DerivedC::RowsAtCompileTime == 4 || + DerivedC::RowsAtCompileTime == Eigen::Dynamic, + "C must have 4 rows."); + // runtime assert that C has 4 rows + assert(C.rows() == 4 && "C must have 4 rows."); + using Scalar = typename DerivedQ::Scalar; + constexpr int ColsAtCompileTime = DerivedQ::ColsAtCompileTime; + const int dim = C.cols(); + using RowVectorSD = Eigen::RowVector; + sqrD.setConstant(Q.rows(),1,std::numeric_limits::infinity()); + S.resize(Q.rows()); + K.resize(Q.rows(),dim); + const int n = Q.rows(); + for(int i = 0;i +void igl::cycodebase::point_cubic_squared_distance( + const Eigen::MatrixBase& q, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& D, + const Eigen::MatrixBase& B, + typename Derivedq::Scalar& sqrD, + typename Derivedq::Scalar& s, + Eigen::PlainObjectBase& k) +{ + using cy::PolynomialRoots; + using Scalar = typename Derivedq::Scalar; + const int dim = C.cols(); + // static assert that C has 4 rows or Dynamic + static_assert( + DerivedC::RowsAtCompileTime == 4 || + DerivedC::RowsAtCompileTime == Eigen::Dynamic, + "C must have 4 rows."); + assert(C.rows() == 4 && "C should be 4 x dim."); + assert(q.cols() == dim && "q and C should have the same number of columns."); + typedef Eigen::Matrix RowVectorSD; + // Fill in coefficients using precomputed data and Q + // => (dC/dt) · (C(t) - Q) = f(t) = 0 + // + // C(t) = (1-t)^3 C0 + 3(1-t)^2 t C1 + 3(1-t) t^2 C2 + t^3 C3 + // C(t) = C0 + t⋅3(C1 - C0) + t^2⋅3(C0 - 2C1 + C2) + t^3⋅(C3 - C0 + 3(C1 - C2)) + // dC/dt = 3(C1 - C0) + t⋅6(C0 - 2C1 + C2) + t^2⋅3(C3 - C0 + 3(C1 - C2)) + Eigen::RowVector coef = B; + for (int j = 0; j < 3; ++j) + { + coef(j) -= D.row(j).dot(q); + } + + sqrD = std::numeric_limits::infinity(); + // f is a quintic polynomial: + constexpr int N = 5; + Scalar r[N]; + int nr = PolynomialRoots(r, coef.data(),Scalar(0),Scalar(1)); + for(int j = 0;j +void igl::cycodebase::point_cubic_squared_distance( + const Eigen::MatrixBase& q, + const Eigen::MatrixBase& C, + typename Derivedq::Scalar& sqrD, + typename Derivedq::Scalar& s, + Eigen::PlainObjectBase& k) +{ + using Scalar = typename Derivedq::Scalar; + constexpr int ColsAtCompileTime = Derivedq::ColsAtCompileTime; + static_assert( + DerivedC::RowsAtCompileTime == 4 || + DerivedC::RowsAtCompileTime == Eigen::Dynamic, + "C must have 4 rows."); + static_assert( + int(Derivedq::ColsAtCompileTime) == int(DerivedC::ColsAtCompileTime), + "q and C must have the same number of columns."); + const int dim = C.cols(); + assert(C.rows() == 4 && "C should be 4 x dim."); + assert(q.cols() == dim && "q and C should have the same number of columns."); + using RowVectorSD = Eigen::RowVector; + using MatrixS4D = Eigen::Matrix; + using MatrixS3D = Eigen::Matrix; + MatrixS4D M; + MatrixS3D D; + Eigen::RowVector B; + cubic_monomial_bases(C, M, D, B); + return point_cubic_squared_distance( q, C, D, B, sqrD, s, k); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +template void igl::cycodebase::point_cubic_squared_distance, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&, Eigen::Matrix::Scalar&, Eigen::PlainObjectBase>&); +template void igl::cycodebase::point_cubic_squared_distance, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&, Eigen::Matrix::Scalar&, Eigen::PlainObjectBase>&); +template void igl::cycodebase::point_cubic_squared_distance, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cycodebase/point_cubic_squared_distance.h b/include/igl/cycodebase/point_cubic_squared_distance.h new file mode 100644 index 000000000..ffe7e4842 --- /dev/null +++ b/include/igl/cycodebase/point_cubic_squared_distance.h @@ -0,0 +1,94 @@ +#ifndef IGL_CYCODEBASE_POINT_CUBIC_SQUARED_DISTANCE_H +#define IGL_CYCODEBASE_POINT_CUBIC_SQUARED_DISTANCE_H + +#include "../igl_inline.h" +#include + +namespace igl { + namespace cycodebase { + /// Compute the squared distance from a set of points to a cubic Bezier + /// curve. + /// + /// @param[in] Q #Q by dim matrix of query points + /// @param[in] C 4 by dim matrix of control points defining the cubic + /// bezier curve + /// @param[out] sqrD #Q vector of squared distances from each query point + /// to the cubic bezier curve + /// @param[out] S #Q vector of parametric locations on the cubic bezier + /// curve of the closest point to each query point + /// @param[out] K #Q by dim matrix of closest point locations + /// + /// \see igl::cubic_monomial_bases, igl::cycodebase::point_spline_squared_distance + template < + typename DerivedQ, + typename DerivedC, + typename DerivedsqrD, + typename DerivedS, + typename DerivedK> + void point_cubic_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K); + /// \brief overload + /// + /// @param[in] D 3 by dim matrix of monomial coefficients for dC/dt + /// @param[in] B 6-vector of inner products of monomial coefficients of C + /// and dCdt + template < + typename DerivedQ, + typename DerivedC, + typename DerivedD, + typename DerivedB, + typename DerivedsqrD, + typename DerivedS, + typename DerivedK> + void point_cubic_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& D, + const Eigen::MatrixBase& B, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K); + /// \brief single point overload + /// + /// @param[in] q dim vector of a query point + template < + typename Derivedq, + typename DerivedC, + typename DerivedD, + typename DerivedB, + typename Derivedk + > + void point_cubic_squared_distance( + const Eigen::MatrixBase& q, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& D, + const Eigen::MatrixBase& B, + typename Derivedq::Scalar& sqrD, + typename Derivedq::Scalar& s, + Eigen::PlainObjectBase& k); + /// \brief single point overload + template < + typename Derivedq, + typename DerivedC, + typename Derivedk + > + void point_cubic_squared_distance( + const Eigen::MatrixBase& q, + const Eigen::MatrixBase& C, + typename Derivedq::Scalar& sqrD, + typename Derivedq::Scalar& s, + Eigen::PlainObjectBase& k); + } +} + +#ifndef IGL_STATIC_LIBRARY + #include "point_cubic_squared_distance.cpp" +#endif + +#endif + + diff --git a/include/igl/cycodebase/point_spline_squared_distance.cpp b/include/igl/cycodebase/point_spline_squared_distance.cpp new file mode 100644 index 000000000..3a5fc65e8 --- /dev/null +++ b/include/igl/cycodebase/point_spline_squared_distance.cpp @@ -0,0 +1,117 @@ +#include "point_spline_squared_distance.h" +#include "point_cubic_squared_distance.h" +#include "spline_eytzinger_aabb.h" +#include "../eytzinger_aabb_sdf.h" +#include "../parallel_for.h" +#include "../cubic_monomial_bases.h" + +template < + typename DerivedQ, + typename DerivedP, + typename DerivedC, + typename DerivedsqrD, + typename DerivedI, + typename DerivedS, + typename DerivedK> +void igl::cycodebase::point_spline_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& I, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K) +{ + static_assert( + DerivedQ::ColsAtCompileTime == DerivedP::ColsAtCompileTime, + "Q and P must have the same number of columns."); + using Scalar = typename DerivedQ::Scalar; + Eigen::Matrix + B1,B2; + Eigen::VectorXi leaf; + spline_eytzinger_aabb(P, C, B1, B2,leaf); + return point_spline_squared_distance( + Q,P,C,B1,B2,leaf,sqrD,I,S,K); +} + +template < + typename DerivedQ, + typename DerivedP, + typename DerivedC, + typename DerivedB, + typename Derivedleaf, + typename DerivedsqrD, + typename DerivedI, + typename DerivedS, + typename DerivedK> +void igl::cycodebase::point_spline_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& B1, + const Eigen::MatrixBase& B2, + const Eigen::MatrixBase& leaf, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& I, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K) +{ + using Scalar = typename DerivedQ::Scalar; + using RowVectorS = Eigen::Matrix; + using MatrixS4D = Eigen::Matrix; + using MatrixS3D = Eigen::Matrix; + + // Cache these. + std::vector C_vec(C.rows()); + std::vector D_vec(C.rows()); + std::vector> B_vec(C.rows()); + { + MatrixS4D M_unused; + for(int j = 0; j < C.rows(); j++) + { + C_vec[j] = P(C.row(j), Eigen::all); + cubic_monomial_bases( C_vec[j], M_unused, D_vec[j], B_vec[j]); + } + } + + sqrD.setConstant(Q.rows(),std::numeric_limits::infinity()); + S.resize(Q.rows()); + I.resize(Q.rows()); + K.resize(Q.rows(),Q.cols()); + + //for(int i = 0; i < Q.rows(); i++) + igl::parallel_for(Q.rows(), [&](const int i) + { + const Eigen::Matrix q = Q.row(i); + const std::function primitive_i = [&](const int j) + { + Scalar sqrD_ij,s_ij; + RowVectorS k_ij; + point_cubic_squared_distance( q, + C_vec[j], D_vec[j], B_vec[j], + sqrD_ij, s_ij, k_ij); + if(sqrD_ij < sqrD(i)) + { + sqrD(i) = sqrD_ij; + S(i) = s_ij; + K.row(i) = k_ij; + I(i) = j; + } + return sqrD_ij; + }; + // This is expected (un)signed **unsquared** _distance. So let it compute + // that in there. We will square it and acculumate our results in the lambda + // anyway. So dist_i is ignored. + Scalar dist_i; + igl::eytzinger_aabb_sdf(q,primitive_i,B1,B2,leaf,dist_i); + },1000); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template void igl::cycodebase::point_spline_squared_distance, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cycodebase::point_spline_squared_distance, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +template void igl::cycodebase::point_spline_squared_distance, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cycodebase/point_spline_squared_distance.h b/include/igl/cycodebase/point_spline_squared_distance.h new file mode 100644 index 000000000..218cfc2c0 --- /dev/null +++ b/include/igl/cycodebase/point_spline_squared_distance.h @@ -0,0 +1,80 @@ +#ifndef IGL_CYCODEBASE_POINT_SPLINE_SQUARED_DISTANCE_H +#define IGL_CYCODEBASE_POINT_SPLINE_SQUARED_DISTANCE_H + +#include "../igl_inline.h" +#include + +namespace igl { + namespace cycodebase { + /// Compute the squared distance from a set of points to a spline defined by + /// cubic Bezier curves. + /// + /// @param[in] Q #Q by dim matrix of query points + /// @param[in] P #P by dim matrix of spline control points + /// @param[in] C #C by 4 matrix of indices into P defining the cubic bezier + /// curves making up the spline + /// @param[out] sqrD #Q vector of squared distances from each query point + /// to the spline + /// @param[out] I #Q vector of indices of the closest cubic bezier curve in + /// the spline to each query point + /// @param[out] S #Q vector of parametric locations on the closest cubic + /// bezier curve of the closest point to each query point + /// @param[out] K #Q by dim matrix of closest point locations + /// + /// \see igl::cycodebase::point_cubic_squared_distance, + /// igl::cycodebase::spline_eytzinger_aabb + template < + typename DerivedQ, + typename DerivedP, + typename DerivedC, + typename DerivedsqrD, + typename DerivedI, + typename DerivedS, + typename DerivedK> + void point_spline_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& I, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K); + /// \brief overload with AABB acceleration structure + /// + /// @param[in] B1 #B by dim matrix of AABB min box corners + /// @param[in] B2 #B by dim matrix of AABB max box corners + /// @param[in] leaf #B by 1 matrix of AABB leaf node indices/flags + /// + /// \see igl::eytzinger_aabb + template < + typename DerivedQ, + typename DerivedP, + typename DerivedC, + typename DerivedB, + typename Derivedleaf, + typename DerivedsqrD, + typename DerivedI, + typename DerivedS, + typename DerivedK> + void point_spline_squared_distance( + const Eigen::MatrixBase& Q, + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& B1, + const Eigen::MatrixBase& B2, + const Eigen::MatrixBase& leaf, + Eigen::PlainObjectBase& sqrD, + Eigen::PlainObjectBase& I, + Eigen::PlainObjectBase& S, + Eigen::PlainObjectBase& K); + } +} + +#ifndef IGL_STATIC_LIBRARY + #include "point_spline_squared_distance.cpp" +#endif + +#endif + + + diff --git a/include/igl/cycodebase/roots.cpp b/include/igl/cycodebase/roots.cpp new file mode 100644 index 000000000..7cc58b993 --- /dev/null +++ b/include/igl/cycodebase/roots.cpp @@ -0,0 +1,83 @@ +#include "roots.h" +#include + +namespace +{ + template + int roots_helper( + const int deg, + const Eigen::MatrixBase& coef, + const typename Derivedcoef::Scalar xmin, + const typename Derivedcoef::Scalar xmax, + Eigen::PlainObjectBase& R) + { + using cy::PolynomialRoots; + if constexpr (MAX_DEG == 0) + { + // Should be impossible to reach here from roots() + return 0; + } + else + { + if (deg == MAX_DEG) + { + using Scalar = typename Derivedcoef::Scalar; + Scalar r[MAX_DEG]; + int nr = PolynomialRoots(r, coef.derived().data(), xmin, xmax); + // will resize if needed. Safe on compile-time sized matrices + R.setConstant(MAX_DEG, std::numeric_limits::quiet_NaN()); + for (int i = 0; i < nr; ++i) + { + R[i] = r[i]; + } + return nr; + } + return roots_helper(deg, coef, xmin, xmax, R); + } + } + +} + +template < + typename Derivedcoef, + typename DerivedR> +IGL_INLINE int igl::cycodebase::roots( + const Eigen::MatrixBase& coef, + const typename Derivedcoef::Scalar xmin, + const typename Derivedcoef::Scalar xmax, + Eigen::PlainObjectBase& R) +{ + using Scalar = typename Derivedcoef::Scalar; + // Static assert that DerivedR::Scalar is same as Derivedcoef::Scalar + static_assert( + std::is_same::value, + "DerivedR::Scalar must be the same as Derivedcoef::Scalar"); + // Check that Derivedcoef and DerivedR are vectors + static_assert( + (Derivedcoef::ColsAtCompileTime == 1 || + Derivedcoef::RowsAtCompileTime == 1) && + (DerivedR::ColsAtCompileTime == 1 || + DerivedR::RowsAtCompileTime == 1), + "Derivedcoef and DerivedR must be vectors"); + + constexpr int MAX_DEG = 16; + const int deg = coef.size() - 1; + + if (deg < 1 || deg > MAX_DEG) + { + throw std::runtime_error("Polynomial degree out of range"); + } + + return ::roots_helper( deg, coef, xmin, xmax, R); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template specialization +template IGL_INLINE int igl::cycodebase::roots( + const Eigen::MatrixBase& coef, + const double xmin, + const double xmax, + Eigen::PlainObjectBase& R); +template int igl::cycodebase::roots, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar, Eigen::Matrix::Scalar, Eigen::PlainObjectBase>&); +template int igl::cycodebase::roots, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::Matrix::Scalar, Eigen::Matrix::Scalar, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cycodebase/roots.h b/include/igl/cycodebase/roots.h new file mode 100644 index 000000000..72344a279 --- /dev/null +++ b/include/igl/cycodebase/roots.h @@ -0,0 +1,33 @@ +#ifndef IGL_CYCODEBASE_ROOTS_H +#define IGL_CYCODEBASE_ROOTS_H + +#include "../igl_inline.h" +#include + +namespace igl { + namespace cycodebase { + /// Compute the real roots of a polynomial with given coefficients within a + /// given interval [xmin, xmax]. + /// + /// @param[in] coef Vector of polynomial coefficients in increasing order + /// of degree (i.e., coef[0] + coef[1]*x + coef[2]*x^2 + ... ) + /// @param[in] xmin Minimum x value of the interval + /// @param[in] xmax Maximum x value of the interval + /// @param[out] R Vector of real roots within [xmin, xmax] + /// @return Number of real roots found within [xmin, xmax] + template < + typename Derivedcoef, + typename DerivedR> + IGL_INLINE int roots( + const Eigen::MatrixBase& coef, + const typename Derivedcoef::Scalar xmin, + const typename Derivedcoef::Scalar xmax, + Eigen::PlainObjectBase& R); + } +} + +#ifndef IGL_STATIC_LIBRARY + #include "roots.cpp" +#endif + +#endif diff --git a/include/igl/cycodebase/spline_eytzinger_aabb.cpp b/include/igl/cycodebase/spline_eytzinger_aabb.cpp new file mode 100644 index 000000000..d4ed501ee --- /dev/null +++ b/include/igl/cycodebase/spline_eytzinger_aabb.cpp @@ -0,0 +1,30 @@ +#include "spline_eytzinger_aabb.h" +#include "box_cubic.h" +#include "../eytzinger_aabb.h" + +template < + typename DerivedP, + typename DerivedC, + typename DerivedB, + typename Derivedleaf> +void igl::cycodebase::spline_eytzinger_aabb( + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& B1, + Eigen::PlainObjectBase& B2, + Eigen::PlainObjectBase& leaf) +{ + using Scalar = typename DerivedP::Scalar; + Eigen::Matrix + PB1,PB2; + box_cubic(P,C,PB1,PB2); + eytzinger_aabb( PB1, PB2, B1, B2,leaf); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +// generated by autoexplicit.sh +template void igl::cycodebase::spline_eytzinger_aabb, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::cycodebase::spline_eytzinger_aabb, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/cycodebase/spline_eytzinger_aabb.h b/include/igl/cycodebase/spline_eytzinger_aabb.h new file mode 100644 index 000000000..a0c319f89 --- /dev/null +++ b/include/igl/cycodebase/spline_eytzinger_aabb.h @@ -0,0 +1,43 @@ +#ifndef IGL_CYCODEBASE_SPLINE_EYTZINGER_AABB_H +#define IGL_CYCODEBASE_SPLINE_EYTZINGER_AABB_H + +#include "../igl_inline.h" +#include + +namespace igl { + namespace cycodebase { + /// Compute an AABB acceleration structure for a spline of cubic Bezier + /// curves in Eytzinger layout. + /// + /// @param[in] P #P by dim matrix of spline control points + /// @param[in] C #C by 4 matrix of indices into P defining the cubic Bézier + /// curves making up the spline + /// @param[out] B1 #B by dim matrix of AABB min box corners + /// @param[out] B2 #B by dim matrix of AABB max box corners + /// @param[out] leaf #B by 1 matrix of AABB leaf node indices/flags + /// + /// \see igl::cycodebase::box_cubic, igl::eytzinger_aabb + /// + template < + typename DerivedP, + typename DerivedC, + typename DerivedB, + typename Derivedleaf> + void spline_eytzinger_aabb( + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + Eigen::PlainObjectBase& B1, + Eigen::PlainObjectBase& B2, + Eigen::PlainObjectBase& leaf); + } +} + +#ifndef IGL_STATIC_LIBRARY + #include "spline_eytzinger_aabb.cpp" +#endif + +#endif + + + + diff --git a/include/igl/eytzinger_aabb.cpp b/include/igl/eytzinger_aabb.cpp index c27af1193..7db6f9b0c 100644 --- a/include/igl/eytzinger_aabb.cpp +++ b/include/igl/eytzinger_aabb.cpp @@ -85,6 +85,12 @@ IGL_INLINE void igl::eytzinger_aabb( #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation // generated by autoexplicit.sh +template void igl::eytzinger_aabb, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::eytzinger_aabb, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh +template void igl::eytzinger_aabb, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); +// generated by autoexplicit.sh template void igl::eytzinger_aabb, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); // generated by autoexplicit.sh template void igl::eytzinger_aabb, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); diff --git a/include/igl/eytzinger_aabb_sdf.cpp b/include/igl/eytzinger_aabb_sdf.cpp index 135a553b6..87630d862 100644 --- a/include/igl/eytzinger_aabb_sdf.cpp +++ b/include/igl/eytzinger_aabb_sdf.cpp @@ -2,6 +2,7 @@ #include "parallel_for.h" template < + bool Squared, typename Derivedp, typename Func, typename DerivedB, @@ -25,9 +26,16 @@ IGL_INLINE void igl::eytzinger_aabb_sdf( // relative position const RowVectorS r = p - 0.5*(b1 + b2); const RowVectorS q = r.cwiseAbs() - b; - const Scalar t1 = q.cwiseMax(0.0).squaredNorm(); - const Scalar t2 = std::min(q.maxCoeff(),0.0); - return std::sqrt(t1 + t2); + if constexpr (Squared) + { + return q.cwiseMax(Scalar(0)).squaredNorm(); + }else + { + const Scalar t1 = q.cwiseMax(0.0).norm(); + const Scalar t2 = std::min(q.maxCoeff(),0.0); + const Scalar sdf = t1 + t2; + return sdf; + } }; // List of current active nodes @@ -63,6 +71,7 @@ IGL_INLINE void igl::eytzinger_aabb_sdf( } template < + bool Squared, typename DerivedP, typename Func, typename DerivedB, @@ -87,16 +96,18 @@ IGL_INLINE void igl::eytzinger_aabb_sdf( { return primitive(p,j); }; - igl::eytzinger_aabb_sdf(p,primitive_i,B1,B2,leaf,S(i)); + igl::eytzinger_aabb_sdf(p,primitive_i,B1,B2,leaf,S(i)); },1000); } #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation // generated by autoexplicit.sh -template void igl::eytzinger_aabb_sdf, std::function, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); +template void igl::eytzinger_aabb_sdf, std::function, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); // generated by autoexplicit.sh -template void igl::eytzinger_aabb_sdf, std::function const&, int)>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, int)> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&); -template void igl::eytzinger_aabb_sdf, std::function, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); -template void igl::eytzinger_aabb_sdf, std::function const&, int)>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, int)> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&); +template void igl::eytzinger_aabb_sdf, std::function, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); +// generated by autoexplicit.sh +template void igl::eytzinger_aabb_sdf, std::function, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); +template void igl::eytzinger_aabb_sdf, std::function const&, int)>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, int)> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&); +template void igl::eytzinger_aabb_sdf, std::function, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, std::function const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); #endif diff --git a/include/igl/eytzinger_aabb_sdf.h b/include/igl/eytzinger_aabb_sdf.h index ff395cd04..b5ebec446 100644 --- a/include/igl/eytzinger_aabb_sdf.h +++ b/include/igl/eytzinger_aabb_sdf.h @@ -21,7 +21,13 @@ namespace igl /// and returns the SDF to that primitive `primitive(i)` /// @param[out] f SDF value at query point x /// + /// eytzinger_aabb_sdf computes the signed distance function and + /// assumes that primitive(i) returns the signed distance to primitive i + /// + /// eytzinger_aabb_sdf computes the squared (unsigned) distance function + /// and assumes that primitive(i) returns the squared distance to primitive i template < + bool Squared, typename Derivedp, typename Func, typename DerivedB, @@ -45,6 +51,7 @@ namespace igl /// `primitive(p,i)` /// @param[out] S #P list of SDF values at query points P template < + bool Squared, typename DerivedP, typename Func, typename DerivedB, diff --git a/include/igl/eytzinger_aabb_winding_number.cpp b/include/igl/eytzinger_aabb_winding_number.cpp index 4ed29a0fd..247cee23e 100644 --- a/include/igl/eytzinger_aabb_winding_number.cpp +++ b/include/igl/eytzinger_aabb_winding_number.cpp @@ -1,4 +1,26 @@ #include "eytzinger_aabb_winding_number.h" +#include "PI.h" + +// signed_angle helper. Not sure why igl::signed_angle is doing all that +// normalization (divisions and square roots) when atan2 doesn't need it. +namespace +{ + template + inline typename MatV::Scalar signed_angle( + const MatV & V, + const VecP & p, + const int e0, + const int e1) + { + const Eigen::Matrix v0 = V.row(e0) - p; + const Eigen::Matrix v1 = V.row(e1) - p; + const typename MatV::Scalar angle = std::atan2( + v0(0)*v1(1) - v0(1)*v1(0), + v0(0)*v1(0) + v0(1)*v1(1)); + return angle; + } +} + template < typename Derivedp, typename DerivedV, @@ -18,16 +40,34 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number( const Eigen::MatrixBase & C, typename DerivedV::Scalar & wn) { - const auto signed_angle = [&V,&p]( const int e0, const int e1)->typename DerivedV::Scalar + using Scalar = typename DerivedV::Scalar; + const auto primitive = [&p,&V,&E](const int e)->typename DerivedV::Scalar { - const Eigen::Matrix v0 = V.row(e0) - p; - const Eigen::Matrix v1 = V.row(e1) - p; - const typename DerivedV::Scalar angle = std::atan2( - v0(0)*v1(1) - v0(1)*v1(0), - v0(0)*v1(0) + v0(1)*v1(1)); - return angle; + return ::signed_angle(V,p,E(e,0),E(e,1)) / Scalar(2.0 * igl::PI); }; + return igl::eytzinger_aabb_winding_number( + p,V,primitive,B1,B2,leaf,I,C,wn); +} +template < + typename Derivedp, + typename DerivedV, + typename DerivedB, + typename Derivedleaf, + typename DerivedI, + typename DerivedC> +IGL_INLINE void igl::eytzinger_aabb_winding_number( + const Eigen::MatrixBase & p, + const Eigen::MatrixBase & V, + const std::function & primitive, + const Eigen::MatrixBase & B1, + const Eigen::MatrixBase & B2, + const Eigen::MatrixBase & leaf, + const Eigen::MatrixBase & I, + const Eigen::MatrixBase & C, + typename DerivedV::Scalar & wn) +{ + using Scalar = typename DerivedV::Scalar; wn = 0; if(leaf.size() == 0) { wn = 0; return; } // I don't think stack or queue matters @@ -39,7 +79,8 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number( stack.pop_back(); if(leaf(r) >= 0) { - wn += signed_angle(E(leaf(r),0),E(leaf(r),1)); + // Default behavior + wn += primitive(leaf(r)); continue; } // otherwise is p outside the box B1.row(r), B2.row(r)? @@ -50,7 +91,7 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number( assert(C(r+1) - C(r) % 2 == 0); for(int i = C(r); i < C(r+1); i+=2) { - wn += signed_angle(I(i),I(i+1)); + wn += ::signed_angle(V,p,I(i),I(i+1)) / Scalar(2.0 * igl::PI); } continue; } @@ -64,6 +105,7 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number( #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation +// generated by autoexplicit.sh +template void igl::eytzinger_aabb_winding_number, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, std::function::Scalar (int)> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); template void igl::eytzinger_aabb_winding_number, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::Matrix::Scalar&); #endif - diff --git a/include/igl/eytzinger_aabb_winding_number.h b/include/igl/eytzinger_aabb_winding_number.h index 0ca1b71fa..7b9e82bf2 100644 --- a/include/igl/eytzinger_aabb_winding_number.h +++ b/include/igl/eytzinger_aabb_winding_number.h @@ -41,6 +41,23 @@ namespace igl const Eigen::MatrixBase & I, const Eigen::MatrixBase & C, typename DerivedV::Scalar & wn); + template < + typename Derivedp, + typename DerivedV, + typename DerivedB, + typename Derivedleaf, + typename DerivedI, + typename DerivedC> + IGL_INLINE void eytzinger_aabb_winding_number( + const Eigen::MatrixBase & p, + const Eigen::MatrixBase & V, + const std::function & primitive, + const Eigen::MatrixBase & B1, + const Eigen::MatrixBase & B2, + const Eigen::MatrixBase & leaf, + const Eigen::MatrixBase & I, + const Eigen::MatrixBase & C, + typename DerivedV::Scalar & wn); } #ifndef IGL_STATIC_LIBRARY diff --git a/include/igl/eytzinger_aabb_winding_number_tree.cpp b/include/igl/eytzinger_aabb_winding_number_tree.cpp index 6cda0ebaa..4b3231e50 100644 --- a/include/igl/eytzinger_aabb_winding_number_tree.cpp +++ b/include/igl/eytzinger_aabb_winding_number_tree.cpp @@ -114,5 +114,7 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number_tree( #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation +// generated by autoexplicit.sh +template void igl::eytzinger_aabb_winding_number_tree, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); template void igl::eytzinger_aabb_winding_number_tree, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&, Eigen::PlainObjectBase>&); #endif diff --git a/include/igl/predicates/cubic_winding_number.cpp b/include/igl/predicates/cubic_winding_number.cpp new file mode 100644 index 000000000..546594816 --- /dev/null +++ b/include/igl/predicates/cubic_winding_number.cpp @@ -0,0 +1,83 @@ +#include "cubic_winding_number.h" +#include "point_in_convex_hull.h" +#include "orient2d.h" +#include "../Orientation.h" +#include "../cubic_is_flat.h" +#include "../cubic_split.h" +#include "../PI.h" +#include "../matlab_format.h" +#include +#include + +namespace +{ + // helper where we assume that the input is not degenerate + template + inline typename DerivedC::Scalar cubic_winding_number_helper( + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& q) + { + using Scalar = typename DerivedC::Scalar; + + if(igl::predicates::point_in_convex_hull( + q, + C.row(0).eval(), + C.row(1).eval(), + C.row(2).eval(), + C.row(3).eval()) == igl::Orientation::NEGATIVE + ) + { + const auto v0 = (C.row(0) - q).eval(); + const auto v3 = (C.row(3) - q).eval(); + const auto w = std::atan2( + v0(0)*v3(1) - v0(1)*v3(0), + v0(0)*v3(0) + v0(1)*v3(1)) / + Scalar(2.0 * igl::PI); + return w; + } + else + { + Eigen::Matrix C1,C2; + igl::cubic_split(C,Scalar(0.5),C1,C2); + const auto w1 = cubic_winding_number_helper(C1,q); + const auto w2 = cubic_winding_number_helper(C2,q); + return w1 + w2; + } + } +} + +template +IGL_INLINE typename DerivedC::Scalar igl::predicates::cubic_winding_number( + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& q) +{ + // Degenerate cases. + using Scalar = typename DerivedC::Scalar; + // check if C are all the same point + if( (C.row(0)-C.row(1)).squaredNorm() == Scalar(0) && + (C.row(0)-C.row(2)).squaredNorm() == Scalar(0) && + (C.row(0)-C.row(3)).squaredNorm() == Scalar(0) ) + { + return Scalar(0); + } + // Check if C are all COLLINEAR + if( + (orient2d(C.row(0),C.row(3),C.row(1)) == Orientation::COLLINEAR) && + (orient2d(C.row(0),C.row(3),C.row(2)) == Orientation::COLLINEAR) ) + { + const auto v0 = (C.row(0) - q).eval(); + const auto v3 = (C.row(3) - q).eval(); + const auto w = std::atan2( + v0(0)*v3(1) - v0(1)*v3(0), + v0(0)*v3(0) + v0(1)*v3(1)) / + Scalar(2.0 * igl::PI); + return w; + } + return cubic_winding_number_helper(C,q); +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +template Eigen::Matrix::Scalar igl::predicates::cubic_winding_number, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +template Eigen::Matrix::Scalar igl::predicates::cubic_winding_number, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +#endif diff --git a/include/igl/predicates/cubic_winding_number.h b/include/igl/predicates/cubic_winding_number.h new file mode 100644 index 000000000..7439e467a --- /dev/null +++ b/include/igl/predicates/cubic_winding_number.h @@ -0,0 +1,30 @@ +#ifndef IGL_PREDICATES_CUBIC_WINDING_NUMBER_H +#define IGL_PREDICATES_CUBIC_WINDING_NUMBER_H +#include "../igl_inline.h" +#include + +namespace igl +{ + namespace predicates + { + /// Computes the (generalized) winding number of a cubic Bézier curve around a query. + /// This implementation is similar to "Robust Containment Queries Over + /// Collections of Rational Parametric Curves via Generalized Winding + /// Numbers" [Spainhour et al. 2024]. + /// + /// @param[in] C 4 by dimensions matrix of control points for a cubic + /// Bézier curve + /// @param[in] q 1 by dimensions query point + template + IGL_INLINE typename DerivedC::Scalar cubic_winding_number( + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& q); + } +} + +#ifndef IGL_STATIC_LIBRARY +# include "cubic_winding_number.cpp" +#endif + +#endif + diff --git a/include/igl/predicates/ear_clipping.cpp b/include/igl/predicates/ear_clipping.cpp index 8c045e6f6..38fd94c26 100644 --- a/include/igl/predicates/ear_clipping.cpp +++ b/include/igl/predicates/ear_clipping.cpp @@ -46,8 +46,8 @@ IGL_INLINE void igl::predicates::ear_clipping( Eigen::Matrix pb = P.row(b); Eigen::Matrix pi = P.row(i); auto r = igl::predicates::orient2d(pa, pi, pb); - if(r == igl::predicates::Orientation::NEGATIVE || - r == igl::predicates::Orientation::COLLINEAR) return false; + if(r == igl::Orientation::NEGATIVE || + r == igl::Orientation::COLLINEAR) return false; // check edge (b, R(b)) Index k = b; @@ -57,8 +57,8 @@ IGL_INLINE void igl::predicates::ear_clipping( auto r1 = igl::predicates::orient2d(pb, pl, pa); auto r2 = igl::predicates::orient2d(pi, pl, pb); if (l == a) return true; // single triangle - if (r1 != igl::predicates::Orientation::POSITIVE && - r2 != igl::predicates::Orientation::POSITIVE) { + if (r1 != igl::Orientation::POSITIVE && + r2 != igl::Orientation::POSITIVE) { return false; } @@ -83,8 +83,8 @@ IGL_INLINE void igl::predicates::ear_clipping( pl = P.row(l); r1 = igl::predicates::orient2d(pb, pk, pl); r2 = igl::predicates::orient2d(pi, pa, pk); - if (r1 != igl::predicates::Orientation::POSITIVE && - r2 != igl::predicates::Orientation::POSITIVE + if (r1 != igl::Orientation::POSITIVE && + r2 != igl::Orientation::POSITIVE ){ return false; } diff --git a/include/igl/predicates/incircle.cpp b/include/igl/predicates/incircle.cpp index 9f726fa9b..d550942fb 100644 --- a/include/igl/predicates/incircle.cpp +++ b/include/igl/predicates/incircle.cpp @@ -41,7 +41,7 @@ IGL_INLINE Orientation incircle( } #ifdef IGL_STATIC_LIBRARY -#define IGL_INCIRCLE(Vector) template igl::predicates::Orientation igl::predicates::incircle(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) +#define IGL_INCIRCLE(Vector) template igl::Orientation igl::predicates::incircle(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) #define IGL_MATRIX(T, R, C) Eigen::Matrix IGL_INCIRCLE(IGL_MATRIX(float, 1, 2)); IGL_INCIRCLE(IGL_MATRIX(float, 2, 1)); diff --git a/include/igl/predicates/incircle.h b/include/igl/predicates/incircle.h index 7d8e0544c..de4dcf236 100644 --- a/include/igl/predicates/incircle.h +++ b/include/igl/predicates/incircle.h @@ -11,7 +11,7 @@ #include "../igl_inline.h" #include -#include "Orientation.h" +#include "../Orientation.h" namespace igl { namespace predicates { diff --git a/include/igl/predicates/insphere.cpp b/include/igl/predicates/insphere.cpp index 73c265515..9f4410042 100644 --- a/include/igl/predicates/insphere.cpp +++ b/include/igl/predicates/insphere.cpp @@ -44,7 +44,7 @@ IGL_INLINE Orientation insphere( } #ifdef IGL_STATIC_LIBRARY -#define IGL_INSPHERE(Vector) template igl::predicates::Orientation igl::predicates::insphere(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) +#define IGL_INSPHERE(Vector) template igl::Orientation igl::predicates::insphere(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) #define IGL_MATRIX(T, R, C) Eigen::Matrix IGL_INSPHERE(IGL_MATRIX(float, 1, 3)); IGL_INSPHERE(IGL_MATRIX(float, 3, 1)); diff --git a/include/igl/predicates/insphere.h b/include/igl/predicates/insphere.h index 134463e07..7ab4dd1d2 100644 --- a/include/igl/predicates/insphere.h +++ b/include/igl/predicates/insphere.h @@ -11,7 +11,7 @@ #include "../igl_inline.h" #include -#include "Orientation.h" +#include "../Orientation.h" namespace igl { namespace predicates { diff --git a/include/igl/predicates/orient2d.cpp b/include/igl/predicates/orient2d.cpp index ac2da5c65..89cff5f5e 100644 --- a/include/igl/predicates/orient2d.cpp +++ b/include/igl/predicates/orient2d.cpp @@ -16,14 +16,38 @@ namespace predicates { using REAL = IGL_PREDICATES_REAL; #include "IGL_PREDICATES_ASSERT_SCALAR.h" -template +template < + typename Derivedpa, + typename Derivedpb, + typename Derivedpc> IGL_INLINE Orientation orient2d( - const Eigen::MatrixBase& pa, - const Eigen::MatrixBase& pb, - const Eigen::MatrixBase& pc) + const Eigen::MatrixBase& pa, + const Eigen::MatrixBase& pb, + const Eigen::MatrixBase& pc) { - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector2D, 2); - IGL_PREDICATES_ASSERT_SCALAR(Vector2D); + static_assert( + (Derivedpa::RowsAtCompileTime == 2 && Derivedpa::ColsAtCompileTime == 1) || + (Derivedpa::RowsAtCompileTime == 1 && Derivedpa::ColsAtCompileTime == 2) || + (Derivedpa::RowsAtCompileTime == Eigen::Dynamic && Derivedpa::ColsAtCompileTime == 1 ) || + (Derivedpa::RowsAtCompileTime == 1 && Derivedpa::ColsAtCompileTime == Eigen::Dynamic ), + "pa must be a 2D point"); + assert(pa.size() == 2 && "pa must be a 2D point"); + static_assert( + (Derivedpb::RowsAtCompileTime == 2 && Derivedpb::ColsAtCompileTime == 1) || + (Derivedpb::RowsAtCompileTime == 1 && Derivedpb::ColsAtCompileTime == 2) || + (Derivedpb::RowsAtCompileTime == Eigen::Dynamic && Derivedpb::ColsAtCompileTime == 1 ) || + (Derivedpb::RowsAtCompileTime == 1 && Derivedpb::ColsAtCompileTime == Eigen::Dynamic ), + "pb must be a 2D point"); + assert(pb.size() == 2 && "pb must be a 2D point"); + static_assert( + (Derivedpc::RowsAtCompileTime == 2 && Derivedpc::ColsAtCompileTime == 1) || + (Derivedpc::RowsAtCompileTime == 1 && Derivedpc::ColsAtCompileTime == 2) || + (Derivedpc::RowsAtCompileTime == Eigen::Dynamic && Derivedpc::ColsAtCompileTime == 1 ) || + (Derivedpc::RowsAtCompileTime == 1 && Derivedpc::ColsAtCompileTime == Eigen::Dynamic ), + "pc must be a 2D point"); + assert(pc.size() == 2 && "pc must be a 2D point"); + + using Point = Eigen::Matrix; Point a{pa[0], pa[1]}; @@ -72,14 +96,11 @@ IGL_INLINE void orient2d( } #ifdef IGL_STATIC_LIBRARY -#define IGL_ORIENT2D(Vector) template igl::predicates::Orientation igl::predicates::orient2d(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) -#define IGL_MATRIX(T, R, C) Eigen::Matrix -IGL_ORIENT2D(IGL_MATRIX(float, 1, 2)); -IGL_ORIENT2D(IGL_MATRIX(float, 2, 1)); -#ifndef LIBIGL_PREDICATES_USE_FLOAT -IGL_ORIENT2D(IGL_MATRIX(double, 1, 2)); -IGL_ORIENT2D(IGL_MATRIX(double, 2, 1)); -#endif -#undef IGL_MATRIX -#undef IGL_ORIENT2D +// Explicit template instantiation +template igl::Orientation igl::predicates::orient2d const, 1, -1, true>, Eigen::Block const, 1, -1, true>, Eigen::Block const, 1, -1, true>>(Eigen::MatrixBase const, 1, -1, true>> const&, Eigen::MatrixBase const, 1, -1, true>> const&, Eigen::MatrixBase const, 1, -1, true>> const&); +template igl::Orientation igl::predicates::orient2d const, 1, 2, false>, Eigen::Block const, 1, 2, false>, Eigen::Block const, 1, 2, false>>(Eigen::MatrixBase const, 1, 2, false>> const&, Eigen::MatrixBase const, 1, 2, false>> const&, Eigen::MatrixBase const, 1, 2, false>> const&); +template igl::Orientation igl::predicates::orient2d, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +template igl::Orientation igl::predicates::orient2d, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +template igl::Orientation igl::predicates::orient2d, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +template igl::Orientation igl::predicates::orient2d, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); #endif diff --git a/include/igl/predicates/orient2d.h b/include/igl/predicates/orient2d.h index 1aadb02d3..432b2b5f1 100644 --- a/include/igl/predicates/orient2d.h +++ b/include/igl/predicates/orient2d.h @@ -11,7 +11,7 @@ #include "../igl_inline.h" #include -#include "Orientation.h" +#include "../Orientation.h" namespace igl { namespace predicates { @@ -25,11 +25,14 @@ namespace igl { /// COLLINEAR if they are collinear. /// /// \fileinfo - template + template< + typename Derivedpa, + typename Derivedpb, + typename Derivedpc> IGL_INLINE Orientation orient2d( - const Eigen::MatrixBase& pa, - const Eigen::MatrixBase& pb, - const Eigen::MatrixBase& pc); + const Eigen::MatrixBase& pa, + const Eigen::MatrixBase& pb, + const Eigen::MatrixBase& pc); /// Compute the orientation of the tetrahedron formed by each 4-tuple of /// points /// diff --git a/include/igl/predicates/orient3d.cpp b/include/igl/predicates/orient3d.cpp index 3c71f906e..afbe1aafe 100644 --- a/include/igl/predicates/orient3d.cpp +++ b/include/igl/predicates/orient3d.cpp @@ -81,7 +81,7 @@ IGL_INLINE void orient3d( // Explicit template instantiation template void igl::predicates::orient3d, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&); -#define IGL_ORIENT3D(Vector) template igl::predicates::Orientation igl::predicates::orient3d(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) +#define IGL_ORIENT3D(Vector) template igl::Orientation igl::predicates::orient3d(const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&, const Eigen::MatrixBase&) #define IGL_MATRIX(T, R, C) Eigen::Matrix IGL_ORIENT3D(IGL_MATRIX(float, 1, 3)); IGL_ORIENT3D(IGL_MATRIX(float, 3, 1)); diff --git a/include/igl/predicates/orient3d.h b/include/igl/predicates/orient3d.h index 0863cd1c9..78c369833 100644 --- a/include/igl/predicates/orient3d.h +++ b/include/igl/predicates/orient3d.h @@ -11,7 +11,7 @@ #include "../igl_inline.h" #include -#include "Orientation.h" +#include "../Orientation.h" namespace igl { namespace predicates { diff --git a/include/igl/predicates/point_in_convex_hull.cpp b/include/igl/predicates/point_in_convex_hull.cpp new file mode 100644 index 000000000..18e062f22 --- /dev/null +++ b/include/igl/predicates/point_in_convex_hull.cpp @@ -0,0 +1,204 @@ +#include "point_in_convex_hull.h" +#include "orient2d.h" +#include +#include + +// Lightly edited LLM code. + +namespace +{ + using igl::Orientation; +template +struct SH +{ + static inline int osign(const Orientation o) + { + using namespace igl::predicates; + switch(o) + { + case Orientation::POSITIVE: return +1; + case Orientation::NEGATIVE: return -1; + case Orientation::COLLINEAR: return 0; + } + return 0; + } + + static inline bool between_1d(const Scalar x, Scalar lo, Scalar hi) + { + if(lo > hi) std::swap(lo, hi); + return (lo <= x) && (x <= hi); + } + + // Assumes q is collinear with segment endpoints (by exact orient2d). + static inline bool on_segment_collinear( + const Scalar qx, const Scalar qy, + const Scalar ax, const Scalar ay, + const Scalar bx, const Scalar by) + { + // Use dominant axis to avoid issues with vertical-ish segments. + const Scalar dx = std::abs(bx - ax); + const Scalar dy = std::abs(by - ay); + if(dx >= dy) return between_1d(qx, ax, bx); + else return between_1d(qy, ay, by); + } + + // Closed triangle membership using only cached signs; handles degenerate triples. + static inline bool in_triangle_cached( + const int o_ab, const int o_bc, const int o_ca, + const Scalar qx, const Scalar qy, + const Scalar ax, const Scalar ay, + const Scalar bx, const Scalar by, + const Scalar cx, const Scalar cy) + { + const bool all_nonneg = (o_ab >= 0 && o_bc >= 0 && o_ca >= 0); + const bool all_nonpos = (o_ab <= 0 && o_bc <= 0 && o_ca <= 0); + if(!(all_nonneg || all_nonpos)) + return false; + + // Fully collinear (or coincident) triple: triangle degenerates to a segment/point. + if(o_ab == 0 && o_bc == 0 && o_ca == 0) + { + const Scalar xmin = std::min(ax, std::min(bx, cx)); + const Scalar xmax = std::max(ax, std::max(bx, cx)); + const Scalar ymin = std::min(ay, std::min(by, cy)); + const Scalar ymax = std::max(ay, std::max(by, cy)); + const Scalar sx = xmax - xmin; + const Scalar sy = ymax - ymin; + + if(sx == 0.0 && sy == 0.0) + return (qx == ax) && (qy == ay); + + if(sx >= sy) return between_1d(qx, xmin, xmax); + else return between_1d(qy, ymin, ymax); + } + + return true; + } + +}; +} + +template < + typename Derivedq, + typename Deriveda, + typename Derivedb, + typename Derivedc, + typename Derivedd> +IGL_INLINE igl::Orientation igl::predicates::point_in_convex_hull( + const Eigen::MatrixBase & q, + const Eigen::MatrixBase & a, + const Eigen::MatrixBase & b, + const Eigen::MatrixBase & c, + const Eigen::MatrixBase & d) +{ + assert(q.size()==2 && a.size()==2 && b.size()==2 && c.size()==2 && d.size()==2); + using Scalar = typename Derivedq::Scalar; + using SH = ::SH; + using igl::Orientation; + + // For cheap coordinate comparisons (between tests) + const Scalar qx = (Scalar)q(0), qy = (Scalar)q(1); + const Scalar px[4] = {(Scalar)a(0), (Scalar)b(0), (Scalar)c(0), (Scalar)d(0)}; + const Scalar py[4] = {(Scalar)a(1), (Scalar)b(1), (Scalar)c(1), (Scalar)d(1)}; + + // Cache the 6 unordered edge orientations with q: sign(orient2d(p_i,p_j,q)) for i Oq = [&](int i, int j)->int + { + // directed sign(orient2d(p_i,p_j,q)) via antisymmetry + if(i == j) return 0; + if(i > j) return -Oq(j,i); + // ibool + { + const int o_ab = Oq(ia,ib); + const int o_bc = Oq(ib,ic); + const int o_ca = Oq(ic,ia); + return SH::in_triangle_cached( + o_ab, o_bc, o_ca, + qx, qy, + px[ia], py[ia], + px[ib], py[ib], + px[ic], py[ic]); + }; + + const bool inside = + tri(0,1,2) || tri(0,1,3) || tri(0,2,3) || tri(1,2,3); + + if(!inside) + return Orientation::NEGATIVE; + + // --- Boundary classification --- + // A segment (i,j) is a supporting segment of conv(P) iff the other two points + // lie in the same closed half-plane of the line through (i,j): + // sign(orient(p_i,p_j,p_k)) and sign(orient(p_i,p_j,p_l)) are NOT opposite. + auto supporting_pair = [&](int i, int j)->bool + { + int other[2], t=0; + for(int k=0;k<4;k++) if(k!=i && k!=j) other[t++] = k; + + // Choose Eigen refs for orient2d calls + auto Pi = [&](int idx)->auto const& { + switch(idx){ + case 0: return a; + case 1: return b; + case 2: return c; + default:return d; + } + }; + + const int s1 = SH::osign(orient2d(Pi(i), Pi(j), Pi(other[0]))); + const int s2 = SH::osign(orient2d(Pi(i), Pi(j), Pi(other[1]))); + + // Opposite strict signs => not supporting + return !((s1 > 0 && s2 < 0) || (s1 < 0 && s2 > 0)); + }; + + // If hull is fully collinear (degenerate), any in-hull point is "COLLINEAR". + const bool nondegenerate = + (orient2d(a,b,c) != Orientation::COLLINEAR) || + (orient2d(a,b,d) != Orientation::COLLINEAR) || + (orient2d(a,c,d) != Orientation::COLLINEAR) || + (orient2d(b,c,d) != Orientation::COLLINEAR); + + if(!nondegenerate) + return Orientation::COLLINEAR; + + // For non-degenerate hull: q is on boundary iff it lies on some supporting segment. + for(int i=0;i<4;i++) + for(int j=i+1;j<4;j++) + { + if(!supporting_pair(i,j)) continue; + + if(Oq(i,j) == 0) // q collinear with (p_i, p_j) by exact orient2d + { + if(SH::on_segment_collinear(qx,qy, px[i],py[i], px[j],py[j])) + return Orientation::COLLINEAR; + } + } + + return Orientation::POSITIVE; // inside and not on boundary +} + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +template igl::Orientation igl::predicates::point_in_convex_hull, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +template igl::Orientation igl::predicates::point_in_convex_hull, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +template igl::Orientation igl::predicates::point_in_convex_hull, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&); +#endif diff --git a/include/igl/predicates/point_in_convex_hull.h b/include/igl/predicates/point_in_convex_hull.h new file mode 100644 index 000000000..d84736b49 --- /dev/null +++ b/include/igl/predicates/point_in_convex_hull.h @@ -0,0 +1,43 @@ +#ifndef IGL_PREDICATES_POINT_IN_CONVEX_POLYGON_H +#define IGL_PREDICATES_POINT_IN_CONVEX_POLYGON_H + +#include "../igl_inline.h" +#include +#include "../Orientation.h" + +namespace igl +{ + namespace predicates + { + /// Determine if the 2D point q is inside, outside, or on the boundary of the + /// convex hull of the 2D points a,b,c,d. The points a,b,c,d may be given in + /// any order and the convex hull may be a point, segment, triangle or + /// quadrilateral. + /// + /// @param[in] q 2D query point + /// @param[in] a 2D point possibly on the convex hull + /// @param[in] b 2D point possibly on the convex hull + /// @param[in] c 2D point possibly on the convex hull + /// @param[in] d 2D point possibly on the convex hull + /// @return Returns +1 if q is inside the convex hull, -1 if outside, and 0 + /// if on the boundary. + template < + typename Derivedq, + typename Deriveda, + typename Derivedb, + typename Derivedc, + typename Derivedd> + IGL_INLINE Orientation point_in_convex_hull( + const Eigen::MatrixBase & q, + const Eigen::MatrixBase & a, + const Eigen::MatrixBase & b, + const Eigen::MatrixBase & c, + const Eigen::MatrixBase & d); + } +} + +#ifndef IGL_STATIC_LIBRARY +#include "point_in_convex_hull.cpp" +#endif + +#endif diff --git a/include/igl/predicates/point_inside_convex_polygon.cpp b/include/igl/predicates/point_inside_convex_polygon.cpp index 41af128fd..16c459e13 100644 --- a/include/igl/predicates/point_inside_convex_polygon.cpp +++ b/include/igl/predicates/point_inside_convex_polygon.cpp @@ -22,8 +22,8 @@ IGL_INLINE bool igl::predicates::point_inside_convex_polygon( Eigen::Matrix a = P.row(i); Eigen::Matrix b = P.row(i_1); auto r = igl::predicates::orient2d(a,b,q); - if(r == igl::predicates::Orientation::COLLINEAR || - r == igl::predicates::Orientation::NEGATIVE) + if(r == igl::Orientation::COLLINEAR || + r == igl::Orientation::NEGATIVE) return false; } return true; diff --git a/include/igl/predicates/segment_segment_intersect.cpp b/include/igl/predicates/segment_segment_intersect.cpp index d3374b59e..d720a1685 100644 --- a/include/igl/predicates/segment_segment_intersect.cpp +++ b/include/igl/predicates/segment_segment_intersect.cpp @@ -36,10 +36,10 @@ IGL_INLINE bool igl::predicates::segment_segment_intersect( }; // colinear case - if((t1 == igl::predicates::Orientation::COLLINEAR && on_segment(a,b,c)) || - (t2 == igl::predicates::Orientation::COLLINEAR && on_segment(c,d,b)) || - (t3 == igl::predicates::Orientation::COLLINEAR && on_segment(a,b,d)) || - (t4 == igl::predicates::Orientation::COLLINEAR && on_segment(c,d,a))) + if((t1 == igl::Orientation::COLLINEAR && on_segment(a,b,c)) || + (t2 == igl::Orientation::COLLINEAR && on_segment(c,d,b)) || + (t3 == igl::Orientation::COLLINEAR && on_segment(a,b,d)) || + (t4 == igl::Orientation::COLLINEAR && on_segment(c,d,a))) return true; // ordinary case diff --git a/include/igl/predicates/spline_winding_number.cpp b/include/igl/predicates/spline_winding_number.cpp new file mode 100644 index 000000000..4082e5879 --- /dev/null +++ b/include/igl/predicates/spline_winding_number.cpp @@ -0,0 +1,51 @@ +#include "spline_winding_number.h" +#include "cubic_winding_number.h" +#include "../parallel_for.h" +#include "../eytzinger_aabb.h" +#include "../eytzinger_aabb.h" +#include "../eytzinger_aabb_winding_number_tree.h" +#include "../eytzinger_aabb_winding_number.h" + +template < + typename DerivedP, + typename DerivedC, + typename DerivedB, + typename Derivedleaf, + typename DerivedQ, + typename DerivedW> +IGL_INLINE void igl::predicates::spline_winding_number( + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& B1, + const Eigen::MatrixBase& B2, + const Eigen::MatrixBase& leaf, + const Eigen::MatrixBase & Q, + Eigen::PlainObjectBase& W) +{ + using Scalar = typename DerivedP::Scalar; + Eigen::Matrix E(C.rows(),2); + E << C.col(0), C.col(3); + Eigen::VectorXi tI,tC; + eytzinger_aabb_winding_number_tree( E, leaf, tI, tC); + + W.resize(Q.rows(),1); + //for(int i = 0;i Mat4; + // point_spline_squared_distance is caching this slice. + const Mat4 Cc = P(C.row(c),Eigen::all); + return igl::predicates::cubic_winding_number(Cc, qi); + }; + igl::eytzinger_aabb_winding_number(qi,P,primitive,B1,B2,leaf,tI,tC,W(i)); + },1000); +} + + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +template void igl::predicates::spline_winding_number, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::MatrixBase> const&, Eigen::PlainObjectBase>&); +#endif diff --git a/include/igl/predicates/spline_winding_number.h b/include/igl/predicates/spline_winding_number.h new file mode 100644 index 000000000..860073b8e --- /dev/null +++ b/include/igl/predicates/spline_winding_number.h @@ -0,0 +1,50 @@ +#ifndef IGL_PREDICATES_SPLINE_WINDING_NUMBER_H +#define IGL_PREDICATES_SPLINE_WINDING_NUMBER_H +#include "../igl_inline.h" +#include + +namespace igl +{ + namespace predicates + { + /// Computes the (generalized) winding number of a spline of cubic Bézier + /// curves around a set of query points. + /// + /// @param[in] q 1 by dim query point + /// @param[in] P #P by dim matrix of spline control points + /// @param[in] C #C by 4 matrix of indices into P defining the cubic Bézier + /// curves making up the spline + /// @param[in] B1 #B by dim matrix of AABB min box corners + /// @param[in] B2 #B by dim matrix of AABB max box corners + /// @param[in] leaf #B by 1 matrix of AABB leaf node indices/flags + /// @param[in] Q #Q by dim matrix of query points + /// @param[out] W #Q by 1 matrix of winding numbers for each query point + /// + /// + /// \see igl::cycodebase::box_cubic, igl::predicates::cubic_winding_number, + /// igl::eytzinger_aabb, igl::cycodebase::spline_eytzinger_aabb + template < + typename DerivedP, + typename DerivedC, + typename DerivedB, + typename Derivedleaf, + typename DerivedQ, + typename DerivedW> + IGL_INLINE void spline_winding_number( + const Eigen::MatrixBase& P, + const Eigen::MatrixBase& C, + const Eigen::MatrixBase& B1, + const Eigen::MatrixBase& B2, + const Eigen::MatrixBase& leaf, + const Eigen::MatrixBase & Q, + Eigen::PlainObjectBase& W); + } +} + +#ifndef IGL_STATIC_LIBRARY +# include "spline_winding_number.cpp" +#endif + +#endif + + diff --git a/include/igl/predicates/triangle_triangle_intersect.cpp b/include/igl/predicates/triangle_triangle_intersect.cpp index 15e09f61f..fa65811e2 100644 --- a/include/igl/predicates/triangle_triangle_intersect.cpp +++ b/include/igl/predicates/triangle_triangle_intersect.cpp @@ -21,9 +21,9 @@ IGL_INLINE bool igl::predicates::triangle_triangle_intersect( exactinit(); using Vector2D = Eigen::Matrix; - constexpr Orientation COPLANAR = igl::predicates::Orientation::COPLANAR; - constexpr Orientation NEGATIVE = igl::predicates::Orientation::NEGATIVE; - constexpr Orientation POSITIVE = igl::predicates::Orientation::POSITIVE; + constexpr Orientation COPLANAR = igl::Orientation::COPLANAR; + constexpr Orientation NEGATIVE = igl::Orientation::NEGATIVE; + constexpr Orientation POSITIVE = igl::Orientation::POSITIVE; // Determine for each vertex of one triangle if it's above, below or on the // plane of the other triangle. diff --git a/include/igl/predicates/triangle_triangle_intersect.h b/include/igl/predicates/triangle_triangle_intersect.h index 45261f3ab..74050df69 100644 --- a/include/igl/predicates/triangle_triangle_intersect.h +++ b/include/igl/predicates/triangle_triangle_intersect.h @@ -6,6 +6,17 @@ namespace igl { namespace predicates { + /// Triangle-triangle intersection test using exact predicates. + /// + /// @param[in] a1 First vertex of triangle A. + /// @param[in] a2 Second vertex of triangle A. + /// @param[in] a3 Third vertex of triangle A. + /// @param[in] b1 First vertex of triangle B. + /// @param[in] b2 Second vertex of triangle B. + /// @param[in] b3 Third vertex of triangle B. + /// @param[out] coplanar True if the triangles are coplanar. + /// @return True if the triangles intersect. + /// template IGL_INLINE bool triangle_triangle_intersect( const Vector3D & a1, diff --git a/include/igl/variable_radius_offset.cpp b/include/igl/variable_radius_offset.cpp index 5dd45ebef..fee686d7e 100644 --- a/include/igl/variable_radius_offset.cpp +++ b/include/igl/variable_radius_offset.cpp @@ -113,12 +113,9 @@ void igl::variable_radius_offset( Eigen::PlainObjectBase & mF) { using RowVector3S = Eigen::Matrix; - IGL_TICTOC_LAMBDA; - tictoc(); Eigen::Matrix B1,B2; Eigen::VectorXi leaf; igl::eytzinger_aabb(PB1,PB2,B1,B2,leaf); - printf("%-20s: %g secs\n","eytzinger_aabb",tictoc()); const std::function sdf = [&](const RowVector3S & p) -> Scalar @@ -128,7 +125,7 @@ void igl::variable_radius_offset( return primitive(p,j); }; Scalar f; - igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f); + igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f); return f; }; const std::function @@ -157,18 +154,14 @@ void igl::variable_radius_offset( }; Eigen::Matrix ijk; igl::lipschitz_octree( origin,h0,max_depth,udf,ijk); - printf("%-20s: %g secs\n","lipschitz_octree",tictoc()); { - tictoc(); // Gather the corners of those leaf cells const Scalar h = h0 / (1 << max_depth); Eigen::Matrix J; Eigen::Matrix unique_ijk; Eigen::Matrix unique_corner_positions; igl::unique_sparse_voxel_corners(origin,h0,max_depth,ijk,unique_ijk,J,unique_corner_positions); - //printf("unique_sparse_voxel_corners: %0.7f seconds\n",tictoc()); - printf("%-20s: %g secs\n","unique_sparse_vo...",tictoc()); /// Evaluate the signed distance function at the corners Eigen::VectorXd S(unique_corner_positions.rows()); //for(int u = 0;u +#include + +TEST_CASE("cubic_is_flat: simple", "[igl]" ) +{ + { + Eigen::Matrix C; + C << 0,0, + 1,1, + 2,-1, + 3,0; + REQUIRE( igl::cubic_is_flat(C,1e-2) == false); + } + { + Eigen::Matrix C; + C << 0,0, + 1,0.1, + 2,-0.1, + 3,0; + REQUIRE( igl::cubic_is_flat(C,1e-2) == true); + } + { + Eigen::Matrix C; + C << 0,0, + 1,0.1, + 2,-0.1, + 3,0; + REQUIRE( igl::cubic_is_flat(C,1e-4) == false); + } +} + +TEST_CASE("cubic_is_flat: degenerate", "[igl]" ) +{ + Eigen::Matrix C; + C << 0,0, + 0,0, + 0,0, + 0,0; + REQUIRE( igl::cubic_is_flat(C,1) == true); +} + diff --git a/tests/include/igl/cubic_split.cpp b/tests/include/igl/cubic_split.cpp new file mode 100644 index 000000000..e81327a30 --- /dev/null +++ b/tests/include/igl/cubic_split.cpp @@ -0,0 +1,36 @@ +#include +#include + +TEST_CASE("cubic_split: simple", "[igl]" ) +{ + { + Eigen::Matrix C; + C << 0,0, + 1,1, + 2,-1, + 3,0; + double t = 0.5; + Eigen::RowVector2d C01,C012,C0123,C123,C23; + igl::cubic_split(C,t,C01,C012,C0123,C123,C23); + // C01 = (0.5,0.5) + // C012 = (1.0,0.25) + // C0123 = (1.5,0.0) + // C123 = (2.0,-0.25) + // C23 = (2.5,-0.5) + REQUIRE( C01.isApprox( Eigen::RowVector2d(0.5,0.5) ) ); + REQUIRE( C012.isApprox( Eigen::RowVector2d(1.0,0.25) ) ); + REQUIRE( C0123.isApprox( Eigen::RowVector2d(1.5,0.0) ) ); + REQUIRE( C123.isApprox( Eigen::RowVector2d(2.0,-0.25) ) ); + REQUIRE( C23.isApprox( Eigen::RowVector2d(2.5,-0.5) ) ); + Eigen::Matrix C1,C2; + igl::cubic_split(C,t,C1,C2); + REQUIRE( C1.row(0).isApprox( C.row(0) ) ); + REQUIRE( C1.row(1).isApprox( C01 ) ); + REQUIRE( C1.row(2).isApprox( C012 ) ); + REQUIRE( C1.row(3).isApprox( C0123 ) ); + REQUIRE( C2.row(0).isApprox( C0123 ) ); + REQUIRE( C2.row(1).isApprox( C123 ) ); + REQUIRE( C2.row(2).isApprox( C23 ) ); + REQUIRE( C2.row(3).isApprox( C.row(3) ) ); + } +} diff --git a/tests/include/igl/cycodebase/box_cubic.cpp b/tests/include/igl/cycodebase/box_cubic.cpp new file mode 100644 index 000000000..1b71168fb --- /dev/null +++ b/tests/include/igl/cycodebase/box_cubic.cpp @@ -0,0 +1,20 @@ +#include +#include + +TEST_CASE("box_cubic: simple", "[igl/cycodebase]" ) +{ + Eigen::MatrixXd C(4,2); + C<<0,0, + 1,2, + 2,-2, + 3,0; + Eigen::RowVectorXd B1,B2; + igl::cycodebase::box_cubic(C,B1,B2); + REQUIRE(B1.size() == 2); + REQUIRE(B2.size() == 2); + REQUIRE(B1[0] == Approx(0.0).margin(1e-12)); + REQUIRE(B1[1] == Approx(-0.57735026918962584).margin(1e-12)); + REQUIRE(B2[0] == Approx(3.0).margin(1e-12)); + REQUIRE(B2[1] == Approx(0.57735026918962584).margin(1e-12)); +} + diff --git a/tests/include/igl/cycodebase/point_cubic_squared_distance.cpp b/tests/include/igl/cycodebase/point_cubic_squared_distance.cpp new file mode 100644 index 000000000..cdc1efd19 --- /dev/null +++ b/tests/include/igl/cycodebase/point_cubic_squared_distance.cpp @@ -0,0 +1,39 @@ +#include +#include + +TEST_CASE("point_cubic_squared_distance: simple2d", "[igl/cycodebase]" ) +{ + Eigen::MatrixXd C(4,2); + C<<0,0, + 1,2, + 2,-2, + 3,0; + Eigen::MatrixXd Q(3,2); + //[1.5 0;2 0.5;2.5 1] + Q<< + 1.5,0, + 2,0.5, + 2.5,1; + Eigen::VectorXd sqrD,S; + Eigen::MatrixXd K; + igl::cycodebase::point_cubic_squared_distance(Q,C,sqrD,S,K); + REQUIRE(sqrD.size() == 3); + REQUIRE(S.size() == 3); + REQUIRE(K.rows() == 3); + REQUIRE(K.cols() == 2); + // sqrd = [0;0.5;1.25] + REQUIRE(sqrD[0] == Approx(0.0).margin(1e-12)); + REQUIRE(sqrD[1] == Approx(0.5).margin(1e-12)); + REQUIRE(sqrD[2] == Approx(1.25).margin(1e-12)); + // S = [0.5;0.5;1.0] + REQUIRE(S[0] == Approx(0.5).margin(1e-12)); + REQUIRE(S[1] == Approx(0.5).margin(1e-12)); + REQUIRE(S[2] == Approx(1.0).margin(1e-12)); + // [1.5 0;1.5 0;3 0] + REQUIRE(K(0,0) == Approx(1.5).margin(1e-12)); + REQUIRE(K(0,1) == Approx(0.0).margin(1e-12)); + REQUIRE(K(1,0) == Approx(1.5).margin(1e-12)); + REQUIRE(K(1,1) == Approx(0.0).margin(1e-12)); + REQUIRE(K(2,0) == Approx(3.0).margin(1e-12)); + REQUIRE(K(2,1) == Approx(0.0).margin(1e-12)); +} diff --git a/tests/include/igl/cycodebase/roots.cpp b/tests/include/igl/cycodebase/roots.cpp new file mode 100644 index 000000000..16a0a5165 --- /dev/null +++ b/tests/include/igl/cycodebase/roots.cpp @@ -0,0 +1,43 @@ +#include +#include + +TEST_CASE("roots: cubic", "[igl/cycodebase]" ) +{ + // t^3 - 6*t^2 + 11*t - 6 + { + Eigen::Vector coef(4); + coef << -6,11,-6,1; + Eigen::Vector R; + const int nr = igl::cycodebase::roots(coef, 0.0, 4.0, R); + REQUIRE(nr == 3); + REQUIRE(R[0] == Approx(1.0).margin(1e-12)); + REQUIRE(R[1] == Approx(2.0).margin(1e-12)); + REQUIRE(R[2] == Approx(3.0).margin(1e-12)); + } + + // With dynamic size + { + Eigen::VectorXd coef(4); + coef << -6,11,-6,1; + { + Eigen::VectorXd R; + const int nr = igl::cycodebase::roots(coef, 0.0, 4.0, R); + REQUIRE(nr == 3); + REQUIRE(R.size() == 3); + REQUIRE(R[0] == Approx(1.0).margin(1e-12)); + REQUIRE(R[1] == Approx(2.0).margin(1e-12)); + REQUIRE(R[2] == Approx(3.0).margin(1e-12)); + } + // Only pick the first root + { + Eigen::VectorXd R; + const int nr = igl::cycodebase::roots(coef, 0.0, 1.5, R); + REQUIRE(nr == 1); + REQUIRE(R.size() == 3); + REQUIRE(R[0] == Approx(1.0).margin(1e-12)); + REQUIRE(std::isnan(R[1])); + REQUIRE(std::isnan(R[2])); + } + } +} + diff --git a/tests/include/igl/predicates/cubic_winding_number.cpp b/tests/include/igl/predicates/cubic_winding_number.cpp new file mode 100644 index 000000000..804ca3b95 --- /dev/null +++ b/tests/include/igl/predicates/cubic_winding_number.cpp @@ -0,0 +1,77 @@ +#include +#include +#include + +TEST_CASE("cubic_winding_number: simple", "[igl/predicates]") +{ + Eigen::Matrix C; + C<<0,0, + 1,1, + 2,-1, + 3,0; + { + Eigen::RowVector2d q(1.1,1.1); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(0.29147615882815 == Approx(w).epsilon(1e-12)); + } + { + Eigen::RowVector2d q(1.45,0.01); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(-0.502124394734 == Approx(w).epsilon(1e-12)); + } + { + Eigen::RowVector2d q(1.45,0); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(-0.5 == Approx(w).epsilon(1e-12)); + } + // Stress test the numerics + { + Eigen::RowVector2d q(1.5+1e-15,0); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(0.5 == Approx(w).epsilon(1e-12)); + } + { + Eigen::RowVector2d q(1.5-1e-15,0); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(-0.5 == Approx(w).epsilon(1e-12)); + } + { + Eigen::RowVector2d q(1.5+1e-16,1e-16); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(0.5 == Approx(w).epsilon(1e-12)); + } + { + Eigen::RowVector2d q(1.5-1e-16,-1e-16); + double w = igl::predicates::cubic_winding_number(C,q); + REQUIRE(-0.5 == Approx(w).epsilon(1e-12)); + } +} + +TEST_CASE("cubic_winding_number: degenerate", "[igl/predicates]") +{ + Eigen::RowVector2d q(1.2,1.0); + { + Eigen::Matrix C; + C<<0,0, + 0,0, + 0,0, + 0,0; + REQUIRE( igl::predicates::cubic_winding_number(C,q) == 0.0); + } + { + Eigen::Matrix C; + C<<0,0, + 1,0, + 2,0, + 3,0; + REQUIRE(0.3087217355796 == Approx(igl::predicates::cubic_winding_number(C,q)).epsilon(1e-12)); + } + { + Eigen::Matrix C; + C<<0,0, + 1,1, + 3,3, + 9,9; + REQUIRE(-0.4835565188700 == Approx(igl::predicates::cubic_winding_number(C,q)).epsilon(1e-12)); + } +} diff --git a/tests/include/igl/predicates/orient3d.cpp b/tests/include/igl/predicates/orient3d.cpp index a5f993cd6..0df8b726e 100644 --- a/tests/include/igl/predicates/orient3d.cpp +++ b/tests/include/igl/predicates/orient3d.cpp @@ -35,9 +35,9 @@ TEST_CASE("orient3d", "[igl][predicates]") Eigen::VectorXi R; igl::predicates::orient3d(A,B,C,D,R); REQUIRE(R.size() == 3); - REQUIRE(R(0) == int(igl::predicates::Orientation::NEGATIVE)); - REQUIRE(R(1) == int(igl::predicates::Orientation::POSITIVE)); - REQUIRE(R(2) == int(igl::predicates::Orientation::COPLANAR)); + REQUIRE(R(0) == int(igl::Orientation::NEGATIVE)); + REQUIRE(R(1) == int(igl::Orientation::POSITIVE)); + REQUIRE(R(2) == int(igl::Orientation::COPLANAR)); } { Eigen::MatrixXd A(1,3); @@ -57,9 +57,9 @@ TEST_CASE("orient3d", "[igl][predicates]") Eigen::VectorXi R; igl::predicates::orient3d(A,B,C,D,R); REQUIRE(R.size() == 3); - REQUIRE(R(0) == int(igl::predicates::Orientation::NEGATIVE)); - REQUIRE(R(1) == int(igl::predicates::Orientation::POSITIVE)); - REQUIRE(R(2) == int(igl::predicates::Orientation::COPLANAR)); + REQUIRE(R(0) == int(igl::Orientation::NEGATIVE)); + REQUIRE(R(1) == int(igl::Orientation::POSITIVE)); + REQUIRE(R(2) == int(igl::Orientation::COPLANAR)); } } diff --git a/tests/include/igl/predicates/point_in_convex_hull.cpp b/tests/include/igl/predicates/point_in_convex_hull.cpp new file mode 100644 index 000000000..8b9615cc3 --- /dev/null +++ b/tests/include/igl/predicates/point_in_convex_hull.cpp @@ -0,0 +1,83 @@ +#include +#include + +TEST_CASE("point_in_convex_hull: simple", "[igl/predicates]") +{ + Eigen::RowVector2d a(0,0); + Eigen::RowVector2d b(1,1); + Eigen::RowVector2d c(2,-1); + Eigen::RowVector2d d(3,0); + + { + Eigen::RowVector2d q(1.5,0); + auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d); + REQUIRE(r == igl::Orientation::POSITIVE); + } + { + Eigen::RowVector2d q(-1,0); + auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d); + REQUIRE(r == igl::Orientation::NEGATIVE); + } + { + Eigen::RowVector2d q(0,0); + auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d); + REQUIRE(r == igl::Orientation::COLLINEAR); + } + { + Eigen::RowVector2d q(0.5,0.5); + auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d); + REQUIRE(r == igl::Orientation::COLLINEAR); + } + { + Eigen::RowVector2d q(1.0,-0.4); + auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d); + REQUIRE(r == igl::Orientation::POSITIVE); + } + { + Eigen::RowVector2d q(2.0,0.6); + auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d); + REQUIRE(r == igl::Orientation::NEGATIVE); + } +} + +TEST_CASE("point_in_convex_hull: degenerate", "[igl/predicates]") +{ + // Triangle + { + Eigen::RowVector2d a(0,0); + Eigen::RowVector2d b(1,0); + Eigen::RowVector2d c(0,1); + Eigen::RowVector2d d(0.1,0.1); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.1,0.1),a,b,c,d) == igl::Orientation::POSITIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.2,0.2),a,b,c,d) == igl::Orientation::POSITIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(-0.2,-0.2),a,b,c,d) == igl::Orientation::NEGATIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.5,0.5),a,b,c,d) == igl::Orientation::COLLINEAR); + } + + // Segment + { + Eigen::RowVector2d a(0,0); + Eigen::RowVector2d b(1,0); + Eigen::RowVector2d c(0.5,0); + Eigen::RowVector2d d(0.25,0); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.25,0),a,b,c,d) == igl::Orientation::COLLINEAR); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.5,0),a,b,c,d) == igl::Orientation::COLLINEAR); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(1.5,0),a,b,c,d) == igl::Orientation::NEGATIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(-0.5,0),a,b,c,d) == igl::Orientation::NEGATIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.5,1),a,b,c,d) == igl::Orientation::NEGATIVE); + } + + // Point + { + Eigen::RowVector2d a(0,0); + Eigen::RowVector2d b(0,0); + Eigen::RowVector2d c(0,0); + Eigen::RowVector2d d(0,0); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0,0),a,b,c,d) == igl::Orientation::COLLINEAR); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(1,0),a,b,c,d) == igl::Orientation::NEGATIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0,1),a,b,c,d) == igl::Orientation::NEGATIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(-1,0),a,b,c,d) == igl::Orientation::NEGATIVE); + REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0,-1),a,b,c,d) == igl::Orientation::NEGATIVE); + } + +} diff --git a/tests/include/igl/predicates/predicates.cpp b/tests/include/igl/predicates/predicates.cpp index 3bd67f1b7..afd119ead 100644 --- a/tests/include/igl/predicates/predicates.cpp +++ b/tests/include/igl/predicates/predicates.cpp @@ -23,17 +23,17 @@ TEST_CASE("predicates", "[igl][predicates]") { e << 0.0, 1.0; f << 0.5, 0.5; - REQUIRE(orient2d(a, b, c) == Orientation::POSITIVE); - REQUIRE(orient2d(a, c, b) == Orientation::NEGATIVE); - REQUIRE(orient2d(a, b, b) == Orientation::COLLINEAR); - REQUIRE(orient2d(a, a, a) == Orientation::COLLINEAR); - REQUIRE(orient2d(a, b, d) == Orientation::COLLINEAR); - REQUIRE(orient2d(a, f, c) == Orientation::COLLINEAR); + REQUIRE(orient2d(a, b, c) == igl::Orientation::POSITIVE); + REQUIRE(orient2d(a, c, b) == igl::Orientation::NEGATIVE); + REQUIRE(orient2d(a, b, b) == igl::Orientation::COLLINEAR); + REQUIRE(orient2d(a, a, a) == igl::Orientation::COLLINEAR); + REQUIRE(orient2d(a, b, d) == igl::Orientation::COLLINEAR); + REQUIRE(orient2d(a, f, c) == igl::Orientation::COLLINEAR); - REQUIRE(incircle(a,b,c,e) == Orientation::COCIRCULAR); - REQUIRE(incircle(a,b,c,a) == Orientation::COCIRCULAR); - REQUIRE(incircle(a,b,c,d) == Orientation::OUTSIDE); - REQUIRE(incircle(a,b,c,f) == Orientation::INSIDE); + REQUIRE(incircle(a,b,c,e) == igl::Orientation::COCIRCULAR); + REQUIRE(incircle(a,b,c,a) == igl::Orientation::COCIRCULAR); + REQUIRE(incircle(a,b,c,d) == igl::Orientation::OUTSIDE); + REQUIRE(incircle(a,b,c,f) == igl::Orientation::INSIDE); } SECTION("3D") { @@ -46,16 +46,16 @@ TEST_CASE("predicates", "[igl][predicates]") { e << 1.0, 1.0, 1.0; f << std::numeric_limits::epsilon(), 0.0, 0.0; - REQUIRE(orient3d(a, b, c, d) == Orientation::NEGATIVE); - REQUIRE(orient3d(a, b, d, c) == Orientation::POSITIVE); - REQUIRE(orient3d(a, b, d, d) == Orientation::COPLANAR); - REQUIRE(orient3d(a, a, a, a) == Orientation::COPLANAR); - REQUIRE(orient3d(a, b, f, c) == Orientation::COPLANAR); + REQUIRE(orient3d(a, b, c, d) == igl::Orientation::NEGATIVE); + REQUIRE(orient3d(a, b, d, c) == igl::Orientation::POSITIVE); + REQUIRE(orient3d(a, b, d, d) == igl::Orientation::COPLANAR); + REQUIRE(orient3d(a, a, a, a) == igl::Orientation::COPLANAR); + REQUIRE(orient3d(a, b, f, c) == igl::Orientation::COPLANAR); - REQUIRE(insphere(a, b, c, d, e) == Orientation::COSPHERICAL); - REQUIRE(insphere(a, b, d, e, c) == Orientation::COSPHERICAL); - REQUIRE(insphere(b, c, e, d, ((a+b)*0.5).eval()) == Orientation::INSIDE); - REQUIRE(insphere(b, c, e, d, (-f).eval()) == Orientation::OUTSIDE); - REQUIRE(insphere(f, b, d, c, e) == Orientation::INSIDE); + REQUIRE(insphere(a, b, c, d, e) == igl::Orientation::COSPHERICAL); + REQUIRE(insphere(a, b, d, e, c) == igl::Orientation::COSPHERICAL); + REQUIRE(insphere(b, c, e, d, ((a+b)*0.5).eval()) == igl::Orientation::INSIDE); + REQUIRE(insphere(b, c, e, d, (-f).eval()) == igl::Orientation::OUTSIDE); + REQUIRE(insphere(f, b, d, c, e) == igl::Orientation::INSIDE); } } diff --git a/tutorial/1002_EytzingerAABB/main.cpp b/tutorial/1002_EytzingerAABB/main.cpp index 58af4fda3..d3e0c64fa 100644 --- a/tutorial/1002_EytzingerAABB/main.cpp +++ b/tutorial/1002_EytzingerAABB/main.cpp @@ -178,7 +178,7 @@ int main(int argc, char * argv[]) if(offset == OffsetType::VARIABLE_RADIUS_OFFSET) { - viewer.data().set_colors(R); + viewer.data().set_data(R); } viewer.append_mesh(); const int m_index = viewer.selected_data_index; @@ -207,7 +207,7 @@ int main(int argc, char * argv[]) GV *= 2; GV.array() -= 1; Eigen::VectorXd S; - igl::eytzinger_aabb_sdf(GV,primitive,B1,B2,leaf,S); + igl::eytzinger_aabb_sdf(GV,primitive,B1,B2,leaf,S); printf("%-20s: %g secs\n","eytzinger_aabb_sdf",tictoc()); tictoc(); @@ -231,7 +231,7 @@ int main(int argc, char * argv[]) return primitive(p,j); }; double f; - igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f); + igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f); return f - isovalue; }; const std::function diff --git a/tutorial/1002_EytzingerAABB/main.cpp.off b/tutorial/1002_EytzingerAABB/main.cpp.off new file mode 100644 index 000000000..05a564417 --- /dev/null +++ b/tutorial/1002_EytzingerAABB/main.cpp.off @@ -0,0 +1,287 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include + +#include +#define sign(x) ((x>0) - (x<0)) + + +inline double dot2(const Eigen::RowVector3d & v) +{ + return v.dot(v); +} + +inline double clamp(const double & x, const double & a, const double & b) +{ + return std::max(a, std::min(b, x)); +} + +// https://iquilezles.org/articles/triangledistance/ +template +typename Tp::Scalar udTriangle( + const T1 & v1, const T2 & v2, const T3 & v3, const Tp & p) +{ + using vec3 = Eigen::Matrix; + // prepare data + vec3 v21 = v2 - v1; vec3 p1 = p - v1; + vec3 v32 = v3 - v2; vec3 p2 = p - v2; + vec3 v13 = v1 - v3; vec3 p3 = p - v3; + vec3 nor = v21.cross( v13 ); + + return sqrt( // inside/outside test + (sign(v21.cross(nor).dot(p1)) + + sign(v32.cross(nor).dot(p2)) + + sign(v13.cross(nor).dot(p3))<2.0) + ? + // 3 edges + std::min( std::min( + dot2(v21*clamp(v21.dot(p1)/dot2(v21),0.0,1.0)-p1), + dot2(v32*clamp(v32.dot(p2)/dot2(v32),0.0,1.0)-p2) ), + dot2(v13*clamp(v13.dot(p3)/dot2(v13),0.0,1.0)-p3) ) + : + // 1 face + nor.dot(p1)*nor.dot(p1)/dot2(nor) ); +} + + + + +template +typename Tp::Scalar udLineSegment(const Tp & p, const Ta & a, const Tb & b) +{ + using vec3 = Eigen::Matrix; + vec3 pa = p - a; + vec3 ba = b - a; + const auto h = clamp( pa.dot(ba)/ba.dot(ba), 0.0, 1.0 ); + return ( pa - ba*h ).norm(); +} + + + +int main(int argc, char * argv[]) +{ + const int ns = 256+1; + Eigen::RowVector3i res(ns,ns,ns); + + IGL_TICTOC_LAMBDA; + // Read mesh + Eigen::Matrix V; + Eigen::Matrix F; + if(!igl::read_triangle_mesh + (argc>1?argv[1]: TUTORIAL_SHARED_PATH "/decimated-knight.off",V,F)) { + std::cout << "Failed to load mesh." << std::endl; + } + // Normalize + V.rowwise() -= 0.5* (V.colwise().maxCoeff() + V.colwise().minCoeff()); + const double scale = 0.8/V.array().abs().maxCoeff(); + V *= scale; + + // if second arg exists read it to dmat + Eigen::VectorXd R; + if(argc>2) + { + igl::readDMAT( argv[2],R); + R *= scale; + }else + { + //Eigen::VectorXd R = (0.1+0.9*(0.5+(Eigen::VectorXd::Random(V.rows())*0.5).array()))*((2.0/(res(0)-1))*2); + R = V.col(1); + R.array() -= R.minCoeff(); + R /= R.maxCoeff(); + R.array() += 1.0; + R = (R.array().pow(4)).eval(); + R *= 0.01; + } + + + tictoc(); + Eigen::Matrix PB1,PB2; + Eigen::Matrix E; + igl::edges(F,E); + + enum class OffsetType + { + EDGE_OFFSET = 0, + TRIANGLE_OFFSET = 1, + VARIABLE_RADIUS_OFFSET = 2 + } offset = OffsetType::VARIABLE_RADIUS_OFFSET; + + double isovalue = 0.16;(2.0/(res(0)-1))*2; + std::function primitive; + std::vector> data; + switch(offset) + { + case OffsetType::EDGE_OFFSET: + { + igl::box_simplices(V,E,PB1,PB2); + primitive = [&](const Eigen::RowVector3d & p, const int i) + { + return udLineSegment(p, V.row(E(i,0)), V.row(E(i,1))); + }; + break; + } + case OffsetType::TRIANGLE_OFFSET: + { + igl::box_simplices(V,F,PB1,PB2); + primitive = [&](const Eigen::RowVector3d & p, const int i) + { + return udTriangle( V.row(F(i,0)), V.row(F(i,1)), V.row(F(i,2)), p); + }; + break; + } + case OffsetType::VARIABLE_RADIUS_OFFSET: + { + isovalue = 0; + igl::variable_radius_offset(V,F,R,PB1,PB2,data,primitive); + break; + } + } + printf("%-20s: %g secs\n","PB",tictoc()); + + tictoc(); + Eigen::Matrix B1,B2; + Eigen::VectorXi leaf; + igl::eytzinger_aabb(PB1,PB2,B1,B2,leaf); + printf("%-20s: %g secs\n","eytzinger_aabb",tictoc()); + //{ + // Eigen::VectorXi U; + // igl::unique(leaf,U); + // printf("%d → %d\n",PB1.rows(),U.size()-2); + //} + + tictoc(); + Eigen::Matrix mV; + Eigen::Matrix mF; + if(ns <= 64+1) + { + Eigen::Matrix GV; + igl::grid(res,GV); + GV *= 2; + GV.array() -= 1; + Eigen::VectorXd S; + igl::eytzinger_aabb_sdf(GV,primitive,B1,B2,leaf,S); + printf("%-20s: %g secs\n","eytzinger_aabb_sdf",tictoc()); + + tictoc(); + igl::marching_cubes(S,GV,res(0),res(1),res(2),isovalue,mV,mF); + printf("%-20s: %g secs\n","marching_cubes",tictoc()); + } + + tictoc(); + const std::function + sdf = [&](const Eigen::RowVector3d & p) -> double + { + const std::function primitive_p = [&](const int j) + { + return primitive(p,j); + }; + double f; + igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f); + return f; + }; + const std::function + udf = [&](const Eigen::RowVector3d & p) -> double + { + return std::abs(sdf(p)); + //// This made performance worse. + //const std::function primitive_p = [&](const int j) + //{ + // const double d = udTriangle( V.row(F(j,0)), V.row(F(j,1)), V.row(F(j,2)), p); + // const Eigen::RowVector3d r(R(F(j,0)),R(F(j,1)),R(F(j,2))); + // const double min_r = r.minCoeff(); + // const double max_r = r.maxCoeff(); + // if(d > max_r) + // { + // return d - max_r; + // }else if(d < min_r) + // { + // return d - min_r; + // } + // return 0.0; + //}; + //double f; + //igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f); + //return f; + }; + Eigen::RowVector3d origin(-1,-1,-1); + const double h0 = 2; + const int max_depth = floor(log2(res(0))); + Eigen::Matrix ijk; + igl::lipschitz_octree( origin,h0,max_depth,udf,ijk); + printf("%-20s: %g secs\n","lipschitz_octree",tictoc()); + + Eigen::Matrix oV; + Eigen::Matrix oF; + { + tictoc(); + // Gather the corners of those leaf cells + const double h = h0 / (1 << max_depth); + Eigen::Matrix J; + Eigen::Matrix unique_ijk; + Eigen::Matrix unique_corner_positions; + igl::unique_sparse_voxel_corners(origin,h0,max_depth,ijk,unique_ijk,J,unique_corner_positions); + //printf("unique_sparse_voxel_corners: %0.7f seconds\n",tictoc()); + printf("%-20s: %g secs\n","unique_sparse_vo...",tictoc()); + /// Evaluate the signed distance function at the corners + Eigen::VectorXd S(unique_corner_positions.rows()); + //for(int u = 0;u(p,primitive_p,B1,B2,leaf,di); } if(quantity == WINDING_NUMBER || quantity == SIGNED_DISTANCE) { diff --git a/tutorial/1004_SplineDistance/main.cpp b/tutorial/1004_SplineDistance/main.cpp new file mode 100644 index 000000000..28391af0c --- /dev/null +++ b/tutorial/1004_SplineDistance/main.cpp @@ -0,0 +1,186 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char * argv[]) +{ + IGL_TICTOC_LAMBDA; + Eigen::MatrixXd P; + Eigen::MatrixXi C; + igl::readDMAT(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/libigl-cursive-P.dmat",P); + igl::readDMAT(argc>1?argv[2]: TUTORIAL_SHARED_PATH "/libigl-cursive-C.dmat",C); + + // Cheezy conversion of (P,C) to polyline using 50 samples per curve + const int ns = 50; + Eigen::MatrixXd V(ns*C.rows(),P.cols()); + Eigen::MatrixXi E((ns-1)*C.rows(),2); + for(int c = 0;c0) + { + E.row(c*(ns-1) + (i-1))<< c*ns + (i-1), c*ns + i; + } + } + } + + + Eigen::MatrixXd GV; + Eigen::MatrixXi GF; + { + Eigen::RowVector2d min_corner = P.colwise().minCoeff(); + Eigen::RowVector2d max_corner = P.colwise().maxCoeff(); + // samples on x-axis + const int sx = 512; + const int sy = sx * (max_corner(1)-min_corner(1)) / (max_corner(0)-min_corner(0)); + Eigen::RowVector2d diag = (max_corner - min_corner).eval(); + const double bbd = diag.norm(); + diag /= diag.norm(); + min_corner -= 0.1*bbd*diag; + max_corner += 0.1*bbd*diag; + + igl::triangulated_grid(sx,sy,GV,GF); + // Scale and translate grid to fit bounding box + GV.col(0) = (GV.col(0).array() * (max_corner(0)-min_corner(0))) + min_corner(0); + GV.col(1) = (GV.col(1).array() * (max_corner(1)-min_corner(1))) + min_corner(1); + } + + + Eigen::VectorXd sqrD,S; + Eigen::VectorXi I; + Eigen::MatrixXd K; + Eigen::VectorXd W; + + Eigen::Matrix B1,B2; + Eigen::VectorXi leaf; + tictoc(); + igl::cycodebase::spline_eytzinger_aabb(P, C, B1, B2,leaf); + printf("AABB build time: %g secs\n",tictoc()); + igl::writeDMAT("P.dmat", Eigen::MatrixXd(P)); + igl::writeDMAT("C.dmat", Eigen::MatrixXi(C)); + igl::writeDMAT("B1.dmat",Eigen::MatrixXd(B1)); + igl::writeDMAT("B2.dmat",Eigen::MatrixXd(B2)); + + { + tictoc(); + igl::cycodebase::point_spline_squared_distance( + GV, P, C, B1,B2,leaf, sqrD, I, S, K); + printf("sqrD: %d points in %g secs\n",GV.rows(),tictoc()); + tictoc(); + igl::predicates::spline_winding_number(P,C,B1,B2,leaf,GV,W); + printf("Wind: %d points in %g secs\n",GV.rows(),tictoc()); + } + Eigen::VectorXd unsigned_D = sqrD.array().sqrt(); + Eigen::VectorXd signed_D = unsigned_D.array() * (0.5 - W.array().abs()) * 2.0; + + // Just leaves + //Eigen::MatrixXd B1,B2; + //igl::cycodebase::box_cubic(P,C,B1,B2); + + const auto nonempty = igl::find((leaf.array() != -2).eval()); + const int nb = nonempty.size(); + Eigen::MatrixXd BV(nb*4,2); + BV<< + B1(nonempty,0), B1(nonempty,1), + B1(nonempty,0), B2(nonempty,1), + B2(nonempty,0), B2(nonempty,1), + B2(nonempty,0), B1(nonempty,1); + Eigen::MatrixXi BE(nb*4,2); + for(int c = 0;cbool + { + switch(key) + { + default: + return false; + case 'W': + case 'w': + quantity = WINDING_NUMBER; + break; + case 'U': + case 'u': + quantity = UNSIGNED_DISTANCE; + break; + case 'S': + case 's': + quantity = SIGNED_DISTANCE; + break; + } + update(); + return true; + }; + std::cout<