diff --git a/src/mcut/CMakeLists.txt b/src/mcut/CMakeLists.txt index 9024986d7..c3e4dc0ac 100644 --- a/src/mcut/CMakeLists.txt +++ b/src/mcut/CMakeLists.txt @@ -61,7 +61,7 @@ set (CMAKE_DEBUG_POSTFIX "d") list (APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") set(MCUT_MAJOR 1) -set(MCUT_MINOR 1) +set(MCUT_MINOR 2) set(MCUT_PATCH 0) set( MCUT_VERSION "${MCUT_MAJOR}.${MCUT_MINOR}.${MCUT_PATCH}" ) diff --git a/src/mcut/include/mcut/internal/frontend.h b/src/mcut/include/mcut/internal/frontend.h index bc27530e9..ccc2a2bab 100644 --- a/src/mcut/include/mcut/internal/frontend.h +++ b/src/mcut/include/mcut/internal/frontend.h @@ -104,6 +104,12 @@ extern "C" void get_info_impl( McVoid* pMem, McSize* pNumBytes) noexcept(false); +extern "C" void bind_state_impl( + const McContext context, + McFlags stateInfo, + McSize bytes, + const McVoid* pMem); + extern "C" void create_user_event_impl(McEvent* event, McContext context); extern "C" void set_user_event_status_impl(McEvent event, McInt32 execution_status); @@ -142,6 +148,20 @@ extern "C" void dispatch_impl( const McEvent* pEventWaitList, McEvent* pEvent) noexcept(false); +extern "C" void dispatch_planar_section_impl( + McContext context, + McFlags flags, + const McVoid* pSrcMeshVertices, + const uint32_t* pSrcMeshFaceIndices, + const uint32_t* pSrcMeshFaceSizes, + uint32_t numSrcMeshVertices, + uint32_t numSrcMeshFaces, + const McDouble* pNormalVector, + const McDouble sectionOffset, + uint32_t numEventsInWaitlist, + const McEvent* pEventWaitList, + McEvent* pEvent) noexcept(false); + extern "C" void get_connected_components_impl( const McContext context, const McConnectedComponentType connectedComponentType, @@ -233,6 +253,8 @@ struct connected_component_t { std::vector face_adjacent_faces_size_cache; bool face_adjacent_faces_size_cache_initialized = false; #endif // #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) + // non-zero if origin source and cut-mesh where perturbed + vec3 perturbation_vector = vec3(0.0); }; // struct representing a fragment @@ -296,8 +318,32 @@ struct event_t { // API command to be able to effectively wait on user events. std::unique_ptr> m_user_API_command_task_emulator; McContext m_context; - event_t() - : m_user_handle(MC_NULL_HANDLE) + + const char* get_cmd_type_str() + { + switch (m_command_type) { + case McCommandType::MC_COMMAND_DISPATCH: { + return "MC_COMMAND_DISPATCH"; + } break; + case McCommandType::MC_COMMAND_GET_CONNECTED_COMPONENT_DATA: { + return "MC_COMMAND_GET_CONNECTED_COMPONENT_DATA"; + } break; + case McCommandType::MC_COMMAND_GET_CONNECTED_COMPONENTS: { + return "MC_COMMAND_GET_CONNECTED_COMPONENTS"; + } break; + case McCommandType::MC_COMMAND_USER: { + return "MC_COMMAND_USER"; + } break; + case McCommandType::MC_COMMAND_UKNOWN: { + return "MC_COMMAND_UKNOWN"; + } break; + default: + fprintf(stderr, "unknown command type value (%d)\n", (int)m_command_type); + return "UNKNOWN VALUE"; + } + } + explicit event_t(McEvent user_handle, McCommandType command_type) + : m_user_handle(user_handle) , m_responsible_thread_id(UINT32_MAX) , m_runtime_exec_status(MC_NO_ERROR) , m_timestamp_submit(0) @@ -305,11 +351,11 @@ struct event_t { , m_timestamp_end(0) , m_command_exec_status(MC_RESULT_MAX_ENUM) , m_profiling_enabled(true) - , m_command_type(MC_COMMAND_UKNOWN) + , m_command_type(command_type) , m_user_API_command_task_emulator(nullptr) , m_context(nullptr) { - log_msg("[MCUT] Create event " << this); + log_msg("[MCUT] Create event (type=" << get_cmd_type_str() <<", handle=" << m_user_handle << ")"); m_callback_info.m_fn_ptr = nullptr; m_callback_info.m_data_ptr = nullptr; @@ -325,7 +371,7 @@ struct event_t { (*(m_callback_info.m_fn_ptr))(m_user_handle, m_callback_info.m_data_ptr); } - log_msg("[MCUT] Destroy event " << this << "(" << m_user_handle << ")"); + log_msg("[MCUT] Destroy event (type=" << get_cmd_type_str() << ", handle=" << m_user_handle << ")"); } inline std::size_t get_time_since_epoch() @@ -431,6 +477,10 @@ private: // The state and flag variable current used to configure the next dispatch call McFlags m_flags = (McFlags)0; + std::atomic m_general_position_enforcement_constant; + std::atomic m_max_num_perturbation_attempts; + std::atomic m_connected_component_winding_order; + void api_thread_main(uint32_t thread_id) { log_msg("[MCUT] Launch API thread " << std::this_thread::get_id() << " (" << thread_id << ")"); @@ -466,6 +516,9 @@ public: : m_done(false) // , m_joiner(m_api_threads) , m_flags(flags) + , m_general_position_enforcement_constant(1e-4) + , m_max_num_perturbation_attempts(1 << 2), + m_connected_component_winding_order(McConnectedComponentFaceWindingOrder::MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_AS_GIVEN) , m_user_handle(handle) , dbgCallbackBitfieldSource(0) , dbgCallbackBitfieldType(0) @@ -505,6 +558,9 @@ public: void shutdown() { m_done = true; + + std::atomic_thread_fence(std::memory_order_acq_rel); + #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) m_compute_threadpool.reset(); #endif @@ -524,6 +580,41 @@ public: return this->m_flags; } + // returns (user controllable) epsilon representing the maximum by which the cut-mesh + // can be perturbed on any axis + McDouble get_general_position_enforcement_constant() const + { + return this->m_general_position_enforcement_constant.load(std::memory_order_acquire); + } + + void set_general_position_enforcement_constant(McDouble new_value) + { + this->m_general_position_enforcement_constant.store(new_value, std::memory_order_release); + } + + // returns (user controllable) maximum number of times by which the cut-mesh + // can be perturbed before giving (input likely need to be preprocessed) + McUint32 get_general_position_enforcement_attempts() const + { + return this->m_max_num_perturbation_attempts.load(std::memory_order_acquire); + } + + void set_general_position_enforcement_attempts(McUint32 new_value) + { + this->m_max_num_perturbation_attempts.store(new_value, std::memory_order_release); + } + + McConnectedComponentFaceWindingOrder get_connected_component_winding_order() const + { + return this->m_connected_component_winding_order.load(std::memory_order_acquire); + } + + void set_connected_component_winding_order(McConnectedComponentFaceWindingOrder new_value) + { + this->m_connected_component_winding_order.store(new_value, std::memory_order_release); + } + + #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) thread_pool& get_shared_compute_threadpool() { @@ -538,15 +629,17 @@ public: // create the event object associated with the enqueued task // - std::shared_ptr event_ptr = std::shared_ptr(new event_t); + std::shared_ptr event_ptr = std::shared_ptr(new event_t( + reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)), + cmdType)); MCUT_ASSERT(event_ptr != nullptr); g_events.push_front(event_ptr); - event_ptr->m_user_handle = reinterpret_cast(g_objects_counter++); + //event_ptr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); event_ptr->m_profiling_enabled = (this->m_flags & MC_PROFILING_ENABLE) != 0; - event_ptr->m_command_type = cmdType; + // event_ptr->m_command_type = cmdType; event_ptr->log_submit_time(); @@ -567,7 +660,9 @@ public: const std::shared_ptr parent_task_event_ptr = g_events.find_first_if([=](std::shared_ptr e) { return e->m_user_handle == parent_task_event_handle; }); - MCUT_ASSERT(parent_task_event_ptr != nullptr); + if (parent_task_event_ptr == nullptr) { + throw std::invalid_argument("invalid event in waitlist"); + } const bool parent_task_is_not_finished = parent_task_event_ptr->m_finished.load() == false; @@ -588,7 +683,7 @@ public: uint32_t thread_with_empty_queue = UINT32_MAX; for (uint32_t i = 0; i < (uint32_t)m_api_threads.size(); ++i) { - if (m_queues[i].empty() == true) { + if (m_queues[(i + 1) % (uint32_t)m_api_threads.size()].empty() == true) { thread_with_empty_queue = i; break; } @@ -670,7 +765,7 @@ public: // McFlags dispatchFlags = (McFlags)0; - // client/user debugging variable + // client/user debugging variables // ------------------------------ // function pointer to user-define callback function for status/erro reporting diff --git a/src/mcut/include/mcut/internal/utils.h b/src/mcut/include/mcut/internal/utils.h index c4f7d1968..4d8988603 100644 --- a/src/mcut/include/mcut/internal/utils.h +++ b/src/mcut/include/mcut/internal/utils.h @@ -99,8 +99,6 @@ typedef char couldnt_parse_cxx_standard[-1]; ///< Error: couldn't parse standard #define SAFE_ACCESS(var, i) var[i] #endif - - static inline int wrap_integer(int x, const int lo, const int hi) { const int range_size = hi - lo + 1; @@ -247,12 +245,12 @@ pair make_pair(const T a, const T b) // Threadsafe logging to console which prevents std::cerr from mixing strings when // concatenating with the operator<< multiple time per string, across multiple // threads. -#define log_msg(msg_str) \ - { \ - std::stringstream ss; \ - ss << msg_str << std::endl; \ - std::cerr << ss.str(); \ - } +#define log_msg(msg_str) \ + { \ + std::stringstream ss; \ + ss << msg_str << std::endl; \ + std::cerr << ss.str() << std::flush; \ + } // used to marked/label unused function parameters to prevent warnings #define UNUSED(x) [&x] {}() diff --git a/src/mcut/include/mcut/mcut.h b/src/mcut/include/mcut/mcut.h index 7e9c1f2d4..76f68f3fd 100644 --- a/src/mcut/include/mcut/mcut.h +++ b/src/mcut/include/mcut/mcut.h @@ -50,6 +50,12 @@ extern "C" { /** MCUT 1.0 version number */ #define MC_API_VERSION_1_0 MC_MAKE_VERSION(1, 0, 0) // Patch version should always be set to 0 +/** MCUT 1.1 version number */ +#define MC_API_VERSION_1_1 MC_MAKE_VERSION(1, 1, 0) // Patch version should always be set to 0 + +/** MCUT 1.2 version number */ +#define MC_API_VERSION_1_2 MC_MAKE_VERSION(1, 2, 0) // Patch version should always be set to 0 + /** Macro to decode MCUT version (MAJOR) from MC_HEADER_VERSION_COMPLETE */ #define MC_VERSION_MAJOR(version) ((uint32_t)(version) >> 22) @@ -63,7 +69,7 @@ extern "C" { #define MC_HEADER_VERSION 100 /** Complete version of this file */ -#define MC_HEADER_VERSION_COMPLETE MC_MAKE_VERSION(1, 0, MC_HEADER_VERSION) +#define MC_HEADER_VERSION_COMPLETE MC_MAKE_VERSION(1, 2, MC_HEADER_VERSION) /** Constant value assigned to null variables and parameters */ #define MC_NULL_HANDLE 0 @@ -275,7 +281,7 @@ typedef enum McConnectedComponentData { // MC_CONNECTED_COMPONENT_DATA_VERTEX_COUNT = (1 << 0), /**< Number of vertices. */ MC_CONNECTED_COMPONENT_DATA_VERTEX_FLOAT = (1 << 1), /**< List of vertex coordinates as an array of 32 bit floating-point numbers. */ MC_CONNECTED_COMPONENT_DATA_VERTEX_DOUBLE = (1 << 2), /**< List of vertex coordinates as an array of 64 bit floating-point numbers. */ - // MC_CONNECTED_COMPONENT_DATA_FACE_COUNT = (1 << 4), /**< Number of faces. */ + MC_CONNECTED_COMPONENT_DATA_DISPATCH_PERTURBATION_VECTOR = (1 << 4), /**< Tuple of three real numbers represneting the vector that was used to nudge/perturb the source- and cut-mesh into general position to produce the respective connected component. This vector will represent the zero-vector if no perturbation was applied prior to the cutting operation that produced this connected component. */ MC_CONNECTED_COMPONENT_DATA_FACE = (1 << 5), /**< List of faces as an array of indices. Each face can also be understood as a "planar straight line graph" (PSLG), which is a collection of vertices and segments that lie on the same plane. Segments are edges whose endpoints are vertices in the PSLG.*/ MC_CONNECTED_COMPONENT_DATA_FACE_SIZE = (1 << 6), /**< List of face sizes (vertices per face) as an array. */ // MC_CONNECTED_COMPONENT_DATA_EDGE_COUNT = (1 << 7), /**< Number of edges. */ @@ -293,7 +299,7 @@ typedef enum McConnectedComponentData { MC_CONNECTED_COMPONENT_DATA_FACE_ADJACENT_FACE_SIZE = (1 << 18), /**< List of adjacent-face-list sizes (number of adjacent faces per face).*/ MC_CONNECTED_COMPONENT_DATA_FACE_TRIANGULATION = (1 << 19), /**< List of 3*N triangulated face indices, where N is the number of triangles that are produced using a [Constrained] Delaunay triangulation. Such a triangulation is similar to a Delaunay triangulation, but each (non-triangulated) face segment is present as a single edge in the triangulation. A constrained Delaunay triangulation is not truly a Delaunay triangulation. Some of its triangles might not be Delaunay, but they are all constrained Delaunay. */ MC_CONNECTED_COMPONENT_DATA_FACE_TRIANGULATION_MAP = (1 << 20), /**< List of a subset of face indices from one of the input meshes (source-mesh or the cut-mesh). Each value will be the index of an input mesh face. This index-value corresponds to the connected-component face at the accessed index. Example: the value at index 0 of the queried array is the index of the face in the original input mesh. Note that all triangulated-faces are mapped to a defined value. In order to clearly distinguish indices of the cut mesh from those of the source mesh, an input-mesh face index value corresponds to a cut-mesh vertex-index if it is great-than-or-equal-to the number of source-mesh faces. The input connected component (source-mesh or cut-mesh) that is referred to must be one stored internally by MCUT (i.e. a connected component queried from the API via ::McInputOrigin), to ensure consistency with any modification done internally by MCUT. */ - // TODO MC_CONNECTED_COMPONENT_DATA_FACE_TRIANGULATION_CDT = (1<<20) /**< List of 3*N triangulated face indices, where N is the number of triangles that are produced using a [Conforming] Delaunay triangulation. A conforming Delaunay triangulation (CDT) of a PSLG (i.e. a face of in the given connected component) is a true Delaunay triangulation in which each PSLG segment/edge may have been subdivided into several edges by the insertion of additional vertices, called Steiner points. Steiner points are necessary to allow the segments to exist in the mesh while maintaining the Delaunay property. This Delaunay property follows from the definition of a "Delaunay triangulation": the Delaunay triangulation is the triangulation such that, if you circumscribe a circle around every triangle, none of those circles will contain any other points. Alternatively, the Delaunay triangulation is the triangulation that “maximizes the minimum angle in all of the triangles.”. Steiner points are not inserted to meet constraints on the minimum angle and maximum triangle area. MCUT computes a CDT of a given connected component starting only from the faces that have more than three vertices. Neighbouring (triangle) faces will also be processed (i.e. refined) if their incident edges are split as a result of performing CDT on the current face. */ + } McConnectedComponentData; /** @@ -391,20 +397,33 @@ typedef enum McDispatchFlags { MC_DISPATCH_FILTER_SEAM_SRCMESH | // MC_DISPATCH_FILTER_SEAM_CUTMESH), /**< Keep all connected components resulting from the dispatched cut. */ /** - * Allow MCUT to perturb the cut-mesh if the inputs are not in general position. + * The following two flags allow MCUT to perturb the cut-mesh if the inputs are found not to be in general position. * - * MCUT is formulated for inputs in general position. Here the notion of general position is defined with - respect to the orientation predicate (as evaluated on the intersecting polygons). Thus, a set of points - is in general position if no three points are collinear and also no four points are coplanar. + * MCUT is formulated for inputs in general position. Here the notion of general position is defined with respect + * to the orientation predicate (as evaluated on the intersecting polygons). Thus, a set of points is in general + * position if no three points are collinear and also no four points are coplanar. + * + * MCUT uses the "MC_DISPATCH_ENFORCE_GENERAL_POSITION.." flags to inform of when to use perturbation (of the + * cut-mesh) so as to bring the input into general position. In such cases, the idea is to solve the cutting + * problem not on the given input, but on a nearby input. The nearby input is obtained by perturbing the given + * input. The perturbed input will then be in general position and, since it is near the original input, + * the result for the perturbed input will hopefully still be useful. This is justified by the fact that + * the task of MCUT is not to decide whether the input is in general position but rather to make perturbation + * on the input (if) necessary within the available precision of the computing device. + * + * HOW GENERAL POSITION IS ENFORCED + * + * When the inputs are found _not_ to be in GP, MCUT will generate a pseudo-random + * 3d vector "p" representing a translation that will be applied to the cut-mesh. + * + * A component "i" of "p" is computed as "p[i] = r() * c", where "r()" is a + * function returning a random variable from a uniform distribution on the + * interval [-1.0, 1.0), and "c" is a scalar computed from the general position + * enforcement constant (see ::MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT). + * */ - MCUT uses the "GENERAL_POSITION_VIOLATION" flag to inform of when to use perturbation (of the - cut-mesh) so as to bring the input into general position. In such cases, the idea is to solve the cutting - problem not on the given input, but on a nearby input. The nearby input is obtained by perturbing the given - input. The perturbed input will then be in general position and, since it is near the original input, - the result for the perturbed input will hopefully still be useful. This is justified by the fact that - the task of MCUT is not to decide whether the input is in general position but rather to make perturbation - on the input (if) necessary within the available precision of the computing device. */ - MC_DISPATCH_ENFORCE_GENERAL_POSITION = (1 << 15) + MC_DISPATCH_ENFORCE_GENERAL_POSITION = (1 << 15), /**< Enforce general position such that the variable "c" (see detailed note above) is computed as the multiplication of the current general position enforcement constant (of current MCUT context) and the diagonal length of the bounding box of the cut-mesh. So this uses a relative perturbation of the cut-mesh based on its scale (see also ::MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT). */ + MC_DISPATCH_ENFORCE_GENERAL_POSITION_ABSOLUTE= (1 << 16), /**< Enforce general position such that the variable "c" (see detailed note above) is the current general position enforcement constant (of current MCUT context). So this uses an absolute perturbation of the cut-mesh based on the stored constant (see also ::MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT). */ } McDispatchFlags; /** @@ -433,6 +452,17 @@ typedef enum McCommandType { MC_COMMAND_UKNOWN } McCommandType; +/** + * \enum McConnectedComponentFaceWindingOrder + * @brief Flags for specifying the state used to determine winding order of queried faces. + * + * This enum structure defines the flags which are used for identifying the winding order that is used to orient the vertex indices defining the faces of connected components. + */ +typedef enum McConnectedComponentFaceWindingOrder { + MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_AS_GIVEN = 1 << 0, /**< Define the order of face-vertex indices using the orientation implied by the input meshes. */ + MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_REVERSED = 1 << 1, /**< Define the order of face-vertex indices using the reversed orientation implied by the input meshes. */ +} McConnectedComponentFaceWindingOrder; + /** * \enum McQueryFlags * @brief Flags for querying fixed API state. @@ -449,7 +479,10 @@ typedef enum McQueryFlags { MC_EVENT_COMMAND_EXECUTION_STATUS = 1 << 6, /**< the execution status of the command identified by event. See also ::McEventCommandExecStatus */ MC_EVENT_CONTEXT = 1 << 7, /**< The context associated with event. */ MC_EVENT_COMMAND_TYPE = 1 << 8, /**< The command associated with event. Can be one of the values in :: */ - MC_MAX_DEBUG_MESSAGE_LENGTH = 1 << 9 + MC_CONTEXT_MAX_DEBUG_MESSAGE_LENGTH = 1 << 9, /**< The maximum length of a single message return from ::mcGetDebugMessageLog */ + MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT = 1 << 10, /**< A constant small real number representing the amount by which to perturb the cut-mesh when two intersecting polygon are found to not be in general position. */ + MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_ATTEMPTS = 1<<11, /**< The number of times that a dispatch operation will attempt to perturb the cut-mesh if the input meshes are found to not be in general position.*/ + MC_CONTEXT_CONNECTED_COMPONENT_FACE_WINDING_ORDER = 1<<12 /**< The winding order that is used when specifying vertex indices that define the faces of connected components. */ } McQueryFlags; /** @@ -641,7 +674,7 @@ extern MCAPI_ATTR McResult MCAPI_CALL mcDebugMessageCallback( * McResult GetFirstNMessages(McContext context, McUint32 numMsgs) * { * McSize maxMsgLen = 0; - * mcGet(MC_MAX_DEBUG_MESSAGE_LENGTH, &maxMsgLen); + * mcGetInfo(MC_CONTEXT_MAX_DEBUG_MESSAGE_LENGTH, &maxMsgLen); * std::vector msgData(numMsgs * maxMsgLen); * std::vector sources(numMsgs); * std::vector types(numMsgs); @@ -919,7 +952,7 @@ extern MCAPI_ATTR McResult MCAPI_CALL mcSetEventCallback( * -# \p numCutMeshFaces is less than one. * -# \p numEventsInWaitlist Number of events in the waitlist. * -# \p pEventWaitList events that need to complete before this particular command can be executed - * -# ::MC_DISPATCH_ENFORCE_GENERAL_POSITION is not set and: 1) Found two intersecting edges between the source-mesh and the cut-mesh and/or 2) An intersection test between a face and an edge failed because an edge vertex only touches (but does not penetrate) the face, and/or 3) One or more source-mesh vertices are colocated with one or more cut-mesh vertices. + * -# ::MC_DISPATCH_ENFORCE_GENERAL_POSITION or ::MC_DISPATCH_ENFORCE_GENERAL_POSITION_ABSOLUTE is not set and: 1) Found two intersecting edges between the source-mesh and the cut-mesh and/or 2) An intersection test between a face and an edge failed because an edge vertex only touches (but does not penetrate) the face, and/or 3) One or more source-mesh vertices are colocated with one or more cut-mesh vertices. * - ::MC_OUT_OF_MEMORY * -# Insufficient memory to perform operation. */ @@ -957,6 +990,20 @@ extern MCAPI_ATTR McResult MCAPI_CALL mcDispatch( uint32_t numCutMeshVertices, uint32_t numCutMeshFaces); +extern MCAPI_ATTR McResult MCAPI_CALL mcEnqueueDispatchPlanarSection( + const McContext context, + McFlags dispatchFlags, + const McVoid* pSrcMeshVertices, + const uint32_t* pSrcMeshFaceIndices, + const uint32_t* pSrcMeshFaceSizes, + uint32_t numSrcMeshVertices, + uint32_t numSrcMeshFaces, + const McDouble* pNormalVector, + const McDouble sectionOffset, + uint32_t numEventsInWaitlist, + const McEvent* pEventWaitList, + McEvent* pEvent); + /** * @brief Return the value of a selected parameter. * @@ -991,8 +1038,6 @@ extern MCAPI_ATTR McResult MCAPI_CALL mcDispatch( * - MC_INVALID_VALUE * -# \p pContext is NULL or \p pContext is not an existing context. * -# \p bytes is greater than the returned size of data type queried - * - * @note Event synchronisation is not implemented. */ extern MCAPI_ATTR McResult MCAPI_CALL mcGetInfo( const McContext context, @@ -1001,6 +1046,42 @@ extern MCAPI_ATTR McResult MCAPI_CALL mcGetInfo( McVoid* pMem, McSize* pNumBytes); +/** + * @brief Set the value of a selected parameter of a context. + * + * @param[in] context The context handle that was created by a previous call to ::mcCreateContext. + * @param[in] info Information being set. ::McQueryFlags + * @param[in] bytes Size in bytes of memory pointed to by \p pMem. + * @param[out] pMem Pointer to memory from where the appropriate result being copied. + * + * This function effectively sets state variables of a context. All API functions using the respective context shall be effected by this state. + * + * An example of usage: + * @code + * McDouble epsilon = 1e-4; + * McResult err = mcBindState(context, MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT, sizeof(McDouble), &epsilon); + * if(err != MC_NO_ERROR) + * { + * // deal with error + * } + * @endcode + * @return Error code. + * + * Error codes + * - MC_NO_ERROR + * -# proper exit + * - MC_INVALID_VALUE + * -# \p pContext is NULL or \p pContext is not an existing context. + * -# \p stateInfo is not an accepted flag. + * -# \p bytes is 0 + * -# \p pMem is NULL + */ +extern MCAPI_ATTR McResult MCAPI_CALL mcBindState( + const McContext context, + McFlags stateInfo, + McSize bytes, + const McVoid* pMem); + /** * @brief Query the connected components available in a context. * diff --git a/src/mcut/source/frontend.cpp b/src/mcut/source/frontend.cpp index da3ec443a..467e11177 100644 --- a/src/mcut/source/frontend.cpp +++ b/src/mcut/source/frontend.cpp @@ -40,7 +40,7 @@ void create_context_impl(McContext* pOutContext, McFlags flags, uint32_t helperT std::call_once(g_objects_counter_init_flag, []() { g_objects_counter.store(0xDECAF); /*any non-ero value*/ }); - const McContext handle = reinterpret_cast(g_objects_counter++); + const McContext handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); g_contexts.push_front(std::shared_ptr( new context_t(handle, flags #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) @@ -193,7 +193,7 @@ void debug_message_control_impl( // for each possible "source" flag for (auto i : { MC_DEBUG_SOURCE_API, MC_DEBUG_SOURCE_KERNEL }) { if (sourceBitfieldParam & i) { // was it set/included by the user (to be enabled/disabled)? - int n = trailing_zeroes(MC_DEBUG_SOURCE_ALL & i); // get position of bit representing current "source" flag + int n = trailing_zeroes(MC_DEBUG_SOURCE_ALL & i); // get coords of bit representing current "source" flag if (enabled) { // does the user want to enabled this information (from being logged in the debug callback function) context_ptr->dbgCallbackBitfieldSource = set_bit(context_ptr->dbgCallbackBitfieldSource, n); } else { // ... user wants to disable this information @@ -264,12 +264,94 @@ void get_info_impl( } break; } - case MC_MAX_DEBUG_MESSAGE_LENGTH: { - McSize sizeMax = 0; - for (McUint32 i = 0; i < (McUint32)context_ptr->m_debug_logs.size(); ++i) { - sizeMax = std::max((McSize)sizeMax, (McSize)context_ptr->m_debug_logs[i].str.size()); + case MC_CONTEXT_MAX_DEBUG_MESSAGE_LENGTH: { + if (pMem == nullptr) { + *pNumBytes = sizeof(McSize); + } else { + std::lock_guard lock(context_ptr->debugCallbackMutex); + McSize sizeMax = 0; + for (McUint32 i = 0; i < (McUint32)context_ptr->m_debug_logs.size(); ++i) { + sizeMax = std::max((McSize)sizeMax, (McSize)context_ptr->m_debug_logs[i].str.size()); + } + memcpy(pMem, reinterpret_cast(&sizeMax), sizeof(McSize)); } - memcpy(pMem, reinterpret_cast(&sizeMax), sizeof(McSize)); + + } break; + case MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT: { + if (pMem == nullptr) { + *pNumBytes = sizeof(McDouble); + } else { + const McDouble gpec = context_ptr->get_general_position_enforcement_constant(); + memcpy(pMem, reinterpret_cast(&gpec), sizeof(McDouble)); + } + } break; + + case MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_ATTEMPTS: { + if (pMem == nullptr) { + *pNumBytes = sizeof(McUint32); + } else { + const McUint32 attempts = context_ptr->get_general_position_enforcement_attempts(); + memcpy(pMem, reinterpret_cast(&attempts), sizeof(McUint32)); + } + } break; + case MC_CONTEXT_CONNECTED_COMPONENT_FACE_WINDING_ORDER: { + if (pMem == nullptr) { + *pNumBytes = sizeof(McConnectedComponentFaceWindingOrder); + } else { + const McConnectedComponentFaceWindingOrder wo = context_ptr->get_connected_component_winding_order(); + memcpy(pMem, reinterpret_cast(&wo), sizeof(McConnectedComponentFaceWindingOrder)); + } + } break; + + default: + throw std::invalid_argument("unknown info parameter"); + break; + } +} + +void bind_state_impl( + const McContext context, + McFlags stateInfo, + McSize bytes, + const McVoid* pMem) +{ + std::shared_ptr context_ptr = g_contexts.find_first_if([=](const std::shared_ptr cptr) { return cptr->m_user_handle == context; }); + + if (context_ptr == nullptr) { + throw std::invalid_argument("invalid context"); + } + + switch (stateInfo) { + case MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT: { + McDouble value; + memcpy(&value, pMem, bytes); + context_ptr->dbg_cb(MC_DEBUG_SOURCE_API, MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_NOTIFICATION, "general coords enforcement constant set to " + std::to_string(value)); + if (value <= 0) { + throw std::invalid_argument("invalid general coords enforcement constant"); + } + context_ptr->set_general_position_enforcement_constant(value); + + } break; + case MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_ATTEMPTS: { + McUint32 value; + memcpy(&value, pMem, bytes); + context_ptr->dbg_cb(MC_DEBUG_SOURCE_API, MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_NOTIFICATION, "general position enforcement attempts set to " + std::to_string(value)); + if (value < 1) { + throw std::invalid_argument("invalid general coords enforcement attempts -> " + std::to_string(value)); + } + context_ptr->set_general_position_enforcement_attempts(value); + + } break; + case MC_CONTEXT_CONNECTED_COMPONENT_FACE_WINDING_ORDER: { + McConnectedComponentFaceWindingOrder value; + memcpy(&value, pMem, bytes); + context_ptr->dbg_cb(MC_DEBUG_SOURCE_API, MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_NOTIFICATION, "winding-order set to " + std::to_string(value)); + + if (value != McConnectedComponentFaceWindingOrder::MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_AS_GIVEN && // + value != McConnectedComponentFaceWindingOrder::MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_REVERSED) { + throw std::invalid_argument("invalid winding-order param value"); + } + context_ptr->set_connected_component_winding_order(value); } break; default: throw std::invalid_argument("unknown info parameter"); @@ -289,15 +371,17 @@ void create_user_event_impl(McEvent* eventHandle, McContext context) // create the event object associated with the enqueued task // - std::shared_ptr user_event_ptr = std::shared_ptr(new event_t); + std::shared_ptr user_event_ptr = std::shared_ptr(new event_t( + reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)), + McCommandType::MC_COMMAND_USER)); MCUT_ASSERT(user_event_ptr != nullptr); g_events.push_front(user_event_ptr); - user_event_ptr->m_user_handle = reinterpret_cast(g_objects_counter++); + // user_event_ptr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); user_event_ptr->m_profiling_enabled = (context_ptr->get_flags() & MC_PROFILING_ENABLE) != 0; - user_event_ptr->m_command_type = McCommandType::MC_COMMAND_USER; + // user_event_ptr->m_command_type = McCommandType::MC_COMMAND_USER; user_event_ptr->m_context = context; user_event_ptr->m_responsible_thread_id = 0; // initialized but unused for user event @@ -558,6 +642,239 @@ void dispatch_impl( *pEvent = event_handle; } +template +T clamp(const T& n, const T& lower, const T& upper) +{ + return std::max(lower, std::min(n, upper)); +} + +void generate_supertriangle_from_mesh_vertices( + std::vector& supertriangle_vertices, + std::vector& supertriangle_indices, + McFlags dispatchFlags, + const McVoid* pMeshVertices, + uint32_t numMeshVertices, + const McDouble* pNormalVector, + const McDouble sectionOffset, const McDouble eps) +{ + // did the user give us an array of doubles? (otherwise floats) + const bool have_double = (dispatchFlags & MC_DISPATCH_VERTEX_ARRAY_DOUBLE) != 0; + // size (number of bytes) of input floating point type + const std::size_t flt_size = (have_double ? sizeof(double) : sizeof(float)); + // normalized input normal vector + const vec3 n = normalize(vec3(pNormalVector[0], pNormalVector[1], pNormalVector[2])); + + // minimum projection of mesh vertices along the normal vector + double proj_min = 1e10; + // maximum projection of mesh vertices along the normal vector + double proj_max = -proj_min; + + // mesh bbox extents + vec3 bbox_min(1e10); + vec3 bbox_max(-1e10); + + vec3 mean(0.0); + + // mesh-vertex with the most-minimum projection onto the normal vector + vec3 most_min_vertex_pos; + // mesh-vertex with the most-maximum projection onto the normal vector + vec3 most_max_vertex_pos; + + for (uint32_t i = 0; i < numMeshVertices; ++i) { // for each vertex + + // input (raw) pointer in bytes + const McChar* vptr = ((McChar*)pMeshVertices) + (i * flt_size * 3); + vec3 coords; // coordinates of current vertex + + for (uint32_t j = 0; j < 3; ++j) { // for each component + + double coord; + const McChar* const srcptr = vptr + (j * flt_size); + + if (have_double) { + memcpy(&coord, srcptr, flt_size); + } else { + float tmp; + memcpy(&tmp, srcptr, flt_size); + coord = tmp; + } + + bbox_min[j] = std::min(bbox_min[j], coord); + bbox_max[j] = std::max(bbox_max[j], coord); + + coords[j] = coord; + } + mean = mean+coords; + const double dot = dot_product(n, coords); + + if (dot < proj_min) { + most_min_vertex_pos = coords; + proj_min = dot; + } + + if (dot > proj_max) { + most_max_vertex_pos = coords; + proj_max = dot; + } + } + mean = mean/numMeshVertices; + // length of bounding box diagonal + const double bbox_diag = length(bbox_max - bbox_min); + // length from vertex with most-minimum projection to vertex with most-maximum projection + const double max_span = length(most_max_vertex_pos - most_min_vertex_pos); + // parameter indicating distance along the span from vertex with most-minimum projection to vertex with most-maximum projection + const double alpha = clamp(sectionOffset, eps, 1.0 - eps); + // actual distance along the span from vertex with most-minimum projection to vertex with most-maximum projection + const double shiftby = (alpha * max_span); + + const vec3 centroid = most_min_vertex_pos + (n * shiftby); + + // absolute value of the largest component of the normal vector + double max_normal_comp_val_abs = -1e10; + // index of the largest component of the normal vector + uint32_t max_normal_comp_idx = 0; + + for (uint32_t i = 0; i < 3; ++i) { + + const double comp = n[i]; + const double comp_abs = std::fabs(comp); + + if (comp_abs > max_normal_comp_val_abs) { + max_normal_comp_idx = i; + max_normal_comp_val_abs = comp_abs; + } + } + + vec3 w(0.0); + w[max_normal_comp_idx] = 1.0; + if (w == n) { + w[max_normal_comp_idx] = 0.0; + w[(max_normal_comp_idx + 1) % 3] = 1.0; + } + + const vec3 u = cross_product(n, w); + const vec3 v = cross_product(n, u); + + const double mean_dot_n = dot_product(mean, n); + vec3 mean_on_plane = mean - n*mean_dot_n; + + vec3 uv_pos = normalize( (u + v) ); + vec3 uv_neg = normalize( (u - v) ); + vec3 vertex0 = mean_on_plane + uv_pos * ( bbox_diag * 2); + vec3 vertex1 = mean_on_plane + uv_neg * ( bbox_diag * 2); + vec3 vertex2 = mean_on_plane - (normalize(uv_pos + uv_neg) * ( bbox_diag * 2));// supertriangle_origin + (u * bbox_diag * 4); + + supertriangle_vertices.resize(9 * flt_size); + + uint32_t counter = 0; + for (uint32_t i = 0; i < 3; ++i) { + void* dst = supertriangle_vertices.data() + (counter * flt_size); + if (have_double) { + memcpy(dst, &vertex0[i], flt_size); + } else { + float tmp = vertex0[i]; + memcpy(dst, &tmp, flt_size); + } + counter++; + // supertriangle_vertices[counter++] = quad_vertex0[i]; + } + + for (uint32_t i = 0; i < 3; ++i) { + void* dst = supertriangle_vertices.data() + (counter * flt_size); + if (have_double) { + memcpy(dst, &vertex1[i], flt_size); + } else { + float tmp = vertex1[i]; + memcpy(dst, &tmp, flt_size); + } + counter++; + // supertriangle_vertices[counter++] = quad_vertex1[i]; + } + + for (uint32_t i = 0; i < 3; ++i) { + void* dst = supertriangle_vertices.data() + (counter * flt_size); + if (have_double) { + memcpy(dst, &vertex2[i], flt_size); + } else { + float tmp = vertex2[i]; + memcpy(dst, &tmp, flt_size); + } + counter++; + // supertriangle_vertices[counter++] = quad_vertex2[i]; + } + + supertriangle_indices.resize(3); + + supertriangle_indices[0] = 0; + supertriangle_indices[1] = 1; + supertriangle_indices[2] = 2; +} + +void dispatch_planar_section_impl( + McContext context, + McFlags flags, + const McVoid* pSrcMeshVertices, + const uint32_t* pSrcMeshFaceIndices, + const uint32_t* pSrcMeshFaceSizes, + uint32_t numSrcMeshVertices, + uint32_t numSrcMeshFaces, + const McDouble* pNormalVector, + const McDouble sectionOffset, + uint32_t numEventsInWaitlist, + const McEvent* pEventWaitList, + McEvent* pEvent) noexcept(false) +{ + std::shared_ptr context_ptr = g_contexts.find_first_if([=](const std::shared_ptr cptr) { return cptr->m_user_handle == context; }); + + if (context_ptr == nullptr) { + throw std::invalid_argument("invalid context"); + } + + std::weak_ptr context_weak_ptr(context_ptr); + + // submit the dispatch call to be executed asynchronously and return the future + // object that will be waited on as an event + const McEvent event_handle = context_ptr->prepare_and_submit_API_task( + MC_COMMAND_DISPATCH, numEventsInWaitlist, pEventWaitList, + [=]() { + if (!context_weak_ptr.expired()) { + std::shared_ptr context = context_weak_ptr.lock(); + if (context) { + std::vector supertriangle_vertices; + std::vector supertriangle_indices; + + generate_supertriangle_from_mesh_vertices( + supertriangle_vertices, + supertriangle_indices, + flags, + pSrcMeshVertices, + numSrcMeshVertices, + pNormalVector, + sectionOffset, + context->get_general_position_enforcement_constant()); + + preproc( + context, + flags, + pSrcMeshVertices, + pSrcMeshFaceIndices, + pSrcMeshFaceSizes, + numSrcMeshVertices, + numSrcMeshFaces, + (McVoid*)supertriangle_vertices.data(), + (McIndex*)&supertriangle_indices[0], + nullptr, + 3, + 1); + } + } + }); + + MCUT_ASSERT(pEvent != nullptr); + + *pEvent = event_handle; +} + void get_connected_components_impl( const McContext contextHandle, const McConnectedComponentType connectedComponentType, @@ -661,7 +978,7 @@ void triangulate_face( // for each vertex in face: get its coordinates for (uint32_t i = 0; i < cc_face_vcount; ++i) { - + cc_face_vcoords2d.clear(); const vertex_descriptor_t cc_face_vertex_descr = SAFE_ACCESS(cc_face_vertices, i); const vec3& coords = cc.vertex(cc_face_vertex_descr); @@ -872,8 +1189,8 @@ void triangulate_face( // Find the duplicates (if any) const cdt::duplicates_info_t duplicates_info_pre = cdt::find_duplicates( - cc_face_vcoords2d.begin(), - cc_face_vcoords2d.end(), + cc_face_vcoords2d.cbegin(), + cc_face_vcoords2d.cend(), cdt::get_x_coord_vec2d, cdt::get_y_coord_vec2d); @@ -915,205 +1232,240 @@ void triangulate_face( // vector along incident edge, pointing from current to next vertex (NOTE: counter-clockwise dir, normal) const vec2 to_next = next_vtx_coords - perturbed_dvertex_coords; - // positive-value if three points are in CCW order (sign_t::ON_POSITIVE_SIDE) - // negative-value if three points are in CW order (sign_t::ON_NEGATIVE_SIDE) - // zero if collinear (sign_t::ON_ORIENTED_BOUNDARY) - const double orient2d_res = orient2d(perturbed_dvertex_coords, next_vtx_coords, prev_vtx_coords); - const sign_t orient2d_sgn = sign(orient2d_res); - - const double to_prev_sqr_len = squared_length(to_prev); - const double to_next_sqr_len = squared_length(to_next); - // - // Now we must determine which side is the perturbation_vector must be - // pointing. i.e. the side of "perturbed_dvertex_coords" or the side - // of its duplicate + // There is a rare case in which MCUT will produce a CC from complete (not partial)! cut where at least + // one face will be defined by a list of vertices such that this list contains + // two vertices with exactly the same coordinates (due to limitation of floating point precision e.g. after shewchuk predicates in kernel) // - // NOTE: this is only really necessary if the partially cut polygon - // Has more that 3 intersection points (i.e. more than the case of - // one tip, and two duplicates) + // In such a case, we "break the tie" by shifting "perturbed_dvertex_id" halfway + // along the vector running from "perturbed_dvertex_id" to "next_vtx_id", + // if "perturbed_dvertex_id" and "next_vtx_id" are the duplicates. Otherwise, we shift + // "perturbed_dvertex_id" halfway + // along the vector running from "perturbed_dvertex_id" to "prev_vtx_id", + // if instead "perturbed_dvertex_id" and "prev_vtx_id" are the duplicates. // - - const int32_t flip = (orient2d_sgn == sign_t::ON_NEGATIVE_SIDE) ? -1 : 1; - + // It remains possible that next_vtx_id+1 (prev_vtx_id-1) may also be duplicates, in which + // case we should just throw our hands up and bail (can be fixed but too hard and rare to justify the effort) // - // Compute the perturbation vector as the average of the two incident edges eminating - // from the current vertex. NOTE: This perturbation vector should generally point in - // the direction of the polygon-interior (i.e. analogous to pushing the polygon at - // the location represented by perturbed_dvertex_coords) to cause a minute dent due to small - // loss of area. - // Normalization happens below - vec2 perturbation_vector = ((to_prev + to_next) / 2.0) * flip; + // These issue can generally be avoided if the input meshes resemble a uniform triangulation + const bool same_as_prev = ((uint32_t)other_dvertex_id == prev_vtx_id); + const bool same_as_next = (next_vtx_id == (uint32_t)other_dvertex_id); + const bool have_adjacent_duplicates = same_as_prev || same_as_next; - // "orient2d()" is exact in the sense that it can depend on computations with numbers - // whose magnitude is lower than the threshold "orient2d_ccwerrboundA". It follows - // that this threshold is too "small" a number for us to be able to reliably compute - // stuff with the result of "orient2d()" that is near this threshold. - const double errbound = 1e-2; + if (have_adjacent_duplicates) { + // const bool same_as_prev = std::abs(idx_dist_to_prev)==1; + const vec2& shiftby = (same_as_prev) ? to_next : to_prev; + // if(same_as_prev) + //{ + // shiftby = to_next; + // } + // else{ /// then we have "std::abs(idx_dist_to_next) ==1" + // shiftby = to_prev; + // } - // We use "errbound", rather than "orient2d_res", to determine if the incident edges - // are parallel to give us sufficient room of numerical-precision to reliably compute - // the perturbation vector. - // In general, if the incident edges are not parallel then the perturbation vector - // is computed as the mean of "to_prev" and "to_next". Thus, being "too close" - // (within some threshold) to the edges being parallel, can induce unpredicatable - // numerical instabilities, where the mean-vector will be too close to the zero-vector - // and can complicate the task of perturbation. - const bool incident_edges_are_parallel = std::fabs(orient2d_res) <= std::fabs(errbound); + perturbed_dvertex_coords = perturbed_dvertex_coords + (shiftby * 0.5); + } else { // case of partial cut + + // positive-value if three points are in CCW order (sign_t::ON_POSITIVE_SIDE) + // negative-value if three points are in CW order (sign_t::ON_NEGATIVE_SIDE) + // zero if collinear (sign_t::ON_ORIENTED_BOUNDARY) + const double orient2d_res = orient2d(perturbed_dvertex_coords, next_vtx_coords, prev_vtx_coords); + const sign_t orient2d_sgn = sign(orient2d_res); + + const double to_prev_sqr_len = squared_length(to_prev); + const double to_next_sqr_len = squared_length(to_next); - if (incident_edges_are_parallel) { // - // pick the shortest of the two incident edges and compute the - // orthogonal perturbation vector as the counter-clockwise rotation - // of this shortest incident edge. + // Now we must determine which side is the perturbation_vector must be + // pointing. i.e. the side of "perturbed_dvertex_coords" or the side + // of its duplicate + // + // NOTE: this is only really necessary if the partially cut polygon + // Has more that 3 intersection points (i.e. more than the case of + // one tip, and two duplicates) // - // flip sign so that the edge is in the CCW dir by pointing from "prev" to "cur" - vec2 edge_vec(-to_prev.x(), -to_prev.y()); + const int32_t flip = (orient2d_sgn == sign_t::ON_NEGATIVE_SIDE) ? -1 : 1; - if (to_prev_sqr_len > to_next_sqr_len) { - edge_vec = to_next; // pick shortest (NOTE: "to_next" is already in CCW dir) - } + // + // Compute the perturbation vector as the average of the two incident edges eminating + // from the current vertex. NOTE: This perturbation vector should generally point in + // the direction of the polygon-interior (i.e. analogous to pushing the polygon at + // the location represented by perturbed_dvertex_coords) to cause a minute dent due to small + // loss of area. + // Normalization happens below + vec2 perturbation_vector = ((to_prev + to_next) / 2.0) * flip; - // rotate the selected edge by 90 degrees - const vec2 edge_vec_rotated90(-edge_vec.y(), edge_vec.x()); + // "orient2d()" is exact in the sense that it can depend on computations with numbers + // whose magnitude is lower than the threshold "orient2d_ccwerrboundA". It follows + // that this threshold is too "small" a number for us to be able to reliably compute + // stuff with the result of "orient2d()" that is near this threshold. + const double errbound = 1e-2; - perturbation_vector = edge_vec_rotated90; - } + // We use "errbound", rather than "orient2d_res", to determine if the incident edges + // are parallel to give us sufficient room of numerical-precision to reliably compute + // the perturbation vector. + // In general, if the incident edges are not parallel then the perturbation vector + // is computed as the mean of "to_prev" and "to_next". Thus, being "too close" + // (within some threshold) to the edges being parallel, can induce unpredicatable + // numerical instabilities, where the mean-vector will be too close to the zero-vector + // and can complicate the task of perturbation. + const bool incident_edges_are_parallel = std::fabs(orient2d_res) <= std::fabs(errbound); - const vec2 perturbation_dir = normalize(perturbation_vector); + if (incident_edges_are_parallel) { + // + // pick the shortest of the two incident edges and compute the + // orthogonal perturbation vector as the counter-clockwise rotation + // of this shortest incident edge. + // - // - // Compute the maximum length between any two vertices in "cc_face_iter" as the - // largest length between any two vertices. - // - // This will be used to scale "perturbation_dir" so that we find the - // closest edge (from "perturbed_dvertex_coords") that is intersected by this ray. - // We will use the resulting information to determine the amount by-which - // "perturbed_dvertex_coords" is to be perturbed. - // + // flip sign so that the edge is in the CCW dir by pointing from "prev" to "cur" + vec2 edge_vec(-to_prev.x(), -to_prev.y()); - // largest squared length between any two vertices in "cc_face_iter" - double largest_sqrd_length = -1.0; - - for (uint32_t i = 0; i < cc_face_vcount; ++i) { - - const vec2& a = SAFE_ACCESS(cc_face_vcoords2d, i); - - for (uint32_t j = 0; j < cc_face_vcount; ++j) { - - if (i == j) { - continue; // skip -> comparison is redundant + if (to_prev_sqr_len > to_next_sqr_len) { + edge_vec = to_next; // pick shortest (NOTE: "to_next" is already in CCW dir) } - const vec2& b = SAFE_ACCESS(cc_face_vcoords2d, j); + // rotate the selected edge by 90 degrees + const vec2 edge_vec_rotated90(-edge_vec.y(), edge_vec.x()); - const double sqrd_length = squared_length(b - a); - largest_sqrd_length = std::max(sqrd_length, largest_sqrd_length); + perturbation_vector = edge_vec_rotated90; } - } - // - // construct the segment with-which will will find the closest - // intersection point from "perturbed_dvertex_coords" to "perturbed_dvertex_coords + perturbation_dir*std::sqrt(largest_sqrd_length)""; - // + const vec2 perturbation_dir = normalize(perturbation_vector); - const double shift_len = std::sqrt(largest_sqrd_length); - const vec2 shift = perturbation_dir * shift_len; + // + // Compute the maximum length between any two vertices in "cc_face_iter" as the + // largest length between any two vertices. + // + // This will be used to scale "perturbation_dir" so that we find the + // closest edge (from "perturbed_dvertex_coords") that is intersected by this ray. + // We will use the resulting information to determine the amount by-which + // "perturbed_dvertex_coords" is to be perturbed. + // - vec2 intersection_point_on_edge = perturbed_dvertex_coords + shift; // some location potentially outside of polygon + // largest squared length between any two vertices in "cc_face_iter" + double largest_sqrd_length = -1.0; - { - struct { - vec2 start; - vec2 end; - } segment; - segment.start = perturbed_dvertex_coords; - segment.end = perturbed_dvertex_coords + shift; + for (uint32_t i = 0; i < cc_face_vcount; ++i) { - // test segment against all edges to find closest intersection point + const vec2& a = SAFE_ACCESS(cc_face_vcoords2d, i); - double segment_min_tval = 1.0; + for (uint32_t j = 0; j < cc_face_vcount; ++j) { - // for each edge of face to be triangulated (number of vertices == number of edges) - for (std::uint32_t i = 0; i < cc_face_vcount; ++i) { - const std::uint32_t edge_start_idx = i; - const std::uint32_t edge_end_idx = (i + 1) % cc_face_vcount; - - if ((edge_start_idx == (uint32_t)perturbed_dvertex_id || edge_end_idx == (uint32_t)perturbed_dvertex_id) || // - (edge_start_idx == (uint32_t)other_dvertex_id || edge_end_idx == (uint32_t)other_dvertex_id)) { - continue; // impossible to properly intersect incident edges - } - - const vec2& edge_start_coords = SAFE_ACCESS(cc_face_vcoords2d, edge_start_idx); - const vec2& edge_end_coords = SAFE_ACCESS(cc_face_vcoords2d, edge_end_idx); - - double segment_tval; // parameter along segment - double edge_tval; // parameter along current edge - vec2 ipoint; // intersection point between segment and current edge - - const char result = compute_segment_intersection( - segment.start, segment.end, edge_start_coords, edge_end_coords, - ipoint, segment_tval, edge_tval); - - if (result == '1' && segment_min_tval > segment_tval) { // we have an clear intersection point - segment_min_tval = segment_tval; - intersection_point_on_edge = ipoint; - } else if ( - // segment and edge are collinear - result == 'e' || - // segment and edge are collinear, or one entity cuts through the vertex of the other - result == 'v') { - // pick the closest vertex of edge and compute "segment_tval" as a ratio of vector length - - // length from segment start to the start of edge - const double sqr_dist_to_edge_start = squared_length(edge_start_coords - segment.start); - // length from segment start to the end of edge - const double sqr_dist_to_edge_end = squared_length(edge_end_coords - segment.start); - - // length from start of segment to either start of edge or end of edge (depending on which is closer) - double sqr_dist_to_closest = sqr_dist_to_edge_start; - const vec2* ipoint_ptr = &edge_start_coords; - - if (sqr_dist_to_edge_start > sqr_dist_to_edge_end) { - sqr_dist_to_closest = sqr_dist_to_edge_end; - ipoint_ptr = &edge_end_coords; + if (i == j) { + continue; // skip -> comparison is redundant } - // ratio along segment - segment_tval = std::sqrt(sqr_dist_to_closest) / shift_len; + const vec2& b = SAFE_ACCESS(cc_face_vcoords2d, j); - if (segment_min_tval > segment_tval) { + const double sqrd_length = squared_length(b - a); + largest_sqrd_length = std::max(sqrd_length, largest_sqrd_length); + } + } + + // + // construct the segment with-which will will find the closest + // intersection point from "perturbed_dvertex_coords" to "perturbed_dvertex_coords + perturbation_dir*std::sqrt(largest_sqrd_length)""; + // + + const double shift_len = std::sqrt(largest_sqrd_length); + const vec2 shift = perturbation_dir * shift_len; + + vec2 intersection_point_on_edge = perturbed_dvertex_coords + shift; // some location potentially outside of polygon + + { + struct { + vec2 start; + vec2 end; + } segment; + segment.start = perturbed_dvertex_coords; + segment.end = perturbed_dvertex_coords + shift; + + // test segment against all edges to find closest intersection point + + double segment_min_tval = 1.0; + + // for each edge of face to be triangulated (number of vertices == number of edges) + for (std::uint32_t i = 0; i < cc_face_vcount; ++i) { + const std::uint32_t edge_start_idx = i; + const std::uint32_t edge_end_idx = (i + 1) % cc_face_vcount; + + if ((edge_start_idx == (uint32_t)perturbed_dvertex_id || edge_end_idx == (uint32_t)perturbed_dvertex_id) || // + (edge_start_idx == (uint32_t)other_dvertex_id || edge_end_idx == (uint32_t)other_dvertex_id)) { + continue; // impossible to properly intersect incident edges + } + + const vec2& edge_start_coords = SAFE_ACCESS(cc_face_vcoords2d, edge_start_idx); + const vec2& edge_end_coords = SAFE_ACCESS(cc_face_vcoords2d, edge_end_idx); + + double segment_tval; // parameter along segment + double edge_tval; // parameter along current edge + vec2 ipoint; // intersection point between segment and current edge + + const char result = compute_segment_intersection( + segment.start, segment.end, edge_start_coords, edge_end_coords, + ipoint, segment_tval, edge_tval); + + if (result == '1' && segment_min_tval > segment_tval) { // we have an clear intersection point segment_min_tval = segment_tval; - intersection_point_on_edge = *ipoint_ptr; // closest point + intersection_point_on_edge = ipoint; + } else if ( + // segment and edge are collinear + result == 'e' || + // segment and edge are collinear, or one entity cuts through the vertex of the other + result == 'v') { + // pick the closest vertex of edge and compute "segment_tval" as a ratio of vector length + + // length from segment start to the start of edge + const double sqr_dist_to_edge_start = squared_length(edge_start_coords - segment.start); + // length from segment start to the end of edge + const double sqr_dist_to_edge_end = squared_length(edge_end_coords - segment.start); + + // length from start of segment to either start of edge or end of edge (depending on which is closer) + double sqr_dist_to_closest = sqr_dist_to_edge_start; + const vec2* ipoint_ptr = &edge_start_coords; + + if (sqr_dist_to_edge_start > sqr_dist_to_edge_end) { + sqr_dist_to_closest = sqr_dist_to_edge_end; + ipoint_ptr = &edge_end_coords; + } + + // ratio along segment + segment_tval = std::sqrt(sqr_dist_to_closest) / shift_len; + + if (segment_min_tval > segment_tval) { + segment_min_tval = segment_tval; + intersection_point_on_edge = *ipoint_ptr; // closest point + } } } + + MCUT_ASSERT(segment_min_tval <= 1.0); // ... because we started from max length between any two vertices } - MCUT_ASSERT(segment_min_tval <= 1.0); // ... because we started from max length between any two vertices - } + // Shortened perturbation vector: shortening from the vector that is as long as the + // max length between any two vertices in "cc_face_iter", to a vector that runs + // from "perturbed_dvertex_coords" and upto the boundary-point of the "cc_face_iter", along + // "perturbation_vector" and passing through the interior of "cc_face_iter") + const vec2 revised_perturbation_vector = (intersection_point_on_edge - perturbed_dvertex_coords); + const double revised_perturbation_len = length(revised_perturbation_vector); - // Shortened perturbation vector: shortening from the vector that is as long as the - // max length between any two vertices in "cc_face_iter", to a vector that runs - // from "perturbed_dvertex_coords" and upto the boundary-point of the "cc_face_iter", along - // "perturbation_vector" and passing through the interior of "cc_face_iter") - const vec2 revised_perturbation_vector = (intersection_point_on_edge - perturbed_dvertex_coords); - const double revised_perturbation_len = length(revised_perturbation_vector); + const double scale = (errbound * revised_perturbation_len); + // The translation by which we perturb "perturbed_dvertex_coords" + // + // NOTE: since "perturbation_vector" was constructed from "to_prev" and "to_next", + // "displacement" is by-default pointing in the positive/CCW direction, which is torward + // the interior of the polygon represented by "cc_face_iter". + // Thus, the cases with "orient2d_sgn == sign_t::ON_POSITIVE_SIDE" and + // "orient2d_sgn == sign_t::ON_ORIENTED_BOUNDARY", result in the same displacement vector + const vec2 displacement = (perturbation_dir * scale); - const double scale = (errbound * revised_perturbation_len); - // The translation by which we perturb "perturbed_dvertex_coords" - // - // NOTE: since "perturbation_vector" was constructed from "to_prev" and "to_next", - // "displacement" is by-default pointing in the positive/CCW direction, which is torward - // the interior of the polygon represented by "cc_face_iter". - // Thus, the cases with "orient2d_sgn == sign_t::ON_POSITIVE_SIDE" and - // "orient2d_sgn == sign_t::ON_ORIENTED_BOUNDARY", result in the same displacement vector - const vec2 displacement = (perturbation_dir * scale); + // perturb + perturbed_dvertex_coords = perturbed_dvertex_coords + displacement; - // perturb - perturbed_dvertex_coords = perturbed_dvertex_coords + displacement; - - //} // for (std::uint32_t dv_iter = 0; dv_iter < 2; ++dv_iter) { + //} // for (std::uint32_t dv_iter = 0; dv_iter < 2; ++dv_iter) { + } // if(have_adjacent_duplicates) } // for (std::vector::const_iterator duplicate_vpair_iter = duplicates_info_pre.duplicates.cbegin(); ... } // if (have_duplicates) { @@ -1126,8 +1478,8 @@ void triangulate_face( // check for duplicate vertices again const cdt::duplicates_info_t duplicates_info_post = cdt::find_duplicates( - cc_face_vcoords2d.begin(), - cc_face_vcoords2d.end(), + cc_face_vcoords2d.cbegin(), + cc_face_vcoords2d.cend(), cdt::get_x_coord_vec2d, cdt::get_y_coord_vec2d); @@ -1628,6 +1980,23 @@ void get_connected_component_data_impl_detail( #endif // #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) } } break; + case MC_CONNECTED_COMPONENT_DATA_DISPATCH_PERTURBATION_VECTOR: { + SCOPED_TIMER("MC_CONNECTED_COMPONENT_DATA_DISPATCH_PERTURBATION_VECTOR"); + if (pMem == nullptr) { + *pNumBytes = sizeof(vec3); + } else { + if (bytes > sizeof(vec3)) { + throw std::invalid_argument("out of bounds memory access"); + } + if (bytes % sizeof(vec3) != 0) { + throw std::invalid_argument("invalid number of bytes"); + } + + ((McDouble*)pMem)[0] = cc_uptr->perturbation_vector[0]; + ((McDouble*)pMem)[1] = cc_uptr->perturbation_vector[1]; + ((McDouble*)pMem)[2] = cc_uptr->perturbation_vector[2]; + } + } break; case MC_CONNECTED_COMPONENT_DATA_FACE: { SCOPED_TIMER("MC_CONNECTED_COMPONENT_DATA_FACE"); if (pMem == nullptr) { // querying for number of bytes @@ -1739,12 +2108,13 @@ void get_connected_component_data_impl_detail( std::vector partial_sum_vec = cc_uptr->face_sizes_cache; // copy parallel_partial_sum(context_ptr->get_shared_compute_threadpool(), partial_sum_vec.begin(), partial_sum_vec.end()); + const bool flip_winding_order = context_ptr->get_connected_component_winding_order() == McConnectedComponentFaceWindingOrder::MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_REVERSED; // // step 3 // - auto fn_face_indices_copy = [&cc_uptr, &partial_sum_vec, &casted_ptr, &num_indices_to_copy](face_array_iterator_t block_start_, face_array_iterator_t block_end_) { + auto fn_face_indices_copy = [flip_winding_order, &cc_uptr, &partial_sum_vec, &casted_ptr, &num_indices_to_copy](face_array_iterator_t block_start_, face_array_iterator_t block_end_) { const uint32_t base_face_offset = std::distance(cc_uptr->kernel_hmesh_data->mesh->faces_begin(), block_start_); MCUT_ASSERT(base_face_offset < (uint32_t)cc_uptr->face_sizes_cache.size()); @@ -1767,14 +2137,27 @@ void get_connected_component_data_impl_detail( MCUT_ASSERT(num_vertices_around_face >= 3u); - // for each vertex in face - for (uint32_t i = 0; i < num_vertices_around_face; ++i) { - const uint32_t vertex_idx = (uint32_t)SAFE_ACCESS(vertices_around_face, i); - *(casted_ptr + index_arr_offset) = vertex_idx; - ++index_arr_offset; + if (flip_winding_order) { + // for each vertex in face + for (int32_t i = (num_vertices_around_face - 1); i >= (int32_t)0; --i) { + const uint32_t vertex_idx = (uint32_t)SAFE_ACCESS(vertices_around_face, i); + *(casted_ptr + index_arr_offset) = vertex_idx; + ++index_arr_offset; - if (index_arr_offset == num_indices_to_copy) { - break; + if (index_arr_offset == num_indices_to_copy) { + break; + } + } + } else { + // for each vertex in face + for (uint32_t i = 0; i < num_vertices_around_face; ++i) { + const uint32_t vertex_idx = (uint32_t)SAFE_ACCESS(vertices_around_face, i); + *(casted_ptr + index_arr_offset) = vertex_idx; + ++index_arr_offset; + + if (index_arr_offset == num_indices_to_copy) { + break; + } } } } @@ -1792,6 +2175,8 @@ void get_connected_component_data_impl_detail( std::vector cc_face_vertices; uint32_t elem_offset = 0; + const bool flip_winding_order = context_ptr->get_connected_component_winding_order() == McConnectedComponentFaceWindingOrder::MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_REVERSED; + for (face_array_iterator_t fiter = cc_uptr->kernel_hmesh_data->mesh->faces_begin(); fiter != cc_uptr->kernel_hmesh_data->mesh->faces_end(); ++fiter) { cc_face_vertices.clear(); @@ -1800,13 +2185,25 @@ void get_connected_component_data_impl_detail( MCUT_ASSERT(num_vertices_around_face >= 3u); - for (uint32_t i = 0; i < num_vertices_around_face; ++i) { - const uint32_t vertex_idx = (uint32_t)cc_face_vertices[i]; - *(casted_ptr + elem_offset) = vertex_idx; - ++elem_offset; + if (flip_winding_order) { + for (int32_t i = (int32_t)(num_vertices_around_face - 1); i >= 0; --i) { + const uint32_t vertex_idx = (uint32_t)cc_face_vertices[i]; + *(casted_ptr + elem_offset) = vertex_idx; + ++elem_offset; - if (elem_offset == num_indices_to_copy) { - break; + if (elem_offset == num_indices_to_copy) { + break; + } + } + } else { + for (uint32_t i = 0; i < num_vertices_around_face; ++i) { + const uint32_t vertex_idx = (uint32_t)cc_face_vertices[i]; + *(casted_ptr + elem_offset) = vertex_idx; + ++elem_offset; + + if (elem_offset == num_indices_to_copy) { + break; + } } } } @@ -2697,9 +3094,13 @@ void get_connected_component_data_impl_detail( triangulate_face(cc_face_triangulation, context_ptr, cc_face_vcount, cc_face_vertices, *(cc.get()), *cc_face_iter); + // NOTE: "cc_face_triangulation" can be empty if the face has near-zero area + for (uint32_t i = 0; i < (uint32_t)cc_face_triangulation.size(); ++i) { const uint32_t local_idx = cc_face_triangulation[i]; // id local within the current face that we are triangulating + MCUT_ASSERT(local_idx < cc_face_vcount); const uint32_t global_idx = (uint32_t)cc_face_vertices[local_idx]; + MCUT_ASSERT(global_idx < (uint32_t)cc->number_of_vertices()); cdt_index_cache_local.push_back(global_idx); @@ -2707,9 +3108,11 @@ void get_connected_component_data_impl_detail( if ((i % 3) == 0) { // every three indices constitute one triangle // map every CDT triangle in "*cc_face_iter" to the index value of "*cc_face_iter" (in the user input mesh) const uint32_t internal_inputmesh_face_idx = (uint32_t)cc_uptr->kernel_hmesh_data->data_maps.face_map[(uint32_t)*cc_face_iter]; + MCUT_ASSERT(internal_inputmesh_face_idx < cc_face_count); const uint32_t user_inputmesh_face_idx = map_internal_inputmesh_face_idx_to_user_inputmesh_face_idx( internal_inputmesh_face_idx, cc_uptr); + MCUT_ASSERT(internal_inputmesh_face_idx < cc_face_count); cdt_face_map_cache_local.push_back(user_inputmesh_face_idx); } } @@ -2815,28 +3218,30 @@ void get_connected_component_data_impl_detail( triangulate_face(cc_face_triangulation, context_ptr, cc_face_vcount, cc_face_vertices, cc.get()[0], *cc_face_iter); - // - // Change local triangle indices to global index values (in CC) and save - // + if (cc_face_triangulation.empty() == false) { + // + // Change local triangle indices to global index values (in CC) and save + // - const uint32_t cc_face_triangulation_index_count = (uint32_t)cc_face_triangulation.size(); - cc_uptr->cdt_index_cache.reserve( - cc_uptr->cdt_index_cache.size() + cc_face_triangulation_index_count); + const uint32_t cc_face_triangulation_index_count = (uint32_t)cc_face_triangulation.size(); + cc_uptr->cdt_index_cache.reserve( + cc_uptr->cdt_index_cache.size() + cc_face_triangulation_index_count); - for (uint32_t i = 0; i < cc_face_triangulation_index_count; ++i) { - const uint32_t local_idx = cc_face_triangulation[i]; // id local within the current face that we are triangulating - const uint32_t global_idx = (uint32_t)cc_face_vertices[local_idx]; + for (uint32_t i = 0; i < cc_face_triangulation_index_count; ++i) { + const uint32_t local_idx = cc_face_triangulation[i]; // id local within the current face that we are triangulating + const uint32_t global_idx = (uint32_t)cc_face_vertices[local_idx]; - cc_uptr->cdt_index_cache.push_back(global_idx); + cc_uptr->cdt_index_cache.push_back(global_idx); - if (user_requested_cdt_face_maps) { - if ((i % 3) == 0) { // every three indices constitute one triangle - // map every CDT triangle in "*cc_face_iter" to the index value of "*cc_face_iter" (in the user input mesh) - const uint32_t internal_inputmesh_face_idx = (uint32_t)cc_uptr->kernel_hmesh_data->data_maps.face_map[(uint32_t)*cc_face_iter]; - const uint32_t user_inputmesh_face_idx = map_internal_inputmesh_face_idx_to_user_inputmesh_face_idx( - internal_inputmesh_face_idx, - cc_uptr); - cc_uptr->cdt_face_map_cache.push_back(user_inputmesh_face_idx); + if (user_requested_cdt_face_maps) { + if ((i % 3) == 0) { // every three indices constitute one triangle + // map every CDT triangle in "*cc_face_iter" to the index value of "*cc_face_iter" (in the user input mesh) + const uint32_t internal_inputmesh_face_idx = (uint32_t)cc_uptr->kernel_hmesh_data->data_maps.face_map[(uint32_t)*cc_face_iter]; + const uint32_t user_inputmesh_face_idx = map_internal_inputmesh_face_idx_to_user_inputmesh_face_idx( + internal_inputmesh_face_idx, + cc_uptr); + cc_uptr->cdt_face_map_cache.push_back(user_inputmesh_face_idx); + } } } } @@ -2886,7 +3291,16 @@ void get_connected_component_data_impl_detail( throw std::invalid_argument("invalid number of bytes"); } - memcpy(pMem, reinterpret_cast(cc_uptr->cdt_index_cache.data()), bytes); + const bool flip_winding_order = context_ptr->get_connected_component_winding_order() == McConnectedComponentFaceWindingOrder::MC_CONNECTED_COMPONENT_FACE_WINDING_ORDER_REVERSED; + + if (flip_winding_order) { + const uint32_t n = (uint32_t)cc_uptr->cdt_index_cache.size(); + for (uint32_t i = 0; i < n; ++i) { + ((uint32_t*)pMem)[n - 1 - i] = cc_uptr->cdt_index_cache[i]; + } + } else { + memcpy(pMem, reinterpret_cast(cc_uptr->cdt_index_cache.data()), bytes); + } } } } break; diff --git a/src/mcut/source/kernel.cpp b/src/mcut/source/kernel.cpp index 8199053e1..4ec72329f 100644 --- a/src/mcut/source/kernel.cpp +++ b/src/mcut/source/kernel.cpp @@ -33,8 +33,8 @@ #include "mcut/internal/hmesh.h" #include "mcut/internal/kernel.h" #include "mcut/internal/math.h" -#include "mcut/internal/utils.h" #include "mcut/internal/timer.h" +#include "mcut/internal/utils.h" #ifndef LICENSE_PURCHASED #define lmsg() printf("NOTE: MCUT is copyrighted and may not be sold or included in commercial products without a license.\n") @@ -656,8 +656,15 @@ hmesh_t extract_connected_components( // structure "mesh" to the (local) connected-component // +//#define EXTRACT_SEAM_HALFEDGES + +#ifdef EXTRACT_SEAM_HALFEDGES + std::vector cc_seam_halfedges; + const std::vector& halfedges_on_face = mesh.get_halfedges_around_face(fd); + bool prev_vertex_belonged_to_seam = false; // previous in face + #endif // for each vertex around the current face - const std::vector vertices_around_face = mesh.get_vertices_around_face(fd); + const std::vector vertices_around_face = mesh.get_vertices_around_face(fd); // order according to "halfedges_on_face" (targets) for (std::vector::const_iterator face_vertex_iter = vertices_around_face.cbegin(); face_vertex_iter != vertices_around_face.cend(); @@ -683,9 +690,31 @@ hmesh_t extract_connected_components( // check if we need to save vertex as being a seam vertex // std::vector::const_iterator fiter = mesh_vertex_to_seam_flag.find(*face_vertex_iter); bool is_seam_vertex = (size_t)(*face_vertex_iter) < mesh_vertex_to_seam_flag.size() && SAFE_ACCESS(mesh_vertex_to_seam_flag, *face_vertex_iter); //(size_t)(*face_vertex_iter) < mesh_vertex_to_seam_flag.size(); //fiter != mesh_vertex_to_seam_flag.cend() && fiter->second == true; + if (is_seam_vertex) { + cc_seam_vertices.push_back(cc_descriptor); +#ifdef EXTRACT_SEAM_HALFEDGES + const uint32_t face_vertex_idx = std::distance(vertices_around_face.cbegin(), face_vertex_iter); + + const bool is_first_face_vertex = (face_vertex_idx == 0); + bool have_seam_halfedge = prev_vertex_belonged_to_seam; + + if (is_first_face_vertex) { + vd_t last_vtx_descr = (*(vertices_around_face.end() - 1)); + bool last_vertex_is_seam_vertex = (size_t)(last_vtx_descr) < mesh_vertex_to_seam_flag.size() && SAFE_ACCESS(mesh_vertex_to_seam_flag, last_vtx_descr); //(size_t)(*face_vertex_iter) < mesh_vertex_to_seam_flag.size(); //fiter != mesh_vertex_to_seam_flag.cend() && fiter->second == true; + have_seam_halfedge = (last_vertex_is_seam_vertex); + } + + if (have_seam_halfedge) { + const halfedge_descriptor_t seam_he = SAFE_ACCESS(halfedges_on_face, face_vertex_idx); // number of halfedge == number of vertices in face + cc_seam_halfedges.push_back(seam_he); + } + #endif } +#ifdef EXTRACT_SEAM_HALFEDGES + prev_vertex_belonged_to_seam = is_seam_vertex; + #endif } } } // for (face_array_iterator_t face_iter = mesh.faces_begin(); face_iter != mesh.faces_end(); ++face_iter) @@ -1453,13 +1482,13 @@ void update_neighouring_ps_iface_m0_edge_list( typedef std::vector traced_polygon_t; bool mesh_is_closed( -#if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) +#if 0 //defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) thread_pool& scheduler, #endif const hmesh_t& mesh) { bool all_halfedges_incident_to_face = true; -#if 0 //defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) +#if 0 // defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) { printf("mesh=%d\n", (int)mesh.number_of_halfedges()); all_halfedges_incident_to_face = parallel_find_if( @@ -1556,7 +1585,7 @@ void dispatch(output_t& output, const input_t& input) TIMESTACK_PUSH("Check source mesh is closed"); const bool sm_is_watertight = mesh_is_closed( -#if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) +#if 0 //defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) *input.scheduler, #endif sm); @@ -1565,7 +1594,7 @@ void dispatch(output_t& output, const input_t& input) TIMESTACK_PUSH("Check cut mesh is closed"); const bool cm_is_watertight = mesh_is_closed( -#if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) +#if 0// defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) *input.scheduler, #endif cs); @@ -1763,7 +1792,7 @@ void dispatch(output_t& output, const input_t& input) /////////////////////////////////////////////////////////////////////////// std::unordered_map> ps_edge_face_intersection_pairs; - + TIMESTACK_PUSH("Prepare edge-to-face pairs"); #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) @@ -1971,7 +2000,7 @@ void dispatch(output_t& output, const input_t& input) // // build bounding boxes for each intersecting edge // - +#if 1 TIMESTACK_PUSH("Build edge bounding boxes"); // http://gamma.cs.unc.edu/RTRI/i3d08_RTRI.pdf @@ -2146,7 +2175,7 @@ void dispatch(output_t& output, const input_t& input) // assuming each edge will produce a new vertex m0.reserve_for_additional_elements((std::uint32_t)ps_edge_face_intersection_pairs.size()); - +#endif TIMESTACK_PUSH("Compute intersecting face properties"); // compute/extract geometry properties of each tested face //-------------------------------------------------------- @@ -2167,6 +2196,8 @@ void dispatch(output_t& output, const input_t& input) OutputStorageTypesTuple; typedef std::map>::const_iterator InputStorageIteratorType; + std::atomic potentially_intersecting_face_with_zero_area(-1); // did any errors occur (e.g. found a face with zero area) + auto fn_compute_intersecting_face_properties = [&](InputStorageIteratorType block_start_, InputStorageIteratorType block_end_) -> OutputStorageTypesTuple { OutputStorageTypesTuple output_res; std::unordered_map& ps_tested_face_to_plane_normal_LOCAL = std::get<0>(output_res); @@ -2196,6 +2227,10 @@ void dispatch(output_t& output, const input_t& input) tested_face_plane_param_d, tested_face_vertices.data(), (int)tested_face_vertices.size()); + + if (squared_length(tested_face_plane_normal) == 0) { + potentially_intersecting_face_with_zero_area.store((int)tested_faces_iter->first, std::memory_order_release); + } } return output_res; }; @@ -2226,6 +2261,10 @@ void dispatch(output_t& output, const input_t& input) OutputStorageTypesTuple future_res = f.get(); + if (potentially_intersecting_face_with_zero_area.load(std::memory_order_acquire) >= 0) { + break; // stop there was a runtime error + } + std::unordered_map& ps_tested_face_to_plane_normal_FUTURE = std::get<0>(future_res); std::unordered_map& ps_tested_face_to_plane_normal_d_param_FUTURE = std::get<1>(future_res); std::unordered_map& ps_tested_face_to_plane_normal_max_comp_FUTURE = std::get<2>(future_res); @@ -2248,6 +2287,19 @@ void dispatch(output_t& output, const input_t& input) ps_tested_face_to_vertices_FUTURE.cend()); } + const int tmp_local = potentially_intersecting_face_with_zero_area.load(std::memory_order_acquire); + + if (tmp_local >= 0) { + const bool is_cutmesh_face = (tmp_local > sm_face_count); + // if "tmp_local" > srcMeshFaceCount then "tmp_local" is a cut-mesh face with id="tmp_local-srcMeshFaceCount" + const std::string msh_name = is_cutmesh_face ? "cut-mesh" : "source-mesh"; + // index/descriptor in the _kernel_ input mesh (note the stress on kernel since frontend might modify user-provided mesh) + const fd_t bad_face_desr = fd_t(is_cutmesh_face ? (tmp_local - sm_face_count) : tmp_local); + lg.set_reason_for_failure("face f" + std::to_string(bad_face_desr) + " of " + msh_name + " is degenerate (has zero area)"); + output.status.store(is_cutmesh_face ? status_t::INVALID_CUT_MESH : status_t::INVALID_SRC_MESH, std::memory_order_release); + return; // stop there was a runtime error + } + } // end of parallel scope #else // for each face that is to be tested for intersection @@ -2277,6 +2329,16 @@ void dispatch(output_t& output, const input_t& input) tested_face_plane_param_d, tested_face_vertices.data(), (int)tested_face_vertices.size()); + + if (squared_length(tested_face_plane_normal) == 0) { + const int tmp_local = (int)tested_faces_iter->first; + const bool is_cutmesh_face = (tmp_local > sm_face_count); + const std::string msh_name = is_cutmesh_face ? "cut-mesh" : "source-mesh"; + const fd_t bad_face_desr = fd_t(is_cutmesh_face ? (tmp_local - sm_face_count) : tmp_local); + lg.set_reason_for_failure("face f" + std::to_string(bad_face_desr) + " of " + msh_name + " is degenerate (has zero area)"); + output.status = (is_cutmesh_face ? status_t::INVALID_CUT_MESH : status_t::INVALID_SRC_MESH); + return; + } } } #endif @@ -3258,7 +3320,7 @@ void dispatch(output_t& output, const input_t& input) std::vector // list of halfedges whose target is the intersection point > ivtx_to_incoming_hlist; -#if 1 // used for debugging colinearity bug, which occur when we have poly with eg. > 3 +#if 0 // used for debugging colinearity bug, which occur when we have poly with eg. > 3 // vertices where at least 3 more-or-less are colinear but exact predicate says no. for (std::map, std::vector>::const_iterator cutpath_edge_creation_info_iter = cutpath_edge_creation_info.cbegin(); cutpath_edge_creation_info_iter != cutpath_edge_creation_info.cend(); @@ -7704,8 +7766,8 @@ void dispatch(output_t& output, const input_t& input) // if (proceed_to_fill_holes == false) { - - return; // exit + printf("[mcut-kernel]: detected a configuration that does not permit filling holes. input-mesh verification advised.\n"); + return; // done } if (false == (input.keep_fragments_below_cutmesh || // diff --git a/src/mcut/source/math.cpp b/src/mcut/source/math.cpp index 8fb42268a..cdb5c70ec 100644 --- a/src/mcut/source/math.cpp +++ b/src/mcut/source/math.cpp @@ -589,9 +589,9 @@ return x; }(); - matrix_t I(3, 3); // 3x3 identity + matrix_t I(3, 3); // 3x3 identity with reflection I(0, 0) = 1.0; - I(1, 1) = 1.0; + I(1, 1) = -1.0; I(2, 2) = 1.0; matrix_t R = I; diff --git a/src/mcut/source/mcut.cpp b/src/mcut/source/mcut.cpp index 8937a0f94..f3b5c1dee 100644 --- a/src/mcut/source/mcut.cpp +++ b/src/mcut/source/mcut.cpp @@ -123,7 +123,7 @@ MCAPI_ATTR McResult MCAPI_CALL mcGetDebugMessageLog( per_thread_api_log_str = "count must be > 0"; } else if (bufSize == 0) { per_thread_api_log_str = "bufSize must be > 0"; - } else if (numFetched == nullptr) { + } else if (numFetched == nullptr) { per_thread_api_log_str = "numFetched undef (NULL)"; } else { try { @@ -198,7 +198,11 @@ MCAPI_ATTR McResult MCAPI_CALL mcGetInfo(const McContext context, McFlags info, per_thread_api_log_str = "context ptr (param0) undef (NULL)"; } else if (bytes != 0 && pMem == nullptr) { per_thread_api_log_str = "invalid specification (param2 & param3)"; - } else if (false == (info == MC_CONTEXT_FLAGS || info == MC_MAX_DEBUG_MESSAGE_LENGTH)) // check all possible values + } else if (false == // + (info == MC_CONTEXT_FLAGS || // + info == MC_CONTEXT_MAX_DEBUG_MESSAGE_LENGTH || // + info == MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT || // + info == MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_ATTEMPTS)) // check all possible values { per_thread_api_log_str = "invalid info flag val (param1)"; } else if ((info == MC_CONTEXT_FLAGS) && (pMem != nullptr && bytes != sizeof(McFlags))) { @@ -221,6 +225,50 @@ MCAPI_ATTR McResult MCAPI_CALL mcGetInfo(const McContext context, McFlags info, return return_value; } +MCAPI_ATTR McResult MCAPI_CALL mcBindState( + const McContext context, + McFlags stateInfo, + McSize bytes, + const McVoid* pMem) +{ + McResult return_value = McResult::MC_NO_ERROR; + per_thread_api_log_str.clear(); + + if (context == nullptr) { + per_thread_api_log_str = "context ptr (param0) undef (NULL)"; + } else if (bytes == 0) { + per_thread_api_log_str = "invalid bytes "; + } else if (pMem == nullptr) { + per_thread_api_log_str = "invalid ptr (pMem)"; + } else if (false == // + (stateInfo == MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT || // + stateInfo == MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_ATTEMPTS ||// + stateInfo == MC_CONTEXT_CONNECTED_COMPONENT_FACE_WINDING_ORDER)) // check all possible values + { + per_thread_api_log_str = "invalid stateInfo "; + } else if ( + ((stateInfo == MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_CONSTANT) && bytes != sizeof(McDouble)) || // + ((stateInfo == MC_CONTEXT_GENERAL_POSITION_ENFORCEMENT_ATTEMPTS) && bytes != sizeof(McUint32))|| // + ((stateInfo == MC_CONTEXT_CONNECTED_COMPONENT_FACE_WINDING_ORDER) && bytes != sizeof(McConnectedComponentFaceWindingOrder))) { + per_thread_api_log_str = "invalid num bytes"; // leads to e.g. "out of bounds" memory access during memcpy + } else { + try { + bind_state_impl(context, stateInfo, bytes, pMem); + } + CATCH_POSSIBLE_EXCEPTIONS(per_thread_api_log_str); + } + + if (!per_thread_api_log_str.empty()) { + std::fprintf(stderr, "%s(...) -> %s\n", __FUNCTION__, per_thread_api_log_str.c_str()); + if (return_value == McResult::MC_NO_ERROR) // i.e. problem with basic local parameter checks + { + return_value = McResult::MC_INVALID_VALUE; + } + } + + return return_value; +} + MCAPI_ATTR McResult MCAPI_CALL mcCreateUserEvent( McEvent* event, McContext context) @@ -396,10 +444,7 @@ MCAPI_ATTR McResult MCAPI_CALL mcEnqueueDispatch( const McEvent* pEventWaitList, McEvent* pEvent) { - TIMESTACK_RESET(); // reset tracking vars - - SCOPED_TIMER(__FUNCTION__); - + McResult return_value = McResult::MC_NO_ERROR; per_thread_api_log_str.clear(); @@ -528,6 +573,89 @@ MCAPI_ATTR McResult MCAPI_CALL mcDispatch( return return_value; } +MCAPI_ATTR McResult MCAPI_CALL mcEnqueueDispatchPlanarSection( + const McContext context, + McFlags dispatchFlags, + const McVoid* pSrcMeshVertices, + const uint32_t* pSrcMeshFaceIndices, + const uint32_t* pSrcMeshFaceSizes, + uint32_t numSrcMeshVertices, + uint32_t numSrcMeshFaces, + const McDouble* pNormalVector, + const McDouble sectionOffset, + uint32_t numEventsInWaitlist, + const McEvent* pEventWaitList, + McEvent* pEvent) + { + McResult return_value = McResult::MC_NO_ERROR; + per_thread_api_log_str.clear(); + + if (context == nullptr) { + per_thread_api_log_str = "context ptr (param0) undef (NULL)"; + } else if (dispatchFlags == 0) { + per_thread_api_log_str = "dispatch flags unspecified"; + } else if ((dispatchFlags & MC_DISPATCH_REQUIRE_THROUGH_CUTS) && // + (dispatchFlags & MC_DISPATCH_FILTER_FRAGMENT_LOCATION_UNDEFINED)) { + // The user states that she does not want a partial cut but yet also states that she + // wants to keep fragments with partial cuts. These two options are mutually exclusive! + per_thread_api_log_str = "use of mutually-exclusive flags: MC_DISPATCH_REQUIRE_THROUGH_CUTS & MC_DISPATCH_FILTER_FRAGMENT_LOCATION_UNDEFINED"; + } else if ((dispatchFlags & MC_DISPATCH_VERTEX_ARRAY_FLOAT) == 0 && (dispatchFlags & MC_DISPATCH_VERTEX_ARRAY_DOUBLE) == 0) { + per_thread_api_log_str = "dispatch vertex aray type unspecified"; + } else if (pSrcMeshVertices == nullptr) { + per_thread_api_log_str = "source-mesh vertex-position array ptr undef (NULL)"; + } else if (numSrcMeshVertices < 3) { + per_thread_api_log_str = "invalid source-mesh vertex count"; + } else if (pSrcMeshFaceIndices == nullptr) { + per_thread_api_log_str = "source-mesh face-index array ptr undef (NULL)"; + } /*else if (pSrcMeshFaceSizes == nullptr) { + per_thread_api_log_str = "source-mesh face-size array ptr undef (NULL)"; + }*/ + else if (numSrcMeshFaces < 1) { + per_thread_api_log_str = "invalid source-mesh vertex count"; + } else if (pNormalVector == nullptr) { + per_thread_api_log_str = "normal vector ptr undef (NULL)"; + } else if (pNormalVector[0] == 0.0 && pNormalVector[1] == 0.0 && pNormalVector[2] == 0.0) { + per_thread_api_log_str = "invalid normal vector (zero vector)"; + }else if (sectionOffset <= 0 && sectionOffset >= 1.0) { + per_thread_api_log_str = "invalid section offset parameter"; + } else if (pEventWaitList == nullptr && numEventsInWaitlist > 0) { + per_thread_api_log_str = "invalid event waitlist ptr (NULL)"; + } else if (pEventWaitList != nullptr && numEventsInWaitlist == 0) { + per_thread_api_log_str = "invalid event waitlist size (zero)"; + } else if (pEventWaitList == nullptr && numEventsInWaitlist == 0 && pEvent == nullptr) { + per_thread_api_log_str = "invalid event ptr (zero)"; + } else { + try { + dispatch_planar_section_impl( + context, + dispatchFlags, + pSrcMeshVertices, + pSrcMeshFaceIndices, + pSrcMeshFaceSizes, + numSrcMeshVertices, + numSrcMeshFaces, + pNormalVector, + sectionOffset, + numEventsInWaitlist, + pEventWaitList, + pEvent); + } + CATCH_POSSIBLE_EXCEPTIONS(per_thread_api_log_str); + } + + if (!per_thread_api_log_str.empty()) { + + std::fprintf(stderr, "%s(...) -> %s\n", __FUNCTION__, per_thread_api_log_str.c_str()); + + if (return_value == McResult::MC_NO_ERROR) // i.e. problem with basic local parameter checks + { + return_value = McResult::MC_INVALID_VALUE; + } + } + + return return_value; + } + MCAPI_ATTR McResult MCAPI_CALL mcEnqueueGetConnectedComponents( const McContext context, const McConnectedComponentType connectedComponentType, diff --git a/src/mcut/source/preproc.cpp b/src/mcut/source/preproc.cpp index fd63ae713..01b21f928 100644 --- a/src/mcut/source/preproc.cpp +++ b/src/mcut/source/preproc.cpp @@ -4,8 +4,8 @@ #include "mcut/internal/hmesh.h" #include "mcut/internal/kernel.h" #include "mcut/internal/math.h" -#include "mcut/internal/utils.h" #include "mcut/internal/timer.h" +#include "mcut/internal/utils.h" #include // std::partial_sum #include @@ -13,8 +13,8 @@ // If the inputs are found to not be in general position, then we perturb the // cut-mesh by this constant (scaled by bbox diag times a random variable [0.1-1.0]). -const double GENERAL_POSITION_ENFORCMENT_CONSTANT = 1e-4; -const int MAX_PERTUBATION_ATTEMPTS = 1 << 3; +// const double GENERAL_POSITION_ENFORCMENT_CONSTANT = 1e-4; +// const int MAX_PERTUBATION_ATTEMPTS = 1 << 3; // this function converts an index array mesh (e.g. as recieved by the dispatch // function) into a halfedge mesh representation for the kernel backend. @@ -111,7 +111,7 @@ bool client_input_arrays_to_hmesh( #if 0 std::partial_sum(partial_sums.begin(), partial_sums.end(), partial_sums.data()); #else - parallel_partial_sum(context_ptr->get_shared_compute_threadpool() , partial_sums.begin(), partial_sums.end()); + parallel_partial_sum(context_ptr->get_shared_compute_threadpool(), partial_sums.begin(), partial_sums.end()); #endif { typedef std::vector::const_iterator InputStorageIteratorType; @@ -131,7 +131,7 @@ bool client_input_arrays_to_hmesh( if (face_vertex_count < 3) { int zero = (int)McResult::MC_NO_ERROR; - bool exchanged = atm_result.compare_exchange_strong(zero, 1); + bool exchanged = atm_result.compare_exchange_strong(zero, 1, std::memory_order_acq_rel); if (exchanged) // first thread to detect error { context_ptr->dbg_cb( // @@ -150,14 +150,9 @@ bool client_input_arrays_to_hmesh( for (int j = 0; j < face_vertex_count; ++j) { uint32_t idx = ((uint32_t*)pFaceIndices)[faceBaseOffset + j]; - MCUT_ASSERT(idx < numVertices); - - const vertex_descriptor_t descr(idx); - const bool isDuplicate = std::find(faceVertices.cbegin(), faceVertices.cend(), descr) != faceVertices.cend(); - - if (isDuplicate) { + if (idx >= numVertices) { int zero = (int)McResult::MC_NO_ERROR; - bool exchanged = atm_result.compare_exchange_strong(zero, 2); + bool exchanged = atm_result.compare_exchange_strong(zero, 2, std::memory_order_acq_rel); if (exchanged) // first thread to detect error { @@ -166,7 +161,26 @@ bool client_input_arrays_to_hmesh( MC_DEBUG_TYPE_ERROR, 0, MC_DEBUG_SEVERITY_HIGH, - "found duplicate vertex in face - " + std::to_string(faceID)); + "vertex index out of range in face - " + std::to_string(faceID)); + } + break; + } + + const vertex_descriptor_t descr(idx); + const bool isDuplicate = std::find(faceVertices.cbegin(), faceVertices.cend(), descr) != faceVertices.cend(); + + if (isDuplicate) { + int zero = (int)McResult::MC_NO_ERROR; + bool exchanged = atm_result.compare_exchange_strong(zero, 2, std::memory_order_acq_rel); + + if (exchanged) // first thread to detect error + { + context_ptr->dbg_cb( + MC_DEBUG_SOURCE_API, + MC_DEBUG_TYPE_ERROR, + 0, + MC_DEBUG_SEVERITY_HIGH, + "duplicate vertex index in face f" + std::to_string(faceID)); } break; } @@ -181,7 +195,7 @@ bool client_input_arrays_to_hmesh( OutputStorageType partial_res; parallel_for( - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), partial_sums.cbegin(), partial_sums.cend(), fn_create_faces, @@ -203,7 +217,7 @@ bool client_input_arrays_to_hmesh( MC_DEBUG_TYPE_ERROR, // 0, // MC_DEBUG_SEVERITY_HIGH, // - "invalid vertices on face - " + std::to_string(faceID)); + "invalid face f" + std::to_string(faceID) + " (potentially contains non-manifold edge)"); return false; } } @@ -216,7 +230,7 @@ bool client_input_arrays_to_hmesh( MCUT_ASSERT(f.valid()); // The behavior is undefined if valid()== false before the call to wait_for OutputStorageType future_res = f.get(); - const int val = atm_result.load(); + const int val = atm_result.load(std::memory_order_acquire); okay = okay && val == 0; if (!okay) { continue; // just go on (to next iteration) in order to at-least wait for all tasks to finish before we return to user @@ -226,10 +240,16 @@ bool client_input_arrays_to_hmesh( okay = okay && result == true; } - if (!okay) { + if (!okay || atm_result.load(std::memory_order_acquire) != 0) { // check if worker threads (or master thread) encountered an error return false; } + // const int val = atm_result.load(); // check if master thread encountered an error + // okay = (val == 0); + // if (!okay) { + // return false; + //} + // add lastly in order to maintain order bool result = add_faces(partial_res.first, partial_res.second); if (!result) { @@ -254,6 +274,18 @@ bool client_input_arrays_to_hmesh( for (int j = 0; j < face_vertex_count; ++j) { uint32_t idx = ((uint32_t*)pFaceIndices)[faceSizeOffset + j]; + + if (idx >= numVertices) { + + context_ptr->dbg_cb( + MC_DEBUG_SOURCE_API, + MC_DEBUG_TYPE_ERROR, + 0, + MC_DEBUG_SEVERITY_HIGH, + "vertex index out of range in face - f" + std::to_string(i)); + return false; + } + const vertex_descriptor_t descr(idx); // = fIter->second; //vmap[*fIter.first]; const bool isDuplicate = std::find(faceVertices.cbegin(), faceVertices.cend(), descr) != faceVertices.cend(); @@ -345,7 +377,7 @@ bool check_input_mesh(std::shared_ptr& context_ptr, const hmesh_t& m) std::vector cc_to_face_count; int n = find_connected_components( #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), #endif fccmap, m, cc_to_vertex_count, cc_to_face_count); @@ -369,7 +401,7 @@ bool check_input_mesh(std::shared_ptr& context_ptr, const hmesh_t& m) MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_NOTIFICATION, - "Vertices (" + std::to_string(fv_count) + ") on face " + std::to_string(*f) + " not coplanar"); + "Vertices (" + std::to_string(fv_count) + ") on face f" + std::to_string(*f) + " are not coplanar"); // No need to return false, simply warn. It is difficult to // know whether the non-coplanarity is severe enough to cause // confusion when computing intersection points between two @@ -397,8 +429,11 @@ McResult convert(const status_t& v) case status_t::INVALID_MESH_INTERSECTION: result = McResult::MC_INVALID_OPERATION; break; + case status_t::INVALID_CUT_MESH: case status_t::INVALID_SRC_MESH: + result = McResult::MC_INVALID_VALUE; + break; default: - std::fprintf(stderr, "[MCUT]: warning - conversion error (McResult=%d)\n", (int)v); + std::fprintf(stderr, "[MCUT]: warning - unknown conversion (McResult=%d)\n", (int)v); } return result; } @@ -1266,7 +1301,7 @@ extern "C" void preproc( input_t kernel_input; // kernel/backend inpout #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - kernel_input.scheduler = &context_ptr->get_shared_compute_threadpool() ; + kernel_input.scheduler = &context_ptr->get_shared_compute_threadpool(); #endif kernel_input.src_mesh = source_hmesh; @@ -1319,7 +1354,7 @@ extern "C" void preproc( kernel_input.keep_cutmesh_seam = true; } - kernel_input.enforce_general_position = (0 != (dispatchFlags & MC_DISPATCH_ENFORCE_GENERAL_POSITION)); + kernel_input.enforce_general_position = (0 != (dispatchFlags & MC_DISPATCH_ENFORCE_GENERAL_POSITION)) || (0 != (dispatchFlags & MC_DISPATCH_ENFORCE_GENERAL_POSITION_ABSOLUTE)); // Construct BVHs // :::::::::::::: @@ -1332,7 +1367,7 @@ extern "C" void preproc( std::vector> source_hmesh_face_aabb_array; build_oibvh( #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), #endif *source_hmesh.get(), source_hmesh_BVH_aabb_array, source_hmesh_BVH_leafdata_array, source_hmesh_face_aabb_array); #else @@ -1404,7 +1439,10 @@ extern "C" void preproc( int cut_mesh_perturbation_count = 0; // number of times we have perturbed the cut mesh int kernel_invocation_counter = -1; // number of times we have called the internal dispatch/intersect function - double numerical_perturbation_constant = 0.0; // = cut_hmesh_aabb_diag * GENERAL_POSITION_ENFORCMENT_CONSTANT; + double relative_perturbation_constant = 0.0; // i.e. relative to the bbox diagonal length + // the (translation) vector to hold the values with which we will + // carry out numerical perturbation of the cutting surface + vec3 perturbation(0.0, 0.0, 0.0); // RESOLVE mesh intersections // :::::::::::::::::::::::::: @@ -1420,6 +1458,8 @@ extern "C" void preproc( // And if floating polygons arise, then we partition the suspected face into two new faces with an edge that is guaranteed to be // severed during the cut. do { + TIMESTACK_RESET(); + kernel_invocation_counter++; // here we check the reason (if any) for entering the loop body. @@ -1440,17 +1480,15 @@ extern "C" void preproc( kernel_output.status = status_t::SUCCESS; #endif - // the (translation) vector to hold the values with which we will - // carry out numerical perturbation of the cutting surface - vec3 perturbation(0.0, 0.0, 0.0); - if (general_position_assumption_was_violated) { // i.e. do we need to perturb the cut-mesh? MCUT_ASSERT(floating_polygon_was_detected == false); // cannot occur at same time! (see kernel) - if (cut_mesh_perturbation_count == MAX_PERTUBATION_ATTEMPTS) { + context_ptr->dbg_cb(MC_DEBUG_SOURCE_KERNEL, MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_HIGH, "general position assumption violated!"); - context_ptr->dbg_cb(MC_DEBUG_SOURCE_KERNEL, MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_MEDIUM, kernel_output.logger.get_reason_for_failure()); + if (cut_mesh_perturbation_count == context_ptr->get_general_position_enforcement_constant()) { + + context_ptr->dbg_cb(MC_DEBUG_SOURCE_KERNEL, MC_DEBUG_TYPE_OTHER, 0, MC_DEBUG_SEVERITY_HIGH, kernel_output.logger.get_reason_for_failure()); throw std::runtime_error("max perturbation iteratons reached"); } @@ -1459,14 +1497,14 @@ extern "C" void preproc( // not intersect at all, which means we need to perturb again. kernel_input.general_position_enforcement_count = cut_mesh_perturbation_count; - MCUT_ASSERT(numerical_perturbation_constant != double(0.0)); + MCUT_ASSERT(relative_perturbation_constant != double(0.0)); static thread_local std::default_random_engine random_engine(1); static thread_local std::mt19937 mersenne_twister_generator(random_engine()); static thread_local std::uniform_real_distribution uniform_distribution(-1.0, 1.0); for (int i = 0; i < 3; ++i) { - perturbation[i] = uniform_distribution(mersenne_twister_generator) * numerical_perturbation_constant; + perturbation[i] = uniform_distribution(mersenne_twister_generator) * relative_perturbation_constant; } cut_mesh_perturbation_count++; @@ -1488,7 +1526,12 @@ extern "C" void preproc( throw std::invalid_argument("invalid cut-mesh arrays"); } - numerical_perturbation_constant = cut_hmesh_aabb_diag * GENERAL_POSITION_ENFORCMENT_CONSTANT; + /*const*/ double perturbation_scalar = cut_hmesh_aabb_diag; + if (dispatchFlags & MC_DISPATCH_ENFORCE_GENERAL_POSITION_ABSOLUTE) { + perturbation_scalar = 1.0; + } + + relative_perturbation_constant = perturbation_scalar * context_ptr->get_general_position_enforcement_constant(); kernel_input.cut_mesh = cut_hmesh; @@ -1498,11 +1541,11 @@ extern "C" void preproc( cut_hmesh_BVH_leafdata_array.clear(); build_oibvh( #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), #endif - *cut_hmesh.get(), cut_hmesh_BVH_aabb_array, cut_hmesh_BVH_leafdata_array, cut_hmesh_face_face_aabb_array, numerical_perturbation_constant); + *cut_hmesh.get(), cut_hmesh_BVH_aabb_array, cut_hmesh_BVH_leafdata_array, cut_hmesh_face_face_aabb_array, relative_perturbation_constant); #else - cut_hmesh_BVH.buildTree(cut_hmesh, numerical_perturbation_constant); + cut_hmesh_BVH.buildTree(cut_hmesh, relative_perturbation_constant); #endif source_or_cut_hmesh_BVH_rebuilt = true; } @@ -1539,7 +1582,7 @@ extern "C" void preproc( source_hmesh_BVH_leafdata_array.clear(); build_oibvh( #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), #endif *source_hmesh.get(), source_hmesh_BVH_aabb_array, @@ -1556,15 +1599,15 @@ extern "C" void preproc( cut_hmesh_BVH_leafdata_array.clear(); build_oibvh( #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), #endif *cut_hmesh.get(), cut_hmesh_BVH_aabb_array, cut_hmesh_BVH_leafdata_array, cut_hmesh_face_face_aabb_array, - numerical_perturbation_constant); + relative_perturbation_constant); #else - cut_hmesh_BVH.buildTree(cut_hmesh, numerical_perturbation_constant); + cut_hmesh_BVH.buildTree(cut_hmesh, relative_perturbation_constant); #endif } @@ -1605,7 +1648,7 @@ extern "C" void preproc( #else BoundingVolumeHierarchy::intersectBVHTrees( #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL) - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), #endif ps_face_to_potentially_intersecting_others, source_hmesh_BVH, @@ -1727,16 +1770,16 @@ extern "C" void preproc( // const std::string cs_patch_loc_str = to_string(j->first); for (std::vector>::const_iterator k = j->second.cbegin(); k != j->second.cend(); ++k) { - + std::shared_ptr cc_ptr = std::shared_ptr(new fragment_cc_t, fn_delete_cc); - + MCUT_ASSERT(cc_ptr != nullptr); - + std::shared_ptr asFragPtr = std::dynamic_pointer_cast(cc_ptr); - + MCUT_ASSERT(asFragPtr != nullptr); - asFragPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + asFragPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); asFragPtr->type = MC_CONNECTED_COMPONENT_TYPE_FRAGMENT; asFragPtr->fragmentLocation = convert(i->first); asFragPtr->patchLocation = convert(j->first); @@ -1757,6 +1800,8 @@ extern "C" void preproc( asFragPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asFragPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asFragPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object } } @@ -1774,15 +1819,15 @@ extern "C" void preproc( for (std::vector>::const_iterator j = i->second.cbegin(); j != i->second.cend(); ++j) { // for each mesh std::shared_ptr cc_ptr = std::shared_ptr(new fragment_cc_t, fn_delete_cc); - + MCUT_ASSERT(cc_ptr != nullptr); - + std::shared_ptr asFragPtr = std::dynamic_pointer_cast(cc_ptr); - + MCUT_ASSERT(asFragPtr != nullptr); - - asFragPtr->m_user_handle = reinterpret_cast(g_objects_counter++); - + + asFragPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); + asFragPtr->type = MC_CONNECTED_COMPONENT_TYPE_FRAGMENT; asFragPtr->fragmentLocation = convert(i->first); asFragPtr->patchLocation = McPatchLocation::MC_PATCH_LOCATION_UNDEFINED; @@ -1800,6 +1845,8 @@ extern "C" void preproc( asFragPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asFragPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asFragPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object } } @@ -1813,20 +1860,20 @@ extern "C" void preproc( it != insidePatches.cend(); ++it) { - std::shared_ptr cc_ptr = std::shared_ptr(new patch_cc_t, fn_delete_cc); - - MCUT_ASSERT(cc_ptr != nullptr); - - std::shared_ptr asPatchPtr = std::dynamic_pointer_cast(cc_ptr); - - MCUT_ASSERT(asPatchPtr != nullptr); - - asPatchPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + std::shared_ptr cc_ptr = std::shared_ptr(new patch_cc_t, fn_delete_cc); + + MCUT_ASSERT(cc_ptr != nullptr); + + std::shared_ptr asPatchPtr = std::dynamic_pointer_cast(cc_ptr); + + MCUT_ASSERT(asPatchPtr != nullptr); + + asPatchPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); #if 0 // std::shared_ptr patchConnComp = std::unique_ptr(new patch_cc_t, fn_delete_cc); // McConnectedComponent clientHandle = reinterpret_cast(patchConnComp.get()); // context_ptr->connected_components.emplace(clientHandle, std::move(patchConnComp)); - const McConnectedComponent handle = reinterpret_cast(g_objects_counter++); + const McConnectedComponent handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); // allocate internal context object (including associated threadpool etc.) context_ptr->connected_components.add_or_update_mapping(handle, std::shared_ptr(new patch_cc_t, fn_delete_cc)); @@ -1851,6 +1898,8 @@ extern "C" void preproc( asPatchPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asPatchPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asPatchPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object } TIMESTACK_POP(); @@ -1860,20 +1909,20 @@ extern "C" void preproc( const std::vector>& outsidePatches = kernel_output.outside_patches[cm_patch_winding_order_t::DEFAULT]; for (std::vector>::const_iterator it = outsidePatches.cbegin(); it != outsidePatches.cend(); ++it) { - std::shared_ptr cc_ptr = std::shared_ptr(new patch_cc_t, fn_delete_cc); - - MCUT_ASSERT(cc_ptr != nullptr); - - std::shared_ptr asPatchPtr = std::dynamic_pointer_cast(cc_ptr); - - MCUT_ASSERT(asPatchPtr != nullptr); - - asPatchPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + std::shared_ptr cc_ptr = std::shared_ptr(new patch_cc_t, fn_delete_cc); + + MCUT_ASSERT(cc_ptr != nullptr); + + std::shared_ptr asPatchPtr = std::dynamic_pointer_cast(cc_ptr); + + MCUT_ASSERT(asPatchPtr != nullptr); + + asPatchPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); #if 0 /// std::shared_ptr patchConnComp = std::unique_ptr(new patch_cc_t, fn_delete_cc); // McConnectedComponent clientHandle = reinterpret_cast(patchConnComp.get()); // context_ptr->connected_components.emplace(clientHandle, std::move(patchConnComp)); - const McConnectedComponent handle = reinterpret_cast(g_objects_counter++); + const McConnectedComponent handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); // allocate internal context object (including associated threadpool etc.) context_ptr->connected_components.add_or_update_mapping(handle, std::shared_ptr(new patch_cc_t, fn_delete_cc)); @@ -1897,6 +1946,8 @@ extern "C" void preproc( asPatchPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asPatchPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asPatchPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object } TIMESTACK_POP(); @@ -1910,22 +1961,22 @@ extern "C" void preproc( if (kernel_output.seamed_src_mesh != nullptr && kernel_output.seamed_src_mesh->mesh->number_of_faces() > 0) { TIMESTACK_PUSH("store source-mesh seam"); - + std::shared_ptr cc_ptr = std::shared_ptr(new seam_cc_t, fn_delete_cc); - + MCUT_ASSERT(cc_ptr != nullptr); - + std::shared_ptr asSrcMeshSeamPtr = std::dynamic_pointer_cast(cc_ptr); - + MCUT_ASSERT(asSrcMeshSeamPtr != nullptr); - - asSrcMeshSeamPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + + asSrcMeshSeamPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); #if 0 // std::shared_ptr srcMeshSeam = std::unique_ptr(new seam_cc_t, fn_delete_cc); // McConnectedComponent clientHandle = reinterpret_cast(srcMeshSeam.get()); // context_ptr->connected_components.emplace(clientHandle, std::move(srcMeshSeam)); // allocate internal context object (including associated threadpool etc.) - const McConnectedComponent handle = reinterpret_cast(g_objects_counter++); + const McConnectedComponent handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); context_ptr->connected_components.add_or_update_mapping(handle, std::shared_ptr(new seam_cc_t, fn_delete_cc)); @@ -1950,6 +2001,8 @@ extern "C" void preproc( asSrcMeshSeamPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asSrcMeshSeamPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asSrcMeshSeamPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object TIMESTACK_POP(); @@ -1959,21 +2012,21 @@ extern "C" void preproc( if (kernel_output.seamed_cut_mesh != nullptr && kernel_output.seamed_cut_mesh->mesh->number_of_faces() > 0) { TIMESTACK_PUSH("store cut-mesh seam"); - + std::shared_ptr cc_ptr = std::shared_ptr(new seam_cc_t, fn_delete_cc); - + MCUT_ASSERT(cc_ptr != nullptr); - + std::shared_ptr asCutMeshSeamPtr = std::dynamic_pointer_cast(cc_ptr); - + MCUT_ASSERT(asCutMeshSeamPtr != nullptr); - - asCutMeshSeamPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + + asCutMeshSeamPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); #if 0 // std::shared_ptr cutMeshSeam = std::unique_ptr(new seam_cc_t, fn_delete_cc); // McConnectedComponent clientHandle = reinterpret_cast(cutMeshSeam.get()); // context_ptr->connected_components.emplace(clientHandle, std::move(cutMeshSeam)); - const McConnectedComponent handle = reinterpret_cast(g_objects_counter++); + const McConnectedComponent handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); context_ptr->connected_components.add_or_update_mapping(handle, std::shared_ptr(new seam_cc_t, fn_delete_cc)); @@ -1997,6 +2050,8 @@ extern "C" void preproc( asCutMeshSeamPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asCutMeshSeamPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asCutMeshSeamPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object TIMESTACK_POP(); @@ -2010,19 +2065,19 @@ extern "C" void preproc( TIMESTACK_PUSH("store original cut-mesh"); std::shared_ptr cc_ptr = std::shared_ptr(new input_cc_t, fn_delete_cc); - + MCUT_ASSERT(cc_ptr != nullptr); - + std::shared_ptr asCutMeshInputPtr = std::dynamic_pointer_cast(cc_ptr); - + MCUT_ASSERT(asCutMeshInputPtr != nullptr); - - asCutMeshInputPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + + asCutMeshInputPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); #if 0 // std::shared_ptr internalCutMesh = std::unique_ptr(new input_cc_t, fn_delete_cc); // McConnectedComponent clientHandle = reinterpret_cast(internalCutMesh.get()); // context_ptr->connected_components.emplace(clientHandle, std::move(internalCutMesh)); - const McConnectedComponent handle = reinterpret_cast(g_objects_counter++); + const McConnectedComponent handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); context_ptr->connected_components.add_or_update_mapping(handle, std::shared_ptr(new input_cc_t, fn_delete_cc)); @@ -2051,7 +2106,7 @@ extern "C" void preproc( }; parallel_for( - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), cut_hmesh->vertices_begin(), cut_hmesh->vertices_end(), fn_fill_vertex_map); @@ -2060,7 +2115,6 @@ extern "C" void preproc( omi->data_maps.vertex_map[*i] = vd_t((*i) + source_hmesh->number_of_vertices()); // apply offset like kernel does } #endif - } if (kernel_input.populate_face_maps) { @@ -2073,7 +2127,7 @@ extern "C" void preproc( }; parallel_for( - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), cut_hmesh->faces_begin(), cut_hmesh->faces_end(), fn_fill_face_map); @@ -2098,6 +2152,8 @@ extern "C" void preproc( asCutMeshInputPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asCutMeshInputPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asCutMeshInputPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object TIMESTACK_POP(); @@ -2108,19 +2164,19 @@ extern "C" void preproc( TIMESTACK_PUSH("store original src-mesh"); std::shared_ptr cc_ptr = std::shared_ptr(new input_cc_t, fn_delete_cc); - + MCUT_ASSERT(cc_ptr != nullptr); - + std::shared_ptr asSrcMeshInputPtr = std::dynamic_pointer_cast(cc_ptr); - + MCUT_ASSERT(asSrcMeshInputPtr != nullptr); - - asSrcMeshInputPtr->m_user_handle = reinterpret_cast(g_objects_counter++); + + asSrcMeshInputPtr->m_user_handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); #if 0 // std::shared_ptr internalSrcMesh = std::unique_ptr(new input_cc_t, fn_delete_cc); // McConnectedComponent clientHandle = reinterpret_cast(internalSrcMesh.get()); // context_ptr->connected_components.emplace(clientHandle, std::move(internalSrcMesh)); - const McConnectedComponent handle = reinterpret_cast(g_objects_counter++); + const McConnectedComponent handle = reinterpret_cast(g_objects_counter.fetch_add(1, std::memory_order_relaxed)); context_ptr->connected_components.add_or_update_mapping(handle, std::shared_ptr(new input_cc_t, fn_delete_cc)); @@ -2145,7 +2201,7 @@ extern "C" void preproc( }; parallel_for( - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), source_hmesh->vertices_begin(), source_hmesh->vertices_end(), fn_fill_vertex_map); @@ -2166,7 +2222,7 @@ extern "C" void preproc( }; parallel_for( - context_ptr->get_shared_compute_threadpool() , + context_ptr->get_shared_compute_threadpool(), source_hmesh->faces_begin(), source_hmesh->faces_end(), fn_fill_face_map); @@ -2191,6 +2247,8 @@ extern "C" void preproc( asSrcMeshInputPtr->internal_sourcemesh_face_count = source_hmesh->number_of_faces(); asSrcMeshInputPtr->client_sourcemesh_face_count = numSrcMeshFaces; // or source_hmesh_face_count + asSrcMeshInputPtr->perturbation_vector = perturbation; + context_ptr->connected_components.push_front(cc_ptr); // copy the connected component ptr into the context object TIMESTACK_POP();