peremato · Moelf · Sep 22, 2022 · Sep 22, 2022 · Sep 22, 2022 · Sep 22, 2022
diff --git a/examples/XRay.jl b/examples/XRay.jl
@@ -7,7 +7,7 @@ const idx = [2 3 1; 1 3 2; 1 2 3]
 const rdx = [3 1 2; 1 3 2; 1 2 3]
 
 # Generate a X-Ray of a geometry----------------------------------------------
-function generateXRay(nav::AbstractNavigator, vol::Volume{T}, npoints::Number, view::Int=1) where T<:AbstractFloat
+function generateXRay(nav::AbstractNavigator, vol::VolumeP{T, S}, npoints::Number, view::Int=1) where {T<:AbstractFloat, S}
     world = getWorld(vol)
     # Setup plane of points and results
     lower, upper = extent(world.shape)
@@ -52,36 +52,17 @@ function generateXRay(nav::AbstractNavigator, vol::Volume{T}, npoints::Number, v
     return xrange, yrange, result
 end
 
-using GLMakie
+const full = processGDML("examples/trackML.gdml", Float64)
+const volume = full[2,1]
+const world = getWorld(volume)
+const nav = BVHNavigator(world)
 
-#if abspath(PROGRAM_FILE) == @__FILE__
-    #-----build detector---------------------------------------
-    #full = processGDML("examples/trackML.gdml", Float64)
-    #volume = full[2,1]
-
-    full = processGDML("examples/cms2018.gdml", Float64)
-    volume = full[1,7,1]
-
-    world = getWorld(volume)
-    nav = BVHNavigator(world)
-    #nav = TrivialNavigator(world)
+const full2 = processGDML("examples/cms2018.gdml", Float64)
+const volume2 = full2[1,7,1]
+const world2 = getWorld(volume2)
+const nav2 = BVHNavigator(world2)
 
-    #-----Generate and Plot results----------------------------
-    @printf "generating x-projection\n"
-    rx = generateXRay(nav, world, 1e6, 1)
-    @printf "generating y-projection\n"
-    ry = generateXRay(nav, world, 1e6, 2)
-    @printf "generating z-projection\n"
-    rz = generateXRay(nav, world, 1e6, 3)
-    limits = (min(minimum(rx[3]), minimum(ry[3]), minimum(rz[3])), max(maximum(rx[3]), maximum(ry[3]), maximum(rz[3])))
-
-    #----Plot the results--------------------------------------
-    fig = Figure(resolution = (2000, 2000))
-    @printf "ploting the results\n"
-    heatmap!(Axis(fig[1, 1], title = "X direction"), rx..., colormap=:grayC, colorrange=limits)
-    heatmap!(Axis(fig[2, 1], title = "Y direction"), ry..., colormap=:grayC, colorrange=limits)
-    heatmap!(Axis(fig[1, 2], title = "Z direction"), rz..., colormap=:grayC, colorrange=limits)
-    draw!(LScene(fig[2, 2]), volume; wireframe=true, maxlevel=2)
-    #display(fig)
-    save("cms2018.png", fig)
-#end
+function mainbench()
+    @time generateXRay(nav, world, 1e4, 1);
+    @time generateXRay(nav2, world2, 1e4, 1);
+end
diff --git a/src/Boolean.jl b/src/Boolean.jl
@@ -80,69 +80,6 @@ function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T
 end
 
 
-function inside(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
-    (; left, right, transformation) = shape
-    lpoint = transformation * point
-
-    positionA = inside(left, point)
-    positionA == kInside && return kInside
-
-    positionB = inside(right, lpoint)
-    positionB == kInside && return kInside
-
-    if positionA == kSurface && positionB == kSurface
-        normalA = normal(left, point)
-        normalB = normal(right, lpoint) * transformation
-        if dot(normalA, normalB) < 0
-            return kInside    # touching solids -)(-
-        else 
-            return kSurface   # overlapping solids =))
-        end
-    elseif positionA == kSurface || positionB == kSurface
-        return kSurface
-    else
-        return kOutside
-    end
-end
-
-function inside(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::Int64  where {T,SL,SR}
-    (; left, right, transformation) = shape
-    lpoint = transformation * point
-
-    positionA = inside(left, point)
-    positionA == kOutside && return kOutside
-
-    positionB = inside(right, lpoint)
-    if positionA == kInside && positionB == kInside
-        return kInside
-    elseif  positionA == kInside && positionB == kSurface ||
-            positionA == kSurface && positionB == kInside ||
-            positionA == kSurface && positionB == kSurface
-        return kSurface
-    else
-        return kOutside
-    end 
-end
-
-function inside(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::Int64  where {T,SL,SR}
-    (; left, right, transformation) = shape
-    lpoint = transformation * point
-
-    positionA = inside(left, point)
-    positionA == kOutside && return kOutside
-
-    positionB = inside(right, lpoint)
-
-    if positionA == kInside && positionB == kOutside
-        return kInside;
-    elseif positionA == kInside && positionB == kSurface ||
-           positionB == kOutside && positionA == kSurface
-        return kSurface
-    else
-        return kOutside
-    end
-end
-
 function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR}
     return 0
 end
@@ -162,173 +99,3 @@ end
 function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR}
     return 0
 end
-
-function distanceToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR}
-    (; left, right, transformation) = shape
-    dist = T(0)
-    positionA = inside(left, point)
-    if positionA != kOutside  # point inside A
-        while(true)
-            distA = distanceToOut(left, point, dir)
-            dist += (distA > 0 && distA < Inf) ? distA : 0
-            dist += kPushTolerance(T)
-            npoint = point + dist * dir
-            lpoint = transformation * npoint
-            if inside(right, lpoint) != kOutside # B could be overlapping with A -- and/or connecting A to another part of A
-                ldir = transformation * dir
-                distB = distanceToOut(right, lpoint, ldir)
-                dist += (distB > 0 && distB < Inf) ? distB : 0
-                dist += kPushTolerance(T)
-                npoint = point + dist * dir
-                if inside(left, npoint) == kOutside
-                    break
-                end
-            else
-                break
-            end
-        end
-        return dist - kPushTolerance(T)
-    else
-        lpoint = transformation * point
-        positionB = inside(right, lpoint)
-        if positionB != kOutside  # point inside B
-            ldir = transformation * dir
-            while(true)
-                distB = distanceToOut(right, lpoint, ldir)
-                dist += (distB > 0 && distB < Inf) ? distB : 0
-                dist += kPushTolerance(T) # Give a push
-                npoint = point + dist * dir
-                if inside(left, npoint) != kOutside # A could be overlapping with B -- and/or connecting B to another part of B
-                    ldir = transformation * dir
-                    distA = distanceToOut(left, npoint, dir)
-                    dist += (distA > 0 && distA < Inf) ? distA : 0
-                    dist += kPushTolerance(T)
-                    npoint = point + dist * dir
-                    lpoint = transformation * npoint
-                    if inside(right, lpoint) == kOutside
-                        break
-                    end
-                else
-                    break
-                end
-            end
-            return dist - kPushTolerance(T)
-        else
-            return T(-1)
-        end
-    end 
-end
-
-function distanceToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR}
-    (; left, right, transformation) = shape
-    distA = distanceToOut(left, point, dir)
-    distB = distanceToOut(right, transformation * point, transformation * dir)
-    return min(distA, distB)
-end
-
-function distanceToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR}
-    (; left, right, transformation) = shape
-    distA = distanceToOut(left, point, dir)
-    distB = distanceToIn(right, transformation * point, transformation * dir)
-    return min(distA, distB)
-end
-
-function distanceToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR}
-    (; left, right, transformation) = shape
-    distA = distanceToIn(left, point, dir)
-    distB = distanceToIn(right, transformation * point, transformation * dir)
-    return min(distA, distB)
-end
-
-function distanceToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR}
-    (; left, right, transformation) = shape
-
-    positionA = inside(left, point)
-    lpoint = transformation * point
-    positionB = inside(right, lpoint)
-
-    inLeft = positionA == kInside
-    inRight = positionB == kInside
-
-    inLeft && inRight && return T(-1)
-
-    dist = T(0)
-    npoint = point
-    ldir = transformation * dir
-    #  main loop
-    while true
-        d1 = d2 = 0
-        if !inLeft
-            d1 = distanceToIn(left, npoint, dir)
-            d1 = max(d1, kTolerance(T))
-            d1 == T(Inf) && return T(Inf)
-        end
-        if !inRight
-            d2 = distanceToIn(right, lpoint, ldir)
-            d2 = max(d2, kTolerance(T))
-            d2 == T(Inf) && return T(Inf)
-        end
-        if d1 > d2
-            # propagate to left shape
-            dist += d1
-            inleft = true
-            npoint += d1 * dir
-            lpoint = transformation * npoint
-            # check if propagated point is inside right shape, check is done with a little push
-            inRight = inside(right, lpoint + kTolerance(T) * ldir) == kInside
-            inRight && return dist
-            # here inleft=true, inright=false
-        else
-            # propagate to right shape
-            dist += d2
-            inright = true
-            # check if propagated point is inside left shape, check is done with a little push
-            npoint += d2 * dir
-            lpoint = transformation * npoint
-            inLeft = inside(left, npoint + kTolerance(T) * dir) == kInside
-            inLeft && return dist
-        end
-    end
-    return dist
-end
-
-function distanceToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR}
-    (; left, right, transformation) = shape
-
-    lpoint = transformation * point
-    ldir = transformation * dir
-    positionB = inside(left, lpoint)
-    inRight = positionB == kInside
-
-    npoint = point
-    dist = T(0)
-
-    while true
-        if inRight
-            # propagate to outside of '- / RightShape'
-            d1 = distanceToOut(right, lpoint, ldir)
-            dist += (d1 >= 0 && d1 < Inf) ? d1 + kPushTolerance(T) : 0
-            npoint = point + (dist + kPushTolerance(T)) * dir
-            lpoint = transformation * npoint
-            # now master outside 'B'; check if inside 'A'
-            inside(left, npoint) == kInside && distanceToOut(left, npoint, dir) > kPushTolerance(T) && return dist
-        end
-
-        # if outside of both we do a max operation master outside '-' and outside '+' ;  find distances to both
-        d2 = distanceToIn(left, npoint, dir)
-        d2 = max(d2, 0)
-        d2 == T(Inf) && return T(Inf)
-
-        d1 = distanceToIn(right, lpoint, ldir)
-        if d2 < d1 - kTolerance(T)
-          dist += d2 + kPushTolerance(T)
-          return dist
-        end
-
-        #   propagate to '-'
-        dist += (d1 >= 0 && d1 < Inf) ? d1 : 0
-        npoint = point + (dist + kPushTolerance(T)) * dir
-        lpoint = transformation * npoint
-        inRight = true
-    end
-end
diff --git a/src/Box.jl b/src/Box.jl
@@ -28,72 +28,6 @@ function normal(box::Box{T}, point::Point3{T}) where T
     end
 end
 
-function inside(box::Box{T}, point::Point3{T})::Int64  where T<:AbstractFloat
-    dist = -Inf
-    for i in 1:3
-        d = abs(point[i]) - box.fDimensions[i]
-        if d  > dist 
-            dist = d 
-        end
-    end
-    abs(dist) <= kTolerance(T)/2 ? kSurface : dist < 0.0 ? kInside : kOutside
-    #dist = maximum(abs.(point) - box.fDimensions)
-    #isapprox(dist, 0.0, atol = kTolerance(T)/2) ? kSurface : dist < 0.0 ? kInside : kOutside
-end
-
-function distanceToOut(box::Box{T}, point::Point3{T}, direction::Vector3{T})::T where T<:AbstractFloat
-    safety = -Inf
-    for i in 1:3
-        d = abs(point[i]) - box.fDimensions[i]
-        if d  > safety 
-            safety = d 
-        end
-    end
-    if safety > kTolerance(T)/2
-        return -1.0
-    end 
-    dist = Inf
-    for i in 1:3
-        d = (copysign(box.fDimensions[i], direction[i]) - point[i])/direction[i]
-        if d < dist
-            dist = d
-        end
-    end
-    return dist
-    #safety = maximum(abs.(point) - box.fDimensions)
-    #distance = minimum((copysign.(box.fDimensions, direction) - point) * (1.0 ./ direction))
-    #safety > kTolerance(T)/2 ? -1.0 : distance
-end
-
-function distanceToIn(box::Box{T}, point::Point3{T}, direction::Vector3{T})::T where T<:AbstractFloat
-    distsurf = Inf
-    distance = -Inf
-    distout = Inf
-    for i in 1:3
-        din  = (-copysign(box.fDimensions[i],direction[i]) - point[i])/direction[i]
-        tout =   copysign(box.fDimensions[i],direction[i]) - point[i]
-        dout = tout/direction[i]
-        dsur = copysign(tout, direction[i])
-        if din > distance 
-            distance = din 
-        end
-        if dout < distout
-            distout = dout
-        end
-        if dsur < distsurf
-            distsurf = dsur
-        end 
-    end
-    (distance >= distout || distout <= kTolerance(T)/2 || abs(distsurf) <= kTolerance(T)/2) ? Inf : distance
-    #invdir = 1.0 ./ direction
-    #tempIn  = -copysign.(box.fDimensions, direction) - point
-    #tempOut =  copysign.(box.fDimensions, direction) - point
-    #distance = maximum(tempIn * invdir)
-    #distout  = minimum(tempOut * invdir)
-    #distsurf = abs(minimum(copysign.(tempOut, direction)))
-    #(distance >= distout || distout <= kTolerance(T)/2 || distsurf <= kTolerance(T)/2) ? Inf : distance
-end
-
 function safetyToOut(box::Box{T}, point::Point3{T}) where T<:AbstractFloat
     minimum(box.fDimensions - abs.(point))
 end