resolve numerical rounding bug that occurred when calculating polygon normals; fix equality operator bug in numeric.h

Author: chitalu

include/mcut/internal/hmesh.h

@@ -455,7 +455,7 @@ typedef halfedge_descriptor_t hd_t;
 typedef edge_descriptor_t ed_t;
 typedef face_descriptor_t fd_t;
 
-void write_off(const char* fpath, const hmesh_t& mesh);
+void write_off(const char* fpath, const hmesh_t& mesh, const double multiplier);
 void read_off(hmesh_t& mesh, const char* fpath);
 
 template <typename V = face_array_t>

include/mcut/internal/math.h

@@ -50,6 +50,19 @@
 
 #include "nfg/numerics.h" // Indirect_Predicates
 
+// Shewchuk predicates : shewchuk.c
+extern "C"
+{
+	// void exactinit();
+	double orient2d(const double* pa, const double* pb, const double* pc);
+	double orient3d(const double* pa, const double* pb, const double* pc, const double* pd);
+	double orient3dfast(const double* pa, const double* pb, const double* pc, const double* pd);
+	double incircle(const double* pa, const double* pb, const double* pc, const double* pd);
+	double insphere(
+		const double* pa, const double* pb, const double* pc, const double* pd, const double* pe);
+}
+
+
 class rational_number : public bigrational
 {
 public:
@@ -786,6 +799,21 @@ void project_to_2d(std::vector<vec2>& out, const std::vector<vec3>& polygon_vert
 bool coplaner(const vec3& pa, const vec3& pb, const vec3& pc,
     const vec3& pd);
 
+static bool collinear(const vec2_<double>& a,
+			   const vec2_<double>& b,
+			   const vec2_<double>& c,
+			   double& predResult)
+{
+	const double pa_[2] = {a.x(), a.y()};
+	const double pb_[2] = {b.x(), b.y()};
+	const double pc_[2] = {c.x(), c.y()};
+
+	predResult = orient2d(pa_, pb_, pc_); // shewchuk predicate
+
+	//predResult = orient2d(a, b, c);
+	return predResult == double(0.);
+}
+
 bool collinear(const vec2& a, const vec2& b, const vec2& c, scalar_t& predResult);
 
 bool collinear(const vec2& a, const vec2& b, const vec2& c);
@@ -909,15 +937,6 @@ void make_bbox(bounding_box_t<vector_type>& bbox, const vector_type* vertices, c
     }
 }
 
-// Shewchuk predicates : shewchuk.c
-extern "C" {
-// void exactinit();
-double orient2d(const double* pa, const double* pb, const double* pc);
-double orient3d(const double* pa, const double* pb, const double* pc, const double* pd);
-double orient3dfast(const double* pa, const double* pb, const double* pc, const double* pd);
-double incircle(const double* pa, const double* pb, const double* pc, const double* pd);
-double insphere(const double* pa, const double* pb, const double* pc, const double* pd, const double* pe);
-}
 
 #	ifdef MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
 static scalar_t orient2d(const scalar_t* pa, const scalar_t* pb, const scalar_t* pc)

include/mcut/internal/nfg/numerics.hpp

@@ -732,14 +732,14 @@ inline void bignatural::operator>>=(uint32_t n) {
 
 inline bool bignatural::operator==(const bignatural& b) const {
 	if (size() != b.size()) return false;
-	for (uint32_t i = 0; i < size(); i++) if (digits[i] == b.digits[i]) return false;
+	for (uint32_t i = 0; i < size(); i++) if (digits[i] != b.digits[i]) return false;
 	return true;
 }
 
 inline bool bignatural::operator!=(const bignatural& b) const {
 	if (size() != b.size()) return true;
-	for (uint32_t i = 0; i < size(); i++) if (digits[i] == b.digits[i]) return true;
-	return false;
+	for (uint32_t i = 0; i < size(); i++) if (digits[i] == b.digits[i]) return false;
+	return true;
 }
 
 inline bool bignatural::operator>=(const bignatural& b) const {

source/hmesh.cpp

@@ -1263,7 +1263,7 @@ const face_array_iterator_t hmesh_t::elements_begin_(id_<array_iterator_t<face_a
     return faces_begin(account_for_removed_elems);
 }
 
-void write_off(const char* fpath, const hmesh_t& mesh)
+void write_off(const char* fpath, const hmesh_t& mesh, const double multiplier)
 {
 
     std::ofstream outfile(fpath);
@@ -1289,7 +1289,10 @@ void write_off(const char* fpath, const hmesh_t& mesh)
     for (vertex_array_iterator_t iter = mesh.vertices_begin(); iter != mesh.vertices_end(); ++iter) {
         // const vertex_data_t& vdata = iter.second;
         const vec3& point = mesh.vertex(*iter);
-        outfile << (double)point.x() << " " << (double)point.y() << " " << (double)point.z() << "\n";
+		outfile << (double)scalar_t::dequantize(point.x(), multiplier) << " "
+				<< (double)scalar_t::dequantize(point.y(), multiplier) << " "
+				<< (double)scalar_t::dequantize(point.z(), multiplier)
+				<< "\n";
     }
 
     //

source/kernel.cpp

@@ -157,7 +157,7 @@ bool inline ps_is_cutmesh_face(const fd_t& ps_fd, const int sm_face_count)
     return ((int)ps_fd) >= sm_face_count;
 }
 
-void dump_mesh(const hmesh_t& mesh, const char* fbasename)
+void dump_mesh(const hmesh_t& mesh, const char* fbasename, const double multiplier)
 {
     const std::string name = std::string(fbasename) + ".off";
 
@@ -184,7 +184,7 @@ void dump_mesh(const hmesh_t& mesh, const char* fbasename)
         }
     }
 
-    write_off(name.c_str(), mesh);
+    write_off(name.c_str(), mesh, multiplier);
 }
 
 #if 0
@@ -1576,8 +1576,8 @@ void dispatch(output_t& output, const input_t& input)
     const hmesh_t& cs = (*input.cut_mesh);
 
     if (input.verbose) {
-        dump_mesh(sm, "src-mesh");
-        dump_mesh(cs, "cut-mesh");
+		dump_mesh(sm, "src-mesh", input.multiplier);
+		dump_mesh(cs, "cut-mesh", input.multiplier);
     }
 
     const int sm_vtx_cnt = sm.number_of_vertices();
@@ -1736,7 +1736,7 @@ void dispatch(output_t& output, const input_t& input)
     // cs_to_ps_vtx.clear();
 
     if (input.verbose) {
-        dump_mesh(ps, "polygon-soup");
+		dump_mesh(ps, "polygon-soup", input.multiplier);
     }
 
     const int ps_vtx_cnt = ps.number_of_vertices();
@@ -3236,7 +3236,11 @@ void dispatch(output_t& output, const input_t& input)
     }
 
     if (input.verbose) {
-        dump_mesh(m0, "m0.v"); // containing only vertices (polygon soup vertices and newly computed intersection points)
+		dump_mesh(
+			m0,
+			"m0.v",
+			input
+				.multiplier); // containing only vertices (polygon soup vertices and newly computed intersection points)
     }
 
     if (partial_cut_detected) {
@@ -3594,7 +3598,7 @@ void dispatch(output_t& output, const input_t& input)
     // intersecting faces in the polygon-soup ("ps").
 
     if (input.verbose) {
-        dump_mesh(m0, "m0.v.e"); // containing only vertices & edges
+		dump_mesh(m0, "m0.v.e", input.multiplier); // containing only vertices & edges
     }
 
     const uint32_t m0_num_cutpath_edges = (uint32_t)m0_cutpath_edges.size();
@@ -6332,7 +6336,9 @@ void dispatch(output_t& output, const input_t& input)
             output.seamed_src_mesh->data_maps = std::move(separated_src_mesh_fragments.begin()->second.front().second.data_maps);
 
             if (input.verbose) {
-                dump_mesh(output.seamed_src_mesh->mesh.get()[0], "src-mesh-traced-poly");
+				dump_mesh(output.seamed_src_mesh->mesh.get()[0],
+						  "src-mesh-traced-poly",
+						  input.multiplier);
             }
         }
     } // if (input.include_seam_srcmesh) {
@@ -6382,7 +6388,8 @@ void dispatch(output_t& output, const input_t& input)
                 output.seamed_cut_mesh->data_maps = std::move(separated_cut_mesh_fragments.begin()->second.front().second.data_maps);
 
                 if (input.verbose) {
-                    dump_mesh(output.seamed_cut_mesh->mesh.get()[0], "cut-mesh-traced-poly");
+					dump_mesh(
+						output.seamed_cut_mesh->mesh.get()[0], "cut-mesh-traced-poly", input.multiplier);
                 }
             }
         }
@@ -7755,7 +7762,11 @@ void dispatch(output_t& output, const input_t& input)
             // is empty before calling "extract_connected_components"
             MCUT_ASSERT(mesh_data.size() == 1);
             if (input.verbose) {
-                dump_mesh(mesh_data.front().first.get()[0], ("fragment.unsealed." + std::to_string(cc_id) + "." + to_string(mesh_data.front().second.location)).c_str());
+				dump_mesh(mesh_data.front().first.get()[0],
+						  ("fragment.unsealed." + std::to_string(cc_id) + "." +
+						   to_string(mesh_data.front().second.location))
+							  .c_str(),
+						  input.multiplier);
             }
             std::pair<std::shared_ptr<hmesh_t>, connected_component_info_t>& md = mesh_data.front();
             std::shared_ptr<output_mesh_info_t> omi = std::shared_ptr<output_mesh_info_t>(new output_mesh_info_t);
@@ -9167,7 +9178,11 @@ void dispatch(output_t& output, const input_t& input)
 #endif // #if defined(MCUT_WITH_COMPUTE_HELPER_THREADPOOL)
 
             if (input.verbose) {
-                dump_mesh(patch_mesh.get()[0], ("patch" + std::to_string(cur_patch_idx) + "." + to_string(patch_location) + "." + cs_patch_descriptor_str).c_str());
+				dump_mesh(patch_mesh.get()[0],
+						  ("patch" + std::to_string(cur_patch_idx) + "." +
+						   to_string(patch_location) + "." + cs_patch_descriptor_str)
+							  .c_str(),
+						  input.multiplier);
             }
 
             std::shared_ptr<output_mesh_info_t> omi = std::shared_ptr<output_mesh_info_t>(new output_mesh_info_t);
@@ -10700,7 +10715,12 @@ void dispatch(output_t& output, const input_t& input)
 
                 if (input.verbose) {
                     // const int idx = (int)std::distance(cc_instances.begin(), cc_instance_iter);
-                    dump_mesh(cc_instance.first.get()[0], (std::string("cc") + std::to_string(idx++) + "." + to_string(cc_instance.second.location) + "." + to_string(patchLocation)).c_str());
+					dump_mesh(cc_instance.first.get()[0],
+							  (std::string("cc") + std::to_string(idx++) + "." +
+							   to_string(cc_instance.second.location) + "." +
+							   to_string(patchLocation))
+								  .c_str(),
+							  input.multiplier);
                 }
 
                 std::shared_ptr<output_mesh_info_t> omi = std::shared_ptr<output_mesh_info_t>(new output_mesh_info_t);

source/math.cpp

@@ -151,7 +151,7 @@
 									scalar_t::dequantize(vec1[2], multiplier));
 
 				auto cross_ = cross_product(vec0_, vec1_);
-				if(squared_length(cross_) > 0.)
+				if(squared_length(cross_) > 1e-7)
 				{
 					auto cross_n = normalize(cross_);
 					vec3 cross = vec3(scalar_t::quantize(cross_n[0], multiplier),
@@ -782,12 +782,13 @@
             const vec3& x = polygon_vertices[i];
             out[i] = P * x; // vertex in xz plane
         }
-#else // This code is not reliable because it shadow-projects a polygons which can skew computations 
+#else // This code is not reliable because it shadow-projects a polygon
         for (int i = 0; i < polygon_vertex_count; ++i) { // for each vertex
             vec2& Tp = out[i];
             int k = 0;
             for (int j = 0; j < 3; j++) { // for each component
-                if (j != polygon_plane_normal_largest_component) { /* skip largest coordinate */
+				if(j != polygon_normal_largest_component)
+				{ /* skip largest coordinate */
 
                     Tp[k] = polygon_vertices[i][j];
                     k++;
@@ -969,7 +970,8 @@
         num = a[0] * (d[1] - c[1]) + //
             c[0] * (a[1] - d[1]) + //
             d[0] * (c[1] - a[1]);
-
+		std::cout << "test\n";
+        std::cout << num.get_str() << std::endl << denom.get_str() << std::endl;
         if ((num == scalar_t(0.0)) || (num == denom)) {
             code = 'v';
         }

source/preproc.cpp

@@ -1079,16 +1079,14 @@ void resolve_floating_polygons(
 						it != polyVerts.cend();
 						++it)
 					{
-						// TODO: call "collinear" with vec3_<double> not vec3. That is the only way to 
-						// get a reliable measure because with double everything is in the same space/units 
-						// (lengths areas etc).
+
 						bool are_collinear = collinear(segStart, segEnd, (*it), predResult);
 						// last ditch attempt to prevent the possibility of creating a partitioning
 						// edge that more-or-less passes through a vertex (of origin-face or the floatig poly itself)
 						// see: test41
 						const double epsilon = 1e-6;
 #if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
-						if(are_collinear || epsilon > scalar_t::dequantize(absolute_value(predResult),multiplier))
+						if(are_collinear || (scalar_t::quantize(epsilon, multiplier) > absolute_value(predResult)))
 #else
 						if(are_collinear || (!are_collinear && epsilon > std::fabs(predResult)))
 #endif