Skip to content

Commit

Permalink
partial update to integration rational arithmetic
Browse files Browse the repository at this point in the history
  • Loading branch information
@chitalu
chitalu committed Nov 22, 2024
1 parent ad23303 commit dc6d85f
Show file tree
Hide file tree
Showing 7 changed files with 3,053 additions and 2,355 deletions.
7 changes: 5 additions & 2 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
# MCUT_BUILD_TUTORIALS [default=OFF] - Build tutorials
# MCUT_BUILD_WITH_COMPUTE_HELPER_THREADPOOL [default=ON] - Build as configurable multi-threaded library
# MCUT_BUILD_WITH_API_EVENT_LOGGING [default=OFF] - Build with logging functionality which dumps event creation and destruction notices to the console
# MCUT_WITH_ARBITRARY_PRECISION_NUMBERS [default=OFF] - Build with exact number representations, which means that mcut can recieve, return and internally compute intyersections with an exact number representation.
# MCUT_WITH_ARBITRARY_PRECISION_NUMBERS [default=ON] - Build with exact number representations, which means that mcut can recieve, return and internally compute intyersections with an exact number representation.
#
# This script will define the following CMake cache variables:
#
Expand Down Expand Up @@ -98,7 +98,7 @@ if (MCUT_TOPLEVEL_PROJECT AND NOT MCUT_BUILD_TUTORIALS)
option(MCUT_BUILD_TUTORIALS "Configure to build MCUT tutorials" ON)
endif()
option(MCUT_BUILD_WITH_API_EVENT_LOGGING "Configure to build MCUT with event logging to console" OFF)

option(MCUT_WITH_ARBITRARY_PRECISION_NUMBERS "Configure to build with arbitrary precision numbers" ON)
#
# machine-precision-numbers library targets
#
Expand All @@ -116,6 +116,9 @@ list(APPEND compilation_flags "")
if(MCUT_BUILD_AS_SHARED_LIB)
list(APPEND preprocessor_defs -DMCUT_SHARED_LIB=1)
endif()
if(MCUT_WITH_ARBITRARY_PRECISION_NUMBERS)
list(APPEND preprocessor_defs -DMCUT_WITH_ARBITRARY_PRECISION_NUMBERS=1)
endif()

find_package(Threads REQUIRED)

Expand Down
22 changes: 16 additions & 6 deletions include/mcut/internal/math.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,16 @@ class rational_number : public bigrational

virtual ~rational_number() { }

inline static rational_number zero()
{
return rational_number::zero();
}

inline static rational_number one()
{
return rational_number(1.0);
}

explicit operator double() const
{
return this->get_d();
Expand Down Expand Up @@ -115,7 +125,7 @@ class rational_number : public bigrational

static rational_number abs(rational_number _a)
{
if(_a < rational_number(0))
if(_a < rational_number::zero())
{
auto copy = _a;
copy.negate();
Expand Down Expand Up @@ -143,12 +153,12 @@ class rational_number : public bigrational
static rational_number quantize(const double& d /*double prec value*/,
const double& m /*multiplier*/)
{
assert(d<=m);
assert(m != 0);
MCUT_ASSERT(d<=m);
MCUT_ASSERT(m != 0);

if(d == 0)
{
return rational_number(0);
return rational_number::zero();
}

// map all into normalized range [-1, 1]^3.
Expand All @@ -162,7 +172,7 @@ class rational_number : public bigrational
static double dequantize(const rational_number& i /*rational*/,
const double& m /*multiplier*/)
{
if(i == 0)
if(i == rational_number::zero())
{
return (0.0);
}
Expand Down Expand Up @@ -513,7 +523,7 @@ class matrix_t {
std::vector<T> m_entries;
};

extern scalar_t square_root(const scalar_t& number);
extern scalar_t square_root(const scalar_t& number, double multiplier=1);
extern scalar_t absolute_value(const scalar_t& number);
extern sign_t sign(const scalar_t& number);
extern std::ostream& operator<<(std::ostream& os, const vec3& v);
Expand Down
16 changes: 8 additions & 8 deletions source/bvh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,14 +77,14 @@ unsigned int clz(unsigned int x) // stub
// next power of two from x
int next_power_of_two(int x)
{
x--;
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
x |= x >> 32;
return x;
x--;
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
x++;
return x;
}

// check if "x" is a power of two
Expand Down
97 changes: 76 additions & 21 deletions source/frontend.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -794,8 +794,11 @@ void generate_supertriangle_from_mesh_vertices(
for (uint32_t i = 0; i < 3; ++i) {

const scalar_t comp = n[i];
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
const scalar_t comp_abs = absolute_value(comp);
#else
const scalar_t comp_abs = std::fabs(comp);

#endif
if (comp_abs > max_normal_comp_val_abs) {
max_normal_comp_idx = i;
max_normal_comp_val_abs = comp_abs;
Expand All @@ -817,19 +820,30 @@ void generate_supertriangle_from_mesh_vertices(

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);
vec3 vertex0 = mean_on_plane + uv_pos * (bbox_diag * scalar_t(2.0));
vec3 vertex1 = mean_on_plane + uv_neg * (bbox_diag * scalar_t(2.0));
vec3 vertex2 = mean_on_plane -
(normalize(uv_pos + uv_neg) *
(bbox_diag * scalar_t(2.0))); // 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);
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
double tmp = (double)vertex0[i].get_d();
#else
double tmp = (double)vertex0[i];
#endif
memcpy(dst, &tmp, flt_size);
} else {
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
float tmp = (float)vertex0[i].get_d();
#else
float tmp = (float)vertex0[i];
#endif
memcpy(dst, &tmp, flt_size);
}
counter++;
Expand All @@ -839,9 +853,18 @@ void generate_supertriangle_from_mesh_vertices(
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);
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
double tmp = (double)vertex1[i].get_d();
#else
double tmp = (double)vertex1[i];
#endif
memcpy(dst, &tmp, flt_size);
} else {
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
float tmp = (float)vertex1[i].get_d();
#else
float tmp = (float)vertex1[i];
#endif
memcpy(dst, &tmp, flt_size);
}
counter++;
Expand All @@ -851,9 +874,18 @@ void generate_supertriangle_from_mesh_vertices(
for (uint32_t i = 0; i < 3; ++i) {
void* dst = supertriangle_vertices.data() + (counter * flt_size);
if (have_double) {
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
double tmp = (double)vertex2[i].get_d();
#else
double tmp = (double)vertex2[i];
#endif
memcpy(dst, &vertex2[i], flt_size);
} else {
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
float tmp = (float)vertex2[i].get_d();
#else
float tmp = (float)vertex2[i];
#endif
memcpy(dst, &tmp, flt_size);
}
counter++;
Expand Down Expand Up @@ -1359,7 +1391,7 @@ void triangulate_face(
// 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 scalar_t errbound = 1e-2;
const double errbound = 1e-2;

// 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
Expand All @@ -1369,8 +1401,13 @@ void triangulate_face(
// (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);

#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS // I think this can be made exact (i.e. without "errbound")
const bool incident_edges_are_parallel =
absolute_value(orient2d_res) <= absolute_value(scalar_t::quantize(errbound, multiplier));
#else
const bool incident_edges_are_parallel = std::fabs(orient2d_res) <= std::fabs(errbound);
#endif
if (incident_edges_are_parallel) {
//
// pick the shortest of the two incident edges and compute the
Expand Down Expand Up @@ -1427,8 +1464,11 @@ void triangulate_face(
// 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)"";
//

#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
const scalar_t shift_len = square_root(largest_sqrd_length, multiplier);
#else
const scalar_t shift_len = std::sqrt(largest_sqrd_length);
#endif
const vec2 shift = perturbation_dir * shift_len;

vec2 intersection_point_on_edge = perturbed_dvertex_coords + shift; // some location potentially outside of polygon
Expand All @@ -1443,7 +1483,7 @@ void triangulate_face(

// test segment against all edges to find closest intersection point

double segment_min_tval = 1.0;
scalar_t 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) {
Expand All @@ -1458,8 +1498,8 @@ void triangulate_face(
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
scalar_t segment_tval; // parameter along segment
scalar_t edge_tval; // parameter along current edge
vec2 ipoint; // intersection point between segment and current edge

const char result = compute_segment_intersection(
Expand All @@ -1477,12 +1517,14 @@ void triangulate_face(
// 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);
const scalar_t 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);
const scalar_t 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;
scalar_t 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) {
Expand All @@ -1491,26 +1533,34 @@ void triangulate_face(
}

// ratio along segment
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
segment_tval = square_root(sqr_dist_to_closest, multiplier) / shift_len;
#else
segment_tval = std::sqrt(sqr_dist_to_closest) / shift_len;

#endif
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 <=
scalar_t::one()); // ... 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);

const double scale = (errbound * revised_perturbation_len);
const scalar_t revised_perturbation_len = length(revised_perturbation_vector);
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
const scalar_t scale = (scalar_t::quantize(errbound,multiplier) * revised_perturbation_len);
#else
const scalar_t scale = (errbound * revised_perturbation_len);
#endif
// The translation by which we perturb "perturbed_dvertex_coords"
//
// NOTE: since "perturbation_vector" was constructed from "to_prev" and "to_next",
Expand Down Expand Up @@ -1929,7 +1979,12 @@ void get_connected_component_data_impl_detail(

// for each component of coordinate
for (int i = 0; i < 3; ++i) {
const float val = static_cast<float>(coords[i]);
#if MCUT_WITH_ARBITRARY_PRECISION_NUMBERS
const float val = static_cast<float>( scalar_t::dequantize(coords[i],multiplier));
#else
const float val = static_cast<float>(coords[i]);
#endif

*(casted_ptr + elem_offset) = val;
elem_offset += 1;
}
Expand Down
Loading

0 comments on commit dc6d85f

Please sign in to comment.