Commit 72784c14 authored by gsell's avatar gsell

Merge branch '611-fix-floating-point-comparisons-in-boundarygeometry' into 'master'

Resolve "fix floating point comparisons in BoundaryGeometry"

Closes #611

See merge request !447
parents 16134d76 e886dc76
......@@ -23,8 +23,10 @@
#include <cmath>
#include <fstream>
#include <string>
#include <algorithm>
#include "H5hut.h"
#include <cfloat>
#include "AbstractObjects/OpalData.h"
#include "Algorithms/PartBunchBase.h"
......@@ -40,8 +42,230 @@ extern Inform* gmsg;
#define PointID(triangle_id, vertex_id) Triangles_m[triangle_id][vertex_id]
#define Point(triangle_id, vertex_id) Points_m[Triangles_m[triangle_id][vertex_id]]
#define EPS 10e-10
/*
In the following namespaces various approximately floating point
comparisons are implemented. The used implementation is selected
via
namespaces cmp = IMPLEMENTATION;
*/
/*
First we define some macros for function common in all namespaces.
*/
#define FUNC_EQ(x, y) inline bool eq(double x, double y) { \
return almost_eq(x, y); \
}
#define FUNC_EQ_ZERO(x) inline bool eq_zero(double x) { \
return almost_eq_zero(x); \
}
#define FUNC_LE(x, y) inline bool le(double x, double y) { \
if (almost_eq(x, y)) { \
return true; \
} \
return x < y; \
}
#define FUNC_LE_ZERO(x) inline bool le_zero(double x) { \
if (almost_eq_zero(x)) { \
return true; \
} \
return x < 0.0; \
}
#define FUNC_LT(x, y) inline bool lt(double x, double y) { \
if (almost_eq(x, y)) { \
return false; \
} \
return x < y; \
}
#define FUNC_LT_ZERO(x) inline bool lt_zero(double x) { \
if (almost_eq_zero(x)) { \
return false; \
} \
return x < 0.0; \
}
#define FUNC_GE(x, y) inline bool ge(double x, double y) { \
if (almost_eq(x, y)) { \
return true; \
} \
return x > y; \
}
#define FUNC_GE_ZERO(x) inline bool ge_zero(double x) { \
if (almost_eq_zero(x)) { \
return true; \
} \
return x > 0.0; \
}
#define FUNC_GT(x, y) inline bool gt(double x, double y) { \
if (almost_eq(x, y)) { \
return false; \
} \
return x > y; \
}
#define FUNC_GT_ZERO(x) inline bool gt_zero(double x) { \
if (almost_eq_zero(x)) { \
return false; \
} \
return x > 0.0; \
}
namespace cmp_diff {
/*
Link:
https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
*/
inline bool almost_eq(double A, double B, double maxDiff = 1e-15, double maxRelDiff = DBL_EPSILON) {
// Check if the numbers are really close -- needed
// when comparing numbers near zero.
const double diff = std::abs(A - B);
if (diff <= maxDiff)
return true;
A = std::abs(A);
B = std::abs(B);
const double largest = (B > A) ? B : A;
if (diff <= largest * maxRelDiff)
return true;
return false;
}
inline bool almost_eq_zero(double A, double maxDiff = 1e-15) {
const double diff = std::abs(A);
return (diff <= maxDiff);
}
FUNC_EQ(x, y);
FUNC_EQ_ZERO(x);
FUNC_LE(x, y);
FUNC_LE_ZERO(x);
FUNC_LT(x, y);
FUNC_LT_ZERO(x);
FUNC_GE(x, y);
FUNC_GE_ZERO(x);
FUNC_GT(x, y);
FUNC_GT_ZERO(x);
}
namespace cmp_ulp_obsolete {
/*
See:
https://www.cygnus-software.com/papers/comparingfloats/comparing_floating_point_numbers_obsolete.htm
*/
inline bool almost_eq(double A, double B, double maxDiff = 1e-20, int maxUlps = 1000) {
// Make sure maxUlps is non-negative and small enough that the
// default NAN won't compare as equal to anything.
// assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
// handle NaN's
// Note: comparing something with a NaN is always false!
if (std::isnan(A) || std::isnan(B)) {
return false;
}
if (std::abs (A - B) <= maxDiff) {
return true;
}
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
auto aInt = *(int64_t*)&A;
#pragma GCC diagnostic pop
// Make aInt lexicographically ordered as a twos-complement int
if (aInt < 0) {
aInt = 0x8000000000000000 - aInt;
}
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
auto bInt = *(int64_t*)&B;
#pragma GCC diagnostic pop
// Make bInt lexicographically ordered as a twos-complement int
if (bInt < 0) {
bInt = 0x8000000000000000 - bInt;
}
if (std::abs (aInt - bInt) <= maxUlps) {
return true;
}
return false;
}
inline bool almost_eq_zero(double A, double maxDiff = 1e-15) {
// no need to handle NaN's!
return (std::abs(A) <= maxDiff);
}
FUNC_EQ(x, y);
FUNC_EQ_ZERO(x);
FUNC_LE(x, y);
FUNC_LE_ZERO(x);
FUNC_LT(x, y);
FUNC_LT_ZERO(x);
FUNC_GE(x, y);
FUNC_GE_ZERO(x);
FUNC_GT(x, y);
FUNC_GT_ZERO(x);
}
namespace cmp_ulp {
/*
See:
https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
*/
inline bool almost_eq (double A, double B, double maxDiff = 1e-20, int maxUlps = 1000) {
// handle NaN's
if (std::isnan (A) || std::isnan (B)) {
return false;
}
// Check if the numbers are really close -- needed
// when comparing numbers near zero.
if (std::abs (A - B) <= maxDiff)
return true;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
auto aInt = *(int64_t*)&A;
auto bInt = *(int64_t*)&B;
#pragma GCC diagnostic pop
// Different signs means they do not match.
// Note: a negative floating point number is also negative as integer.
if (std::signbit (aInt) != std::signbit (bInt))
return false;
// Find the difference in ULPs.
return (std::abs (aInt - bInt) <= maxUlps);
}
inline bool almost_eq_zero(double A, double maxDiff = 1e-15) {
return (std::abs (A) <= maxDiff);
}
FUNC_EQ(x, y);
FUNC_EQ_ZERO(x);
FUNC_LE(x, y);
FUNC_LE_ZERO(x);
FUNC_LT(x, y);
FUNC_LT_ZERO(x);
FUNC_GE(x, y);
FUNC_GE_ZERO(x);
FUNC_GT(x, y);
FUNC_GT_ZERO(x);
}
namespace cmp = cmp_ulp;
/*
Some
......@@ -57,19 +281,19 @@ extern Inform* gmsg;
namespace {
struct VectorLessX {
bool operator() (Vector_t x1, Vector_t x2) {
return x1 (0) < x2 (0);
return cmp::lt (x1(0), x2(0));
}
};
struct VectorLessY {
bool operator() (Vector_t x1, Vector_t x2) {
return x1 (1) < x2 (1);
return cmp::lt(x1(1), x2 (1));
}
};
struct VectorLessZ {
bool operator() (Vector_t x1, Vector_t x2) {
return x1 (2) < x2 (2);
return cmp::lt(x1(2), x2(2));
}
};
......@@ -83,7 +307,7 @@ Vector_t get_max_extent (std::vector<Vector_t>& coords) {
coords.begin (), coords.end (), VectorLessY ());
const Vector_t z = *max_element (
coords.begin (), coords.end (), VectorLessZ ());
return Vector_t (x (0), y (1), z (2));
return Vector_t (x(0), y(1), z(2));
}
......@@ -97,7 +321,7 @@ Vector_t get_min_extent (std::vector<Vector_t>& coords) {
coords.begin (), coords.end (), VectorLessY ());
const Vector_t z = *min_element (
coords.begin (), coords.end (), VectorLessZ ());
return Vector_t (x (0), y (1), z (2));
return Vector_t (x(0), y(1), z(2));
}
/*
......@@ -183,11 +407,6 @@ static void write_voxel_mesh (
*/
#define LERP( A, B, C) ((B)+(A)*((C)-(B)))
#define MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define MAX2(a,b) (((a) > (b)) ? (a) : (b))
#define MIN3(a,b,c) ((((a)<(b))&&((a)<(c))) ? (a) : (((b)<(c)) ? (b) : (c)))
#define MAX3(a,b,c) ((((a)>(b))&&((a)>(c))) ? (a) : (((b)>(c)) ? (b) : (c)))
#define INSIDE 0
#define OUTSIDE 1
......@@ -252,12 +471,12 @@ face_plane (
) {
int outcode_fcmp = 0;
if (gsl_fcmp (p[0], 0.5, EPS) > 0) outcode_fcmp |= 0x01;
if (gsl_fcmp (p[0],-0.5, EPS) < 0) outcode_fcmp |= 0x02;
if (gsl_fcmp (p[1], 0.5, EPS) > 0) outcode_fcmp |= 0x04;
if (gsl_fcmp (p[1],-0.5, EPS) < 0) outcode_fcmp |= 0x08;
if (gsl_fcmp (p[2], 0.5, EPS) > 0) outcode_fcmp |= 0x10;
if (gsl_fcmp (p[2],-0.5, EPS) < 0) outcode_fcmp |= 0x20;
if (cmp::gt(p[0], 0.5)) outcode_fcmp |= 0x01;
if (cmp::lt(p[0], -0.5)) outcode_fcmp |= 0x02;
if (cmp::gt(p[1], 0.5)) outcode_fcmp |= 0x04;
if (cmp::lt(p[1], -0.5)) outcode_fcmp |= 0x08;
if (cmp::gt(p[2], 0.5)) outcode_fcmp |= 0x10;
if (cmp::lt(p[2], -0.5)) outcode_fcmp |= 0x20;
return(outcode_fcmp);
}
......@@ -272,18 +491,18 @@ bevel_2d (
) {
int outcode_fcmp = 0;
if (gsl_fcmp( p[0] + p[1], 1.0, EPS) > 0) outcode_fcmp |= 0x001;
if (gsl_fcmp( p[0] - p[1], 1.0, EPS) > 0) outcode_fcmp |= 0x002;
if (gsl_fcmp(-p[0] + p[1], 1.0, EPS) > 0) outcode_fcmp |= 0x004;
if (gsl_fcmp(-p[0] - p[1], 1.0, EPS) > 0) outcode_fcmp |= 0x008;
if (gsl_fcmp( p[0] + p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x010;
if (gsl_fcmp( p[0] - p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x020;
if (gsl_fcmp(-p[0] + p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x040;
if (gsl_fcmp(-p[0] - p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x080;
if (gsl_fcmp( p[1] + p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x100;
if (gsl_fcmp( p[1] - p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x200;
if (gsl_fcmp(-p[1] + p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x400;
if (gsl_fcmp(-p[1] - p[2], 1.0, EPS) > 0) outcode_fcmp |= 0x800;
if (cmp::gt( p[0] + p[1], 1.0)) outcode_fcmp |= 0x001;
if (cmp::gt( p[0] - p[1], 1.0)) outcode_fcmp |= 0x002;
if (cmp::gt(-p[0] + p[1], 1.0)) outcode_fcmp |= 0x004;
if (cmp::gt(-p[0] - p[1], 1.0)) outcode_fcmp |= 0x008;
if (cmp::gt( p[0] + p[2], 1.0)) outcode_fcmp |= 0x010;
if (cmp::gt( p[0] - p[2], 1.0)) outcode_fcmp |= 0x020;
if (cmp::gt(-p[0] + p[2], 1.0)) outcode_fcmp |= 0x040;
if (cmp::gt(-p[0] - p[2], 1.0)) outcode_fcmp |= 0x080;
if (cmp::gt( p[1] + p[2], 1.0)) outcode_fcmp |= 0x100;
if (cmp::gt( p[1] - p[2], 1.0)) outcode_fcmp |= 0x200;
if (cmp::gt(-p[1] + p[2], 1.0)) outcode_fcmp |= 0x400;
if (cmp::gt(-p[1] - p[2], 1.0)) outcode_fcmp |= 0x800;
return(outcode_fcmp);
}
......@@ -298,14 +517,14 @@ bevel_3d (
) {
int outcode_fcmp = 0;
if (gsl_fcmp( p[0] + p[1] + p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x01;
if (gsl_fcmp( p[0] + p[1] - p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x02;
if (gsl_fcmp( p[0] - p[1] + p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x04;
if (gsl_fcmp( p[0] - p[1] - p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x08;
if (gsl_fcmp(-p[0] + p[1] + p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x10;
if (gsl_fcmp(-p[0] + p[1] - p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x20;
if (gsl_fcmp(-p[0] - p[1] + p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x40;
if (gsl_fcmp(-p[0] - p[1] - p[2], 1.5, EPS) > 0) outcode_fcmp |= 0x80;
if (cmp::gt( p[0] + p[1] + p[2], 1.5)) outcode_fcmp |= 0x01;
if (cmp::gt( p[0] + p[1] - p[2], 1.5)) outcode_fcmp |= 0x02;
if (cmp::gt( p[0] - p[1] + p[2], 1.5)) outcode_fcmp |= 0x04;
if (cmp::gt( p[0] - p[1] - p[2], 1.5)) outcode_fcmp |= 0x08;
if (cmp::gt(-p[0] + p[1] + p[2], 1.5)) outcode_fcmp |= 0x10;
if (cmp::gt(-p[0] + p[1] - p[2], 1.5)) outcode_fcmp |= 0x20;
if (cmp::gt(-p[0] - p[1] + p[2], 1.5)) outcode_fcmp |= 0x40;
if (cmp::gt(-p[0] - p[1] - p[2], 1.5)) outcode_fcmp |= 0x80;
return(outcode_fcmp);
}
......@@ -326,9 +545,12 @@ check_point (
) {
Vector_t plane_point;
plane_point[0] = LERP(alpha, p1[0], p2[0]);
plane_point[1] = LERP(alpha, p1[1], p2[1]);
plane_point[2] = LERP(alpha, p1[2], p2[2]);
#define LERP(a, b, t) (a + t * (b - a))
// with C++20 we can use: std::lerp(a, b, t)
plane_point[0] = LERP(p1[0], p2[0], alpha);
plane_point[1] = LERP(p1[1], p2[1], alpha);
plane_point[2] = LERP(p1[2], p2[2], alpha);
#undef LERP
return(face_plane(plane_point) & mask);
}
......@@ -364,7 +586,7 @@ check_line (
Test if 3D point is inside 3D triangle
*/
constexpr double EPS = 10e-15;
static inline int
SIGN3 (
Vector_t A
......@@ -383,12 +605,12 @@ point_triangle_intersection (
First, a quick bounding-box test:
If P is outside triangle bbox, there cannot be an intersection.
*/
if (gsl_fcmp (p[0], MAX3(t.v1(0), t.v2(0), t.v3(0)), EPS) > 0) return(OUTSIDE);
if (gsl_fcmp (p[1], MAX3(t.v1(1), t.v2(1), t.v3(1)), EPS) > 0) return(OUTSIDE);
if (gsl_fcmp (p[2], MAX3(t.v1(2), t.v2(2), t.v3(2)), EPS) > 0) return(OUTSIDE);
if (gsl_fcmp (p[0], MIN3(t.v1(0), t.v2(0), t.v3(0)), EPS) < 0) return(OUTSIDE);
if (gsl_fcmp (p[1], MIN3(t.v1(1), t.v2(1), t.v3(1)), EPS) < 0) return(OUTSIDE);
if (gsl_fcmp (p[2], MIN3(t.v1(2), t.v2(2), t.v3(2)), EPS) < 0) return(OUTSIDE);
if (cmp::gt(p[0], std::max({t.v1(0), t.v2(0), t.v3(0)}))) return(OUTSIDE);
if (cmp::gt(p[1], std::max({t.v1(1), t.v2(1), t.v3(1)}))) return(OUTSIDE);
if (cmp::gt(p[2], std::max({t.v1(2), t.v2(2), t.v3(2)}))) return(OUTSIDE);
if (cmp::lt(p[0], std::min({t.v1(0), t.v2(0), t.v3(0)}))) return(OUTSIDE);
if (cmp::lt(p[1], std::min({t.v1(1), t.v2(1), t.v3(1)}))) return(OUTSIDE);
if (cmp::lt(p[2], std::min({t.v1(2), t.v2(2), t.v3(2)}))) return(OUTSIDE);
/*
For each triangle side, make a vector out of it by subtracting vertexes;
......@@ -515,30 +737,37 @@ triangle_intersects_cube (
if one of the diagonals is parallel to the plane, the other will
intersect the plane
*/
double denom;
if(std::abs(denom=(norm[0] + norm[1] + norm[2]))>EPS) {
double denom = norm[0] + norm[1] + norm[2];
if (cmp::eq_zero(std::abs(denom)) == false) {
/* skip parallel diagonals to the plane; division by 0 can occure */
Vector_t hitpp = d / denom;
if (std::abs(hitpp[0]) <= 0.5)
if (point_triangle_intersection(hitpp,t) == INSIDE) return(INSIDE);
if (cmp::le(std::abs(hitpp[0]), 0.5))
if (point_triangle_intersection(hitpp,t) == INSIDE)
return(INSIDE);
}
if(std::abs(denom=(norm[0] + norm[1] - norm[2]))>EPS) {
denom = norm[0] + norm[1] - norm[2];
if (cmp::eq_zero(std::abs(denom)) == false) {
Vector_t hitpn;
hitpn[2] = -(hitpn[0] = hitpn[1] = d / denom);
if (std::abs(hitpn[0]) <= 0.5)
if (point_triangle_intersection(hitpn,t) == INSIDE) return(INSIDE);
if (cmp::le(std::abs(hitpn[0]), 0.5))
if (point_triangle_intersection(hitpn,t) == INSIDE)
return(INSIDE);
}
if(std::abs(denom=(norm[0] - norm[1] + norm[2]))>EPS) {
denom = norm[0] - norm[1] + norm[2];
if (cmp::eq_zero(std::abs(denom)) == false) {
Vector_t hitnp;
hitnp[1] = -(hitnp[0] = hitnp[2] = d / denom);
if (std::abs(hitnp[0]) <= 0.5)
if (point_triangle_intersection(hitnp,t) == INSIDE) return(INSIDE);
if (cmp::le(std::abs(hitnp[0]), 0.5))
if (point_triangle_intersection(hitnp,t) == INSIDE)
return(INSIDE);
}
if(std::abs(denom=(norm[0] - norm[1] - norm[2]))>EPS) {
denom = norm[0] - norm[1] - norm[2];
if (cmp::eq_zero(std::abs(denom)) == false) {
Vector_t hitnn;
hitnn[1] = hitnn[2] = -(hitnn[0] = d / denom);
if (std::abs(hitnn[0]) <= 0.5)
if (point_triangle_intersection(hitnn,t) == INSIDE) return(INSIDE);
if (cmp::le(std::abs(hitnn[0]), 0.5))
if (point_triangle_intersection(hitnn,t) == INSIDE)
return(INSIDE);
}
/*
......@@ -622,23 +851,23 @@ public:
double tmax_ = (pts[1-r.sign[0]][0] - r.origin[0]) * r.inv_direction[0];
const double tymin = (pts[r.sign[1]][1] - r.origin[1]) * r.inv_direction[1];
const double tymax = (pts[1-r.sign[1]][1] - r.origin[1]) * r.inv_direction[1];
if ( (tmin_ > tymax) || (tymin > tmax_) )
if ( cmp::gt(tmin_, tymax) || cmp::gt(tymin, tmax_) )
return 0; // no intersection
if (tymin > tmin_)
if (cmp::gt(tymin, tmin_))
tmin_ = tymin;
if (tymax < tmax_)
if (cmp::lt(tymax, tmax_))
tmax_ = tymax;
const double tzmin = (pts[r.sign[2]][2] - r.origin[2]) * r.inv_direction[2];
const double tzmax = (pts[1-r.sign[2]][2] - r.origin[2]) * r.inv_direction[2];
if ( (tmin_ > tzmax) || (tzmin > tmax_) )
if ( cmp::gt(tmin_, tzmax) || cmp::gt(tzmin, tmax_) )
return 0; // no intersection
if (tzmin > tmin_)
if (cmp::gt(tzmin, tmin_))
tmin_ = tzmin;
tmin = tmin_;
if (tzmax < tmax_)
if (cmp::lt(tzmax, tmax_))
tmax_ = tzmax;
tmax = tmax_;
return (tmax >= 0);
return cmp::ge_zero(tmax);
}
inline bool intersect (
......@@ -668,12 +897,12 @@ public:
const Vector_t& P
) const {
return (
P[0] >= pts[0][0] &&
P[1] >= pts[0][1] &&
P[2] >= pts[0][2] &&
P[0] <= pts[1][0] &&
P[1] <= pts[1][1] &&
P[2] <= pts[1][2]);
cmp::ge(P[0], pts[0][0])
&& cmp::ge(P[1], pts[0][1])
&& cmp::ge(P[2], pts[0][2])
&& cmp::le(P[0], pts[1][0])
&& cmp::le(P[1], pts[1][1])
&& cmp::le(P[2], pts[1][2]));
}
Vector_t pts[2];
......@@ -686,7 +915,7 @@ static inline Vector_t normalVector (
) {
const Vector_t N = cross (B - A, C - A);
const double magnitude = std::sqrt (SQR (N (0)) + SQR (N (1)) + SQR (N (2)));
PAssert (gsl_fcmp (magnitude, 0.0, EPS) > 0); // in case we have degenerated triangles
PAssert (cmp::gt_zero(magnitude)); // in case we have degenerated triangles
return N / magnitude;
}
......@@ -958,8 +1187,8 @@ BoundaryGeometry::intersectLineTriangle (
const Vector_t w0 = P0 - V0;
const double a = -dot(n,w0);
const double b = dot(n,dir);
if (gsl_fcmp (b, 0.0, EPS) == 0) { // ray is parallel to triangle plane
if (a == 0) { // ray lies in triangle plane
if (cmp::eq_zero(b)) { // ray is parallel to triangle plane
if (cmp::eq_zero(a)) { // ray lies in triangle plane
return 1;
} else { // ray disjoint from plane
return 0;
......@@ -970,12 +1199,12 @@ BoundaryGeometry::intersectLineTriangle (
const double r = a / b;
switch (kind) {
case RAY:
if (r < 0.0) { // ray goes away from triangle
if (cmp::lt_zero(r)) { // ray goes away from triangle
return 0; // => no intersect
}
break;
case SEGMENT:
if (r < 0 || 1.0 < r) { // intersection on extended
if (cmp::lt_zero(r) || cmp::lt(1.0, r)) { // intersection on extended
return 0; // segment
}
break;
......@@ -995,17 +1224,17 @@ BoundaryGeometry::intersectLineTriangle (
// get and test parametric coords
const double s = (uv * wv - vv * wu) / D;
if (s < 0.0 || s > 1.0) { // I is outside T
if (cmp::lt_zero(s) || cmp::gt(s, 1.0)) { // I is outside T
return 0;
}
const double t = (uv * wu - uu * wv) / D;
if (t < 0.0 || (s + t) > 1.0) { // I is outside T
if (cmp::lt_zero(t) || cmp::gt((s + t), 1.0)) { // I is outside T
return 0;
}
// intersection point is in triangle
if (r < 0.0) { // in extended segment in opposite
if (cmp::lt_zero(r)) { // in extended segment in opposite
return 2; // direction of ray
} else if ((0.0 <= r) && (r <= 1.0)) { // in segment
} else if (cmp::ge_zero(r) && cmp::le(r, 1.0)) { // in segment
return 3;
} else { // in extended segment in
return 4; // direction of ray
......@@ -1033,33 +1262,33 @@ BoundaryGeometry::isInside (
// left boundary of bounding box (x direction)
x = maxExtent_m[0] + 0.01;
if (x - P[0] < distance) {
if (cmp::lt(x - P[0], distance)) {
distance = x - P[0];
ref_pt = {x, P[1], P[2]};
}
// lower boundary of bounding box (y direction)
double y = minExtent_m[1] - 0.01;
if (P[1] - y < distance) {
if (cmp::lt(P[1] - y, distance)) {
distance = P[1] -y;
ref_pt = {P[0], y, P[1]};
}
// upper boundary of bounding box (y direction)
y = maxExtent_m[1] + 0.01;
if (y - P[1] < distance) {
if (cmp::lt(y - P[1], distance)) {
distance = y - P[1];
ref_pt = {P[0], y, P[2]};
}
// front boundary of bounding box (z direction)
double z = minExtent_m[2] - 0.01;
if (P[2] - z < distance) {
if (cmp::lt(P[2] - z, distance)) {
distance = P[2] - z;
ref_pt = {P[0], P[1], z};
}
// back boundary of bounding box (z direction)
z = maxExtent_m[2] + 0.01;
if (z - P[2] < distance) {
if (cmp::lt(z - P[2], distance)) {
ref_pt = {P[0], P[1], z};
}
......@@ -1206,9 +1435,9 @@ BoundaryGeometry::fastIsInside (
}
#endif
const Vector_t v = reference_pt - P;
const int N = std::ceil (magnitude (v) / MIN3 (voxelMesh_m.sizeOfVoxel [0],
const int N = std::ceil (magnitude (v) / std