Improvement in bounding box intersection gives 15% improvement in run time + minor bug fix

[ndinsmore]

First, the bug fix is to make unitize return a vector with a magnitude of 1.

I thought this would just be a micro-optimization but it ends up being fairly significant.

There are two caveats to this improvement:

1. As you will see in the code both the affected methods calculated a term distsurf, which I could not really figure out what it consistently calculated so I was unable to implement it. I don't think that should change the results much.
2. I do get some minor differences in the heatmap from the x-ray, when I take a difference between the two methods. I honestly don't know what is right.

Running the benchmark:

using Revise, Geom4hep
using BenchmarkTools
includet("./XRay.jl")
const full2 = processGDML("examples/cms2018.gdml", Float64)
const volume2 = full2[1,7,1]
const world2 = getWorld(volume2)
const nav2 = BVHNavigator(world2)
@time img= generateXRay(nav2, world2, 1e4, 1)
@time img= generateXRay(nav2, world2, 1e4, 1)

Master:

julia> @time img= generateXRay(nav2, world2, 1e4, 1);
 2.374147 seconds (10 allocations: 80.766 KiB)

This PR:

julia> @time img= generateXRay(nav2, world2, 1e4, 1);
  0.773122 seconds (10 allocations: 80.766 KiB)

[peremato]

Thank-you very much. It looks really interesting.

• fortunately the buggy unitize was not used anywhere
• the unit tests are failing. test/runtests.jl. It is because the conditions with the return of intersectAABoxRay(...) needs to be adjusted.
• passing the reciprocal of direction as argument can save time in some loops. For example in Navigators.jl#L119 where the direction does not change.

[ndinsmore]

I fixed the issue wit the test. There are two left failing because the distout is 2.09e-9 and the failure criteria is < 1e-9/2

That said a large part of the improvement has disappeared. This is still some but it is only like 15%-20% now

[peremato]

The following implementation of distanceToIn(box::Box{T}...) fixes all the tests. I have not observed a nan.

function distanceToIn(box::Box{T}, point::Point3{T}, direction::Vector3{T},rcp_direction::Vector3{T}=inv.(direction))::T where T<:AbstractFloat
    max = Point3{T}(box.fDimensions)
    (distin, distout) = intersectAABoxRay(-max, max, point, direction, rcp_direction)
    (distin >= distout || distin <= -kTolerance(T)) ? Inf : distin
end

However I still have problem generating the X-ray in the z direction.

[ndinsmore]

What was the edge case that the surface checks were trying to fix? I noticed that there was not something similar VecGeom. I believe that the new method is effectively doing what the previous methods were, this is just a streamlined slabs, so i don't know how that could effect the z direction.

Also do you have any opinion on the naming of intersectAABoxRay it was originally intersectAABBRay but I thought that was to close the AABB type.

[ndinsmore]

The isnan is important because it will happen for 0*Inf then when you do a comparison to the nan it returns false. The fixed version I posted that has isnan was important for that reason.

[peremato]

Yes, you are right the handling of NaN is very important. But this way produces the wrong result. See for example a unit Box and a Point in one of the vertices.

Geom4hep.distanceToIn(Box{T}(1,1,1), Point3{T}(-1,-1,-1),Vector3{T}(0,1,0))

The result should be 0 since you are already in the surface of the Box while it returns Inf. The problem is that the functions min() and max() do not handle correctly the NaN case.

[ndinsmore]

Should work now. There is an interesting post about all the floating point logic at: https://tavianator.com/2022/ray_box_boundary.html

[peremato]

Yes, it does work now. Thanks for the pointer to the post. I can obtain the same X-ray images. However the speedup has vanished or almost vanished. Did you measure it again with this latest version?

[ndinsmore]

Yeah I was still getting about 10% 2.2s -> 2s

[ndinsmore]

I would not merge this anymore just merge #6 instead of this it includes it and actually results in a performance benefit

[peremato]

Supplemented by #6