moab
|
Classes | |
class | VolMap |
Class representing a 3-D mapping function (e.g. shape function for volume element) More... | |
class | LinearHexMap |
Shape function for trilinear hexahedron. More... | |
Enumerations | |
enum | intersection_type { NONE, INTERIOR, NODE0, NODE1, NODE2, EDGE0, EDGE1, EDGE2 } |
Plücker test for intersection between a ray and a triangle. More... | |
Functions | |
static void | min_max_3 (double a, double b, double c, double &min, double &max) |
static double | dot_abs (const CartVect &u, const CartVect &v) |
bool | segment_box_intersect (CartVect box_min, CartVect box_max, const CartVect &seg_pt, const CartVect &seg_unit_dir, double &seg_start, double &seg_end) |
bool | first (const CartVect &a, const CartVect &b) |
bool | plucker_ray_tri_intersect (const CartVect vertices[3], const CartVect &origin, const CartVect &direction, double, double &dist_out, const double *nonneg_ray_len, const double *neg_ray_len, const int *orientation, intersection_type *type) |
bool | ray_tri_intersect (const CartVect vertices[3], const CartVect &ray_point, const CartVect &ray_unit_direction, double tolerance, double &t_out, const double *ray_length=0) |
Test for intersection between a ray and a triangle. | |
bool | ray_box_intersect (const CartVect &box_min, const CartVect &box_max, const CartVect &ray_pt, const CartVect &ray_dir, double &t_enter, double &t_exit) |
bool | box_plane_overlap (const CartVect &plane_normal, double plane_coeff, CartVect box_min_corner, CartVect box_max_corner) |
Test if plane intersects axis-aligned box. | |
bool | box_tri_overlap (const CartVect triangle_corners[3], const CartVect &box_center, const CartVect &box_half_dims) |
Test if triangle intersects axis-aligned box. | |
bool | box_tri_overlap (const CartVect triangle_corners[3], const CartVect &box_min_corner, const CartVect &box_max_corner, double tolerance) |
Test if triangle intersects axis-aligned box. | |
bool | box_elem_overlap (const CartVect *elem_corners, EntityType elem_type, const CartVect &box_center, const CartVect &box_half_dims) |
Test if the specified element intersects an axis-aligned box. | |
static CartVect | quad_norm (const CartVect &v1, const CartVect &v2, const CartVect &v3, const CartVect &v4) |
static CartVect | tri_norm (const CartVect &v1, const CartVect &v2, const CartVect &v3) |
bool | box_linear_elem_overlap (const CartVect *elem_corners, EntityType elem_type, const CartVect &box_center, const CartVect &box_half_dims) |
Test if the specified element intersects an axis-aligned box. | |
bool | box_linear_elem_overlap (const CartVect *elem_corners, EntityType elem_type, const CartVect &box_half_dims) |
Test if the specified element intersects an axis-aligned box. | |
bool | box_hex_overlap (const CartVect *elem_corners, const CartVect ¢er, const CartVect &dims) |
static bool | box_tet_overlap_edge (const CartVect &dims, const CartVect &edge, const CartVect &ve, const CartVect &v1, const CartVect &v2) |
bool | box_tet_overlap (const CartVect *corners_in, const CartVect ¢er, const CartVect &dims) |
void | closest_location_on_tri (const CartVect &location, const CartVect *vertices, CartVect &closest_out) |
find closest location on triangle | |
void | closest_location_on_tri (const CartVect &location, const CartVect *vertices, double tolerance, CartVect &closest_out, int &closest_topo) |
find closest topological location on triangle | |
void | closest_location_on_polygon (const CartVect &location, const CartVect *vertices, int num_vertices, CartVect &closest_out) |
find closest location on polygon | |
void | closest_location_on_box (const CartVect &min, const CartVect &max, const CartVect &point, CartVect &closest) |
bool | box_point_overlap (const CartVect &box_min_corner, const CartVect &box_max_corner, const CartVect &point, double tolerance) |
bool | boxes_overlap (const CartVect &box_min1, const CartVect &box_max1, const CartVect &box_min2, const CartVect &box_max2, double tolerance) |
bool | bounding_boxes_overlap (const CartVect *list1, int num1, const CartVect *list2, int num2, double tolerance) |
bool | nat_coords_trilinear_hex (const CartVect *corner_coords, const CartVect &x, CartVect &xi, double tol) |
bool | point_in_trilinear_hex (const CartVect *hex, const CartVect &xyz, double etol) |
bool | box_hex_overlap (const CartVect hexv[8], const CartVect &box_center, const CartVect &box_dims) |
bool | box_tet_overlap (const CartVect tet_corners[4], const CartVect &box_center, const CartVect &box_dims) |
Plücker test for intersection between a ray and a triangle.
vertices | Nodes of the triangle. |
ray_point | The start point of the ray. |
ray_unit_direction | The direction of the ray. Must be a unit vector. |
tolerance | Absolute distance tolerance for point equality |
t_out | Output: The distance along the ray from ray_point in the direction of ray_unit_direction at which the ray intersected the triangle. |
nonneg_ray_length | Optional: If non-null, a maximum length for the ray, converting this function to a segment-tri-intersect test. |
neg_ray_length | Optional: If non-null, a maximum length for the ray behind the origin, converting this function to a segment-tri-intersect test. |
orientation | Optional: Reject intersections without the specified orientation of ray with respect to triangle normal vector. Indicate desired orientation by passing 1 (forward), -1 (reverse), or 0 (no preference). |
int_type | Optional Output: The type of intersection; used to identify edge/node intersections. |
Definition at line 106 of file GeomUtil.hpp.
bool moab::GeomUtil::bounding_boxes_overlap | ( | const CartVect * | list1, |
int | num1, | ||
const CartVect * | list2, | ||
int | num2, | ||
double | tolerance | ||
) |
Definition at line 1400 of file GeomUtil.cpp.
{ assert(num1>=1 && num2>=1); CartVect box_min1=list1[0], box_max1=list1[0]; CartVect box_min2=list2[0], box_max2=list2[0]; for (int i=1; i<num1; i++) { for (int k=0; k<3; k++) { double val=list1[i][k]; if (box_min1[k] > val) box_min1[k] = val; if (box_max1[k] < val) box_max1[k]=val; } } for (int i=1; i<num2; i++) { for (int k=0; k<3; k++) { double val=list2[i][k]; if (box_min2[k] > val) box_min2[k] = val; if (box_max2[k] < val) box_max2[k]=val; } } return boxes_overlap(box_min1, box_max1, box_min2, box_max2, tolerance); }
bool moab::GeomUtil::box_elem_overlap | ( | const CartVect * | elem_corners, |
EntityType | elem_type, | ||
const CartVect & | box_center, | ||
const CartVect & | box_half_dims | ||
) |
Test if the specified element intersects an axis-aligned box.
Test if element intersects axis-aligned box. Use element-specific optimization if available, otherwise call box_general_elem_overlap.
elem_corners | The coordinates of the element vertices |
elem_type | The toplogy of the element. |
box_center | The center of the axis-aligned box |
box_half_dims | Half of the width of the box in each axial direction. |
Definition at line 546 of file GeomUtil.cpp.
{ switch (elem_type) { case MBTRI: return box_tri_overlap( elem_corners, center, dims ); case MBTET: return box_tet_overlap( elem_corners, center, dims ); case MBHEX: return box_hex_overlap( elem_corners, center, dims ); case MBPOLYGON: case MBPOLYHEDRON: assert(false); return false; default: return box_linear_elem_overlap( elem_corners, elem_type, center, dims ); } }
bool moab::GeomUtil::box_hex_overlap | ( | const CartVect | hexv[8], |
const CartVect & | box_center, | ||
const CartVect & | box_dims | ||
) |
bool moab::GeomUtil::box_hex_overlap | ( | const CartVect * | elem_corners, |
const CartVect & | center, | ||
const CartVect & | dims | ||
) |
Definition at line 764 of file GeomUtil.cpp.
{ // Do Separating Axis Theorem: // If the element and the box overlap, then the 1D projections // onto at least one of the axes in the following three sets // must overlap (assuming convex polyhedral element). // 1) The normals of the faces of the box (the principal axes) // 2) The crossproduct of each element edge with each box edge // (crossproduct of each edge with each principal axis) // 3) The normals of the faces of the element unsigned i, e, f; // loop counters double dot, cross[2], tmp; CartVect norm; const CartVect corners[8] = { elem_corners[0] - center, elem_corners[1] - center, elem_corners[2] - center, elem_corners[3] - center, elem_corners[4] - center, elem_corners[5] - center, elem_corners[6] - center, elem_corners[7] - center }; // test box face normals (principal axes) int not_less[3] = { 8, 8, 8 }; int not_greater[3] = { 8, 8, 8 }; int not_inside; for (i = 0; i < 8; ++i) { // for each element corner not_inside = 3; if (corners[i][0] < -dims[0]) --not_less[0]; else if (corners[i][0] > dims[0]) --not_greater[0]; else --not_inside; if (corners[i][1] < -dims[1]) --not_less[1]; else if (corners[i][1] > dims[1]) --not_greater[1]; else --not_inside; if (corners[i][2] < -dims[2]) --not_less[2]; else if (corners[i][2] > dims[2]) --not_greater[2]; else --not_inside; if (!not_inside) return true; } // If all points less than min_x of box, then // not_less[0] == 0, and therfore // the following product is zero. if (not_greater[0] * not_greater[1] * not_greater[2] * not_less[0] * not_less[1] * not_less[2] == 0) return false; // Test edge-edge crossproducts const unsigned edges[12][2] = { { 0, 1 }, { 0, 4 }, { 0, 3 }, { 2, 3 }, { 2, 1 }, { 2, 6 }, { 5, 6 }, { 5, 1 }, { 5, 4 }, { 7, 4 }, { 7, 3 }, { 7, 6 } }; // Edge directions for box are principal axis, so // for each element edge, check along the cross-product // of that edge with each of the tree principal axes. for (e = 0; e < 12; ++e) { // for each element edge // get which element vertices bound the edge const CartVect& v0 = corners[ edges[e][0] ]; const CartVect& v1 = corners[ edges[e][1] ]; // X-Axis // calculate crossproduct: axis x (v1 - v0), // where v1 and v0 are edge vertices. cross[0] = v0[2] - v1[2]; cross[1] = v1[1] - v0[1]; // skip if parallel if ((cross[0]*cross[0] + cross[1]*cross[1]) >= std::numeric_limits<double>::epsilon()) { dot = fabs(cross[0] * dims[1]) + fabs(cross[1] * dims[2]); not_less[0] = not_greater[0] = 7; for (i = (edges[e][0]+1)%8; i != edges[e][0]; i = (i+1)%8) { // for each element corner tmp = cross[0] * corners[i][1] + cross[1] * corners[i][2]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } // Y-Axis // calculate crossproduct: axis x (v1 - v0), // where v1 and v0 are edge vertices. cross[0] = v0[0] - v1[0]; cross[1] = v1[2] - v0[2]; // skip if parallel if ((cross[0]*cross[0] + cross[1]*cross[1]) >= std::numeric_limits<double>::epsilon()) { dot = fabs(cross[0] * dims[2]) + fabs(cross[1] * dims[0]); not_less[0] = not_greater[0] = 7; for (i = (edges[e][0]+1)%8; i != edges[e][0]; i = (i+1)%8) { // for each element corner tmp = cross[0] * corners[i][2] + cross[1] * corners[i][0]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } // Z-Axis // calculate crossproduct: axis x (v1 - v0), // where v1 and v0 are edge vertices. cross[0] = v0[1] - v1[1]; cross[1] = v1[0] - v0[0]; // skip if parallel if ((cross[0]*cross[0] + cross[1]*cross[1]) >= std::numeric_limits<double>::epsilon()) { dot = fabs(cross[0] * dims[0]) + fabs(cross[1] * dims[1]); not_less[0] = not_greater[0] = 7; for (i = (edges[e][0]+1)%8; i != edges[e][0]; i = (i+1)%8) { // for each element corner tmp = cross[0] * corners[i][0] + cross[1] * corners[i][1]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } } // test element face normals const unsigned faces[6][4] = { { 0, 1, 2, 3 }, { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 6, 7, 3 }, { 3, 7, 4, 0 }, { 7, 4, 5, 6 } }; for (f = 0; f < 6; ++f) { norm = quad_norm( corners[faces[f][0]], corners[faces[f][1]], corners[faces[f][2]], corners[faces[f][3]] ); dot = dot_abs(norm, dims); // for each element vertex not_less[0] = not_greater[0] = 8; for (i = 0; i < 8; ++i) { tmp = norm % corners[i]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } // Overlap on all tested axes. return true; }
bool moab::GeomUtil::box_linear_elem_overlap | ( | const CartVect * | elem_corners, |
EntityType | elem_type, | ||
const CartVect & | box_center, | ||
const CartVect & | box_half_dims | ||
) |
Test if the specified element intersects an axis-aligned box.
Uses MBCN and separating axis theorem for general algorithm that works for all fixed-size elements (not poly*).
elem_corners | The coordinates of the element vertices |
elem_type | The toplogy of the element. |
box_center | The center of the axis-aligned box |
box_half_dims | Half of the width of the box in each axial direction. |
Definition at line 580 of file GeomUtil.cpp.
{ CartVect corners[8]; const unsigned num_corner = CN::VerticesPerEntity( type ); assert( num_corner <= sizeof(corners)/sizeof(corners[0]) ); for (unsigned i = 0; i < num_corner; ++i) corners[i] = elem_corners[i] - box_center; return box_linear_elem_overlap( corners, type, box_halfdims ); }
bool moab::GeomUtil::box_linear_elem_overlap | ( | const CartVect * | elem_corners, |
EntityType | elem_type, | ||
const CartVect & | box_half_dims | ||
) |
Test if the specified element intersects an axis-aligned box.
Uses MBCN and separating axis theorem for general algorithm that works for all fixed-size elements (not poly*). Box and element vertices must be translated such that box center is at origin.
elem_corners | The coordinates of the element vertices, in local coordinate system of box. |
elem_type | The toplogy of the element. |
box_half_dims | Half of the width of the box in each axial direction. |
Definition at line 594 of file GeomUtil.cpp.
{ // Do Separating Axis Theorem: // If the element and the box overlap, then the 1D projections // onto at least one of the axes in the following three sets // must overlap (assuming convex polyhedral element). // 1) The normals of the faces of the box (the principal axes) // 2) The crossproduct of each element edge with each box edge // (crossproduct of each edge with each principal axis) // 3) The normals of the faces of the element unsigned e, f; // loop counters int i; double dot, cross[2], tmp; CartVect norm; int indices[4]; // element edge/face vertex indices // test box face normals (principal axes) const int num_corner = CN::VerticesPerEntity( type ); int not_less[3] = { num_corner, num_corner, num_corner }; int not_greater[3] = { num_corner, num_corner, num_corner }; int not_inside; for (i = 0; i < num_corner; ++i) { // for each element corner not_inside = 3; if (elem_corners[i][0] < -dims[0]) --not_less[0]; else if (elem_corners[i][0] > dims[0]) --not_greater[0]; else --not_inside; if (elem_corners[i][1] < -dims[1]) --not_less[1]; else if (elem_corners[i][1] > dims[1]) --not_greater[1]; else --not_inside; if (elem_corners[i][2] < -dims[2]) --not_less[2]; else if (elem_corners[i][2] > dims[2]) --not_greater[2]; else --not_inside; if (!not_inside) return true; } // If all points less than min_x of box, then // not_less[0] == 0, and therfore // the following product is zero. if (not_greater[0] * not_greater[1] * not_greater[2] * not_less[0] * not_less[1] * not_less[2] == 0) return false; // Test edge-edge crossproducts // Edge directions for box are principal axis, so // for each element edge, check along the cross-product // of that edge with each of the tree principal axes. const unsigned num_edge = CN::NumSubEntities( type, 1 ); for (e = 0; e < num_edge; ++e) { // for each element edge // get which element vertices bound the edge CN::SubEntityVertexIndices( type, 1, e, indices ); // X-Axis // calculate crossproduct: axis x (v1 - v0), // where v1 and v0 are edge vertices. cross[0] = elem_corners[indices[0]][2] - elem_corners[indices[1]][2]; cross[1] = elem_corners[indices[1]][1] - elem_corners[indices[0]][1]; // skip if parallel if ((cross[0]*cross[0] + cross[1]*cross[1]) >= std::numeric_limits<double>::epsilon()) { dot = fabs(cross[0] * dims[1]) + fabs(cross[1] * dims[2]); not_less[0] = not_greater[0] = num_corner - 1; for (i = (indices[0]+1)%num_corner; i != indices[0]; i = (i+1)%num_corner) { // for each element corner tmp = cross[0] * elem_corners[i][1] + cross[1] * elem_corners[i][2]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } // Y-Axis // calculate crossproduct: axis x (v1 - v0), // where v1 and v0 are edge vertices. cross[0] = elem_corners[indices[0]][0] - elem_corners[indices[1]][0]; cross[1] = elem_corners[indices[1]][2] - elem_corners[indices[0]][2]; // skip if parallel if ((cross[0]*cross[0] + cross[1]*cross[1]) >= std::numeric_limits<double>::epsilon()) { dot = fabs(cross[0] * dims[2]) + fabs(cross[1] * dims[0]); not_less[0] = not_greater[0] = num_corner - 1; for (i = (indices[0]+1)%num_corner; i != indices[0]; i = (i+1)%num_corner) { // for each element corner tmp = cross[0] * elem_corners[i][2] + cross[1] * elem_corners[i][0]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } // Z-Axis // calculate crossproduct: axis x (v1 - v0), // where v1 and v0 are edge vertices. cross[0] = elem_corners[indices[0]][1] - elem_corners[indices[1]][1]; cross[1] = elem_corners[indices[1]][0] - elem_corners[indices[0]][0]; // skip if parallel if ((cross[0]*cross[0] + cross[1]*cross[1]) >= std::numeric_limits<double>::epsilon()) { dot = fabs(cross[0] * dims[0]) + fabs(cross[1] * dims[1]); not_less[0] = not_greater[0] = num_corner - 1; for (i = (indices[0]+1)%num_corner; i != indices[0]; i = (i+1)%num_corner) { // for each element corner tmp = cross[0] * elem_corners[i][0] + cross[1] * elem_corners[i][1]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } } // test element face normals const unsigned num_face = CN::NumSubEntities( type, 2 ); for (f = 0; f < num_face; ++f) { CN::SubEntityVertexIndices( type, 2, f, indices ); switch (CN::SubEntityType( type, 2, f )) { case MBTRI: norm = tri_norm( elem_corners[indices[0]], elem_corners[indices[1]], elem_corners[indices[2]] ); break; case MBQUAD: norm = quad_norm( elem_corners[indices[0]], elem_corners[indices[1]], elem_corners[indices[2]], elem_corners[indices[3]] ); break; default: assert(false); continue; } dot = dot_abs(norm, dims); // for each element vertex not_less[0] = not_greater[0] = num_corner; for (i = 0; i < num_corner; ++i) { tmp = norm % elem_corners[i]; not_less[0] -= (tmp < -dot); not_greater[0] -= (tmp > dot); } if (not_less[0] * not_greater[0] == 0) return false; } // Overlap on all tested axes. return true; }
bool moab::GeomUtil::box_plane_overlap | ( | const CartVect & | plane_normal, |
double | plane_coeff, | ||
CartVect | box_min_corner, | ||
CartVect | box_max_corner | ||
) |
Test if plane intersects axis-aligned box.
Test for intersection between an unbounded plane and an axis-aligned box.
plane_normal | Vector in plane normal direction (need *not* be a unit vector). The N in the plane equation: N . X + D = 0 |
plane_coeff | The scalar 'D' term in the plane equation: N . X + D = 0 |
box_min_corner | The smallest coordinates of the box along each axis. The corner of the box for which all three coordinate values are smaller than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the -X, -Y, and -Z sides respectively. |
box_max_corner | The largest coordinates of the box along each axis. The corner of the box for which all three coordinate values are larger than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the +X, +Y, and +Z sides respectively. |
Definition at line 408 of file GeomUtil.cpp.
{ if (normal[0] < 0.0) std::swap( min[0], max[0] ); if (normal[1] < 0.0) std::swap( min[1], max[1] ); if (normal[2] < 0.0) std::swap( min[2], max[2] ); return (normal % min <= -d) && (normal % max >= -d); }
bool moab::GeomUtil::box_point_overlap | ( | const CartVect & | box_min_corner, |
const CartVect & | box_max_corner, | ||
const CartVect & | point, | ||
double | tolerance | ||
) |
Definition at line 1372 of file GeomUtil.cpp.
{ CartVect closest; closest_location_on_box( box_min_corner, box_max_corner, point, closest ); closest -= point; return closest % closest < tolerance * tolerance; }
bool moab::GeomUtil::box_tet_overlap | ( | const CartVect | tet_corners[4], |
const CartVect & | box_center, | ||
const CartVect & | box_dims | ||
) |
bool moab::GeomUtil::box_tet_overlap | ( | const CartVect * | corners_in, |
const CartVect & | center, | ||
const CartVect & | dims | ||
) |
Definition at line 979 of file GeomUtil.cpp.
{ // Do Separating Axis Theorem: // If the element and the box overlap, then the 1D projections // onto at least one of the axes in the following three sets // must overlap (assuming convex polyhedral element). // 1) The normals of the faces of the box (the principal axes) // 2) The crossproduct of each element edge with each box edge // (crossproduct of each edge with each principal axis) // 3) The normals of the faces of the element // Translate problem such that box center is at origin. const CartVect corners[4] = { corners_in[0] - center, corners_in[1] - center, corners_in[2] - center, corners_in[3] - center }; // 0) Check if any vertex is within the box if (fabs(corners[0][0]) <= dims[0] && fabs(corners[0][1]) <= dims[1] && fabs(corners[0][2]) <= dims[2]) return true; if (fabs(corners[1][0]) <= dims[0] && fabs(corners[1][1]) <= dims[1] && fabs(corners[1][2]) <= dims[2]) return true; if (fabs(corners[2][0]) <= dims[0] && fabs(corners[2][1]) <= dims[1] && fabs(corners[2][2]) <= dims[2]) return true; if (fabs(corners[3][0]) <= dims[0] && fabs(corners[3][1]) <= dims[1] && fabs(corners[3][2]) <= dims[2]) return true; // 1) Check for overlap on each principal axis (box face normal) // X if (corners[0][0] < -dims[0] && corners[1][0] < -dims[0] && corners[2][0] < -dims[0] && corners[3][0] < -dims[0]) return false; if (corners[0][0] > dims[0] && corners[1][0] > dims[0] && corners[2][0] > dims[0] && corners[3][0] > dims[0]) return false; // Y if (corners[0][1] < -dims[1] && corners[1][1] < -dims[1] && corners[2][1] < -dims[1] && corners[3][1] < -dims[1]) return false; if (corners[0][1] > dims[1] && corners[1][1] > dims[1] && corners[2][1] > dims[1] && corners[3][1] > dims[1]) return false; // Z if (corners[0][2] < -dims[2] && corners[1][2] < -dims[2] && corners[2][2] < -dims[2] && corners[3][2] < -dims[2]) return false; if (corners[0][2] > dims[2] && corners[1][2] > dims[2] && corners[2][2] > dims[2] && corners[3][2] > dims[2]) return false; // 3) test element face normals CartVect norm; double dot, dot1, dot2; const CartVect v01 = corners[1] - corners[0]; const CartVect v02 = corners[2] - corners[0]; norm = v01 * v02; dot = dot_abs(norm, dims); dot1 = norm % corners[0]; dot2 = norm % corners[3]; if (dot1 > dot2) std::swap(dot1, dot2); if (dot2 < -dot || dot1 > dot) return false; const CartVect v03 = corners[3] - corners[0]; norm = v03 * v01; dot = dot_abs(norm, dims); dot1 = norm % corners[0]; dot2 = norm % corners[2]; if (dot1 > dot2) std::swap(dot1, dot2); if (dot2 < -dot || dot1 > dot) return false; norm = v02 * v03; dot = dot_abs(norm, dims); dot1 = norm % corners[0]; dot2 = norm % corners[1]; if (dot1 > dot2) std::swap(dot1, dot2); if (dot2 < -dot || dot1 > dot) return false; const CartVect v12 = corners[2] - corners[1]; const CartVect v13 = corners[3] - corners[1]; norm = v13 * v12; dot = dot_abs(norm, dims); dot1 = norm % corners[0]; dot2 = norm % corners[1]; if (dot1 > dot2) std::swap(dot1, dot2); if (dot2 < -dot || dot1 > dot) return false; // 2) test edge-edge cross products const CartVect v23 = corners[3] - corners[2]; return box_tet_overlap_edge( dims, v01, corners[0], corners[2], corners[3] ) && box_tet_overlap_edge( dims, v02, corners[0], corners[1], corners[3] ) && box_tet_overlap_edge( dims, v03, corners[0], corners[1], corners[2] ) && box_tet_overlap_edge( dims, v12, corners[1], corners[0], corners[3] ) && box_tet_overlap_edge( dims, v13, corners[1], corners[0], corners[2] ) && box_tet_overlap_edge( dims, v23, corners[2], corners[0], corners[1] ); }
static bool moab::GeomUtil::box_tet_overlap_edge | ( | const CartVect & | dims, |
const CartVect & | edge, | ||
const CartVect & | ve, | ||
const CartVect & | v1, | ||
const CartVect & | v2 | ||
) | [inline, static] |
Definition at line 935 of file GeomUtil.cpp.
{ double dot, dot1, dot2, dot3, min, max; // edge x X if (fabs(edge[1]*edge[2]) > std::numeric_limits<double>::epsilon()) { dot = fabs(edge[2]) * dims[1] + fabs(edge[1]) * dims[2]; dot1 = edge[2] * ve[1] - edge[1] * ve[2]; dot2 = edge[2] * v1[1] - edge[1] * v1[2]; dot3 = edge[2] * v2[1] - edge[1] * v2[2]; min_max_3( dot1, dot2, dot3, min, max ); if (max < -dot || min > dot) return false; } // edge x Y if (fabs(edge[1]*edge[2]) > std::numeric_limits<double>::epsilon()) { dot = fabs(edge[2]) * dims[0] + fabs(edge[0]) * dims[2]; dot1 = -edge[2] * ve[0] + edge[0] * ve[2]; dot2 = -edge[2] * v1[0] + edge[0] * v1[2]; dot3 = -edge[2] * v2[0] + edge[0] * v2[2]; min_max_3( dot1, dot2, dot3, min, max ); if (max < -dot || min > dot) return false; } // edge x Z if (fabs(edge[1]*edge[2]) > std::numeric_limits<double>::epsilon()) { dot = fabs(edge[1]) * dims[0] + fabs(edge[0]) * dims[1]; dot1 = edge[1] * ve[0] - edge[0] * ve[1]; dot2 = edge[1] * v1[0] - edge[0] * v1[1]; dot3 = edge[1] * v2[0] - edge[0] * v2[1]; min_max_3( dot1, dot2, dot3, min, max ); if (max < -dot || min > dot) return false; } return true; }
bool moab::GeomUtil::box_tri_overlap | ( | const CartVect | triangle_corners[3], |
const CartVect & | box_center, | ||
const CartVect & | box_half_dims | ||
) |
Test if triangle intersects axis-aligned box.
Test if a triangle intersects an axis-aligned box.
triangle_corners | The corners of the triangle. |
box_center | The center of the box. |
box_hanf_dims | The distance along each axis, respectively, from the box_center to the boundary of the box. |
Definition at line 441 of file GeomUtil.cpp.
{ // translate everthing such that box is centered at origin const CartVect v0( vertices[0] - box_center ); const CartVect v1( vertices[1] - box_center ); const CartVect v2( vertices[2] - box_center ); // do case 1) tests if (v0[0] > box_dims[0] && v1[0] > box_dims[0] && v2[0] > box_dims[0]) return false; if (v0[1] > box_dims[1] && v1[1] > box_dims[1] && v2[1] > box_dims[1]) return false; if (v0[2] > box_dims[2] && v1[2] > box_dims[2] && v2[2] > box_dims[2]) return false; if (v0[0] < -box_dims[0] && v1[0] < -box_dims[0] && v2[0] < -box_dims[0]) return false; if (v0[1] < -box_dims[1] && v1[1] < -box_dims[1] && v2[1] < -box_dims[1]) return false; if (v0[2] < -box_dims[2] && v1[2] < -box_dims[2] && v2[2] < -box_dims[2]) return false; // compute triangle edge vectors const CartVect e0( vertices[1] - vertices[0] ); const CartVect e1( vertices[2] - vertices[1] ); const CartVect e2( vertices[0] - vertices[2] ); // do case 3) tests double fex, fey, fez, p0, p1, p2, rad; fex = fabs(e0[0]); fey = fabs(e0[1]); fez = fabs(e0[2]); p0 = e0[2]*v0[1] - e0[1]*v0[2]; p2 = e0[2]*v2[1] - e0[1]*v2[2]; rad = fez * box_dims[1] + fey * box_dims[2]; CHECK_RANGE( p0, p2, rad ); p0 = -e0[2]*v0[0] + e0[0]*v0[2]; p2 = -e0[2]*v2[0] + e0[0]*v2[2]; rad = fez * box_dims[0] + fex * box_dims[2]; CHECK_RANGE( p0, p2, rad ); p1 = e0[1]*v1[0] - e0[0]*v1[1]; p2 = e0[1]*v2[0] - e0[0]*v2[1]; rad = fey * box_dims[0] + fex * box_dims[1]; CHECK_RANGE( p1, p2, rad ); fex = fabs(e1[0]); fey = fabs(e1[1]); fez = fabs(e1[2]); p0 = e1[2]*v0[1] - e1[1]*v0[2]; p2 = e1[2]*v2[1] - e1[1]*v2[2]; rad = fez * box_dims[1] + fey * box_dims[2]; CHECK_RANGE( p0, p2, rad ); p0 = -e1[2]*v0[0] + e1[0]*v0[2]; p2 = -e1[2]*v2[0] + e1[0]*v2[2]; rad = fez * box_dims[0] + fex * box_dims[2]; CHECK_RANGE( p0, p2, rad ); p0 = e1[1]*v0[0] - e1[0]*v0[1]; p1 = e1[1]*v1[0] - e1[0]*v1[1]; rad = fey * box_dims[0] + fex * box_dims[1]; CHECK_RANGE( p0, p1, rad ); fex = fabs(e2[0]); fey = fabs(e2[1]); fez = fabs(e2[2]); p0 = e2[2]*v0[1] - e2[1]*v0[2]; p1 = e2[2]*v1[1] - e2[1]*v1[2]; rad = fez * box_dims[1] + fey * box_dims[2]; CHECK_RANGE( p0, p1, rad ); p0 = -e2[2]*v0[0] + e2[0]*v0[2]; p1 = -e2[2]*v1[0] + e2[0]*v1[2]; rad = fez * box_dims[0] + fex * box_dims[2]; CHECK_RANGE( p0, p1, rad ); p1 = e2[1]*v1[0] - e2[0]*v1[1]; p2 = e2[1]*v2[0] - e2[0]*v2[1]; rad = fey * box_dims[0] + fex * box_dims[1]; CHECK_RANGE( p1, p2, rad ); // do case 2) test CartVect n = e0 * e1; return box_plane_overlap( n, -(n % v0), -box_dims, box_dims ); }
bool moab::GeomUtil::box_tri_overlap | ( | const CartVect | triangle_corners[3], |
const CartVect & | box_min_corner, | ||
const CartVect & | box_max_corner, | ||
double | tolerance | ||
) |
Test if triangle intersects axis-aligned box.
Test if a triangle intersects an axis-aligned box.
triangle_corners | The corners of the triangle. |
box_min_corner | The smallest coordinates of the box along each axis. The corner of the box for which all three coordinate values are smaller than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the -X, -Y, and -Z sides respectively. |
box_max_corner | The largest coordinates of the box along each axis. The corner of the box for which all three coordinate values are larger than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the +X, +Y, and +Z sides respectively. |
tolerance | The tolerance used in the intersection test. The box size is increased by this amount before the intersection test. |
Definition at line 534 of file GeomUtil.cpp.
{ const CartVect box_center = 0.5 * (box_max_corner + box_min_corner); const CartVect box_hf_dim = 0.5 * (box_max_corner - box_min_corner); return box_tri_overlap( triangle_corners, box_center, box_hf_dim + CartVect(tolerance) ); }
bool moab::GeomUtil::boxes_overlap | ( | const CartVect & | box_min1, |
const CartVect & | box_max1, | ||
const CartVect & | box_min2, | ||
const CartVect & | box_max2, | ||
double | tolerance | ||
) |
Definition at line 1383 of file GeomUtil.cpp.
void moab::GeomUtil::closest_location_on_box | ( | const CartVect & | min, |
const CartVect & | max, | ||
const CartVect & | point, | ||
CartVect & | closest | ||
) |
Definition at line 1362 of file GeomUtil.cpp.
{ closest[0] = point[0] < min[0] ? min[0] : point[0] > max[0] ? max[0] : point[0]; closest[1] = point[1] < min[1] ? min[1] : point[1] > max[1] ? max[1] : point[1]; closest[2] = point[2] < min[2] ? min[2] : point[2] > max[2] ? max[2] : point[2]; }
void moab::GeomUtil::closest_location_on_polygon | ( | const CartVect & | location, |
const CartVect * | vertices, | ||
int | num_vertices, | ||
CartVect & | closest_out | ||
) |
find closest location on polygon
Find closest location on polygon
location | Input position to evaluate from |
vertices | Array of corner vertex coordinates. |
num_vertices | Length of 'vertices' array. |
closest_out | Result position |
Definition at line 1296 of file GeomUtil.cpp.
{ const int n = num_vertices; CartVect d, p, v; double shortest_sqr, dist_sqr, t_closest, t; int i, e; // Find closest edge of polygon. e = n - 1; v = vertices[0] - vertices[e]; t_closest = (v % (location - vertices[e])) / (v % v); if (t_closest < 0.0) d = location - vertices[e]; else if (t_closest > 1.0) d = location - vertices[0]; else d = location - vertices[e] - t_closest * v; shortest_sqr = d % d; for (i = 0; i < n - 1; ++i) { v = vertices[i+1] - vertices[i]; t = (v % (location - vertices[i])) / (v % v); if (t < 0.0) d = location - vertices[i]; else if (t > 1.0) d = location - vertices[i+1]; else d = location - vertices[i] - t * v; dist_sqr = d % d; if (dist_sqr < shortest_sqr) { e = i; shortest_sqr = dist_sqr; t_closest = t; } } // If we are beyond the bounds of the edge, then // the point is outside and closest to a vertex if (t_closest <= 0.0) { closest_out = vertices[e]; return; } else if (t_closest >= 1.0) { closest_out = vertices[(e+1)%n]; return; } // Now check which side of the edge we are one const CartVect v0 = vertices[e] - vertices[(e+n-1)%n]; const CartVect v1 = vertices[(e+1)%n] - vertices[e]; const CartVect v2 = vertices[(e+2)%n] - vertices[(e+1)%n]; const CartVect norm = (1.0 - t_closest) * (v0 * v1) + t_closest * (v1 * v2); // if on outside of edge, result is closest point on edge if ((norm % ((vertices[e] - location) * v1)) <= 0.0) { closest_out = vertices[e] + t_closest * v1; return; } // Inside. Project to plane defined by point and normal at // closest edge const double D = -(norm % (vertices[e] + t_closest * v1)); closest_out = (location - (norm % location + D) * norm)/(norm % norm); }
void moab::GeomUtil::closest_location_on_tri | ( | const CartVect & | location, |
const CartVect * | vertices, | ||
CartVect & | closest_out | ||
) |
find closest location on triangle
Find closest location on linear triangle.
location | Input position to evaluate from |
vertices | Array of three corner vertex coordinates. |
closest_out | Result position |
Definition at line 1133 of file GeomUtil.cpp.
{ // ops comparisons const CartVect sv( vertices[1] - vertices[0] ); // +3 = 3 const CartVect tv( vertices[2] - vertices[0] ); // +3 = 6 const CartVect pv( vertices[0] - location ); // +3 = 9 const double ss = sv % sv; // +5 = 14 const double st = sv % tv; // +5 = 19 const double tt = tv % tv; // +5 = 24 const double sp = sv % pv; // +5 = 29 const double tp = tv % pv; // +5 = 34 const double det = ss*tt - st*st; // +3 = 37 double s = st*tp - tt*sp; // +3 = 40 double t = st*sp - ss*tp; // +3 = 43 if (s+t < det) { // +1 = 44, +1 = 1 if (s < 0) { // +1 = 2 if (t < 0) { // +1 = 3 // region 4 if (sp < 0) { // +1 = 4 if (-sp > ss) // +1 = 5 closest_out = vertices[1]; // 44 5 else closest_out = vertices[0] - (sp/ss) * sv; // +7 = 51, 5 } else if (tp < 0) { // +1 = 5 if (-tp > tt) // +1 = 6 closest_out = vertices[2]; // 44, 6 else closest_out = vertices[0] - (tp/tt) * tv; // +7 = 51, 6 } else { closest_out = vertices[0]; // 44, 5 } } else { // region 3 if (tp >= 0) // +1 = 4 closest_out = vertices[0]; // 44, 4 else if (-tp >= tt) // +1 = 5 closest_out = vertices[2]; // 44, 5 else closest_out = vertices[0] - (tp/tt) * tv; // +7 = 51, 5 } } else if (t < 0) { // +1 = 3 // region 5; if (sp >= 0.0) // +1 = 4 closest_out = vertices[0]; // 44, 4 else if (-sp >= ss) // +1 = 5 closest_out = vertices[1]; // 44 5 else closest_out = vertices[0] - (sp/ss) * sv; // +7 = 51, 5 } else { // region 0 const double inv_det = 1.0 / det; // +1 = 45 s *= inv_det; // +1 = 46 t *= inv_det; // +1 = 47 closest_out = vertices[0] + s*sv + t*tv; //+12 = 59, 3 } } else { if (s < 0) { // +1 = 2 // region 2 s = st + sp; // +1 = 45 t = tt + tp; // +1 = 46 if (t > s) { // +1 = 3 const double num = t - s; // +1 = 47 const double den = ss - 2*st + tt; // +3 = 50 if (num > den) // +1 = 4 closest_out = vertices[1]; // 50, 4 else { s = num/den; // +1 = 51 t = 1 - s; // +1 = 52 closest_out = s*vertices[1] + t*vertices[2];// +9 = 61, 4 } } else if (t <= 0) // +1 = 4 closest_out = vertices[2]; // 46, 4 else if (tp >= 0) // +1 = 5 closest_out = vertices[0]; // 46, 5 else closest_out = vertices[0] - (tp/tt) * tv; // +7 = 53, 5 } else if (t < 0) { // +1 = 3 // region 6 t = st + tp; // +1 = 45 s = ss + sp; // +1 = 46 if (s > t) { // +1 = 4 const double num = t - s; // +1 = 47 const double den = tt - 2*st + ss; // +3 = 50 if (num > den) // +1 = 5 closest_out = vertices[2]; // 50, 5 else { t = num/den; // +1 = 51 s = 1 - t; // +1 = 52 closest_out = s*vertices[1] + t*vertices[2];// +9 = 61, 5 } } else if (s <= 0) // +1 = 5 closest_out = vertices[1]; // 46, 5 else if (sp >= 0) // +1 = 6 closest_out = vertices[0]; // 46, 6 else closest_out = vertices[0] - (sp/ss) * sv; // +7 = 53, 6 } else { // region 1 const double num = tt + tp - st - sp; // +3 = 47 if (num <= 0) { // +1 = 4 closest_out = vertices[2]; // 47, 4 } else { const double den = ss - 2*st + tt; // +3 = 50 if (num >= den) // +1 = 5 closest_out = vertices[1]; // 50, 5 else { s = num/den; // +1 = 51 t = 1 - s; // +1 = 52 closest_out = s*vertices[1] + t*vertices[2];// +9 = 61, 5 } } } } }
void moab::GeomUtil::closest_location_on_tri | ( | const CartVect & | location, |
const CartVect * | vertices, | ||
double | tolerance, | ||
CartVect & | closest_out, | ||
int & | closest_topo | ||
) |
find closest topological location on triangle
Find closest location on linear triangle.
location | Input position to evaluate from |
vertices | Array of three corner vertex coordinates. |
tolerance | Tolerance to use when comparing to corners and edges |
closest_out | Result position |
closest_topo | Closest topological entity 0-2 : vertex index 3-5 : edge beginning at closest_topo - 3 6 : triangle interior |
Definition at line 1260 of file GeomUtil.cpp.
{ const double tsqr = tolerance*tolerance; int i; CartVect pv[3], ev, ep; double t; closest_location_on_tri( location, vertices, closest_out ); for (i = 0; i < 3; ++i) { pv[i] = vertices[i] - closest_out; if ((pv[i] % pv[i]) <= tsqr) { closest_topo = i; return; } } for (i = 0; i < 3; ++i) { ev = vertices[(i+1)%3] - vertices[i]; t = (ev % pv[i]) / (ev % ev); ep = closest_out - (vertices[i] + t * ev); if ((ep % ep) <= tsqr) { closest_topo = i+3; return; } } closest_topo = 6; }
static double moab::GeomUtil::dot_abs | ( | const CartVect & | u, |
const CartVect & | v | ||
) | [inline, static] |
Definition at line 64 of file GeomUtil.cpp.
{ return fabs(u[0]*v[0]) + fabs(u[1]*v[1]) + fabs(u[2]*v[2]); }
bool moab::GeomUtil::first | ( | const CartVect & | a, |
const CartVect & | b | ||
) | [inline] |
Definition at line 111 of file GeomUtil.cpp.
{ if(a[0] < b[0]) { return true; } else if(a[0] == b[0]) { if(a[1] < b[1]) { return true; } else if(a[1] == b[1]) { if(a[2] < b[2]) { return true; } else { return false; } } else { return false; } } else { return false; } }
static void moab::GeomUtil::min_max_3 | ( | double | a, |
double | b, | ||
double | c, | ||
double & | min, | ||
double & | max | ||
) | [inline, static] |
Definition at line 41 of file GeomUtil.cpp.
{ if (a < b) { if (a < c) { min = a; max = b > c ? b : c; } else { min = c; max = b; } } else if (b < c) { min = b; max = a > c ? a : c; } else { min = c; max = a; } }
bool moab::GeomUtil::nat_coords_trilinear_hex | ( | const CartVect * | corner_coords, |
const CartVect & | x, | ||
CartVect & | xi, | ||
double | tol | ||
) |
Definition at line 1521 of file GeomUtil.cpp.
{
return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
}
bool moab::GeomUtil::plucker_ray_tri_intersect | ( | const CartVect | vertices[3], |
const CartVect & | origin, | ||
const CartVect & | direction, | ||
double | , | ||
double & | dist_out, | ||
const double * | nonneg_ray_len, | ||
const double * | neg_ray_len, | ||
const int * | orientation, | ||
intersection_type * | type | ||
) |
Definition at line 144 of file GeomUtil.cpp.
{ const CartVect raya = direction; const CartVect rayb = direction*origin; // edge 0 double pip0; if(first(vertices[0],vertices[1])) { const CartVect edge0a = vertices[1]-vertices[0]; const CartVect edge0b = edge0a*vertices[0]; pip0 = raya % edge0b + rayb % edge0a; } else { const CartVect edge0a = vertices[0]-vertices[1]; const CartVect edge0b = edge0a*vertices[1]; pip0 = raya % edge0b + rayb % edge0a; pip0 = -pip0; } // try to exit early if(orientation && (*orientation)*pip0 > 0) { if(type) *type = NONE; return false; } // edge 1 double pip1; if(first(vertices[1],vertices[2])) { const CartVect edge1a = vertices[2]-vertices[1]; const CartVect edge1b = edge1a*vertices[1]; pip1 = raya % edge1b + rayb % edge1a; } else { const CartVect edge1a = vertices[1]-vertices[2]; const CartVect edge1b = edge1a*vertices[2]; pip1 = raya % edge1b + rayb % edge1a; pip1 = -pip1; } // try to exit early if(orientation) { if( (*orientation)*pip1 > 0) { if(type) *type = NONE; return false; } // If the orientation is not specified, all pips must be the same sign or zero. } else if( (0.0<pip0 && 0.0>pip1) || (0.0>pip0 && 0.0<pip1) ) { if(type) *type = NONE; return false; } // edge 2 double pip2; if(first(vertices[2],vertices[0])) { const CartVect edge2a = vertices[0]-vertices[2]; const CartVect edge2b = edge2a*vertices[2]; pip2 = raya % edge2b + rayb % edge2a; } else { const CartVect edge2a = vertices[2]-vertices[0]; const CartVect edge2b = edge2a*vertices[0]; pip2 = raya % edge2b + rayb % edge2a; pip2 = -pip2; } // try to exit early if(orientation) { if( (*orientation)*pip2 > 0) { if(type) *type = NONE; return false; } // If the orientation is not specified, all pips must be the same sign or zero. } else if( (0.0<pip1 && 0.0>pip2) || (0.0>pip1 && 0.0<pip2) || (0.0<pip0 && 0.0>pip2) || (0.0>pip0 && 0.0<pip2) ) { if(type) *type = NONE; return false; } // check for coplanar case to avoid dividing by zero if(0==pip0 && 0==pip1 && 0==pip2) { //std::cout << "plucker: coplanar" << std::endl; if(type) *type = NONE; return false; } // get the distance to intersection const double inverse_sum = 1.0/(pip0+pip1+pip2); assert(0.0 != inverse_sum); const CartVect intersection(pip0*inverse_sum*vertices[2]+ pip1*inverse_sum*vertices[0]+ pip2*inverse_sum*vertices[1]); // To minimize numerical error, get index of largest magnitude direction. int idx = 0; double max_abs_dir = 0; for(unsigned int i=0; i<3; ++i) { if( fabs(direction[i]) > max_abs_dir ) { idx = i; max_abs_dir = fabs(direction[i]); } } const double dist = (intersection[idx]-origin[idx])/direction[idx]; // is the intersection within distance limits? if(nonneg_ray_len && *nonneg_ray_len<dist) { if(type) *type = NONE; return false; } if(neg_ray_len && *neg_ray_len>=dist) { if(type) *type = NONE; return false; // Unless a neg_ray_len is used, don't return negative distances } else if (!neg_ray_len && 0>dist) { if(type) *type = NONE; return false; } dist_out = dist; // check for special cases if(0==pip0 || 0==pip1 || 0==pip2) { if (0==pip0 && 0==pip1) { //std::cout << "plucker: node1" << std::endl; if(type) *type = NODE1; return true; } else if(0==pip1 && 0==pip2) { //std::cout << "plucker: node2" << std::endl; if(type) *type = NODE2; return true; } else if(0==pip2 && 0==pip0) { //std::cout << "plucker: node0" << std::endl; if(type) *type = NODE0; return true; } else if(0==pip0) { //std::cout << "plucker: edge0" << std::endl; if(type) *type = EDGE0; return true; } else if(0==pip1) { //std::cout << "plucker: edge1" << std::endl; if(type) *type = EDGE1; return true; } else if(0==pip2) { //std::cout << "plucker: edge2" << std::endl; if(type) *type = EDGE2; return true; } } // if here, ray intersects interior of tri if(type) *type = INTERIOR; return true; }
bool moab::GeomUtil::point_in_trilinear_hex | ( | const CartVect * | hex, |
const CartVect & | xyz, | ||
double | etol | ||
) |
Definition at line 1530 of file GeomUtil.cpp.
{ CartVect xi; return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && fabs(xi[0])-1 < etol && fabs(xi[1])-1 < etol && fabs(xi[2])-1 < etol; }
static CartVect moab::GeomUtil::quad_norm | ( | const CartVect & | v1, |
const CartVect & | v2, | ||
const CartVect & | v3, | ||
const CartVect & | v4 | ||
) | [inline, static] |
Definition at line 568 of file GeomUtil.cpp.
{ return (-v1+v2+v3-v4) * (-v1-v2+v3+v4); }
bool moab::GeomUtil::ray_box_intersect | ( | const CartVect & | box_min, |
const CartVect & | box_max, | ||
const CartVect & | ray_pt, | ||
const CartVect & | ray_dir, | ||
double & | t_enter, | ||
double & | t_exit | ||
) |
Find range of overlap between ray and axis-aligned box.
box_min | Box corner with minimum coordinate values |
box_max | Box corner with minimum coordinate values |
ray_pt | Coordinates of start point of ray |
ray_dir | Directionion vector for ray such that the ray is defined by r(t) = ray_point + t * ray_vect for t > 0. |
t_enter | Output: if return value is true, this value is the parameter location along the ray at which the ray entered the leaf. If return value is false, then this value is undefined. |
t_exit | Output: if return value is true, this value is the parameter location along the ray at which the ray exited the leaf. If return value is false, then this value is undefined. |
Definition at line 355 of file GeomUtil.cpp.
{ const double epsilon = 1e-12; double t1, t2; // Use 'slabs' method from 13.6.1 of Akenine-Moller t_enter = 0.0; t_exit = std::numeric_limits<double>::infinity(); // Intersect with each pair of axis-aligned planes bounding // opposite faces of the leaf box bool ray_is_valid = false; // is ray direction vector zero? for (int axis = 0; axis < 3; ++axis) { if (fabs(ray_dir[axis]) < epsilon) { // ray parallel to planes if (ray_pt[axis] >= box_min[axis] && ray_pt[axis] <= box_max[axis]) continue; else return false; } // find t values at which ray intersects each plane ray_is_valid = true; t1 = (box_min[axis] - ray_pt[axis]) / ray_dir[axis]; t2 = (box_max[axis] - ray_pt[axis]) / ray_dir[axis]; // t_enter = max( t_enter_x, t_enter_y, t_enter_z ) // t_exit = min( t_exit_x, t_exit_y, t_exit_z ) // where // t_enter_x = min( t1_x, t2_x ); // t_exit_x = max( t1_x, t2_x ) if (t1 < t2) { if (t_enter < t1) t_enter = t1; if (t_exit > t2) t_exit = t2; } else { if (t_enter < t2) t_enter = t2; if (t_exit > t1) t_exit = t1; } } return ray_is_valid && (t_enter <= t_exit); }
bool moab::GeomUtil::ray_tri_intersect | ( | const CartVect | vertices[3], |
const CartVect & | ray_point, | ||
const CartVect & | ray_unit_direction, | ||
double | tolerance, | ||
double & | t_out, | ||
const double * | ray_length = 0 |
||
) |
Test for intersection between a ray and a triangle.
ray_point | The start point of the ray. |
ray_unit_direciton | The direction of the ray. Must be a unit vector. |
tolerance | Absolute distance tolerance for point equality |
t_out | Output: The distance along the ray from ray_point in the direction of ray_unit_direction at which the ray itersected the triangle. |
ray_length | Optional: If non-null, a pointer a maximum length for the ray, converting this function to a segment-tri- intersect test. |
Definition at line 303 of file GeomUtil.cpp.
{ const CartVect p0 = vertices[0] - vertices[1]; // abc const CartVect p1 = vertices[0] - vertices[2]; // def // ghi<-v const CartVect p = vertices[0] - b; // jkl const CartVect c = p1 * v; // eiMinushf,gfMinusdi,dhMinuseg const double mP = p0 % c; const double betaP = p % c; if (mP > 0) { if (betaP < 0) return false; } else if (mP < 0) { if (betaP > 0) return false; } else { return false; } const CartVect d = p0 * p; // jcMinusal,blMinuskc,akMinusjb double gammaP = v % d; if (mP > 0) { if (gammaP < 0 || betaP + gammaP > mP) return false; } else if (betaP + gammaP < mP || gammaP > 0) return false; const double tP = p1 % d; const double m = 1.0 / mP; const double beta = betaP * m; const double gamma = gammaP * m; const double t = -tP * m; if (ray_length && t > *ray_length) return false; if (beta < 0 || gamma < 0 || beta + gamma > 1 || t < 0.0) return false; t_out = t; return true; }
bool moab::GeomUtil::segment_box_intersect | ( | CartVect | box_min, |
CartVect | box_max, | ||
const CartVect & | seg_pt, | ||
const CartVect & | seg_unit_dir, | ||
double & | seg_start, | ||
double & | seg_end | ||
) |
Given a line segment and an axis-aligned box, return the sub-segment of the line segment that itersects the box.
Can be used to intersect ray with box by passing seg_end as HUGE_VAL or std::numeric_limits<double>::maximum().
box_min | Minimum corner of axis-aligned box |
box_max | Maximum corner of axis-aligned box |
seg_pt | A point in the line containing the segement |
seg_unit_dir | A unit vector in the direction of the line containing the semgent. |
seg_start | The distance from seg_pt in the direction of seg_unit_dir at which the segment begins. As input, the start of the original segment, as output, the start of the sub-segment intersecting the box. Note: seg_start must be less than seg_end |
seg_end | The distance from seg_pt in the direction of seg_unit_dir at which the segment ends. As input, the end of the original segment, as output, the end of the sub-segment intersecting the box. Note: seg_start must be less than seg_end |
Definition at line 67 of file GeomUtil.cpp.
{ // translate so that seg_pt is at origin box_min -= seg_pt; box_max -= seg_pt; for (unsigned i = 0; i < 3; ++i) { // X, Y, and Z slabs // intersect line with slab planes const double t_min = box_min[i] / seg_unit_dir[i]; const double t_max = box_max[i] / seg_unit_dir[i]; // check if line is parallel to planes if (!finite(t_min)) { if (box_min[i] > 0.0 || box_max[i] < 0.0) return false; continue; } if (seg_unit_dir[i] < 0) { if (t_min < seg_end) seg_end = t_min; if (t_max > seg_start) seg_start = t_max; } else { // seg_unit_dir[i] > 0 if (t_min > seg_start) seg_start = t_min; if (t_max < seg_end) seg_end = t_max; } } return seg_start <= seg_end; }
static CartVect moab::GeomUtil::tri_norm | ( | const CartVect & | v1, |
const CartVect & | v2, | ||
const CartVect & | v3 | ||
) | [inline, static] |
Definition at line 574 of file GeomUtil.cpp.
{ return (v2-v1) * (v3-v1); }