Skip to content

Commit

Permalink
Add closest surface point search capabilities
Browse files Browse the repository at this point in the history
  • Loading branch information
@nmoehrle
nmoehrle committed Jul 5, 2016
1 parent 97b12b7 commit dcf3b3e
Show file tree
Hide file tree
Showing 3 changed files with 185 additions and 40 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ libacc
-------------------------------------------------------------------------------
This project offers a header only acceleration structure library including
implementations for a BVH- and KD-Tree. Applications may include ray
intersection tests or nearest neighbor searches.
intersection tests, closest surface point or nearest neighbor searches.

Requirements
-------------------------------------------------------------------------------
Expand Down
103 changes: 91 additions & 12 deletions bvh_tree.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@ class BVHTree {
std::atomic<int> * num_threads);

bool intersect(Ray const & ray, typename Node::ID node_id, Hit * hit) const;
Vec3fType closest_point(Vec3fType vertex, typename Node::ID node_id) const;

public:
static
Expand All @@ -109,6 +110,7 @@ class BVHTree {
int max_threads = std::thread::hardware_concurrency());

bool intersect(Ray ray, Hit * hit_ptr = nullptr) const;
Vec3fType closest_point(Vec3fType vertex, float max_dist = inf) const;
};


Expand Down Expand Up @@ -166,7 +168,7 @@ BVHTree<IdxType, Vec3fType>::bsplit(typename Node::ID node_id,
std::array<AABB, NUM_BINS> right_aabbs;
std::vector<IdxType> bin(n);

float min_cost = std::numeric_limits<float>::infinity();
float min_cost = inf;
std::pair<IdxType, char> split;
for (char d = 0; d < 3; ++d) {
float min = node.aabb.min[d];
Expand Down Expand Up @@ -258,7 +260,7 @@ BVHTree<IdxType, Vec3fType>::ssplit(typename Node::ID node_id, std::vector<AABB>
Node & node = nodes[node_id];
IdxType n = node.last - node.first;

float min_cost = std::numeric_limits<float>::infinity();
float min_cost = inf;
std::pair<char, IdxType> split;
std::vector<AABB> right_aabbs(n);
for (char d = 0; d < 3; ++d) {
Expand Down Expand Up @@ -364,15 +366,15 @@ BVHTree<IdxType, Vec3fType>::intersect(Ray const & ray, typename Node::ID node_i
}
return ret;
}

template <typename IdxType, typename Vec3fType> bool
BVHTree<IdxType, Vec3fType>::intersect(Ray ray, Hit * hit_ptr) const {
Hit hit;
hit.t = std::numeric_limits<float>::infinity();
std::stack<typename Node::ID> s;
hit.t = inf;

s.push(0);
while (!s.empty()) {
typename Node::ID node_id = s.top(); s.pop();
typename Node::ID node_id = 0;
std::stack<typename Node::ID> s;
while (true) {
Node const & node = nodes[node_id];
if (node.left != NAI && node.right != NAI) {
float tmin_left, tmin_right;
Expand All @@ -381,23 +383,31 @@ BVHTree<IdxType, Vec3fType>::intersect(Ray ray, Hit * hit_ptr) const {
if (left && right) {
if (tmin_left < tmin_right) {
s.push(node.right);
s.push(node.left);
node_id = node.left;
} else {
s.push(node.left);
s.push(node.right);
node_id = node.right;
}
} else {
if (right) s.push(node.right);
if (left) s.push(node.left);
if (right) node_id = node.right;
if (left) node_id = node.left;
}

if (!left && !right) {
if (s.empty()) break;
node_id = s.top(); s.pop();
}
} else {
if (intersect(ray, node_id, &hit)) {
ray.tmax = hit.t;
}

if (s.empty()) break;
node_id = s.top(); s.pop();
}
}

if (hit.t < std::numeric_limits<float>::infinity()) {
if (hit.t < inf) {
if (hit_ptr != nullptr) {
*hit_ptr = hit;
}
Expand All @@ -407,6 +417,75 @@ BVHTree<IdxType, Vec3fType>::intersect(Ray ray, Hit * hit_ptr) const {
}
}

template <typename IdxType, typename Vec3fType> Vec3fType
BVHTree<IdxType, Vec3fType>::closest_point(Vec3fType vertex, typename Node::ID node_id) const {
Node const & node = nodes[node_id];

Vec3fType closest;
float dist = inf;

for (std::size_t i = node.first; i < node.last; ++i) {
Vec3fType closest_tri = acc::closest_point(vertex, tris[i]);
float dist_tri = (closest_tri - vertex).square_norm();
if (dist_tri < dist) {
closest = closest_tri;
dist = dist_tri;
}
}

return closest;
}

template <typename IdxType, typename Vec3fType> Vec3fType
BVHTree<IdxType, Vec3fType>::closest_point(Vec3fType vertex, float max_dist) const {

float dist = max_dist * max_dist;
Vec3fType closest;

typename Node::ID node_id = 0;
std::stack<typename Node::ID> s;
while (true) {
Node const & node = nodes[node_id];
if (node.left != NAI && node.right != NAI) {
Vec3fType closest_left = acc::closest_point(vertex, nodes[node.left].aabb);
Vec3fType closest_right = acc::closest_point(vertex, nodes[node.right].aabb);
float dmin_left = (closest_left - vertex).square_norm();
float dmin_right = (closest_right - vertex).square_norm();
bool left = dmin_left < dist;
bool right = dmin_right < dist;
if (left && right) {
if (dmin_left < dmin_right) {
s.push(node.right);
node_id = node.left;
} else {
s.push(node.left);
node_id = node.right;
}
} else {
if (right) node_id = node.right;
if (left) node_id = node.left;
}

if (!left && !right) {
if (s.empty()) break;
node_id = s.top(); s.pop();
}
} else {
Vec3fType closest_leaf = closest_point(vertex, node_id);
float dist_leaf = (closest_leaf - vertex).square_norm();
if (dist_leaf < dist) {
dist = dist_leaf;
closest = closest_leaf;
}

if (s.empty()) break;
node_id = s.top(); s.pop();
}
}

return closest;
}

ACC_NAMESPACE_END

#endif /* ACC_BVHTREE_HEADER */
120 changes: 93 additions & 27 deletions primitives.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@ struct Ray {
};

template <typename Vec3fType> inline
AABB<Vec3fType> operator+(AABB<Vec3fType> const & a, AABB<Vec3fType> const & b) {
AABB<Vec3fType> operator+(AABB<Vec3fType> const & a, AABB<Vec3fType> const & b)
{
AABB<Vec3fType> aabb;
for (std::size_t i = 0; i < 3; ++i) {
aabb.min[i] = std::min(a.min[i], b.min[i]);
Expand All @@ -49,42 +50,47 @@ AABB<Vec3fType> operator+(AABB<Vec3fType> const & a, AABB<Vec3fType> const & b)
}

template <typename Vec3fType> inline
void operator+=(AABB<Vec3fType> & a, AABB<Vec3fType> const & b) {
void operator+=(AABB<Vec3fType> & a, AABB<Vec3fType> const & b)
{
for (int i = 0; i < 3; ++i) {
a.min[i] = std::min(a.min[i], b.min[i]);
a.max[i] = std::max(a.max[i], b.max[i]);
}
}

template <typename Vec3fType> inline
void calculate_aabb(Tri<Vec3fType> const & tri, AABB<Vec3fType> * aabb) {
void calculate_aabb(Tri<Vec3fType> const & tri, AABB<Vec3fType> * aabb)
{
for (int i = 0; i < 3; ++i) {
aabb->min[i] = std::min(tri.a[i], std::min(tri.b[i], tri.c[i]));
aabb->max[i] = std::max(tri.a[i], std::max(tri.b[i], tri.c[i]));
}
}

template <typename Vec3fType> inline
float surface_area(AABB<Vec3fType> const & aabb) {
float surface_area(AABB<Vec3fType> const & aabb)
{
float e0 = aabb.max[0] - aabb.min[0];
float e1 = aabb.max[1] - aabb.min[1];
float e2 = aabb.max[2] - aabb.min[2];
return 2.0f * (e0 * e1 + e1 * e2 + e2 * e0);
}

template <typename Vec3fType> inline
float mid(AABB<Vec3fType> const & aabb, std::size_t d) {
float mid(AABB<Vec3fType> const & aabb, std::size_t d)
{
return (aabb.min[d] + aabb.max[d]) / 2.0f;
}

constexpr float inf = std::numeric_limits<float>::infinity();
constexpr float eps = 1e-3;//std::numeric_limits<float>::epsilon();
constexpr float flt_eps = std::numeric_limits<float>::epsilon();

/* Derived form Tavian Barnes implementation posted in
* http://tavianator.com/fast-branchless-raybounding-box-intersections-part-2-nans/
* on 23rd March 2015 */
template <typename Vec3fType> inline
bool intersect(Ray<Vec3fType> const & ray, AABB<Vec3fType> const & aabb, float * tmin_ptr) {
bool intersect(Ray<Vec3fType> const & ray, AABB<Vec3fType> const & aabb, float * tmin_ptr)
{
float tmin = ray.tmin, tmax = ray.tmax;
for (int i = 0; i < 3; ++i) {
float t1 = (aabb.min[i] - ray.origin[i]) / ray.dir[i];
Expand All @@ -98,21 +104,9 @@ bool intersect(Ray<Vec3fType> const & ray, AABB<Vec3fType> const & aabb, float *
}

template <typename Vec3fType> inline
bool intersect(Ray<Vec3fType> const & ray, Tri<Vec3fType> const & tri, float * t_ptr, Vec3fType * bcoords_ptr) {
Vec3fType v0 = tri.b - tri.a;
Vec3fType v1 = tri.c - tri.a;
Vec3fType normal = v1.cross(v0);
if (normal.norm() < std::numeric_limits<float>::epsilon()) return false;
normal.normalize();

double cosine = normal.dot(ray.dir);
if (std::abs(cosine) < std::numeric_limits<double>::epsilon()) return false;

float t = -normal.dot(ray.origin - tri.a) / cosine;

if (t < ray.tmin || ray.tmax < t) return false;
Vec3fType v2 = (ray.origin - tri.a) + t * ray.dir;

Vec3fType barycentric_coordinates(Vec3fType const & v0, Vec3fType const & v1,
Vec3fType const & v2)
{
/* Derived from the book "Real-Time Collision Detection"
* by Christer Ericson published by Morgan Kaufmann in 2005 */
float d00 = v0.dot(v0);
Expand All @@ -127,12 +121,84 @@ bool intersect(Ray<Vec3fType> const & ray, Tri<Vec3fType> const & tri, float * t
bcoords[2] = (d00 * d21 - d01 * d20) / denom;
bcoords[0] = 1.0f - bcoords[1] - bcoords[2];

*t_ptr = t;
*bcoords_ptr = bcoords;
return bcoords;
}

template <typename Vec3fType> inline
Vec3fType closest_point(Vec3fType const & v, AABB<Vec3fType> const & aabb)
{
Vec3fType ret;
for (int i = 0; i < 3; ++i) {
ret[i] = std::max(aabb.min[i], std::min(v[i], aabb.max[i]));
}
return ret;
}

template <typename Vec3fType> inline
Vec3fType closest_point(Vec3fType const & vertex, Tri<Vec3fType> const & tri)
{
Vec3fType v0 = tri.b - tri.a;
Vec3fType v1 = tri.c - tri.a;
Vec3fType normal = v1.cross(v0).normalized();

Vec3fType p = vertex - normal.dot(vertex - tri.a) * normal;
Vec3fType v2 = p - tri.a;

Vec3fType bcoords = barycentric_coordinates(v0, v1, v2);

if (bcoords[0] < 0.0f) {
Vec3fType v12 = tri.c - tri.b;
float n = v12.norm();
float t = std::max(0.0f, std::min(v12.dot(p - tri.b) / n, n));
return tri.b + t / n * v12;
}

if (bcoords[1] < 0.0f) {
Vec3fType v20 = tri.a - tri.c;
float n = v20.norm();
float t = std::max(0.0f, std::min(v20.dot(p - tri.c) / n, n));
return tri.c + t / n * v20;
}

if (bcoords[2] < 0.0f) {
Vec3fType v01 = tri.b - tri.a;
float n = v01.norm();
float t = std::max(0.0f, std::min(v01.dot(p - tri.a) / n, n));
return tri.a + t / n * v01;
}

return tri.a * bcoords[0] + tri.b * bcoords[1] + tri.c * bcoords[2];
}

template <typename Vec3fType> inline
bool intersect(Ray<Vec3fType> const & ray, Tri<Vec3fType> const & tri,
float * t_ptr, Vec3fType * bcoords_ptr)
{
Vec3fType v0 = tri.b - tri.a;
Vec3fType v1 = tri.c - tri.a;
Vec3fType normal = v1.cross(v0);
if (normal.norm() < flt_eps) return false;
normal.normalize();

double cosine = normal.dot(ray.dir);
if (std::abs(cosine) < flt_eps) return false;

float t = -normal.dot(ray.origin - tri.a) / cosine;

if (t < ray.tmin || ray.tmax < t) return false;
Vec3fType v2 = (ray.origin - tri.a) + t * ray.dir;

Vec3fType bcoords = barycentric_coordinates(v0, v1, v2);

constexpr float eps = 1e-3f;
if (-eps > bcoords[0] || bcoords[0] > 1.0f + eps) return false;
if (-eps > bcoords[1] || bcoords[1] > 1.0f + eps) return false;
if (-eps > bcoords[2] || bcoords[2] > 1.0f + eps) return false;

if (t_ptr != nullptr) *t_ptr = t;
if (bcoords_ptr != nullptr) *bcoords_ptr = bcoords;

return -eps <= bcoords[0] && bcoords[0] <= 1.0f + eps
&& -eps <= bcoords[1] && bcoords[1] <= 1.0f + eps
&& -eps <= bcoords[2] && bcoords[2] <= 1.0f + eps;
return true;
}

ACC_NAMESPACE_END
Expand Down

0 comments on commit dcf3b3e

Please sign in to comment.