madmann91 / bvh Goto Github PK

Bit field packing is implementation dependent. Unfortunately, you could argue MSVC is not as smart as GCC (or whatever you developed for) in terms of packing: https://docs.microsoft.com/en-us/cpp/cpp/cpp-bit-fields?redirectedfrom=MSDN&view=msvc-160

For instance, sizeof(Bvh::Node) ends up at 72 instead of 64 bytes because the packed in bool effectively becomes an unpacked IndexType. Switching to a uniform type across the bit fields will fix the issue:

IndexType is_leaf : 1;
IndexType primitive_count : sizeof(IndexType) * CHAR_BIT - 1;
IndexType first_child_or_primitive;

Obviously has the draw back of losing the explicit bool type. There may be another way to get the desired packing out of MSVC, but I am not aware of it.

Not sure if/how you want to handle this, so didn't do a PR.

How hard would it be to modify the traversal algorithm to include all hits?

Are they not supported? I use indexed meshes because it's faster applying transformations to them. The difficulty is that this library seems to require that each primitive is an instanced class.

Specifically, it's arranged like this.

template <typename Float, typename Index>
struct Mesh final {
  Eigen::Matrix<Float, 3, Eigen::Dynamic> vertices;
  Eigen::Matrix<Float, 3, Eigen::Dynamic> normals;
  Eigen::Matrix<Float, 2, Eigen::Dynamic> tex_coords;

  /* +3 vertices, +3 normals, +3 texture coordinates */
  Eigen::Matrix<Index, 9, Eigen::Dynamic> faces;
};

I don't think this is too uncommon, it's how libigl seems to handle meshes as well.

Embree has a progress monitor for build operations, that also allows you to cancel the build operation:
https://github.com/embree/embree/blob/master/README.md#description-description-21

It would be great to have something like this for this library as well. The ability to cancel in particular can be very handy for large scenes where building the bvh can take a very long time.

Here I'm using 32bit floats. The test code generates a depth image using SingleRayTraverser. Seems that a few pixels missed intersections at triangle edges. I'm not very sure the is issue is with bvh. Just trying to ask and see if there's any suggestion. Thanks.

1. Is there a way to deal with motion blur or deforming primitives?

2. For spatial splits builder, what is the general approach to splitting an arbitrary primitive? Just split bbox
   at split point?

Currently, if I'm reading the builders correctly, they always allocate a new bvh instance when building. It would be nice to have the option of re-using allocations.

Hi,
I'm a little further in my ray tracing application and now want to use the normals of the triangles to calculate the angle of incidence of a ray.

auto dirDotN = dot(rayDirection, triangle.n);
auto incidenceAngle = std::acos(dirDotN);

I discovered that the normal inside the bvh::Triangle is not normalized after it has been calculated.
I was under the impression that when a vector is used to express a direction (and a triangle's normal is such a vector), it should always be normalized...
I tried to add a normalize() call within the Triangle's constructor:

n = normalize(LeftHandedNormal ? cross(e1, e2) : cross(e2, e1));

but then the triangle's intersection routine doesn't work anymore.
(I used the same test program as I did in #18, and it doesn't report any hits).

Of course, I can add another member that contains the normalized normal, but I would like to understand why the intersection routine apparently cannot cope with a normalized normal...

Could you please enlighten me on this subject?

simple_example and benchmark need pthread to compile on linux with clang

Hi,

I'm still exploring the library so apologies if I'm doing something stupid. I'm primarily interested in fast "picking" and frustum intersection for box selection. So I was looking into using less memory (at the cost of more cpu) by setting config.min_leaf_size = 8, but it lead to an assertion error.

Thanks

Can you clarify how the flags array is initialized in the bottom up algorithm?
It seems to me as if its uninitialized in the constructor and first accessed here. The leaf collapser also does not seem to touch the flags array.

Right now I see that if I cast rays where ray.x = 0 or ray.z = 0, the BVH tree does not return a triangle. This could affect ray.y but I have not really tried that yet. Adding epsilon in those corner cases fixes the problem but there's probably a better solution.

BVH tree
no BVH tree, just raycast against all triangles
intersection test always returns true, closest intersector
intersection test always returns true, any intersector

In the benchmark, if I force dir.x = 0 (then normalize afterwards), nothing at all renders to the screen:

benchmark.cpp:121

            Ray ray(camera.eye, bvh::normalize(image_u * u + image_v * v + dir));
            ray.direction.values[0] = 0.0f;
            ray.direction = bvh::normalize(ray.direction);

without modification
setting value to 0.000001f
setting value to 0.0f

This seems to vary by each model, models with high triangle density don't seem to suffer from this issue.

The below function processes from node of leaf in order.
However, it looks like it is continuing when is_leaf is true.
Is this correct?

bottom_up_algorithm.hpp
void traverse_in_parallel
...
// Only process leaves
if (MaintainChildIndices ? children[i] != 0 : bvh.nodes[i].is_leaf)
continue;
...

Each node index can be packed into a single Node::Index::Type like so:

using Node_t = decltype(node);

//Serialize:
Node_t::Index::Type serialized = (node.index.first_id << Node_t::prim_count_bits) |
                                 node.index.prim_count;

//Deserialize:
node.index.first_id = serialized >> Node_t::prim_count_bits;
node.index.prim_count = serialized & ((1 << Node_t::prim_count_bits) - 1);

If you're worried about backwards compatibility or users who prefer speed over size, you can hide the optimization behind a bool template parameter.

I'm writing a JavaScript BVH library and using your implementation as a reference. Not for ray tracing, but for picking and collisions. The library would be intended to run in a Worker (a separate long-lived process, much heavier than a thread and only one "CPU" for it)

I'd rather just use WASM and your library directly, but lack of easy shared memory between WASM and JS means copying all the vertex data over from JS to WASM (since the WebGL stuff is in JS) which is tricky because browsers are more memory-constrained than desktop software.

Do you have a suggestion for which construction method to use as a reference? Construction speed is a priority. I'm thinking bvh::LocallyOrderedClusteringBuilder, but what do you think?

Using the Coderrect Scanner coderrect.com I discovered a data race in bvh.hpp between line 42 and itself.

 40|
 41| BoundingBoxProxy& operator = (const BoundingBox<Scalar>& bbox) {
>42|     node.bounds[0] = bbox.min[0];
 43|     node.bounds[1] = bbox.max[0];
 44|     node.bounds[2] = bbox.min[1];

Hi,
there seems to be a problem when there is a -0.0 value in the direction vector. I slightly modified your example from simple_example.cpp to reproduce this error:

I just copied the first triangle and moved it a bit further away (in the z-direction) and also changed the origin of the ray so that there should be a clear intersection:

    std::vector<Triangle> triangles;
    triangles.emplace_back(
        Vector3( 1.0, -1.0, 1.0),
        Vector3( 1.0,  1.0, 1.0),
        Vector3(-1.0,  1.0, 1.0)
    );
    triangles.emplace_back(
        Vector3( 1.0, -1.0, 2.1),
        Vector3( 1.0,  1.0, 2.1),
        Vector3(-1.0,  1.0, 2.1)
    );

    // Version 1: hit: triangle:0, distance:1, u:0.25, v:0.375
    Ray ray(
        Vector3(0.25, 0.25, 0.0),
        Vector3(0.0, 0.0, 1.0), 
        0.0,
        100.0
    );

    // Version 2: no hit
    Ray ray2(
        Vector3(0.25, 0.25, 0.0),
        Vector3(0.0, -0.0, 1.0), 
        0.0,
        100.0
    );

But there is still a hit in the second version when the second triangle is at z=2.0. I also tried RobustNodeIntersector, but with no luck :/

Only want to let you know.

I'd like to use this library for something more akin to collision/nearest neighbors detection rather than ray-tracing. But this would require a bit more dynamism. I propose adding at least a DefaultAdder, DefaultRemover, and maybe even DefaultUpdater, possibly with alternative versions; in the style of DefaultBuilder and its alternatives. Adding and removing elements dynamically while maintaining a good tree is a hard problem I'm sure, but even a simplistic implementation would be nicer than having to roll a bad implementation myself.

Are there face culling options in bvh? It seems to me that a ray can intersect a triangle no matter which direction it comes from.
I hope there can be options to skip triangles based on their winding directions. Is that easy to do?
Thanks.

I ran your benchmark using a sponza model and compared it to the latest Fast-BVH after I PR I made there. I think results are a bit different at this point. Fast-BVH is now a template library and built the sponza model I had in 70 ms. In your library, it built the same model in 133 ms. Maybe you could verify the results. In my patch to Fast-BVH, I also included an .obj file renderer.

Hello!
How to use library with instances?
I need BLAS and TLAS.
Can you create example for it?
Thank you

This was on VS 16.9.5 not sure if it's harmless or not, I silenced it with a #pragma.

When I updated my code to use v2, I feel that build() in debug build got slower.
Is there anything I could do to make build() faster in debug builds? I am already using low quality setting for the default builder. Thanks!

This library can easily be used for 2D contexts as well. The only additional feature it would need is a traverser that takes advantage of the simpler AABB test. Here's the AABB test, where x and y are the pixel coordinates.

float intersect(float x, float y, const Bvh<float>::Node& node)
{
  if ((x < node.bounds[0])
   || (x > node.bounds[1])
   || (y < node.bounds[2])
   || (y > node.bounds[3])) {
    return std::numeric_limits<float>::infinity();
  }
  return std::min(node.bounds[4], node.bounds[5]);
}

This is how I currently use this library for a 2D renderer that I've been working on. Would you be interested in a PR for this? The only thing I'd need to do is add support for the Intersector template parameter and Statistics parameter, the rest is already written and quite simple.

As it stands now, there's only one traversal algorithm. I'd like to implement the following algorithms:

• A frustum-based ray traversal algorithm,
• A aligned/oriented bounding-box intersection algorithm,
• A BVH/BVH intersection algorithm.
• Closest point detection.

It would be great to have this project available via vcpkg. Given that it is header-only, this should be relatively straight forward.

Relevant tutorials/docs:

1. Create and publish packages
2. Add a port to the vcpkg open-source registry

The default usage of #pragma STDC FP_CONTRACT ON in utilities.hpp causes warnings in MSVC 2022:

bvh\include\bvh/utilities.hpp(41,13): error C2220: the following warning is treated as an error
bvh\include\bvh/utilities.hpp(41,13): warning C4068: unknown pragma 'STDC'
bvh\include\bvh/utilities.hpp(50,13): warning C4068: unknown pragma 'STDC'

Hello!
There is an undefined behavior in file bvh/include/bvh/morton.hpp, at this function

template <typename Morton>
Morton morton_split(Morton x) {
    constexpr size_t log_bits = round_up_log2(sizeof(Morton) * CHAR_BIT);
    auto mask = Morton(-1);
    for (size_t i = log_bits, n = 1 << log_bits; i > 0; --i, n >>= 1) {
        mask = (mask | (mask << n)) & ~(mask << (n / 2));
        x = (x | (x << n)) & mask;
    }
    return x;
}

mask << n and x << n is a shift of n-bit operand on n bits. This is an undefined behavior, because of

if the value of the right operand is negative or is greater or equal to the number of bits in the promoted left operand, the behavior is undefined.
(https://en.cppreference.com/w/cpp/language/operator_arithmetic)

I found it when i tried to run this code with a clang. All my morton codes were 0.

Unfortunately the LocallyOrderedClustering does not terminate with the lucy 3D model, maybe some race condition?! I am on Linux with gcc 10.2. I'll will debug this a bit more, but in case you want to try to reproduce, I have uploaded a minimal reproducer here.

Just out of curiosity, is there any metrics that this library can calculate that indicate the quality of the generated BVHs aside from the traversal performance? Perhaps like the accumulative amount of empty space between parent and sub nodes? I don't think I've seen this measured before in other projects, but it seems like it would be a useful metric for comparisons between algorithms.

I am trying to enable use of precompiled headers in my project which uses bvh. Using PCH seems to cause MSVC compiler have issues with identifiers near and far in bvh.h. I wonder if these could be modified in order to avoid the issue? Thanks.

• on vs2017
• should include header file

I get very good performance out of your library for general ray tracing usage. But when ray tracing with vertical rays to get height of a mesh at given X,Y position, performance is terrible. Instead of visiting just a few leafs algorithm often visits almost all leafs.

Do you have any suggestions how to fix or mitigate this?

MSVC gives this warning from bvh/vector.hpp:

warning C4141: 'inline': used more than once when it encounters bvh_always_inline inline.

To me bvh_always_inline inline seems redundant, would it be okay to make it just bvh_always_inline?

Great lib!

However I get very low performance when casting axis-aligned rays.
It takes 7 seconds or so, instead of ~10ms (~30k rays on a ~3M meshes, the mesh has axis aligned triangles as well).

When I add a small epsilon to the ray to make sure all components are non null, things become fast again.

I didn't dig the issue further.
I can just say there is no division by 0 here

I believe there's a potential optimization here. If you assume the ray direction is normalized, the length of the vector will always equal 1. Since sqrt(1) is always 1, you can probably just leave it as:

const Scalar a = 1;

Embree still takes the dot product, perhaps because they don't assume the ray direction is normalized. Since there's some template magic at work here, you could probably do something like:

template <typename Scalar, bool NormalizedRayDir>
struct Sphere final {
  /* ... */
  static Scalar get_ray_direction_squared_length(const Ray<Scalar> &ray);
};

template <typename Scalar>
Scalar Sphere<Scalar, true>::get_ray_direction_squared_length(const Ray<Scalar>&)
{
  return 1;
}

template <typename Scalar>
Scalar Sphere<Scalar, false>::get_ray_direction_squared_length(const Ray<Scalar>& ray)
{
  return dot(ray.direction, ray.direction);
}

Maybe it would also do well as a template parameter for Ray instead, hard to say.

Hi, is there any way to adapt this library to do closest point / distance queries in addition to pure ray casting?

OpenMP requires separate redistributable for MSVC, an option to remove the dependency or use C++11 threads would be appreciated (even at the cost of slower BVH building).

Hey, so ever since commit d2f5e2f , which added the global bounding box, my renderer does not seem to render anything under the same scene/camera settings. I'm using bvh::compute_bounding_boxes_union, yet despite this, I can't get it to properly intersect with any triangle.

  auto globalBbox =
    bvh::compute_bounding_boxes_union(
      boundingBoxes.data()
    , self->triangles.size()
    );

  { // -- build tree
    bvh::BinnedSahBuilder<decltype(self->boundingVolume), 32>
      boundingVolumeBuilder { self->boundingVolume };

    boundingVolumeBuilder.build(
      globalBbox, boundingBoxes.data(), centers.data(), self->triangles.size()
    );
  }

This code works (without globalBbox references) before the aforementioned commit. Let me know if I'm doing something wrong here.

Hello, sorry for polluting the other thread: I just subscribed to github or this very issue and didn't realize I was replying to another issue :D

Here my original message.

Hello, great library!
I just got back to my old raytracing project (abandoned 10 years ago) and decided that it would have been nice to work on the performances, horrible at best. I've started replacing my quick BIH implementation with a BVH but decided to switch to an external lib and found this one. I could make it work, but performances were exactly the same as with my BIH:

bvh::LocallyOrderedClusteringBuilder<bvh::Bvh, int> builder(bvh);
builder.build(bvh::BoundingBox(MasterBVHPrimitive::toVector3(worldbb.min), MasterBVHPrimitive::toVector3(worldbb.max)), bboxes, centers, shapes.size());
bvh::LeafCollapserbvh::Bvh collapser(bvh);
collapser.collapse();

I've started profiling the execution and found that 35% of the time was being spent inside std::fill, and I could track it back to
struct Vector {
Scalar values[N];

Vector() = default;
bvh__always_inline__ Vector(Scalar s) { std::fill(values, values + N, s); }
As I'm working with 3-dimensional vectors, I then replaced the call to fill with:
values[0] = values[1] = values[2] = s;

My test scene rendering time went from 8 seconds to 3!

Am I using the lib the wrong way? Is there an option that would disable/improve that initialization? Would my change break something somewhere?

Thank you!

Hey,

first, great work with this lib, really like using it!

I think it would be great to have a benchmark comparing with fcpw,
which claims to be a very fast BVH implementation. It claims to achieve very fast query times,
by (presumably) vectorizing the slab test with enoki. My feeling is that the tree quality is much more
important, since queries should be latency bound.

Hello, me again!
In the test scene I'm currently using there are a few triangles apparently missing. I was sure it was a problem with the scene file (exported from blender), with Assimp (the lib I use to load the scene) or my code. I've started to play a bit with the mesh in blender (simplify, remove vertices, triangulate etc), then by exporting to different file formats. Nothing changed. Then I tried to render the scene with my BIH, and, to my surprise, the scene was fine.
That's how it looks with my BIH

and with the BVH:

(sorry for the bad quality, hope will do anyway). Please note the fireplace, or the corner opposite to the point of view. There are a couple more, but you can't see it in the that image (I can of course provide you with better ones).
I understand that this is not much to work with, but thought to report it should you be interested in analyzing the issue. Let mw know if you want the collada/gltf/blend file of the scene.

That's how I create the bvh:

	void Preprocess()
	{
		//bvh::SweepSahBuilder<bvh::Bvh<real>> builder(bvh);
		bvh::LocallyOrderedClusteringBuilder<bvh::Bvh<real>, int>  builder(bvh);

		//Adapt data format
		bvh::BoundingBox<real>* bboxes = new bvh::BoundingBox<real>[shapes.size()];
		bvh::Vector3<real>* centers = new bvh::Vector3<real>[shapes.size()];

		qlaabox<real> worldbb;
		for (int i = 0; i < shapes.size(); ++i) {
			qlaabox<real> bb = shapes[i].origShape->GetBoundingBox();
			bboxes[i] = bvh::BoundingBox<real>(MasterBVHPrimitive::toVector3(bb.min), MasterBVHPrimitive::toVector3(bb.max));
			worldbb.Union(bb);
			centers[i] = MasterBVHPrimitive::toVector3(bb.Center());
		}
		
		builder.build(bvh::BoundingBox<real>(MasterBVHPrimitive::toVector3(worldbb.min), MasterBVHPrimitive::toVector3(worldbb.max)), bboxes, centers, shapes.size());
		bvh::LeafCollapser<bvh::Bvh<real>> collapser(bvh);
		collapser.collapse();
		std::cout << "Radius: " << builder.search_radius << std::endl;
		std::cout << "Shapes count : " << shapes.size() << std::endl;
		//std::cout << "Tree depth : " << TreeDepth(root) << std::endl;
		std::cout << "Node count : " << bvh.node_count << std::endl;
		//std::cout << "Leaf count : " << LeafCount(root) << std::endl;
		/*std::cout << "Average primitives/leaf : " << ((real)ShapesCount(root)) / LeafCount(root) << std::endl;
		std::cout << "Empty leaves : " << EmptyLeaves(root) << std::endl;
		std::cout << "Max shapes : " << MaxShapes(root) << std::endl;*/
	}

And the traversal code:

	bool Intersect(const ray<vector4<real>>& r, SurfacePoint& hitdesc) const
	{
		// Intersect a ray with the data structure
		Ray ray = MasterBVHPrimitive::toRay(r);
		
		bvh::ClosestPrimitiveIntersector<bvh::Bvh<real>, MasterBVHPrimitive> primitive_intersector(bvh, shapes.data());
		bvh::SingleRayTraverser<bvh::Bvh<real>> traverser(bvh);

		auto hit = traverser.traverse(ray, primitive_intersector);
		if (hit)
		{
			hitdesc = hit->intersection.hitdesc;
			return true;
		}
		else
			return false;
	}

Thank you and regards,
Alessandro

Hi, I'm wondering if the BVH created using this library can insert primitives after construction, or whether the BVH would need to be reconstructed from scratch if new primitives are added to the scene?

Additionally, I'm curious whether it's possible to query the BVH with an AABB to find which nodes overlap with the AABB (I think this may be in the plans in #30)

I am currently workin on a raytrcer in opengl, do ou have a api to get the raw bvh node data (and any other essential data) to be uploaded to the gpu as a ssbo ?

Any interest in supporting acceleration structures for height fields?

There's a paper that goes into describing the process. I'm currently using this library for calculating ambient occlusion in terrains, and this approach may lead to a performance boost.

Even so much as a dedicated BVH builder for height fields would be really cool to see.

I noticed that the CMake version got bumped up again to version 3.13.5
This breaks the build for anyone who has the standard distribution of CMake on the latest Ubuntu (19.10), which is CMake version 3.13.4.

Was this bumped up because of target_link_options? That's supported by 3.13.4, even though the documentation doesn't specify it.

It might also make the Travis build simpler, since you won't have to install a new CMake version.

Hello! I was wondering whether it would be possible to either include as part of the data, or to recover based on the tree, which access a node was split on?

Hi,
I'm trying to understand why the normal of bvh::Triangle is calculated differently than the 'standard' way...
By 'standard', I am referring to the way of calculating the normal of a trial that I found described on many internet sources...
If a triangle is given by vertices p0, p1, p2, the standard way states that the normal is calculated as follows:
normal = cross(p1 - p0, p2 - p0)
If I understand correctly, in bvh::Triangle the edges are stored e1 = (p0 - p1) and e2 = (p2 - p0).
In the constructor n is calculated as cross(e1,e2), which makes the normal point in the opposite direction than the normal of the 'standard' case.

Could you please enlighten me on this subject?