From 13d2fabf77e620977d114469bda571df2897c7f5 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 16:45:04 -0400 Subject: [PATCH 01/16] good state --- examples/XRay.jl | 43 ++-- src/Boolean.jl | 77 -------- src/Box.jl | 53 ----- src/Distance.jl | 46 +++++ src/DistanceIn.jl | 478 +++++++++++++++++++++++++++++++++++++++++++++ src/DistanceOut.jl | 443 +++++++++++++++++++++++++++++++++++++++++ src/Geom4hep.jl | 3 + src/Polycone.jl | 2 +- src/Trd.jl | 2 +- src/Tube.jl | 167 ---------------- 10 files changed, 984 insertions(+), 330 deletions(-) create mode 100644 src/Distance.jl create mode 100644 src/DistanceIn.jl create mode 100644 src/DistanceOut.jl diff --git a/examples/XRay.jl b/examples/XRay.jl index ff332cd..7fd574f 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -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 index 976452d..d853501 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -163,83 +163,6 @@ function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T 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 diff --git a/src/Box.jl b/src/Box.jl index d1f0805..b4a9552 100644 --- a/src/Box.jl +++ b/src/Box.jl @@ -41,59 +41,6 @@ function inside(box::Box{T}, point::Point3{T})::Int64 where T<:AbstractFloat #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 diff --git a/src/Distance.jl b/src/Distance.jl new file mode 100644 index 0000000..4c5842c --- /dev/null +++ b/src/Distance.jl @@ -0,0 +1,46 @@ +function distanceToIn(shape, point, dir) + if shape isa Trap + distanceToIn_trap(shape::Trap, point, dir) + elseif shape isa Trd + distanceToIn_trd(shape::Trd, point, dir) + elseif shape isa Cone + distanceToIn_cone(shape::Cone, point, dir) + elseif shape isa Box + distanceToIn_box(shape::Box, point, dir) + elseif shape isa Tube + distanceToIn_tube(shape::Tube, point, dir) + elseif shape isa Polycone + distanceToIn_polycone(shape::Polycone, point, dir) + elseif shape isa CutTube + distanceToIn_cuttube(shape::CutTube, point, dir) + elseif shape isa BooleanUnion + distanceToIn_booleanunion(shape::BooleanUnion, point, dir) + elseif shape isa BooleanSubtraction + distanceToIn_booleansubtraction(shape::BooleanSubtraction, point, dir) + elseif shape isa BooleanIntersection + distanceToIn_booleanintersection(shape::BooleanIntersection, point, dir) + end +end +function distanceToOut(shape, point, dir) + if shape isa Trap + distanceToOut_trap(shape::Trap, point, dir) + elseif shape isa Trd + distanceToOut_trd(shape::Trd, point, dir) + elseif shape isa Cone + distanceToOut_cone(shape::Cone, point, dir) + elseif shape isa Box + distanceToOut_box(shape::Box, point, dir) + elseif shape isa Tube + distanceToOut_tube(shape::Tube, point, dir) + elseif shape isa Polycone + distanceToOut_polycone(shape::Polycone, point, dir) + elseif shape isa CutTube + distanceToOut_cuttube(shape::CutTube, point, dir) + elseif shape isa BooleanUnion + distanceToOut_booleanunion(shape::BooleanUnion, point, dir) + elseif shape isa BooleanSubtraction + distanceToOut_booleansubtraction(shape::BooleanSubtraction, point, dir) + elseif shape isa BooleanIntersection + distanceToOut_booleanintersection(shape::BooleanIntersection, point, dir) + end +end diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl new file mode 100644 index 0000000..b59ddb9 --- /dev/null +++ b/src/DistanceIn.jl @@ -0,0 +1,478 @@ +## Boolean +function distanceToIn_booleanunion(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_booleanintersection(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_booleansubtraction(shape, point::Point3{T}, dir) where T + (; 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 + +# Plane +function distanceToIn_plane(plane::Plane{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + # The function returns a negative distance for points already inside or + # direction going outwards (along the normal) + d = T(-Inf) + ndd = dot(normalize(dir), plane.normal) + saf = safety(plane, point) + ndd < 0 && saf > -kTolerance(T) && (d = -saf/ndd) + return d +end + +# Trd +function distanceToIn_trd(trd::Trd{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + x, y, z = point + dx, dy, dz = dir + distance::T = Inf + + inz = abs(z) < (trd.z - kTolerance(T)/2) + distx = trd.halfx1plusx2 - trd.fx * z + inx = (distx - abs(x)) * trd.calfx > kTolerance(T)/2 + disty = trd.halfy1plusy2 - trd.fy * z + iny = (disty - abs(y)) * trd.calfy > kTolerance(T)/2 + inside = inx && iny && inz + if inside + distance = -1. + end + done = inside + done && return distance + + okz = z * dz < 0 + okz &= !inz + if okz + distz = (abs(z) - trd.z)/abs(dz) + hitx = abs(x + distz * dx) + hity = abs(y + distz * dy) + okzt = z > (trd.z - kTolerance(T)/2) && hitx <= trd.x2 && hity <= trd.y2 + okzb = z < (-trd.z + kTolerance(T)/2) && hitx <= trd.x1 && hity <= trd.y1 + okz &= ( okzt || okzb ) + if okz + distance = distz + end + end + done |= okz + if done + return abs(distance) < kTolerance(T)/2 ? 0. : distance + end + + # hitting X faces? + okx = false + if !inx + okx, distx = faceIntersection(trd, point, dir, false, false, true) + if okx distance = distx end + okx, distx = faceIntersection(trd, point, dir, false, true, true) + if okx distance = distx end + end + done |= okx + if done + return abs(distance) < kTolerance(T)/2 ? 0. : distance + end + if !iny + oky, disty = faceIntersection(trd, point, dir, true, false, true) + if oky distance = disty end + oky, disty = faceIntersection(trd, point, dir, true, true, true) + if oky distance = disty end + end + return abs(distance) < kTolerance(T)/2 ? 0. : distance +end + +# Cone +function distanceToIn_cone(cone::Cone{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + x, y, z = point + dx, dy, dz = dir + (; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge, + outerSlope, outerOffset, innerSlope, innerOffset) = cone + + distance = T(Inf) + done = false + + # outside of Z range and going away? + distz = abs(z) - cone.z + done |= (distz > kTolerance(T)/2 && z * dz >= 0) || + (abs(distz) < kTolerance(T)/2 && z * dz > 0) + + # outside of outer tube and going away? + rsq = x * x + y * y + rmax = rmax1 == rmax2 ? rmax1 : outerOffset + outerSlope * z + crmax = rsq - rmax * rmax + rdotn = dx * x + dy * y + done |= crmax > kTolerance(T) * rmax && rdotn >= 0. + done && return distance + + # Next, check all dimensions of the tube, whether points are inside --> return -1 + distance = T(-1) + + # For points inside z-range, return -1 + inside = distz < -kTolerance(T)/2 + inside &= crmax < -kTolerance(T) * rmax + if rmin1 > 0. || rmin2 > 0. + rmin = rmin1 == rmin2 ? rmin1 : innerOffset + innerSlope * z + crmin = rsq - rmin * rmin + inside &= crmin > kTolerance(T) * rmin + end + if Δϕ < 2π + inside &= isInside(ϕWedge, x, y) + end + done |= inside + done && return distance + + # Next step: check if z-plane is the right entry point (both r,phi + # should be valid at z-plane crossing) + distance = T(Inf) + distz /= nonzero(abs(dz)) + hitx = x + distz * dx + hity = y + distz * dy + r2 = (hitx * hitx) + (hity * hity) + isHittingTopPlane = z >= cone.z - kTolerance(T)/2 && r2 <= rmax2 * rmax2 + kTolerance(T)/2 + isHittingBottomPlane = z <= -cone.z + kTolerance(T)/2 && r2 <= rmax1 * rmax1 + kTolerance(T)/2 + okz = isHittingTopPlane || isHittingBottomPlane + if rmin1 > 0 || rmin2 > 0 + isHittingTopPlane &= r2 >= rmin2 * rmin2 - kTolerance(T)/2 + isHittingBottomPlane &= r2 >= rmin1 * rmin1 - kTolerance(T)/2 + okz &= isHittingTopPlane || isHittingBottomPlane + end + if Δϕ < 2π + okz &= isInside(ϕWedge, hitx, hity) + end + + !done && okz && (distance = distz) + done |= okz + + # Next step: intersection of the trajectories with the two conical surfaces + distOuter = distanceToConicalSurface(cone, point, dir; distToIn=true, innerSurface=false) + distOuter < distance && (distance = distOuter) + if rmin1 > 0 || rmin2 > 0 + distInner = distanceToConicalSurface(cone, point, dir; distToIn=true, innerSurface=true) + distInner < distance && (distance = distInner) + end + # Find the intersection of trajectories with Phi planes + if Δϕ < 2π + p = Point2{T}(x,y) + d = Vector2{T}(dx,dy) + distPhi = phiPlaneIntersection(p, d, ϕWedge.along1, ϕWedge.normal1, toIn=true) + hit = point + distPhi * dir + !isnan(distPhi) && abs(z + distPhi * dz) <= cone.z && isInsideR(cone, hit) && isInside(ϕWedge, hit[1],hit[2], tol=kTolerance(T)) && distPhi < distance && (distance = distPhi) + if Δϕ != π + distPhi = phiPlaneIntersection(p, d, ϕWedge.along2, ϕWedge.normal2, toIn=true) + hit = point + distPhi * dir + !isnan(distPhi) && abs(z + distPhi * dz) <= cone.z && isInsideR(cone, hit) && isInside(ϕWedge, hit[1],hit[2], tol=kTolerance(T)) && distPhi < distance && (distance = distPhi) + end + end + return distance +end + +# Box +function distanceToIn_box(box::Box{T}, point::Point3{T}, direction::Vector3{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 + +#Tube +function distanceToIn_tube(tub::Tube{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + x, y, z = point + dx, dy, dz = dir + w = tub.ϕWedge + distance::T = Inf + done = false + + # outside of Z range and going away? + distz = abs(z) - tub.z + done |= distz > kTolerance(T)/2 && z * dz >= 0 + + # outside of outer tube and going away? + rsq = x*x + y*y + rdotn = dx*x + dy*y + done |= rsq > tub.tolIrmax2 && rdotn >= 0. + done && return distance + + # Next, check all dimensions of the tube, whether points are inside --> return -1 + distance = T(-1.) + + # For points inside z-range, return -1 + inside = distz < -kTolerance(T)/2 + inside &= rsq < tub.tolIrmax2 + if tub.rmin > 0. + inside &= rsq > tub.tolIrmin2 + end + if tub.Δϕ < 2π + inside &= isInside(w, x, y, tol=kTolerance(T)/2) + end + done |= inside + done && return distance + + # Next step: check if z-plane is the right entry point (both r,phi) + # should be valid at z-plane crossing) + distance = Inf + distz /= nonzero(abs(dz)) + hitx = x + distz * dx + hity = y + distz * dy + r2 = hitx * hitx + hity * hity; # radius of intersection with z-plane + okz = distz > -kTolerance(T)/2 && z * dz < 0 + okz &= (r2 <= tub.rmax2) + if tub.rmin > 0. + okz &= tub.rmin2 <= r2 + end + if tub.Δϕ < 2π + okz &= isInside(w, hitx, hity, tol=kTolerance(T)/2) + end + !done && okz && (distance = distz) + done |= okz + + # Next step: is in surface and going inside + onsurface = rsq >= tub.tolIrmax2 && rsq <= tub.tolOrmax2 && abs(z) < tub.z + kTolerance(T) && x*dx + y*dy <= 0. + if tub.rmin > 0. + onsurface |= rsq >= tub.tolOrmin2 && rsq <= tub.tolIrmin2 && abs(z) < tub.z + kTolerance(T) && x*dx + y*dy >= 0. + end + if tub.Δϕ < 2π + insector = isInside(w, x, y) + onsurface &= insector + end + !done && onsurface && (distance = 0.) + done |= onsurface + done && return distance + + # Next step: intersection of the trajectories with the two circles + p = Point2{T}(x,y) + d = Vector2{T}(dx,dy) + dist_rmax = circleIntersection(p, d, tub.rmax2) + dist_rmax >= kTolerance(T)/2 && abs(z + dist_rmax * dz) <= tub.z && isInside(w, p + dist_rmax * d) && dist_rmax < distance && (distance = dist_rmax) + if tub.rmin > 0. + dist_rmin = circleIntersection(p, d, tub.rmin2, largest=true ) + dist_rmin >= kTolerance(T)/2 && abs(z + dist_rmin * dz) <= tub.z && isInside(w, p + dist_rmin * d) && dist_rmin < distance && (distance = dist_rmin) + end + + # Calculate intersection between trajectory and the two phi planes + if tub.Δϕ < 2π + dist_phi = phiPlaneIntersection(p, d, w.along1, w.normal1, toIn=true) + hit = p + dist_phi * d + !isnan(dist_phi) && abs(z + dist_phi * dz) <= tub.z && isInsideR(tub, hit) && isInside(w, hit) && dist_phi < distance && (distance = dist_phi) + if tub.Δϕ != π + dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2, toIn=true) + hit = p + dist_phi * d + !isnan(dist_phi) && abs(z + dist_phi * dz) <= tub.z && isInsideR(tub, hit) && isInside(w, hit) && dist_phi < distance && (distance = dist_phi) + end + end + return distance +end + +#Volume +function distanceToIn_volume(agg::Aggregate{T}, point::Point3{T}, dir::Vector3{T}) where T<:AbstractFloat + distance = T(Inf) + for pvol in agg.pvolumes + lpoint = pvol.transformation * point + inout = inside(pvol.volume.shape, lpoint) + inout == kInside && return T(-1) + ldir = pvol.transformation * dir + inout == kSurface && normal(pvol.volume.shape, lpoint) ⋅ ldir < 0 && return T(0) + dist = distanceToIn(pvol.volume.shape, lpoint, ldir) + dist < distance && (distance = dist) + end + return distance +end + +#Polycone +function distanceToIn_polycone(pcone::Polycone{T,N}, point::Point3{T}, dir::Vector3{T})::T where {T,N} + z = point[3] + dz = dir[3] + distance = Inf + increment = dz > 0 ? 1 : -1 + abs(dz) < kTolerance(T) && (increment = 0) + index = getSectionIndex(pcone, z) + index == -1 && (index = 1) + index == -2 && (index = N) + while true + distance = distanceToIn(pcone.sections[index], point - Vector3{T}(0,0,pcone.zᵢ[index]), dir) + (distance < Inf || increment == 0) && break + index += increment + (index < 1 || index > N) && break + end + return distance +end + +#CutTube +function distanceToIn_cuttube(ctube::CutTube{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + bot, top = ctube.planes + + distance = T(Inf) + # Points already inside have to return negative distance + pinside = inside(ctube, point) + pinside == kInside && return T(-1) + + # Compute distance to cut planes + d0 = dot(dir, bot.normal) > 0 ? T(-Inf) : distanceToIn(bot, point, dir) + d1 = dot(dir, top.normal) > 0 ? T(-Inf) : distanceToIn(top, point, dir) + dplanes = max(d0, d1) + splanes = max(safety(top,point), safety(bot, point)) + + # Mark tracks hitting the planes + hitplanes = dplanes < T(Inf) && splanes > -kTolerance(T) + #!hitplanes && pinside != kInside && return distance + + # Propagate with dplanes only the particles that are hitting + !hitplanes && (dplanes = 0) + propagated = point + dplanes * dir + + # Check now which of the propagated points already entering the solid + intube = inside(ctube.tube, propagated) + done = hitplanes && intube != kOutside + done && return abs(dplanes) < kTolerance(T) ? 0 : dplanes + + # The limit distance for tube crossing cannot exceed the distance to exiting the cut planes + dexit = min(distanceToOut(bot, propagated, dir), distanceToOut(top, propagated, dir)) + + # Compute distance to tube + dtube = distanceToIn(ctube.tube, propagated, dir) + dtube < 0 && (dtube = T(0)) + dexit < dtube && (dtube = T(Inf)) + distance = dtube + dplanes + distance < kTolerance(T) && (distance = T(0)) + return distance +end + +#Trap +function distanceToIn_trap(trap::Trap{T}, point::Point3{T}, dir::Vector3{T}) where T<:AbstractFloat + (; z, planes) = trap + dz = dir[3] + z = point[3] + + # step 1.a: input particle is moving away --> return infinity + signZdir = copysign(1,dz) + max = signZdir * trap.z - z + signZdir * max < kTolerance(T)/2 && return Inf + + # step 1.b General case: + smax = max/dz + smin = -(signZdir * trap.z + z)/dz + + # Step 2: find distances for intersections with side planes. + done = false + for i in 1:4 + ndir = dot(dir, planes[i].normal) + safe = safety(planes[i], point) + dist = -safe/ndir + + done |= (safe > kTolerance(T)/2 && ndir >= 0) || (safe > -kTolerance(T) && ndir > 0) + + posPoint = safe > -kTolerance(T)/2 + posDir = ndir > 0 + !posPoint && posDir && dist < smax && (smax = dist) + posPoint && !posDir && dist > smin && (smin = dist) + end + done && return T(Inf) + distance = smin <= smax ? smin : Inf + distance < -kTolerance(T)/2 && return T(-1) + return distance +end diff --git a/src/DistanceOut.jl b/src/DistanceOut.jl new file mode 100644 index 0000000..0e6e0a0 --- /dev/null +++ b/src/DistanceOut.jl @@ -0,0 +1,443 @@ +#Boolean +function distanceToOut_booleanunion(shape, point, dir::Vector3{T})::T where T + (; 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_booleanintersection(shape, point, dir::Vector3{T})::T where T + (; left, right, transformation) = shape + distA = distanceToOut(left, point, dir) + distB = distanceToOut(right, transformation * point, transformation * dir) + return min(distA, distB) +end +function distanceToOut_booleansubtraction(shape, point, dir::Vector3{T})::T where T + (; left, right, transformation) = shape + distA = distanceToOut(left, point, dir) + distB = distanceToIn(right, transformation * point, transformation * dir) + return min(distA, distB) +end + +#Plane +function distanceToOut_plane(plane::Plane{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + #The function returns infinity if the plane is not hit from inside, negative + #if the point is outside + d = T(Inf) + ndd = dot(normalize(dir), plane.normal) + saf = safety(plane, point) + saf > kTolerance(T) && (d = -T(Inf)) + ndd > 0 && saf < kTolerance(T) && (d = -saf/ndd) + return d +end + +#Trd +function distanceToOut_trd(trd::Trd{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + x, y, z = point + dx, dy, dz = dir + distance = 0.0 + + # hit top Z face? + trd.calfx + trd.calfy + safz = trd.z - abs(z) + out = safz < -kTolerance(T)/2 + distx = trd.halfx1plusx2 - trd.fx * z + out |= (distx - abs(x)) * trd.calfx < -kTolerance(T)/2 + disty = trd.halfy1plusy2 - trd.fy * z + out |= (disty - abs(y)) * trd.calfy < -kTolerance(T)/2 + out && return -1. + if dz > 0. + distz = (trd.z - z)/abs(dz) + hitx = abs(x + distz * dx) + hity = abs(y + distz * dy) + if hitx <= trd.x2 && hity <= trd.y2 + distance = distz + abs(distance) < -kTolerance(T)/2 && return 0. + end + end + # hit bottom Z face? + if dz < 0. + distz = (trd.z + z)/abs(dz) + hitx = abs(x + distz * dx) + hity = abs(y + distz * dy) + if hitx <= trd.x1 && hity <= trd.y1 + distance = distz + abs(distance) < -kTolerance(T)/2 && return 0. + end + end + + # hitting X faces? + okx, distx = faceIntersection(trd, point, dir, false, false, false) + if okx + distance = distx + return distx < kTolerance(T)/2 ? 0.0 : distx + end + okx, distx = faceIntersection(trd, point, dir, false, true, false) + if okx + distance = distx + return distx < kTolerance(T)/2 ? 0.0 : distx + end + + # hitting Y faces? + oky, disty = faceIntersection(trd, point, dir, true, false, false) + if oky + distance = disty + return disty < kTolerance(T)/2 ? 0.0 : disty + end + oky, disty = faceIntersection(trd, point, dir, true, true, false) + if oky + distance = disty + return disty < kTolerance(T)/2 ? 0.0 : disty + end + return distance +end +function distanceToOut_trd(trd::TTrd{T}, point::Point3{T}, direction::Vector3{T}) where T<:AbstractFloat + for triangle in trd.triangles + dist, ok = intersect(point, direction, triangle) + ok && return dist + end + -1.0 +end + +#Cone +function distanceToOut_cone(cone::Cone{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + x, y, z = point + dx, dy, dz = dir + (; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge, secRMax, secRMin, + outerSlope, outerOffset, innerSlope, innerOffset) = cone + + distance::T = -1. + done = false + + #---First we check all dimensions of the cone, whether points are outside (return -1) + distz = cone.z - abs(z) + done |= distz < -kTolerance(T)/2 + + r2 = x * x + y * y + rmax = rmax1 == rmax2 ? rmax1 : outerOffset + outerSlope * z + crmax = r2 - rmax * rmax + + # if outside of Rmax, return -1 + done |= crmax > kTolerance(T) * rmax + done && return distance + + if rmin1 > 0. || rmin2 > 0. + rmin = rmin1 == rmin2 ? rmin1 : innerOffset + innerSlope * z + crmin = r2 - rmin * rmin + done |= crmin < -kTolerance(T) * rmin + done && return distance + end + + if Δϕ < 2π + done |= isOutside(ϕWedge, x, y) + done && return distance + end + + # OK, since we're here, then distance must be non-negative, and the + # smallest of possible intersections + !done && (distance = Inf) + + invdirz = 1. / nonzero(dz) + distz = dz >= 0 ? (cone.z - z) * invdirz : (-cone.z - z) * invdirz + dz != 0 && distz < distance && (distance = distz) + + # Find the intersection of the trajectories with the conical surfaces + distOuter = distanceToConicalSurface(cone, point, dir; distToIn=false, innerSurface=false) + distOuter < distance && (distance = distOuter) + if rmin1 > 0. || rmin2 > 0. + distInner = distanceToConicalSurface(cone, point, dir; distToIn=false, innerSurface=true) + distInner < distance && (distance = distInner) + end + + # Find the intersection of trajectories with Phi planes + if Δϕ < 2π + p = Point2{T}(x,y) + d = Vector2{T}(dx,dy) + isOnStartPhi = isOnSurface(ϕWedge.along1, ϕWedge.normal1, x, y) + isOnEndPhi = isOnSurface(ϕWedge.along2, ϕWedge.normal2, x, y) + cond = (isOnStartPhi && d ⋅ ϕWedge.normal1 < 0.) || (isOnEndPhi && d ⋅ ϕWedge.normal2 < 0.) + cond && (distance = 0.) + if Δϕ < π + dist_phi = phiPlaneIntersection(p, d, ϕWedge.along1, ϕWedge.normal1) + dist_phi < distance && (distance = dist_phi) + dist_phi = phiPlaneIntersection(p, d, ϕWedge.along2, ϕWedge.normal2) + dist_phi < distance && (distance = dist_phi) + elseif Δϕ == π + dist_phi = phiPlaneIntersection(p, d, ϕWedge.along2, ϕWedge.normal2) + dist_phi < distance && (distance = dist_phi) + else + dist_phi = phiPlaneIntersection(p, d, ϕWedge.along1, ϕWedge.normal1, ϕgtπ=true) + dist_phi < distance && (distance = dist_phi) + dist_phi = phiPlaneIntersection(p, d, ϕWedge.along2, ϕWedge.normal2, ϕgtπ=true) + dist_phi < distance && (distance = dist_phi) + end + end + return distance +end + +#Box +function distanceToOut_box(box::Box{T}, point::Point3{T}, direction::Vector3{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 + +#Trap +function distanceToOut_trap(trap::Trap{T}, point::Point3{T}, dir::Vector3{T}) where T<:AbstractFloat + (; z, planes) = trap + dz = dir[3] + + #step 0: if point is outside any plane --> return -1, otherwise initialize at Infinity + outside = abs(point[3]) > z + kTolerance(T)/2 + outside && return T(-1) + distance = T(Inf) + + #Step 1: find range of distances along dir between Z-planes (smin, smax) + dz != 0 && (distance = (copysign(z,dz) - point[3]) / dz) + + #Step 2: find distances for intersections with side planes. + for i in 1:4 + dist = T(Inf) + ndir = dot(dir, planes[i].normal) + safe = safety(planes[i], point) + + safe > kTolerance(T) && return T(-1) + ndir > 0 && safe < kTolerance(T) && (dist = -safe/ndir) + dist < distance && (distance = dist) + end + return distance +end + +#Tube +function distanceToOut_tube(tub::Tube{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + x, y, z = point + dx, dy, dz = dir + w = tub.ϕWedge + distance::T = -1. + done = false + + #First we check all dimensions of the tube, whether points are outside (return -1) + distz = tub.z - abs(z) + done |= distz < -kTolerance(T)/2 + + rsq = x*x + y*y + rdotn = dx*x + dy*y + crmax = rsq - tub.rmax2 + + #if outside of Rmax, return -1 + done |= crmax > kTolerance(T) * tub.rmax + done && return distance + + if tub.rmin > 0. + crmin = rsq - tub.rmin2 + done |= crmin < -kTolerance(T) * tub.rmin + done && return distance + end + if tub.Δϕ < 2π + done |= isOutside(w, x, y) + done && return distance + end + + #OK, since we're here, then distance must be non-negative, and the + #smallest of possible intersections + !done && (distance = Inf) + + invdirz = 1. / nonzero(dz) + distz = dz >= 0 ? (tub.z - z) * invdirz : (-tub.z - z) * invdirz + !done && dz != 0 && distz < distance && (distance = distz) + + #Find the intersection of the trajectories with the two circles. + #Here I compute values used in both rmin and rmax calculations. + + invnsq = 1. / nonzero(1. - dz*dz) + b = invnsq * rdotn + + #rmin + if tub.rmin > 0. + crmin *= invnsq + delta = b * b - crmin + dist_rmin = -b + (delta > 0. ? -sqrt(delta) : 0.) + delta > 0. && dist_rmin >= -kTolerance(T)/2 && dist_rmin < distance && (distance = dist_rmin) + end + + #rmax + crmax *= invnsq + delta = b * b - crmax + dist_rmax = -b + (delta >= 0. ? sqrt(delta) : 0.) + dist_rmax >= -kTolerance(T)/2 && dist_rmax < distance && (distance = dist_rmax) + + #Phi planes + + if tub.Δϕ < 2π + w = tub.ϕWedge + p = Point2{T}(x,y) + d = Vector2{T}(dx,dy) + if tub.Δϕ < π + dist_phi = phiPlaneIntersection(p, d, w.along1, w.normal1) + dist_phi < distance && (distance = dist_phi) + dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2) + dist_phi < distance && (distance = dist_phi) + elseif tub.Δϕ == π + dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2) + dist_phi < distance && (distance = dist_phi) + else + dist_phi = phiPlaneIntersection(p, d, w.along1, w.normal1, ϕgtπ=true) + dist_phi < distance && (distance = dist_phi) + dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2, ϕgtπ=true) + dist_phi < distance && (distance = dist_phi) + end + end + return distance +end + +#Polycone +function distanceToOut_polycone(pcone::Polycone{T,N}, point::Point3{T}, dir::Vector3{T})::T where {T,N} + z = point[3] + dz = dir[3] + distance = Inf + if N == 1 + return distanceToOut(pcone.sections[1], point - Vector3{T}(0,0,pcone.zᵢ[1]), dir) + end + indexLow = getSectionIndex(pcone, z - kTolerance(T)) + indexHigh = getSectionIndex(pcone, z + kTolerance(T)) + if indexLow < 0 && indexHigh < 0 + return -1. + elseif indexLow < 0 && indexHigh > 0 + index = indexHigh + inside(pcone.sections[index], point - Vector3{T}(0,0,pcone.zᵢ[index])) == kOutside && return -1. + elseif indexLow != indexHigh && indexLow > 0 + isin = inside(pcone.sections[indexLow], point - Vector3{T}(0,0,pcone.zᵢ[indexLow])) + index = isin == kOutside ? indexHigh : indexLow + else + index = indexLow + end + + if index < 0 + return 0. + else + inside(pcone.sections[index], point - Vector3{T}(0,0,pcone.zᵢ[index])) == kOutside && return -1. + end + + distance = T(0) + increment = dz > 0 ? 1 : -1 + abs(dz) < kTolerance(T) && (increment = 0) + pn = point + istep = 0 + while true + #section surface case + if distance != 0. || istep < 2 + pn = point + distance * dir - Vector3{T}(0,0,pcone.zᵢ[index]) + inside(pcone.sections[index], pn ) == kOutside && break + else + pn -= Vector3{T}(0,0,pcone.zᵢ[index]) + end + istep += 1 + dist = distanceToOut(pcone.sections[index], pn, dir) + dist == -1. && return distance + if abs(dist) < kTolerance(T)/2 + index1 = index + if index > 1 && index < N + index1 += increment + else + index == 1 && increment > 0 && (index1 += increment) + index == N && increment < 0 && (index1 += increment) + end + pte = point + (distance + dist) * dir - Vector3{T}(0,0,pcone.zᵢ[index1]) + (inside(pcone.sections[index], pte ) == kOutside || increment == 0) && break + end + distance += dist + index += increment + (increment == 0 || index < 1 || index > N ) && break + end + return distance +end + +#CutTube +function distanceToOut_cuttube(ctube::CutTube{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat + bot, top = ctube.planes + #Compute distance to cut planes + distance = min(distanceToOut(bot, point, dir), distanceToOut(top, point, dir)) + #Compute distance to tube + dtube = distanceToOut(ctube.tube, point, dir) + dtube < distance && (distance = dtube) + return distance +end + + +#Volume +function distanceToOut_volume(agg::Aggregate{T}, point::Point3{T}, dir::Vector3{T}) where T<:AbstractFloat + for pvol in agg.pvolumes + lpoint = pvol.transformation * point + inout = inside(pvol.volume.shape, lpoint) + inout == kSurface && return T(0) + inout == kInside && return distanceToOut(pvol.volume.shape, lpoint, pvol.transformation * dir) + end + return T(-1) +end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 8cc19b1..1ac2cd1 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -34,6 +34,9 @@ include("Volume.jl") include("BVH.jl") include("Navigators.jl") include("GDML.jl") +include("./DistanceIn.jl") +include("./DistanceOut.jl") +include("./Distance.jl") include("Benchmark.jl") #include("CuGeom.jl") # PackageCompiler fails? diff --git a/src/Polycone.jl b/src/Polycone.jl index d8b1579..5faef36 100644 --- a/src/Polycone.jl +++ b/src/Polycone.jl @@ -241,4 +241,4 @@ function distanceToOut(pcone::Polycone{T,N}, point::Point3{T}, dir::Vector3{T}): (increment == 0 || index < 1 || index > N ) && break end return distance -end \ No newline at end of file +end diff --git a/src/Trd.jl b/src/Trd.jl index b112223..655019a 100644 --- a/src/Trd.jl +++ b/src/Trd.jl @@ -363,4 +363,4 @@ function extent(trd::TTrd{T})::Tuple{Point3{T},Point3{T}} where T<:AbstractFloat p = Point3{T}( max(trd.x1, trd.x2), max(trd.y1, trd.y2), trd.z) ( -p, p ) end -=# \ No newline at end of file +=# diff --git a/src/Tube.jl b/src/Tube.jl index 938034a..140f36b 100644 --- a/src/Tube.jl +++ b/src/Tube.jl @@ -295,173 +295,6 @@ function circleIntersection(point::Point2{T}, dir::Vector2{T}, R²::T; largest:: return dist end -function distanceToOut(tub::Tube{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat - x, y, z = point - dx, dy, dz = dir - w = tub.ϕWedge - distance::T = -1. - done = false - - #---First we check all dimensions of the tube, whether points are outside (return -1) - distz = tub.z - abs(z) - done |= distz < -kTolerance(T)/2 - - rsq = x*x + y*y - rdotn = dx*x + dy*y - crmax = rsq - tub.rmax2 - - # if outside of Rmax, return -1 - done |= crmax > kTolerance(T) * tub.rmax - done && return distance - - if tub.rmin > 0. - crmin = rsq - tub.rmin2 - done |= crmin < -kTolerance(T) * tub.rmin - done && return distance - end - if tub.Δϕ < 2π - done |= isOutside(w, x, y) - done && return distance - end - - # OK, since we're here, then distance must be non-negative, and the - # smallest of possible intersections - !done && (distance = Inf) - - invdirz = 1. / nonzero(dz) - distz = dz >= 0 ? (tub.z - z) * invdirz : (-tub.z - z) * invdirz - !done && dz != 0 && distz < distance && (distance = distz) - - # Find the intersection of the trajectories with the two circles. - # Here I compute values used in both rmin and rmax calculations. - - invnsq = 1. / nonzero(1. - dz*dz) - b = invnsq * rdotn - - # rmin - if tub.rmin > 0. - crmin *= invnsq - delta = b * b - crmin - dist_rmin = -b + (delta > 0. ? -sqrt(delta) : 0.) - delta > 0. && dist_rmin >= -kTolerance(T)/2 && dist_rmin < distance && (distance = dist_rmin) - end - - # rmax - crmax *= invnsq - delta = b * b - crmax - dist_rmax = -b + (delta >= 0. ? sqrt(delta) : 0.) - dist_rmax >= -kTolerance(T)/2 && dist_rmax < distance && (distance = dist_rmax) - - # Phi planes - - if tub.Δϕ < 2π - w = tub.ϕWedge - p = Point2{T}(x,y) - d = Vector2{T}(dx,dy) - if tub.Δϕ < π - dist_phi = phiPlaneIntersection(p, d, w.along1, w.normal1) - dist_phi < distance && (distance = dist_phi) - dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2) - dist_phi < distance && (distance = dist_phi) - elseif tub.Δϕ == π - dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2) - dist_phi < distance && (distance = dist_phi) - else - dist_phi = phiPlaneIntersection(p, d, w.along1, w.normal1, ϕgtπ=true) - dist_phi < distance && (distance = dist_phi) - dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2, ϕgtπ=true) - dist_phi < distance && (distance = dist_phi) - end - end - return distance -end - -function distanceToIn(tub::Tube{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat - x, y, z = point - dx, dy, dz = dir - w = tub.ϕWedge - distance::T = Inf - done = false - - # outside of Z range and going away? - distz = abs(z) - tub.z - done |= distz > kTolerance(T)/2 && z * dz >= 0 - - # outside of outer tube and going away? - rsq = x*x + y*y - rdotn = dx*x + dy*y - done |= rsq > tub.tolIrmax2 && rdotn >= 0. - done && return distance - - # Next, check all dimensions of the tube, whether points are inside --> return -1 - distance = T(-1.) - - # For points inside z-range, return -1 - inside = distz < -kTolerance(T)/2 - inside &= rsq < tub.tolIrmax2 - if tub.rmin > 0. - inside &= rsq > tub.tolIrmin2 - end - if tub.Δϕ < 2π - inside &= isInside(w, x, y, tol=kTolerance(T)/2) - end - done |= inside - done && return distance - - # Next step: check if z-plane is the right entry point (both r,phi) - # should be valid at z-plane crossing) - distance = Inf - distz /= nonzero(abs(dz)) - hitx = x + distz * dx - hity = y + distz * dy - r2 = hitx * hitx + hity * hity; # radius of intersection with z-plane - okz = distz > -kTolerance(T)/2 && z * dz < 0 - okz &= (r2 <= tub.rmax2) - if tub.rmin > 0. - okz &= tub.rmin2 <= r2 - end - if tub.Δϕ < 2π - okz &= isInside(w, hitx, hity, tol=kTolerance(T)/2) - end - !done && okz && (distance = distz) - done |= okz - - # Next step: is in surface and going inside - onsurface = rsq >= tub.tolIrmax2 && rsq <= tub.tolOrmax2 && abs(z) < tub.z + kTolerance(T) && x*dx + y*dy <= 0. - if tub.rmin > 0. - onsurface |= rsq >= tub.tolOrmin2 && rsq <= tub.tolIrmin2 && abs(z) < tub.z + kTolerance(T) && x*dx + y*dy >= 0. - end - if tub.Δϕ < 2π - insector = isInside(w, x, y) - onsurface &= insector - end - !done && onsurface && (distance = 0.) - done |= onsurface - done && return distance - - # Next step: intersection of the trajectories with the two circles - p = Point2{T}(x,y) - d = Vector2{T}(dx,dy) - dist_rmax = circleIntersection(p, d, tub.rmax2) - dist_rmax >= kTolerance(T)/2 && abs(z + dist_rmax * dz) <= tub.z && isInside(w, p + dist_rmax * d) && dist_rmax < distance && (distance = dist_rmax) - if tub.rmin > 0. - dist_rmin = circleIntersection(p, d, tub.rmin2, largest=true ) - dist_rmin >= kTolerance(T)/2 && abs(z + dist_rmin * dz) <= tub.z && isInside(w, p + dist_rmin * d) && dist_rmin < distance && (distance = dist_rmin) - end - - # Calculate intersection between trajectory and the two phi planes - if tub.Δϕ < 2π - dist_phi = phiPlaneIntersection(p, d, w.along1, w.normal1, toIn=true) - hit = p + dist_phi * d - !isnan(dist_phi) && abs(z + dist_phi * dz) <= tub.z && isInsideR(tub, hit) && isInside(w, hit) && dist_phi < distance && (distance = dist_phi) - if tub.Δϕ != π - dist_phi = phiPlaneIntersection(p, d, w.along2, w.normal2, toIn=true) - hit = p + dist_phi * d - !isnan(dist_phi) && abs(z + dist_phi * dz) <= tub.z && isInsideR(tub, hit) && isInside(w, hit) && dist_phi < distance && (distance = dist_phi) - end - end - return distance -end function GeometryBasics.coordinates(tub::Tube{T}, facets=36) where {T<:AbstractFloat} issector = tub.Δϕ < 2π From 2f2a458d8e3cb7c409f917e43a8c38319d246eee Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 16:54:03 -0400 Subject: [PATCH 02/16] retain original performance --- src/Boolean.jl | 93 ---------------------------------------- src/DistanceIn.jl | 104 ++++++++++++++++++++++----------------------- src/DistanceOut.jl | 6 +-- 3 files changed, 55 insertions(+), 148 deletions(-) diff --git a/src/Boolean.jl b/src/Boolean.jl index d853501..79ec626 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -162,96 +162,3 @@ end function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 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/DistanceIn.jl b/src/DistanceIn.jl index b59ddb9..a1c8c94 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -5,58 +5,58 @@ function distanceToIn_booleanunion(shape::BooleanUnion{T, SL, SR}, point::Point3 distB = distanceToIn(right, transformation * point, transformation * dir) return min(distA, distB) end -# function distanceToIn_booleanintersection(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_booleansubtraction(shape, point::Point3{T}, dir) where T +function distanceToIn_booleanintersection(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_booleansubtraction(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR} (; left, right, transformation) = shape lpoint = transformation * point diff --git a/src/DistanceOut.jl b/src/DistanceOut.jl index 0e6e0a0..6ab14cf 100644 --- a/src/DistanceOut.jl +++ b/src/DistanceOut.jl @@ -1,5 +1,5 @@ #Boolean -function distanceToOut_booleanunion(shape, point, dir::Vector3{T})::T where T +function distanceToOut_booleanunion(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) @@ -54,13 +54,13 @@ function distanceToOut_booleanunion(shape, point, dir::Vector3{T})::T where T end end end -function distanceToOut_booleanintersection(shape, point, dir::Vector3{T})::T where T +function distanceToOut_booleanintersection(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_booleansubtraction(shape, point, dir::Vector3{T})::T where T +function distanceToOut_booleansubtraction(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) From 9de89eb3a6959aedb2180a505220e25ae89ed870 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 17:02:58 -0400 Subject: [PATCH 03/16] recover original performance --- src/Boolean.jl | 63 ------------- src/Box.jl | 13 --- src/Cone.jl | 32 ------- src/CutTube.jl | 17 ---- src/Geom4hep.jl | 1 + src/Inside.jl | 228 ++++++++++++++++++++++++++++++++++++++++++++++++ src/Polycone.jl | 17 ---- src/Trap.jl | 14 --- src/Trd.jl | 20 ----- src/Tube.jl | 30 ------- src/Volume.jl | 15 ++-- 11 files changed, 237 insertions(+), 213 deletions(-) create mode 100644 src/Inside.jl diff --git a/src/Boolean.jl b/src/Boolean.jl index 79ec626..538903e 100644 --- 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 diff --git a/src/Box.jl b/src/Box.jl index b4a9552..ca17cb2 100644 --- a/src/Box.jl +++ b/src/Box.jl @@ -28,19 +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 safetyToOut(box::Box{T}, point::Point3{T}) where T<:AbstractFloat minimum(box.fDimensions - abs.(point)) end diff --git a/src/Cone.jl b/src/Cone.jl index 3caac5e..8746a45 100644 --- a/src/Cone.jl +++ b/src/Cone.jl @@ -93,38 +93,6 @@ function extent(cone::Cone{T})::Tuple{Point3{T},Point3{T}} where T<:AbstractFloa extent(Tube{T}(min(rmin1,rmin2), max(rmax1,rmax2),z, ϕ₀, Δϕ)) end -function inside(cone::Cone{T}, point::Point3{T})::Int64 where T<:AbstractFloat - x, y, z = point - (; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge, - outerSlope, outerOffset, innerSlope, innerOffset) = cone - - # Check Z - outside = abs(z) > cone.z + kTolerance(T)/2 - outside && return kOutside - cinside = abs(z) < cone.z - kTolerance(T)/2 - - # Check on RMax - r2 = x * x + y * y - rmax = rmax1 == rmax2 ? rmax1 : outerOffset + outerSlope * z - outside |= r2 > rmax * rmax + kTolerance(T) * rmax - outside && return kOutside - cinside &= r2 < rmax * rmax - kTolerance(T) * rmax - - # Check on RMin - if rmin1 > 0. || rmin2 > 0. - rmin = rmin1 == rmin2 ? rmin1 : innerOffset + innerSlope * z - outside |= r2 <= rmin * rmin - kTolerance(T) * rmin - outside && return kOutside - cinside &= r2 > rmin * rmin + kTolerance(T) * rmin - end - - # Check on Phi - if Δϕ < 2π - outside |= isOutside(ϕWedge, x, y) - cinside &= isInside(ϕWedge, x, y) - end - return outside ? kOutside : cinside ? kInside : kSurface -end function safetyToIn(cone::Cone{T}, point::Point3{T}) where T<:AbstractFloat x,y,z = point diff --git a/src/CutTube.jl b/src/CutTube.jl index ff6b730..6d98ef2 100644 --- a/src/CutTube.jl +++ b/src/CutTube.jl @@ -95,23 +95,6 @@ function extent(ctube::CutTube{T})::Tuple{Point3{T},Point3{T}} where T<:Abstract return (Point3{T}(aMin), Point3{T}(aMax)) end -function inside(ctube::CutTube{T}, point::Point3{T})::Int64 where T<:AbstractFloat - bot, top = ctube.planes - tub = ctube.tube - x, y, z = point - - # Check the cut planes first - a = safety(bot, point) - b = safety(top, point) - outside = a > kTolerance(T)/2 || b > kTolerance(T)/2 - outside && return kOutside - cinside = a < -kTolerance(T)/2 && b < -kTolerance(T)/2 - inplanes = cinside ? kInside : kSurface - # Check the tube - intube = inside(tub, point) - return inplanes == kSurface && intube != kOutside ? inplanes : intube -end - function safetyToIn(ctube::CutTube{T}, point::Point3{T}) where T<:AbstractFloat end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 1ac2cd1..41984e6 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -37,6 +37,7 @@ include("GDML.jl") include("./DistanceIn.jl") include("./DistanceOut.jl") include("./Distance.jl") +include("./Inside.jl") include("Benchmark.jl") #include("CuGeom.jl") # PackageCompiler fails? diff --git a/src/Inside.jl b/src/Inside.jl new file mode 100644 index 0000000..11ca99d --- /dev/null +++ b/src/Inside.jl @@ -0,0 +1,228 @@ +function inside(shape, point) + @nospecialize shape + if shape isa Trap + inside_trap(shape, point) + elseif shape isa Trd + inside_trd(shape, point) + elseif shape isa Cone + inside_cone(shape, point) + elseif shape isa Box + inside_box(shape, point) + elseif shape isa Tube + inside_tube(shape, point) + elseif shape isa Polycone + inside_polycone(shape, point) + elseif shape isa CutTube + inside_cuttube(shape, point) + elseif shape isa BooleanUnion + inside_booleanunion(shape, point) + elseif shape isa BooleanSubtraction + inside_booleansubtraction(shape, point) + elseif shape isa BooleanIntersection + inside_booleanintersection(shape, point) + end +end +function inside_booleanunion(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_booleanintersection(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_booleansubtraction(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 inside_trap(trap::Trap{T}, point::Point3{T}) where T<:AbstractFloat + z = point[3] + planes = trap.planes + # inside z? + outside = abs(z) > trap.z + kTolerance(T)/2 + cinside = abs(z) < trap.z - kTolerance(T)/2 + + for i in 1:4 + s = safety(planes[i], point) # positive if the point is on the outside halfspace, negative otherwise. + outside |= s > kTolerance(T)/2 + cinside &= s < -kTolerance(T)/2 + end + cinside ? kInside : outside ? kOutside : kSurface +end +function inside_cone(cone::Cone{T}, point::Point3{T}) where T<:AbstractFloat + x, y, z = point + (; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge, + outerSlope, outerOffset, innerSlope, innerOffset) = cone + + # Check Z + outside = abs(z) > cone.z + kTolerance(T)/2 + outside && return kOutside + cinside = abs(z) < cone.z - kTolerance(T)/2 + + # Check on RMax + r2 = x * x + y * y + rmax = rmax1 == rmax2 ? rmax1 : outerOffset + outerSlope * z + outside |= r2 > rmax * rmax + kTolerance(T) * rmax + outside && return kOutside + cinside &= r2 < rmax * rmax - kTolerance(T) * rmax + + # Check on RMin + if rmin1 > 0. || rmin2 > 0. + rmin = rmin1 == rmin2 ? rmin1 : innerOffset + innerSlope * z + outside |= r2 <= rmin * rmin - kTolerance(T) * rmin + outside && return kOutside + cinside &= r2 > rmin * rmin + kTolerance(T) * rmin + end + + # Check on Phi + if Δϕ < 2π + outside |= isOutside(ϕWedge, x, y) + cinside &= isInside(ϕWedge, x, y) + end + return outside ? kOutside : cinside ? kInside : kSurface +end +function inside_polycone(pcone::Polycone{T,N}, point::Point3{T}) where {T,N} + x, y, z = point + indexLow = getSectionIndex(pcone, z - kTolerance(T)) + indexHigh = getSectionIndex(pcone, z + kTolerance(T)) + indexLow < 0 && indexHigh < 0 && return kOutside + indexLow < 0 && indexHigh == 1 && return inside(pcone.sections[1], point-Vector3{T}(0,0,pcone.zᵢ[1])) + indexHigh < 0 && indexLow == N && return inside(pcone.sections[N], point-Vector3{T}(0,0,pcone.zᵢ[N])) + if indexLow == indexHigh + return inside(pcone.sections[indexLow], point-Vector3{T}(0,0,pcone.zᵢ[indexLow])) + else + insideLow = inside(pcone.sections[indexLow], point-Vector3{T}(0,0,pcone.zᵢ[indexLow])) + insideHigh = inside(pcone.sections[indexHigh], point-Vector3{T}(0,0,pcone.zᵢ[indexHigh])) + insideLow == kSurface && insideHigh == kSurface && return kInside + insideLow == kOutside && insideHigh == kOutside && return kOutside + return kSurface + end +end +function inside_box(box::Box{T}, point::Point3{T}) 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 inside_trd(trd::Trd{T}, point::Point3{T}) where T<:AbstractFloat + x, y, z = point + + # inside z? + outside = abs(z) > (trd.z + kTolerance(T)/2) + inside = abs(z) < (trd.z - kTolerance(T)/2) + + # inside x? + c = cross(abs(x)-trd.x1, z+trd.z, trd.x2-trd.x1, 2*trd.z) + outside |= (c < -kTolerance(T)/2) + inside &= (c > kTolerance(T)/2) + + # inside y? + c = cross(abs(y)-trd.y1, z+trd.z, trd.y2-trd.y1, 2*trd.z) + outside |= (c < -kTolerance(T)/2) + inside &= (c > kTolerance(T)/2) + + # return + inside ? kInside : outside ? kOutside : kSurface +end +function inside_cuttube(ctube::CutTube{T}, point::Point3{T}) where T<:AbstractFloat + bot, top = ctube.planes + tub = ctube.tube + x, y, z = point + + # Check the cut planes first + a = safety(bot, point) + b = safety(top, point) + outside = a > kTolerance(T)/2 || b > kTolerance(T)/2 + outside && return kOutside + cinside = a < -kTolerance(T)/2 && b < -kTolerance(T)/2 + inplanes = cinside ? kInside : kSurface + # Check the tube + intube = inside(tub, point) + return inplanes == kSurface && intube != kOutside ? inplanes : intube +end +function inside_tube(tub::Tube{T}, point::Point3{T}) where T<:AbstractFloat + x, y, z = point + + # Check Z + outside = abs(z) > tub.z + kTolerance(T)/2 + outside && return kOutside + cinside = abs(z) < tub.z - kTolerance(T)/2 + + # Check on RMax + r2 = x * x + y * y + outside |= r2 > tub.rmax2 + kTolerance(T) * tub.rmax + outside && return kOutside + cinside &= r2 < tub.rmax2 - kTolerance(T) * tub.rmax + + # Check on RMin + if tub.rmin > 0. + outside |= r2 <= tub.rmin2 - kTolerance(T) * tub.rmin + outside && return kOutside + cinside &= r2 > tub.rmin2 + kTolerance(T) * tub.rmin + end + + # Check on Phi + if tub.Δϕ < 2π + outside |= isOutside(tub.ϕWedge, x, y) + cinside &= isInside(tub.ϕWedge, x, y) + end + + return outside ? kOutside : cinside ? kInside : kSurface +end diff --git a/src/Polycone.jl b/src/Polycone.jl index 5faef36..131c4cf 100644 --- a/src/Polycone.jl +++ b/src/Polycone.jl @@ -81,23 +81,6 @@ function extent(pcone::Polycone{T,N})::Tuple{Point3{T},Point3{T}} where {T,N} return (low + Vector3{T}(0, 0, pcone.zᵢ[1] - pcone.sections[1].z), high + Vector3{T}(0, 0, pcone.zᵢ[N] + pcone.sections[N].z)) end -function inside(pcone::Polycone{T,N}, point::Point3{T})::Int64 where {T,N} - x, y, z = point - indexLow = getSectionIndex(pcone, z - kTolerance(T)) - indexHigh = getSectionIndex(pcone, z + kTolerance(T)) - indexLow < 0 && indexHigh < 0 && return kOutside - indexLow < 0 && indexHigh == 1 && return inside(pcone.sections[1], point-Vector3{T}(0,0,pcone.zᵢ[1])) - indexHigh < 0 && indexLow == N && return inside(pcone.sections[N], point-Vector3{T}(0,0,pcone.zᵢ[N])) - if indexLow == indexHigh - return inside(pcone.sections[indexLow], point-Vector3{T}(0,0,pcone.zᵢ[indexLow])) - else - insideLow = inside(pcone.sections[indexLow], point-Vector3{T}(0,0,pcone.zᵢ[indexLow])) - insideHigh = inside(pcone.sections[indexHigh], point-Vector3{T}(0,0,pcone.zᵢ[indexHigh])) - insideLow == kSurface && insideHigh == kSurface && return kInside - insideLow == kOutside && insideHigh == kOutside && return kOutside - return kSurface - end -end Base.contains(pcone::Polycone{T,N}, p::Point3{T}) where {T,N} = inside(pcone, p) == kInside function distanceToIn(pcone::Polycone{T,N}, point::Point3{T}, dir::Vector3{T})::T where {T,N} diff --git a/src/Trap.jl b/src/Trap.jl index accd72f..6c6a512 100644 --- a/src/Trap.jl +++ b/src/Trap.jl @@ -75,20 +75,6 @@ function extent(trap::Trap{T})::Tuple{Point3{T},Point3{T}} where T<:AbstractFloa ( -p, p ) end -function inside(trap::Trap{T}, point::Point3{T})::Int64 where T<:AbstractFloat - z = point[3] - planes = trap.planes - # inside z? - outside = abs(z) > trap.z + kTolerance(T)/2 - cinside = abs(z) < trap.z - kTolerance(T)/2 - - for i in 1:4 - s = safety(planes[i], point) # positive if the point is on the outside halfspace, negative otherwise. - outside |= s > kTolerance(T)/2 - cinside &= s < -kTolerance(T)/2 - end - cinside ? kInside : outside ? kOutside : kSurface -end function safetyToOut(trap::Trap{T}, point::Point3{T}) where T<:AbstractFloat (; z, planes) = trap diff --git a/src/Trd.jl b/src/Trd.jl index 655019a..07d9f94 100644 --- a/src/Trd.jl +++ b/src/Trd.jl @@ -53,26 +53,6 @@ end cross(px, py, vx, vy) = vx * py - vy * px -function inside(trd::Trd{T}, point::Point3{T})::Int64 where T<:AbstractFloat - x, y, z = point - - # inside z? - outside = abs(z) > (trd.z + kTolerance(T)/2) - inside = abs(z) < (trd.z - kTolerance(T)/2) - - # inside x? - c = cross(abs(x)-trd.x1, z+trd.z, trd.x2-trd.x1, 2*trd.z) - outside |= (c < -kTolerance(T)/2) - inside &= (c > kTolerance(T)/2) - - # inside y? - c = cross(abs(y)-trd.y1, z+trd.z, trd.y2-trd.y1, 2*trd.z) - outside |= (c < -kTolerance(T)/2) - inside &= (c > kTolerance(T)/2) - - # return - inside ? kInside : outside ? kOutside : kSurface -end function normal(trd::Trd{T}, point::Point3{T}) where T<:AbstractFloat x, y, z = point diff --git a/src/Tube.jl b/src/Tube.jl index 140f36b..0fd6b20 100644 --- a/src/Tube.jl +++ b/src/Tube.jl @@ -115,36 +115,6 @@ function extent(tub::Tube{T})::Tuple{Point3{T},Point3{T}} where T<:AbstractFloat return (Point3{T}(aMin), Point3{T}(aMax)) end -function inside(tub::Tube{T}, point::Point3{T})::Int64 where T<:AbstractFloat - x, y, z = point - - # Check Z - outside = abs(z) > tub.z + kTolerance(T)/2 - outside && return kOutside - cinside = abs(z) < tub.z - kTolerance(T)/2 - - # Check on RMax - r2 = x * x + y * y - outside |= r2 > tub.rmax2 + kTolerance(T) * tub.rmax - outside && return kOutside - cinside &= r2 < tub.rmax2 - kTolerance(T) * tub.rmax - - # Check on RMin - if tub.rmin > 0. - outside |= r2 <= tub.rmin2 - kTolerance(T) * tub.rmin - outside && return kOutside - cinside &= r2 > tub.rmin2 + kTolerance(T) * tub.rmin - end - - # Check on Phi - if tub.Δϕ < 2π - outside |= isOutside(tub.ϕWedge, x, y) - cinside &= isInside(tub.ϕWedge, x, y) - end - - return outside ? kOutside : cinside ? kInside : kSurface -end - function normal(tub::Tube{T}, point::Point3{T}) where T<:AbstractFloat x, y, z = point norm = zeros(T,3) diff --git a/src/Volume.jl b/src/Volume.jl index 1813153..8345579 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -205,13 +205,14 @@ function extent(agg::Aggregate{T})::Tuple{Point3{T},Point3{T}} where T<:Abstract max((extent(pvol.volume.shape)[2] * pvol.transformation for pvol in agg.pvolumes)...)) end -function inside(agg::Aggregate{T}, point::Point3{T})::Int64 where T<:AbstractFloat - for pvol in agg.pvolumes - inout = inside(pvol.volume.shape, pvol.transformation * point) - (inout == kInside || inout == kSurface) && return inout - end - return kOutside -end +# FIXME never used +# function inside(agg::Aggregate{T}, point::Point3{T})::Int64 where T<:AbstractFloat +# for pvol in agg.pvolumes +# inout = inside(pvol.volume.shape, pvol.transformation * point) +# (inout == kInside || inout == kSurface) && return inout +# end +# return kOutside +# end function safetyToOut(agg::Aggregate{T}, point::Point3{T}) where T<:AbstractFloat end From 05e5a1f1030ecd92fccd629acae4299fdd1a94f6 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 17:13:31 -0400 Subject: [PATCH 04/16] less useless crap --- src/Distance.jl | 40 ++++++++++++++++++++-------------------- src/Inside.jl | 1 - 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/src/Distance.jl b/src/Distance.jl index 4c5842c..8db5836 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -1,46 +1,46 @@ function distanceToIn(shape, point, dir) if shape isa Trap - distanceToIn_trap(shape::Trap, point, dir) + distanceToIn_trap(shape, point, dir) elseif shape isa Trd - distanceToIn_trd(shape::Trd, point, dir) + distanceToIn_trd(shape, point, dir) elseif shape isa Cone - distanceToIn_cone(shape::Cone, point, dir) + distanceToIn_cone(shape, point, dir) elseif shape isa Box - distanceToIn_box(shape::Box, point, dir) + distanceToIn_box(shape, point, dir) elseif shape isa Tube - distanceToIn_tube(shape::Tube, point, dir) + distanceToIn_tube(shape, point, dir) elseif shape isa Polycone - distanceToIn_polycone(shape::Polycone, point, dir) + distanceToIn_polycone(shape, point, dir) elseif shape isa CutTube - distanceToIn_cuttube(shape::CutTube, point, dir) + distanceToIn_cuttube(shape, point, dir) elseif shape isa BooleanUnion - distanceToIn_booleanunion(shape::BooleanUnion, point, dir) + distanceToIn_booleanunion(shape, point, dir) elseif shape isa BooleanSubtraction - distanceToIn_booleansubtraction(shape::BooleanSubtraction, point, dir) + distanceToIn_booleansubtraction(shape, point, dir) elseif shape isa BooleanIntersection - distanceToIn_booleanintersection(shape::BooleanIntersection, point, dir) + distanceToIn_booleanintersection(shape, point, dir) end end function distanceToOut(shape, point, dir) if shape isa Trap - distanceToOut_trap(shape::Trap, point, dir) + distanceToOut_trap(shape, point, dir) elseif shape isa Trd - distanceToOut_trd(shape::Trd, point, dir) + distanceToOut_trd(shape, point, dir) elseif shape isa Cone - distanceToOut_cone(shape::Cone, point, dir) + distanceToOut_cone(shape, point, dir) elseif shape isa Box - distanceToOut_box(shape::Box, point, dir) + distanceToOut_box(shape, point, dir) elseif shape isa Tube - distanceToOut_tube(shape::Tube, point, dir) + distanceToOut_tube(shape, point, dir) elseif shape isa Polycone - distanceToOut_polycone(shape::Polycone, point, dir) + distanceToOut_polycone(shape, point, dir) elseif shape isa CutTube - distanceToOut_cuttube(shape::CutTube, point, dir) + distanceToOut_cuttube(shape, point, dir) elseif shape isa BooleanUnion - distanceToOut_booleanunion(shape::BooleanUnion, point, dir) + distanceToOut_booleanunion(shape, point, dir) elseif shape isa BooleanSubtraction - distanceToOut_booleansubtraction(shape::BooleanSubtraction, point, dir) + distanceToOut_booleansubtraction(shape, point, dir) elseif shape isa BooleanIntersection - distanceToOut_booleanintersection(shape::BooleanIntersection, point, dir) + distanceToOut_booleanintersection(shape, point, dir) end end diff --git a/src/Inside.jl b/src/Inside.jl index 11ca99d..200dcfe 100644 --- a/src/Inside.jl +++ b/src/Inside.jl @@ -1,5 +1,4 @@ function inside(shape, point) - @nospecialize shape if shape isa Trap inside_trap(shape, point) elseif shape isa Trd From 4399af80286b4545551f7f47ed544bc451e1c4ad Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 17:39:11 -0400 Subject: [PATCH 05/16] good case --- src/Plane.jl | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/Plane.jl b/src/Plane.jl index a763118..bb630ab 100644 --- a/src/Plane.jl +++ b/src/Plane.jl @@ -28,26 +28,6 @@ function safety(plane::Plane{T}, point::Point3{T})::T where T<:AbstractFloat return dot(plane.normal, point) + plane.distance end -function distanceToIn(plane::Plane{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat - # The function returns a negative distance for points already inside or - # direction going outwards (along the normal) - d = T(-Inf) - ndd = dot(normalize(dir), plane.normal) - saf = safety(plane, point) - ndd < 0 && saf > -kTolerance(T) && (d = -saf/ndd) - return d -end - -function distanceToOut(plane::Plane{T}, point::Point3{T}, dir::Vector3{T})::T where T<:AbstractFloat - # The function returns infinity if the plane is not hit from inside, negative - # if the point is outside - d = T(Inf) - ndd = dot(normalize(dir), plane.normal) - saf = safety(plane, point) - saf > kTolerance(T) && (d = -T(Inf)) - ndd > 0 && saf < kTolerance(T) && (d = -saf/ndd) - return d -end function zLimit(plane::Plane{T}, r::T, ϕ::T)::T where T<:AbstractFloat (; normal, distance) = plane From eafe81e0c735a4c6c7daa0d7a3937e074638613e Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 17:39:30 -0400 Subject: [PATCH 06/16] bad case --- src/Distance.jl | 2 ++ src/DistanceIn.jl | 3 +++ src/Volume.jl | 3 --- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Distance.jl b/src/Distance.jl index 8db5836..a88de2a 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -19,6 +19,8 @@ function distanceToIn(shape, point, dir) distanceToIn_booleansubtraction(shape, point, dir) elseif shape isa BooleanIntersection distanceToIn_booleanintersection(shape, point, dir) + elseif shape isa PlacedVolume + distanceToIn_placedvolume(shape, point, dir) end end function distanceToOut(shape, point, dir) diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index a1c8c94..1bc5268 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -1,3 +1,6 @@ +function distanceToIn_placedvolume(pvol::PlacedVolume{T}, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat + distanceToIn(pvol.volume.shape, pvol.transformation * p, pvol.transformation * d) +end ## Boolean function distanceToIn_booleanunion(shape::BooleanUnion{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR} (; left, right, transformation) = shape diff --git a/src/Volume.jl b/src/Volume.jl index 8345579..678bda7 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -64,9 +64,6 @@ end function contains(vol::Volume{T}, p::Point3{T})::Bool where T<:AbstractFloat inside(vol.shape, p) == kInside end -function distanceToIn(pvol::PlacedVolume{T}, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat - distanceToIn(pvol.volume.shape, pvol.transformation * p, pvol.transformation * d) -end function placeDaughter!(volume::Volume{T}, placement::Transformation3D{T}, subvol::Volume{T}) where T<:AbstractFloat From 01565e14ada0753167ac8ab596b7cab1d1932dd4 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Thu, 22 Sep 2022 18:20:41 -0400 Subject: [PATCH 07/16] recover again --- src/Distance.jl | 3 ++- src/DistanceIn.jl | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Distance.jl b/src/Distance.jl index a88de2a..146ec50 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -20,7 +20,8 @@ function distanceToIn(shape, point, dir) elseif shape isa BooleanIntersection distanceToIn_booleanintersection(shape, point, dir) elseif shape isa PlacedVolume - distanceToIn_placedvolume(shape, point, dir) + xf = shape.transformation + distanceToIn_placedvolume(shape.volume.shape, xf*point, xf*dir) end end function distanceToOut(shape, point, dir) diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index 1bc5268..4401214 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -1,5 +1,5 @@ -function distanceToIn_placedvolume(pvol::PlacedVolume{T}, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat - distanceToIn(pvol.volume.shape, pvol.transformation * p, pvol.transformation * d) +function distanceToIn_placedvolume(pvol, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat + distanceToIn(pvol, p, d) end ## Boolean function distanceToIn_booleanunion(shape::BooleanUnion{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR} From a90b0a1b9cd8ba8ee129365db1623f9c10726980 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 14:49:51 -0400 Subject: [PATCH 08/16] Unityper does nothing --- Project.toml | 1 + examples/XRay.jl | 2 +- src/Boolean.jl | 92 +++++++++++----------- src/Distance.jl | 23 ++---- src/DistanceIn.jl | 192 +++++++++++++++++++++++---------------------- src/DistanceOut.jl | 118 ++++++++++++++-------------- src/GDML.jl | 4 +- src/Geom4hep.jl | 1 + src/Navigators.jl | 5 +- src/Volume.jl | 4 +- 10 files changed, 218 insertions(+), 224 deletions(-) diff --git a/Project.toml b/Project.toml index 66ba0bb..d1964ae 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnionArrays = "d6dd79e4-993b-11e9-1366-0de1c9fe1122" +Unityper = "a7c27f48-0311-42f6-a7f8-2c11e75eb415" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/XRay.jl b/examples/XRay.jl index 7fd574f..e72db49 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -8,7 +8,7 @@ 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 - world = getWorld(vol) + world = vol # Setup plane of points and results lower, upper = extent(world.shape) ix, iy, iz = idx[view, :] diff --git a/src/Boolean.jl b/src/Boolean.jl index 538903e..75f61f5 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -1,31 +1,41 @@ -#---Boolean Types---------------------------------------------------------------------------------- -struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} - left::SL # the mother (or left) volume A in unplaced form - right::SR # (or right) volume B in placed form, acting on A with a boolean operation - transformation::Transformation3D{T} # placement of "right" with respect of "left" -end -struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} - left::SL # the mother (or left) volume A in unplaced form - right::SR # (or right) volume B in placed form, acting on A with a boolean operation - transformation::Transformation3D{T} # placement of "right" with respect of "left" -end -struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} - left::SL # the mother (or left) volume A in unplaced form - right::SR # (or right) volume B in placed form, acting on A with a boolean operation - transformation::Transformation3D{T} # placement of "right" with respect of "left" +@compactify begin + @abstract struct AbstractBoolean{T, SL, SR} <: AbstractShape{T} + left::SL = NoShape() + right::SR = NoShape() + transformation::Transformation3D{T} = Transformation3D{T}(0, 0, 0) + end + struct BooleanUnion{T, SL, SR} <: AbstractBoolean{T, SL, SR} end + struct BooleanSubtraction{T, SL, SR} <: AbstractBoolean{T, SL, SR} end + struct BooleanIntersection{T, SL, SR} <: AbstractBoolean{T, SL, SR} end end - #---Constructor------------------------------------------------------------------------------------ function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanUnion{T,typeof(left),typeof(right)}(left, right, place) + BooleanUnion{T,typeof(left),typeof(right)}(; left, right, transformation = place) end function BooleanSubtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanSubtraction{T,typeof(left),typeof(right)}(left, right, place) + BooleanSubtraction{T,typeof(left),typeof(right)}(; left, right, transformation = place) end function BooleanIntersection(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanIntersection{T,typeof(left),typeof(right)}(left, right, place) + BooleanIntersection{T,typeof(left),typeof(right)}(; left, right, transformation = place) end +#---Boolean Types---------------------------------------------------------------------------------- +# struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +# left::SL # the mother (or left) volume A in unplaced form +# right::SR # (or right) volume B in placed form, acting on A with a boolean operation +# transformation::Transformation3D{T} # placement of "right" with respect of "left" +# end +# struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +# left::SL # the mother (or left) volume A in unplaced form +# right::SR # (or right) volume B in placed form, acting on A with a boolean operation +# transformation::Transformation3D{T} # placement of "right" with respect of "left" +# end +# struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +# left::SL # the mother (or left) volume A in unplaced form +# right::SR # (or right) volume B in placed form, acting on A with a boolean operation +# transformation::Transformation3D{T} # placement of "right" with respect of "left" +# end + #---Utilities--------------------------------------------------------------------------------------- #---Printing and Plotting--------------------------------------------------------------------------- @@ -63,39 +73,27 @@ function GeometryBasics.mesh(shape::BooleanIntersection{T, SL, SR}) where {T,SL, end #---Basic functions--------------------------------------------------------------------------------- -function extent(shape::BooleanUnion{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} - (; left, right, transformation) = shape - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (min.(minLeft, minRight), max.(maxLeft, maxRight)) -end -function extent(shape::BooleanSubtraction{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} - extent(shape.left) -end -function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} +function extent(shape::AbstractBoolean) (; left, right, transformation) = shape - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (max.(minLeft, minRight), min.(maxLeft, maxRight)) + @compactified shape::AbstractBoolean begin + BooleanUnion => begin + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (min.(minLeft, minRight), max.(maxLeft, maxRight)) + end + BooleanIntersection => begin + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (max.(minLeft, minRight), min.(maxLeft, maxRight)) + end + BooleanSubtraction => extent(shape.left) + end end -function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end - -function safetyToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} +function safetyToOut(shape::AbstractBoolean) return 0 end -function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} +function safetyToIn(shape::AbstractBoolean) return 0 end diff --git a/src/Distance.jl b/src/Distance.jl index 146ec50..37e11f0 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -1,4 +1,4 @@ -function distanceToIn(shape, point, dir) +function distanceToIn(shape, point::Point3{T}, dir)::T where T if shape isa Trap distanceToIn_trap(shape, point, dir) elseif shape isa Trd @@ -13,18 +13,11 @@ function distanceToIn(shape, point, dir) distanceToIn_polycone(shape, point, dir) elseif shape isa CutTube distanceToIn_cuttube(shape, point, dir) - elseif shape isa BooleanUnion - distanceToIn_booleanunion(shape, point, dir) - elseif shape isa BooleanSubtraction - distanceToIn_booleansubtraction(shape, point, dir) - elseif shape isa BooleanIntersection - distanceToIn_booleanintersection(shape, point, dir) - elseif shape isa PlacedVolume - xf = shape.transformation - distanceToIn_placedvolume(shape.volume.shape, xf*point, xf*dir) + elseif shape isa AbstractBoolean + distanceToIn_boolean(shape, point, dir) end end -function distanceToOut(shape, point, dir) +function distanceToOut(shape, point::Point3{T}, dir)::T where T if shape isa Trap distanceToOut_trap(shape, point, dir) elseif shape isa Trd @@ -39,11 +32,7 @@ function distanceToOut(shape, point, dir) distanceToOut_polycone(shape, point, dir) elseif shape isa CutTube distanceToOut_cuttube(shape, point, dir) - elseif shape isa BooleanUnion - distanceToOut_booleanunion(shape, point, dir) - elseif shape isa BooleanSubtraction - distanceToOut_booleansubtraction(shape, point, dir) - elseif shape isa BooleanIntersection - distanceToOut_booleanintersection(shape, point, dir) + elseif shape isa AbstractBoolean + distanceToOut_boolean(shape, point, dir) end end diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index 4401214..0e652f6 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -1,102 +1,104 @@ -function distanceToIn_placedvolume(pvol, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat - distanceToIn(pvol, p, d) -end ## Boolean -function distanceToIn_booleanunion(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_booleanintersection(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 +function distanceToIn_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} + @compactified shape::AbstractBoolean begin + BooleanUnion => begin + (; left, right, transformation) = shape + distA = distanceToIn(left, point, dir) + distB = distanceToIn(right, transformation * point, transformation * dir) + return min(distA, distB) end - end - return dist -end -function distanceToIn_booleansubtraction(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 + BooleanIntersection => begin + (; 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 - - # 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 + BooleanSubtraction => begin + (; 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 + 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/DistanceOut.jl b/src/DistanceOut.jl index 6ab14cf..6205b60 100644 --- a/src/DistanceOut.jl +++ b/src/DistanceOut.jl @@ -1,70 +1,74 @@ #Boolean -function distanceToOut_booleanunion(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) +function distanceToOut_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} + @compactified shape::AbstractBoolean begin + BooleanUnion => begin + (; 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 + 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 - break + return T(-1) end - end - return dist - kPushTolerance(T) - else - return T(-1) + end end - end -end -function distanceToOut_booleanintersection(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_booleansubtraction(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) + BooleanIntersection => begin + (; left, right, transformation) = shape + distA = distanceToOut(left, point, dir) + distB = distanceToOut(right, transformation * point, transformation * dir) + return min(distA, distB) + end + BooleanSubtraction => begin + (; left, right, transformation) = shape + distA = distanceToOut(left, point, dir) + distB = distanceToIn(right, transformation * point, transformation * dir) + return min(distA, distB) + end + end end #Plane diff --git a/src/GDML.jl b/src/GDML.jl index 73afee1..0f36d69 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -225,7 +225,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac volname = attrs["name"] shape = nothing material = nothing - pvols = PlacedVolume{T}[] + pvols = [] for cc in child_nodes(e) #--- loop over ref to solid, material and daughters if is_elementnode(cc) aa = attributes_dict(XMLElement(cc)) @@ -244,7 +244,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - push!(pvols, PlacedVolume{T}(-1, transformation, daughter)) + push!(pvols, PlacedVolume(-1, transformation, daughter)) end end end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 41984e6..6f2e4cc 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -27,6 +27,7 @@ include("Tube.jl") include("Cone.jl") include("Polycone.jl") include("CutTube.jl") +using Unityper include("Boolean.jl") include("Trap.jl") include("Materials.jl") diff --git a/src/Navigators.jl b/src/Navigators.jl index b90ac38..ab1e1c3 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -158,10 +158,11 @@ function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir:: step = step_limit candidate = 0 #---Linear loop over the daughters------------------------------------------------- - for i in intersectedDaughters(nav, volume, point, dir) + @inbounds for i in intersectedDaughters(nav, volume, point, dir) pvol = volume.daughters[i] #---Assuming that it is not yet inside the daughter (otherwise it returns -1.) - dist = distanceToIn(pvol, point, dir) + xf = pvol.transformation + dist = distanceToIn(pvol.volume.shape, xf*point, xf*dir) dist < 0. && return (zero(T), pvol.idx ) if dist > 0. && dist != Inf && dist < step step = dist diff --git a/src/Volume.jl b/src/Volume.jl index 678bda7..00ad235 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -13,9 +13,7 @@ const Shape{T} = Union{NoShape{T}, Cone{T}, Polycone{T}, CutTube{T}, - BooleanUnion{T}, - BooleanIntersection{T}, - BooleanSubtraction{T}} where T<:AbstractFloat + AbstractBoolean{T}} where T<:AbstractFloat #---Volume------------------------------------------------------------------------- struct VolumeP{T<:AbstractFloat,PV} From 6cb6ee22021bb00e38fcbc4397f45b1223736ada Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 15:08:50 -0400 Subject: [PATCH 09/16] Revert "Unityper does nothing" This reverts commit a90b0a1b9cd8ba8ee129365db1623f9c10726980. --- Project.toml | 1 - examples/XRay.jl | 2 +- src/Boolean.jl | 92 +++++++++++----------- src/Distance.jl | 23 ++++-- src/DistanceIn.jl | 192 ++++++++++++++++++++++----------------------- src/DistanceOut.jl | 118 ++++++++++++++-------------- src/GDML.jl | 4 +- src/Geom4hep.jl | 1 - src/Navigators.jl | 5 +- src/Volume.jl | 4 +- 10 files changed, 224 insertions(+), 218 deletions(-) diff --git a/Project.toml b/Project.toml index d1964ae..66ba0bb 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,6 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnionArrays = "d6dd79e4-993b-11e9-1366-0de1c9fe1122" -Unityper = "a7c27f48-0311-42f6-a7f8-2c11e75eb415" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/XRay.jl b/examples/XRay.jl index e72db49..7fd574f 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -8,7 +8,7 @@ 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 - world = vol + world = getWorld(vol) # Setup plane of points and results lower, upper = extent(world.shape) ix, iy, iz = idx[view, :] diff --git a/src/Boolean.jl b/src/Boolean.jl index 75f61f5..538903e 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -1,41 +1,31 @@ -@compactify begin - @abstract struct AbstractBoolean{T, SL, SR} <: AbstractShape{T} - left::SL = NoShape() - right::SR = NoShape() - transformation::Transformation3D{T} = Transformation3D{T}(0, 0, 0) - end - struct BooleanUnion{T, SL, SR} <: AbstractBoolean{T, SL, SR} end - struct BooleanSubtraction{T, SL, SR} <: AbstractBoolean{T, SL, SR} end - struct BooleanIntersection{T, SL, SR} <: AbstractBoolean{T, SL, SR} end +#---Boolean Types---------------------------------------------------------------------------------- +struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} + left::SL # the mother (or left) volume A in unplaced form + right::SR # (or right) volume B in placed form, acting on A with a boolean operation + transformation::Transformation3D{T} # placement of "right" with respect of "left" +end +struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} + left::SL # the mother (or left) volume A in unplaced form + right::SR # (or right) volume B in placed form, acting on A with a boolean operation + transformation::Transformation3D{T} # placement of "right" with respect of "left" +end +struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} + left::SL # the mother (or left) volume A in unplaced form + right::SR # (or right) volume B in placed form, acting on A with a boolean operation + transformation::Transformation3D{T} # placement of "right" with respect of "left" end + #---Constructor------------------------------------------------------------------------------------ function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanUnion{T,typeof(left),typeof(right)}(; left, right, transformation = place) + BooleanUnion{T,typeof(left),typeof(right)}(left, right, place) end function BooleanSubtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanSubtraction{T,typeof(left),typeof(right)}(; left, right, transformation = place) + BooleanSubtraction{T,typeof(left),typeof(right)}(left, right, place) end function BooleanIntersection(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanIntersection{T,typeof(left),typeof(right)}(; left, right, transformation = place) + BooleanIntersection{T,typeof(left),typeof(right)}(left, right, place) end -#---Boolean Types---------------------------------------------------------------------------------- -# struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} -# left::SL # the mother (or left) volume A in unplaced form -# right::SR # (or right) volume B in placed form, acting on A with a boolean operation -# transformation::Transformation3D{T} # placement of "right" with respect of "left" -# end -# struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} -# left::SL # the mother (or left) volume A in unplaced form -# right::SR # (or right) volume B in placed form, acting on A with a boolean operation -# transformation::Transformation3D{T} # placement of "right" with respect of "left" -# end -# struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} -# left::SL # the mother (or left) volume A in unplaced form -# right::SR # (or right) volume B in placed form, acting on A with a boolean operation -# transformation::Transformation3D{T} # placement of "right" with respect of "left" -# end - #---Utilities--------------------------------------------------------------------------------------- #---Printing and Plotting--------------------------------------------------------------------------- @@ -73,27 +63,39 @@ function GeometryBasics.mesh(shape::BooleanIntersection{T, SL, SR}) where {T,SL, end #---Basic functions--------------------------------------------------------------------------------- -function extent(shape::AbstractBoolean) +function extent(shape::BooleanUnion{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} + (; left, right, transformation) = shape + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (min.(minLeft, minRight), max.(maxLeft, maxRight)) +end +function extent(shape::BooleanSubtraction{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} + extent(shape.left) +end +function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} (; left, right, transformation) = shape - @compactified shape::AbstractBoolean begin - BooleanUnion => begin - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (min.(minLeft, minRight), max.(maxLeft, maxRight)) - end - BooleanIntersection => begin - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (max.(minLeft, minRight), min.(maxLeft, maxRight)) - end - BooleanSubtraction => extent(shape.left) - end + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (max.(minLeft, minRight), min.(maxLeft, maxRight)) end -function safetyToOut(shape::AbstractBoolean) +function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end + +function safetyToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 end -function safetyToIn(shape::AbstractBoolean) +function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 end diff --git a/src/Distance.jl b/src/Distance.jl index 37e11f0..146ec50 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -1,4 +1,4 @@ -function distanceToIn(shape, point::Point3{T}, dir)::T where T +function distanceToIn(shape, point, dir) if shape isa Trap distanceToIn_trap(shape, point, dir) elseif shape isa Trd @@ -13,11 +13,18 @@ function distanceToIn(shape, point::Point3{T}, dir)::T where T distanceToIn_polycone(shape, point, dir) elseif shape isa CutTube distanceToIn_cuttube(shape, point, dir) - elseif shape isa AbstractBoolean - distanceToIn_boolean(shape, point, dir) + elseif shape isa BooleanUnion + distanceToIn_booleanunion(shape, point, dir) + elseif shape isa BooleanSubtraction + distanceToIn_booleansubtraction(shape, point, dir) + elseif shape isa BooleanIntersection + distanceToIn_booleanintersection(shape, point, dir) + elseif shape isa PlacedVolume + xf = shape.transformation + distanceToIn_placedvolume(shape.volume.shape, xf*point, xf*dir) end end -function distanceToOut(shape, point::Point3{T}, dir)::T where T +function distanceToOut(shape, point, dir) if shape isa Trap distanceToOut_trap(shape, point, dir) elseif shape isa Trd @@ -32,7 +39,11 @@ function distanceToOut(shape, point::Point3{T}, dir)::T where T distanceToOut_polycone(shape, point, dir) elseif shape isa CutTube distanceToOut_cuttube(shape, point, dir) - elseif shape isa AbstractBoolean - distanceToOut_boolean(shape, point, dir) + elseif shape isa BooleanUnion + distanceToOut_booleanunion(shape, point, dir) + elseif shape isa BooleanSubtraction + distanceToOut_booleansubtraction(shape, point, dir) + elseif shape isa BooleanIntersection + distanceToOut_booleanintersection(shape, point, dir) end end diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index 0e652f6..4401214 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -1,104 +1,102 @@ +function distanceToIn_placedvolume(pvol, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat + distanceToIn(pvol, p, d) +end ## Boolean -function distanceToIn_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} - @compactified shape::AbstractBoolean begin - BooleanUnion => begin - (; left, right, transformation) = shape - distA = distanceToIn(left, point, dir) - distB = distanceToIn(right, transformation * point, transformation * dir) - return min(distA, distB) +function distanceToIn_booleanunion(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_booleanintersection(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 - BooleanIntersection => begin - (; 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 + return dist +end +function distanceToIn_booleansubtraction(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 - BooleanSubtraction => begin - (; 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 - return dist + + # 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/DistanceOut.jl b/src/DistanceOut.jl index 6205b60..6ab14cf 100644 --- a/src/DistanceOut.jl +++ b/src/DistanceOut.jl @@ -1,74 +1,70 @@ #Boolean -function distanceToOut_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} - @compactified shape::AbstractBoolean begin - BooleanUnion => begin - (; left, right, transformation) = shape - dist = T(0) - positionA = inside(left, point) - if positionA != kOutside # point inside A - while(true) - distA = distanceToOut(left, point, dir) +function distanceToOut_booleanunion(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 # 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 + if inside(right, lpoint) == kOutside 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) + break end - end - end - BooleanIntersection => begin - (; left, right, transformation) = shape - distA = distanceToOut(left, point, dir) - distB = distanceToOut(right, transformation * point, transformation * dir) - return min(distA, distB) - end - BooleanSubtraction => begin - (; left, right, transformation) = shape - distA = distanceToOut(left, point, dir) - distB = distanceToIn(right, transformation * point, transformation * dir) - return min(distA, distB) + end + return dist - kPushTolerance(T) + else + return T(-1) end - end + end +end +function distanceToOut_booleanintersection(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_booleansubtraction(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 #Plane diff --git a/src/GDML.jl b/src/GDML.jl index 0f36d69..73afee1 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -225,7 +225,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac volname = attrs["name"] shape = nothing material = nothing - pvols = [] + pvols = PlacedVolume{T}[] for cc in child_nodes(e) #--- loop over ref to solid, material and daughters if is_elementnode(cc) aa = attributes_dict(XMLElement(cc)) @@ -244,7 +244,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - push!(pvols, PlacedVolume(-1, transformation, daughter)) + push!(pvols, PlacedVolume{T}(-1, transformation, daughter)) end end end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 6f2e4cc..41984e6 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -27,7 +27,6 @@ include("Tube.jl") include("Cone.jl") include("Polycone.jl") include("CutTube.jl") -using Unityper include("Boolean.jl") include("Trap.jl") include("Materials.jl") diff --git a/src/Navigators.jl b/src/Navigators.jl index ab1e1c3..b90ac38 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -158,11 +158,10 @@ function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir:: step = step_limit candidate = 0 #---Linear loop over the daughters------------------------------------------------- - @inbounds for i in intersectedDaughters(nav, volume, point, dir) + for i in intersectedDaughters(nav, volume, point, dir) pvol = volume.daughters[i] #---Assuming that it is not yet inside the daughter (otherwise it returns -1.) - xf = pvol.transformation - dist = distanceToIn(pvol.volume.shape, xf*point, xf*dir) + dist = distanceToIn(pvol, point, dir) dist < 0. && return (zero(T), pvol.idx ) if dist > 0. && dist != Inf && dist < step step = dist diff --git a/src/Volume.jl b/src/Volume.jl index 00ad235..678bda7 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -13,7 +13,9 @@ const Shape{T} = Union{NoShape{T}, Cone{T}, Polycone{T}, CutTube{T}, - AbstractBoolean{T}} where T<:AbstractFloat + BooleanUnion{T}, + BooleanIntersection{T}, + BooleanSubtraction{T}} where T<:AbstractFloat #---Volume------------------------------------------------------------------------- struct VolumeP{T<:AbstractFloat,PV} From cc820071aa1f5d0b782c7b0f1d4f6568c8756263 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 15:57:06 -0400 Subject: [PATCH 10/16] try this --- Project.toml | 1 + examples/XRay.jl | 2 +- src/Boolean.jl | 45 ++++++--------------------------------------- src/GDML.jl | 4 ++-- src/Geom4hep.jl | 2 +- src/Navigators.jl | 26 +++++++++++++------------- src/Volume.jl | 40 ++++++++++++++++++++-------------------- 7 files changed, 44 insertions(+), 76 deletions(-) diff --git a/Project.toml b/Project.toml index 66ba0bb..5dce908 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" diff --git a/examples/XRay.jl b/examples/XRay.jl index 7fd574f..6b0c1cb 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -20,7 +20,7 @@ function generateXRay(nav::AbstractNavigator, vol::Volume{T}, npoints::Number, v z = lower[iz] + kTolerance(T) result = zeros(T, nx,ny) - state = NavigatorState{T}(world, nav) + state = NavigatorState(world, nav) _dir = (0, 0, 1) dir = Vector3{T}(_dir[xi], _dir[yi], _dir[zi]) diff --git a/src/Boolean.jl b/src/Boolean.jl index 538903e..07a00b7 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -1,5 +1,7 @@ #---Boolean Types---------------------------------------------------------------------------------- -struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +abstract type AbstractBoolean{T, SL, SR} <: AbstractShape{T} end + +struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractBoolean{T, SL, SR} left::SL # the mother (or left) volume A in unplaced form right::SR # (or right) volume B in placed form, acting on A with a boolean operation transformation::Transformation3D{T} # placement of "right" with respect of "left" @@ -15,18 +17,6 @@ struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShap transformation::Transformation3D{T} # placement of "right" with respect of "left" end -#---Constructor------------------------------------------------------------------------------------ -function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanUnion{T,typeof(left),typeof(right)}(left, right, place) -end -function BooleanSubtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanSubtraction{T,typeof(left),typeof(right)}(left, right, place) -end -function BooleanIntersection(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanIntersection{T,typeof(left),typeof(right)}(left, right, place) -end - -#---Utilities--------------------------------------------------------------------------------------- #---Printing and Plotting--------------------------------------------------------------------------- function Base.show(io::IO, shape::BooleanUnion{T, SL, SR}) where {T,SL,SR} @@ -46,17 +36,7 @@ function Base.show(io::IO, shape::BooleanIntersection{T, SL, SR}) where {T,SL,SR end -function GeometryBasics.mesh(shape::BooleanUnion{T, SL, SR}) where {T,SL,SR} - (; left, right, transformation) = shape - rmesh = mesh(right) - merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))]) -end -function GeometryBasics.mesh(shape::BooleanSubtraction{T, SL, SR}) where {T,SL,SR} - (; left, right, transformation) = shape - rmesh = mesh(right) - merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))]) -end -function GeometryBasics.mesh(shape::BooleanIntersection{T, SL, SR}) where {T,SL,SR} +function GeometryBasics.mesh(shape::AbstractBoolean{T, SL, SR}) where {T,SL,SR} (; left, right, transformation) = shape rmesh = mesh(right) merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))]) @@ -80,22 +60,9 @@ function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T end -function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end - -function safetyToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} +function safetyToOut(shape::AbstractBoolean) return 0 end -function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} +function safetyToIn(shape::AbstractBoolean) return 0 end diff --git a/src/GDML.jl b/src/GDML.jl index 73afee1..52a00e1 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -244,11 +244,11 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - push!(pvols, PlacedVolume{T}(-1, transformation, daughter)) + push!(pvols, PlacedVolume(-1, transformation, daughter)) end end end - volume = elemname == "volume" ? Volume{T}(volname, shape, material) : Assembly{T}(volname) + volume = elemname == "volume" ? Volume(volname, shape, material) : Assembly{T}(volname) for pvol in pvols placeDaughter!(volume, pvol.transformation, pvol.volume) end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 41984e6..0541735 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -16,7 +16,7 @@ export Triangle, Intersection, intersect, distanceToPlane export Tesselation, coordinates, faces, normals, mesh using Requires -using StaticArrays, GeometryBasics, LinearAlgebra, Rotations, AbstractTrees +using StaticArrays, GeometryBasics, LinearAlgebra, Rotations, AbstractTrees, Accessors include("BasicTypes.jl") include("Transformation3D.jl") include("TriangleIntersect.jl") diff --git a/src/Navigators.jl b/src/Navigators.jl index b90ac38..c8fb541 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -49,18 +49,18 @@ function Base.show(io::IO, nav::BVHNavigator{T}) where T end #---NavigatorState---------------------------------------------------------------------------------- -mutable struct NavigatorState{T<:AbstractFloat, NAV<:AbstractNavigator} <: AbstractNavigatorState +struct NavigatorState{T<:AbstractFloat, VT1, VT2, NAV<:AbstractNavigator} <: AbstractNavigatorState navigator::NAV # Navigator to be used - topvol::Volume{T} # Typically the unplaced world - currvol::Volume{T} # the current volume + topvol::VT1 + currvol::VT2 isinworld::Bool # inside world volume volstack::Vector{Int64} # keep the indexes of all daughters up to the current one tolocal::Vector{Transformation3D{T}} # Stack of transformations end #---Constructors------------------------------------------------------------------------------------ -function NavigatorState{T}(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T - x = NavigatorState{T,typeof(nav)}(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) +function NavigatorState(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T + x = NavigatorState(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) sizehint!(x.volstack,16) sizehint!(x.tolocal,16) return x @@ -69,8 +69,8 @@ end @inline function reset!(state::NavigatorState{T}) where T<:AbstractFloat empty!(state.volstack) empty!(state.tolocal) - state.currvol = state.topvol - state.isinworld = false + @set state.currvol = state.topvol + @set state.isinworld = false end @inline function currentVolume(state::NavigatorState{T}) where T<:AbstractFloat @@ -85,17 +85,17 @@ end for idx in state.volstack vol = vol.daughters[idx].volume end - state.currvol = vol + @set state.currvol = vol else - state.currvol = state.topvol - state.isinworld = false + @set state.currvol = state.topvol + @set state.isinworld = false end end @inline function pushIn!(state::NavigatorState{T}, pvol::PlacedVolume{T}) where T<:AbstractFloat push!(state.volstack, pvol.idx) push!(state.tolocal, pvol.transformation) - state.currvol = pvol.volume + @set state.currvol = pvol.volume end #---Basic loops implemented using acceleration structures------------------------------------------- @@ -129,8 +129,8 @@ function locateGlobalPoint!(state::NavigatorState{T,NAV}, point::Point3{T}) wher reset!(state) vol = state.topvol if contains(vol, point) - state.currvol = vol - state.isinworld = true + @set state.currvol = vol + @set state.isinworld = true end # check the daughters in a recursive manner isinside = true diff --git a/src/Volume.jl b/src/Volume.jl index 678bda7..0399239 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -5,39 +5,39 @@ end NoShape{T}() where T = NoShape{T,Nothing}([]) #---Shape-------------------------------------------------------------------------- -const Shape{T} = Union{NoShape{T}, - Box{T}, - Trd{T}, - Trap{T}, - Tube{T}, - Cone{T}, - Polycone{T}, - CutTube{T}, - BooleanUnion{T}, - BooleanIntersection{T}, - BooleanSubtraction{T}} where T<:AbstractFloat +# const Shape{T} = Union{NoShape{T}, +# Box{T}, +# Trd{T}, +# Trap{T}, +# Tube{T}, +# Cone{T}, +# Polycone{T}, +# CutTube{T}, +# BooleanUnion{T}, +# BooleanIntersection{T}, +# BooleanSubtraction{T}} where T<:AbstractFloat #---Volume------------------------------------------------------------------------- -struct VolumeP{T<:AbstractFloat,PV} +struct Volume{T<:AbstractFloat, S, PV} label::String - shape::Shape{T} # Reference to the actual shape - material::Material{T} # Reference to material + shape::S + material::Material{T} daughters::Vector{PV} end #---PlacedVolume------------------------------------------------------------------- -struct PlacedVolume{T<:AbstractFloat} +struct PlacedVolume{T<:AbstractFloat, VT} idx::Int64 transformation::Transformation3D{T} - volume::VolumeP{T,PlacedVolume{T}} + volume::VT end #---Convenient Alias to simplify signatures--------------------------------------- -const Volume{T} = VolumeP{T,PlacedVolume{T}} where T<:AbstractFloat +# const Volume{T} = VolumeP{T,S,PlacedVolume{T}} where {T<:AbstractFloat, S<:AbstractShape{T}} #---Constructor-------------------------------------------------------------------- -function Volume{T}(label::String, shape::Shape{T}, material::Material{T}) where T<:AbstractFloat - Volume{T}(label, shape, material, Vector{PlacedVolume{T}}()) # call the default constructor +function Volume(label, shape, material::Material{T}) where T<:AbstractFloat + Volume(label, shape, material, Vector{PlacedVolume{T}}()) # call the default constructor end #---Utilities---------------------------------------------------------------------- @@ -85,7 +85,7 @@ function getWorld(vol::Volume{T}) where T<:AbstractFloat box = Box{T}((high - low)/2. .+ 0.1) # increase the bounding box tra = Transformation3D{T}(one(RotMatrix3{T}), -(high + low)/2.) mat = Material{T}("vacuum"; density=0) - world = Volume{T}("world", box, mat) + world = Volume("world", box, mat) placeDaughter!(world, tra, vol) return world end From 772f040aea2a8b32725efcc76a5be9678f210809 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 16:27:23 -0400 Subject: [PATCH 11/16] immutable + ref is no-no --- src/Navigators.jl | 44 +++++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/src/Navigators.jl b/src/Navigators.jl index c8fb541..039d077 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -49,18 +49,18 @@ function Base.show(io::IO, nav::BVHNavigator{T}) where T end #---NavigatorState---------------------------------------------------------------------------------- -struct NavigatorState{T<:AbstractFloat, VT1, VT2, NAV<:AbstractNavigator} <: AbstractNavigatorState +struct NavigatorState{T<:AbstractFloat, NAV<:AbstractNavigator} <: AbstractNavigatorState navigator::NAV # Navigator to be used - topvol::VT1 - currvol::VT2 - isinworld::Bool # inside world volume + topvol:: Base.RefValue{Volume} + currvol::Base.RefValue{Volume} + isinworld::Ref{Bool} # inside world volume volstack::Vector{Int64} # keep the indexes of all daughters up to the current one tolocal::Vector{Transformation3D{T}} # Stack of transformations end #---Constructors------------------------------------------------------------------------------------ function NavigatorState(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T - x = NavigatorState(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) + x = NavigatorState(nav, Ref{Volume}(top), Ref{Volume}(top), Ref(false), Vector{Int64}(), Vector{Transformation3D{T}}()) sizehint!(x.volstack,16) sizehint!(x.tolocal,16) return x @@ -69,33 +69,35 @@ end @inline function reset!(state::NavigatorState{T}) where T<:AbstractFloat empty!(state.volstack) empty!(state.tolocal) - @set state.currvol = state.topvol - @set state.isinworld = false + state.currvol[] = state.topvol[] + state.isinworld[] = false + nothing end @inline function currentVolume(state::NavigatorState{T}) where T<:AbstractFloat - return state.currvol + return state.currvol[] end @inline function popOut!(state::NavigatorState{T}) where T<:AbstractFloat - if !isempty(state.volstack) + res = if !isempty(state.volstack) pop!(state.volstack) pop!(state.tolocal) - vol = state.topvol + vol = state.topvol[] for idx in state.volstack vol = vol.daughters[idx].volume end - @set state.currvol = vol + state.currvol[] = vol else - @set state.currvol = state.topvol - @set state.isinworld = false + state.currvol[] = state.topvol[] + state.isinworld[] = false end + nothing end @inline function pushIn!(state::NavigatorState{T}, pvol::PlacedVolume{T}) where T<:AbstractFloat push!(state.volstack, pvol.idx) push!(state.tolocal, pvol.transformation) - @set state.currvol = pvol.volume + state.currvol[] = pvol.volume end #---Basic loops implemented using acceleration structures------------------------------------------- @@ -127,10 +129,10 @@ end function locateGlobalPoint!(state::NavigatorState{T,NAV}, point::Point3{T}) where {T, NAV} nav = state.navigator reset!(state) - vol = state.topvol + vol = state.topvol[] if contains(vol, point) - @set state.currvol = vol - @set state.isinworld = true + state.currvol[] = vol + state.isinworld[] = true end # check the daughters in a recursive manner isinside = true @@ -147,11 +149,11 @@ function locateGlobalPoint!(state::NavigatorState{T,NAV}, point::Point3{T}) wher end end end - return state.isinworld + nothing end @inline function isInVolume(state::NavigatorState{T}) where T<:AbstractFloat - state.isinworld + state.isinworld[] end function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir::Vector3{T}, step_limit::T ) where {T,NAV} @@ -180,7 +182,7 @@ function computeStep!(state::NavigatorState{T,NAV}, gpoint::Point3{T}, gdir::Vec if idx == 0 step = distanceToOut(volume.shape, lpoint, ldir) if step > -kTolerance(T) - popOut!(state) + state = popOut!(state) elseif step < 0 shape = volume.shape volstack = state.volstack @@ -190,7 +192,7 @@ function computeStep!(state::NavigatorState{T,NAV}, gpoint::Point3{T}, gdir::Vec #---We hit a daughter, push it into the stack step += kPushTolerance(T) # to ensure that do not stay in the surface of the daughter pvol = volume.daughters[idx] - pushIn!(state, pvol) + state = pushIn!(state, pvol) end return step end From 0a02f13fa3faf04d706987b960dda2752a508a66 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 16:34:28 -0400 Subject: [PATCH 12/16] Revert "immutable + ref is no-no" This reverts commit 772f040aea2a8b32725efcc76a5be9678f210809. --- src/Navigators.jl | 44 +++++++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/src/Navigators.jl b/src/Navigators.jl index 039d077..c8fb541 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -49,18 +49,18 @@ function Base.show(io::IO, nav::BVHNavigator{T}) where T end #---NavigatorState---------------------------------------------------------------------------------- -struct NavigatorState{T<:AbstractFloat, NAV<:AbstractNavigator} <: AbstractNavigatorState +struct NavigatorState{T<:AbstractFloat, VT1, VT2, NAV<:AbstractNavigator} <: AbstractNavigatorState navigator::NAV # Navigator to be used - topvol:: Base.RefValue{Volume} - currvol::Base.RefValue{Volume} - isinworld::Ref{Bool} # inside world volume + topvol::VT1 + currvol::VT2 + isinworld::Bool # inside world volume volstack::Vector{Int64} # keep the indexes of all daughters up to the current one tolocal::Vector{Transformation3D{T}} # Stack of transformations end #---Constructors------------------------------------------------------------------------------------ function NavigatorState(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T - x = NavigatorState(nav, Ref{Volume}(top), Ref{Volume}(top), Ref(false), Vector{Int64}(), Vector{Transformation3D{T}}()) + x = NavigatorState(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) sizehint!(x.volstack,16) sizehint!(x.tolocal,16) return x @@ -69,35 +69,33 @@ end @inline function reset!(state::NavigatorState{T}) where T<:AbstractFloat empty!(state.volstack) empty!(state.tolocal) - state.currvol[] = state.topvol[] - state.isinworld[] = false - nothing + @set state.currvol = state.topvol + @set state.isinworld = false end @inline function currentVolume(state::NavigatorState{T}) where T<:AbstractFloat - return state.currvol[] + return state.currvol end @inline function popOut!(state::NavigatorState{T}) where T<:AbstractFloat - res = if !isempty(state.volstack) + if !isempty(state.volstack) pop!(state.volstack) pop!(state.tolocal) - vol = state.topvol[] + vol = state.topvol for idx in state.volstack vol = vol.daughters[idx].volume end - state.currvol[] = vol + @set state.currvol = vol else - state.currvol[] = state.topvol[] - state.isinworld[] = false + @set state.currvol = state.topvol + @set state.isinworld = false end - nothing end @inline function pushIn!(state::NavigatorState{T}, pvol::PlacedVolume{T}) where T<:AbstractFloat push!(state.volstack, pvol.idx) push!(state.tolocal, pvol.transformation) - state.currvol[] = pvol.volume + @set state.currvol = pvol.volume end #---Basic loops implemented using acceleration structures------------------------------------------- @@ -129,10 +127,10 @@ end function locateGlobalPoint!(state::NavigatorState{T,NAV}, point::Point3{T}) where {T, NAV} nav = state.navigator reset!(state) - vol = state.topvol[] + vol = state.topvol if contains(vol, point) - state.currvol[] = vol - state.isinworld[] = true + @set state.currvol = vol + @set state.isinworld = true end # check the daughters in a recursive manner isinside = true @@ -149,11 +147,11 @@ function locateGlobalPoint!(state::NavigatorState{T,NAV}, point::Point3{T}) wher end end end - nothing + return state.isinworld end @inline function isInVolume(state::NavigatorState{T}) where T<:AbstractFloat - state.isinworld[] + state.isinworld end function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir::Vector3{T}, step_limit::T ) where {T,NAV} @@ -182,7 +180,7 @@ function computeStep!(state::NavigatorState{T,NAV}, gpoint::Point3{T}, gdir::Vec if idx == 0 step = distanceToOut(volume.shape, lpoint, ldir) if step > -kTolerance(T) - state = popOut!(state) + popOut!(state) elseif step < 0 shape = volume.shape volstack = state.volstack @@ -192,7 +190,7 @@ function computeStep!(state::NavigatorState{T,NAV}, gpoint::Point3{T}, gdir::Vec #---We hit a daughter, push it into the stack step += kPushTolerance(T) # to ensure that do not stay in the surface of the daughter pvol = volume.daughters[idx] - state = pushIn!(state, pvol) + pushIn!(state, pvol) end return step end From f43bdf486a23deb36b6b65a06b1fa72e8b57c2d1 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 16:34:30 -0400 Subject: [PATCH 13/16] Revert "try this" This reverts commit cc820071aa1f5d0b782c7b0f1d4f6568c8756263. --- Project.toml | 1 - examples/XRay.jl | 2 +- src/Boolean.jl | 45 +++++++++++++++++++++++++++++++++++++++------ src/GDML.jl | 4 ++-- src/Geom4hep.jl | 2 +- src/Navigators.jl | 26 +++++++++++++------------- src/Volume.jl | 40 ++++++++++++++++++++-------------------- 7 files changed, 76 insertions(+), 44 deletions(-) diff --git a/Project.toml b/Project.toml index 5dce908..66ba0bb 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "0.1.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" -Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" diff --git a/examples/XRay.jl b/examples/XRay.jl index 6b0c1cb..7fd574f 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -20,7 +20,7 @@ function generateXRay(nav::AbstractNavigator, vol::Volume{T}, npoints::Number, v z = lower[iz] + kTolerance(T) result = zeros(T, nx,ny) - state = NavigatorState(world, nav) + state = NavigatorState{T}(world, nav) _dir = (0, 0, 1) dir = Vector3{T}(_dir[xi], _dir[yi], _dir[zi]) diff --git a/src/Boolean.jl b/src/Boolean.jl index 07a00b7..538903e 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -1,7 +1,5 @@ #---Boolean Types---------------------------------------------------------------------------------- -abstract type AbstractBoolean{T, SL, SR} <: AbstractShape{T} end - -struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractBoolean{T, SL, SR} +struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} left::SL # the mother (or left) volume A in unplaced form right::SR # (or right) volume B in placed form, acting on A with a boolean operation transformation::Transformation3D{T} # placement of "right" with respect of "left" @@ -17,6 +15,18 @@ struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShap transformation::Transformation3D{T} # placement of "right" with respect of "left" end +#---Constructor------------------------------------------------------------------------------------ +function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat + BooleanUnion{T,typeof(left),typeof(right)}(left, right, place) +end +function BooleanSubtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat + BooleanSubtraction{T,typeof(left),typeof(right)}(left, right, place) +end +function BooleanIntersection(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat + BooleanIntersection{T,typeof(left),typeof(right)}(left, right, place) +end + +#---Utilities--------------------------------------------------------------------------------------- #---Printing and Plotting--------------------------------------------------------------------------- function Base.show(io::IO, shape::BooleanUnion{T, SL, SR}) where {T,SL,SR} @@ -36,7 +46,17 @@ function Base.show(io::IO, shape::BooleanIntersection{T, SL, SR}) where {T,SL,SR end -function GeometryBasics.mesh(shape::AbstractBoolean{T, SL, SR}) where {T,SL,SR} +function GeometryBasics.mesh(shape::BooleanUnion{T, SL, SR}) where {T,SL,SR} + (; left, right, transformation) = shape + rmesh = mesh(right) + merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))]) +end +function GeometryBasics.mesh(shape::BooleanSubtraction{T, SL, SR}) where {T,SL,SR} + (; left, right, transformation) = shape + rmesh = mesh(right) + merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))]) +end +function GeometryBasics.mesh(shape::BooleanIntersection{T, SL, SR}) where {T,SL,SR} (; left, right, transformation) = shape rmesh = mesh(right) merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))]) @@ -60,9 +80,22 @@ function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T end -function safetyToOut(shape::AbstractBoolean) +function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end + +function safetyToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 end -function safetyToIn(shape::AbstractBoolean) +function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 end diff --git a/src/GDML.jl b/src/GDML.jl index 52a00e1..73afee1 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -244,11 +244,11 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - push!(pvols, PlacedVolume(-1, transformation, daughter)) + push!(pvols, PlacedVolume{T}(-1, transformation, daughter)) end end end - volume = elemname == "volume" ? Volume(volname, shape, material) : Assembly{T}(volname) + volume = elemname == "volume" ? Volume{T}(volname, shape, material) : Assembly{T}(volname) for pvol in pvols placeDaughter!(volume, pvol.transformation, pvol.volume) end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 0541735..41984e6 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -16,7 +16,7 @@ export Triangle, Intersection, intersect, distanceToPlane export Tesselation, coordinates, faces, normals, mesh using Requires -using StaticArrays, GeometryBasics, LinearAlgebra, Rotations, AbstractTrees, Accessors +using StaticArrays, GeometryBasics, LinearAlgebra, Rotations, AbstractTrees include("BasicTypes.jl") include("Transformation3D.jl") include("TriangleIntersect.jl") diff --git a/src/Navigators.jl b/src/Navigators.jl index c8fb541..b90ac38 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -49,18 +49,18 @@ function Base.show(io::IO, nav::BVHNavigator{T}) where T end #---NavigatorState---------------------------------------------------------------------------------- -struct NavigatorState{T<:AbstractFloat, VT1, VT2, NAV<:AbstractNavigator} <: AbstractNavigatorState +mutable struct NavigatorState{T<:AbstractFloat, NAV<:AbstractNavigator} <: AbstractNavigatorState navigator::NAV # Navigator to be used - topvol::VT1 - currvol::VT2 + topvol::Volume{T} # Typically the unplaced world + currvol::Volume{T} # the current volume isinworld::Bool # inside world volume volstack::Vector{Int64} # keep the indexes of all daughters up to the current one tolocal::Vector{Transformation3D{T}} # Stack of transformations end #---Constructors------------------------------------------------------------------------------------ -function NavigatorState(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T - x = NavigatorState(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) +function NavigatorState{T}(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T + x = NavigatorState{T,typeof(nav)}(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) sizehint!(x.volstack,16) sizehint!(x.tolocal,16) return x @@ -69,8 +69,8 @@ end @inline function reset!(state::NavigatorState{T}) where T<:AbstractFloat empty!(state.volstack) empty!(state.tolocal) - @set state.currvol = state.topvol - @set state.isinworld = false + state.currvol = state.topvol + state.isinworld = false end @inline function currentVolume(state::NavigatorState{T}) where T<:AbstractFloat @@ -85,17 +85,17 @@ end for idx in state.volstack vol = vol.daughters[idx].volume end - @set state.currvol = vol + state.currvol = vol else - @set state.currvol = state.topvol - @set state.isinworld = false + state.currvol = state.topvol + state.isinworld = false end end @inline function pushIn!(state::NavigatorState{T}, pvol::PlacedVolume{T}) where T<:AbstractFloat push!(state.volstack, pvol.idx) push!(state.tolocal, pvol.transformation) - @set state.currvol = pvol.volume + state.currvol = pvol.volume end #---Basic loops implemented using acceleration structures------------------------------------------- @@ -129,8 +129,8 @@ function locateGlobalPoint!(state::NavigatorState{T,NAV}, point::Point3{T}) wher reset!(state) vol = state.topvol if contains(vol, point) - @set state.currvol = vol - @set state.isinworld = true + state.currvol = vol + state.isinworld = true end # check the daughters in a recursive manner isinside = true diff --git a/src/Volume.jl b/src/Volume.jl index 0399239..678bda7 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -5,39 +5,39 @@ end NoShape{T}() where T = NoShape{T,Nothing}([]) #---Shape-------------------------------------------------------------------------- -# const Shape{T} = Union{NoShape{T}, -# Box{T}, -# Trd{T}, -# Trap{T}, -# Tube{T}, -# Cone{T}, -# Polycone{T}, -# CutTube{T}, -# BooleanUnion{T}, -# BooleanIntersection{T}, -# BooleanSubtraction{T}} where T<:AbstractFloat +const Shape{T} = Union{NoShape{T}, + Box{T}, + Trd{T}, + Trap{T}, + Tube{T}, + Cone{T}, + Polycone{T}, + CutTube{T}, + BooleanUnion{T}, + BooleanIntersection{T}, + BooleanSubtraction{T}} where T<:AbstractFloat #---Volume------------------------------------------------------------------------- -struct Volume{T<:AbstractFloat, S, PV} +struct VolumeP{T<:AbstractFloat,PV} label::String - shape::S - material::Material{T} + shape::Shape{T} # Reference to the actual shape + material::Material{T} # Reference to material daughters::Vector{PV} end #---PlacedVolume------------------------------------------------------------------- -struct PlacedVolume{T<:AbstractFloat, VT} +struct PlacedVolume{T<:AbstractFloat} idx::Int64 transformation::Transformation3D{T} - volume::VT + volume::VolumeP{T,PlacedVolume{T}} end #---Convenient Alias to simplify signatures--------------------------------------- -# const Volume{T} = VolumeP{T,S,PlacedVolume{T}} where {T<:AbstractFloat, S<:AbstractShape{T}} +const Volume{T} = VolumeP{T,PlacedVolume{T}} where T<:AbstractFloat #---Constructor-------------------------------------------------------------------- -function Volume(label, shape, material::Material{T}) where T<:AbstractFloat - Volume(label, shape, material, Vector{PlacedVolume{T}}()) # call the default constructor +function Volume{T}(label::String, shape::Shape{T}, material::Material{T}) where T<:AbstractFloat + Volume{T}(label, shape, material, Vector{PlacedVolume{T}}()) # call the default constructor end #---Utilities---------------------------------------------------------------------- @@ -85,7 +85,7 @@ function getWorld(vol::Volume{T}) where T<:AbstractFloat box = Box{T}((high - low)/2. .+ 0.1) # increase the bounding box tra = Transformation3D{T}(one(RotMatrix3{T}), -(high + low)/2.) mat = Material{T}("vacuum"; density=0) - world = Volume("world", box, mat) + world = Volume{T}("world", box, mat) placeDaughter!(world, tra, vol) return world end From f87c1fb1bfdfd607af2a166144599d1d56e87166 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 16:34:32 -0400 Subject: [PATCH 14/16] Revert "Revert "Unityper does nothing"" This reverts commit 6cb6ee22021bb00e38fcbc4397f45b1223736ada. --- Project.toml | 1 + examples/XRay.jl | 2 +- src/Boolean.jl | 92 +++++++++++----------- src/Distance.jl | 23 ++---- src/DistanceIn.jl | 192 +++++++++++++++++++++++---------------------- src/DistanceOut.jl | 118 ++++++++++++++-------------- src/GDML.jl | 4 +- src/Geom4hep.jl | 1 + src/Navigators.jl | 5 +- src/Volume.jl | 4 +- 10 files changed, 218 insertions(+), 224 deletions(-) diff --git a/Project.toml b/Project.toml index 66ba0bb..d1964ae 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnionArrays = "d6dd79e4-993b-11e9-1366-0de1c9fe1122" +Unityper = "a7c27f48-0311-42f6-a7f8-2c11e75eb415" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/XRay.jl b/examples/XRay.jl index 7fd574f..e72db49 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -8,7 +8,7 @@ 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 - world = getWorld(vol) + world = vol # Setup plane of points and results lower, upper = extent(world.shape) ix, iy, iz = idx[view, :] diff --git a/src/Boolean.jl b/src/Boolean.jl index 538903e..75f61f5 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -1,31 +1,41 @@ -#---Boolean Types---------------------------------------------------------------------------------- -struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} - left::SL # the mother (or left) volume A in unplaced form - right::SR # (or right) volume B in placed form, acting on A with a boolean operation - transformation::Transformation3D{T} # placement of "right" with respect of "left" -end -struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} - left::SL # the mother (or left) volume A in unplaced form - right::SR # (or right) volume B in placed form, acting on A with a boolean operation - transformation::Transformation3D{T} # placement of "right" with respect of "left" -end -struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} - left::SL # the mother (or left) volume A in unplaced form - right::SR # (or right) volume B in placed form, acting on A with a boolean operation - transformation::Transformation3D{T} # placement of "right" with respect of "left" +@compactify begin + @abstract struct AbstractBoolean{T, SL, SR} <: AbstractShape{T} + left::SL = NoShape() + right::SR = NoShape() + transformation::Transformation3D{T} = Transformation3D{T}(0, 0, 0) + end + struct BooleanUnion{T, SL, SR} <: AbstractBoolean{T, SL, SR} end + struct BooleanSubtraction{T, SL, SR} <: AbstractBoolean{T, SL, SR} end + struct BooleanIntersection{T, SL, SR} <: AbstractBoolean{T, SL, SR} end end - #---Constructor------------------------------------------------------------------------------------ function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanUnion{T,typeof(left),typeof(right)}(left, right, place) + BooleanUnion{T,typeof(left),typeof(right)}(; left, right, transformation = place) end function BooleanSubtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanSubtraction{T,typeof(left),typeof(right)}(left, right, place) + BooleanSubtraction{T,typeof(left),typeof(right)}(; left, right, transformation = place) end function BooleanIntersection(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanIntersection{T,typeof(left),typeof(right)}(left, right, place) + BooleanIntersection{T,typeof(left),typeof(right)}(; left, right, transformation = place) end +#---Boolean Types---------------------------------------------------------------------------------- +# struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +# left::SL # the mother (or left) volume A in unplaced form +# right::SR # (or right) volume B in placed form, acting on A with a boolean operation +# transformation::Transformation3D{T} # placement of "right" with respect of "left" +# end +# struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +# left::SL # the mother (or left) volume A in unplaced form +# right::SR # (or right) volume B in placed form, acting on A with a boolean operation +# transformation::Transformation3D{T} # placement of "right" with respect of "left" +# end +# struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} +# left::SL # the mother (or left) volume A in unplaced form +# right::SR # (or right) volume B in placed form, acting on A with a boolean operation +# transformation::Transformation3D{T} # placement of "right" with respect of "left" +# end + #---Utilities--------------------------------------------------------------------------------------- #---Printing and Plotting--------------------------------------------------------------------------- @@ -63,39 +73,27 @@ function GeometryBasics.mesh(shape::BooleanIntersection{T, SL, SR}) where {T,SL, end #---Basic functions--------------------------------------------------------------------------------- -function extent(shape::BooleanUnion{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} - (; left, right, transformation) = shape - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (min.(minLeft, minRight), max.(maxLeft, maxRight)) -end -function extent(shape::BooleanSubtraction{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} - extent(shape.left) -end -function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} +function extent(shape::AbstractBoolean) (; left, right, transformation) = shape - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (max.(minLeft, minRight), min.(maxLeft, maxRight)) + @compactified shape::AbstractBoolean begin + BooleanUnion => begin + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (min.(minLeft, minRight), max.(maxLeft, maxRight)) + end + BooleanIntersection => begin + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (max.(minLeft, minRight), min.(maxLeft, maxRight)) + end + BooleanSubtraction => extent(shape.left) + end end -function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end - -function safetyToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} - return 0 -end -function safetyToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} +function safetyToOut(shape::AbstractBoolean) return 0 end -function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} +function safetyToIn(shape::AbstractBoolean) return 0 end diff --git a/src/Distance.jl b/src/Distance.jl index 146ec50..37e11f0 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -1,4 +1,4 @@ -function distanceToIn(shape, point, dir) +function distanceToIn(shape, point::Point3{T}, dir)::T where T if shape isa Trap distanceToIn_trap(shape, point, dir) elseif shape isa Trd @@ -13,18 +13,11 @@ function distanceToIn(shape, point, dir) distanceToIn_polycone(shape, point, dir) elseif shape isa CutTube distanceToIn_cuttube(shape, point, dir) - elseif shape isa BooleanUnion - distanceToIn_booleanunion(shape, point, dir) - elseif shape isa BooleanSubtraction - distanceToIn_booleansubtraction(shape, point, dir) - elseif shape isa BooleanIntersection - distanceToIn_booleanintersection(shape, point, dir) - elseif shape isa PlacedVolume - xf = shape.transformation - distanceToIn_placedvolume(shape.volume.shape, xf*point, xf*dir) + elseif shape isa AbstractBoolean + distanceToIn_boolean(shape, point, dir) end end -function distanceToOut(shape, point, dir) +function distanceToOut(shape, point::Point3{T}, dir)::T where T if shape isa Trap distanceToOut_trap(shape, point, dir) elseif shape isa Trd @@ -39,11 +32,7 @@ function distanceToOut(shape, point, dir) distanceToOut_polycone(shape, point, dir) elseif shape isa CutTube distanceToOut_cuttube(shape, point, dir) - elseif shape isa BooleanUnion - distanceToOut_booleanunion(shape, point, dir) - elseif shape isa BooleanSubtraction - distanceToOut_booleansubtraction(shape, point, dir) - elseif shape isa BooleanIntersection - distanceToOut_booleanintersection(shape, point, dir) + elseif shape isa AbstractBoolean + distanceToOut_boolean(shape, point, dir) end end diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index 4401214..0e652f6 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -1,102 +1,104 @@ -function distanceToIn_placedvolume(pvol, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat - distanceToIn(pvol, p, d) -end ## Boolean -function distanceToIn_booleanunion(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_booleanintersection(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 +function distanceToIn_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} + @compactified shape::AbstractBoolean begin + BooleanUnion => begin + (; left, right, transformation) = shape + distA = distanceToIn(left, point, dir) + distB = distanceToIn(right, transformation * point, transformation * dir) + return min(distA, distB) end - end - return dist -end -function distanceToIn_booleansubtraction(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 + BooleanIntersection => begin + (; 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 - - # 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 + BooleanSubtraction => begin + (; 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 + 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/DistanceOut.jl b/src/DistanceOut.jl index 6ab14cf..6205b60 100644 --- a/src/DistanceOut.jl +++ b/src/DistanceOut.jl @@ -1,70 +1,74 @@ #Boolean -function distanceToOut_booleanunion(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) +function distanceToOut_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} + @compactified shape::AbstractBoolean begin + BooleanUnion => begin + (; 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 + 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 - break + return T(-1) end - end - return dist - kPushTolerance(T) - else - return T(-1) + end end - end -end -function distanceToOut_booleanintersection(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_booleansubtraction(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) + BooleanIntersection => begin + (; left, right, transformation) = shape + distA = distanceToOut(left, point, dir) + distB = distanceToOut(right, transformation * point, transformation * dir) + return min(distA, distB) + end + BooleanSubtraction => begin + (; left, right, transformation) = shape + distA = distanceToOut(left, point, dir) + distB = distanceToIn(right, transformation * point, transformation * dir) + return min(distA, distB) + end + end end #Plane diff --git a/src/GDML.jl b/src/GDML.jl index 73afee1..0f36d69 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -225,7 +225,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac volname = attrs["name"] shape = nothing material = nothing - pvols = PlacedVolume{T}[] + pvols = [] for cc in child_nodes(e) #--- loop over ref to solid, material and daughters if is_elementnode(cc) aa = attributes_dict(XMLElement(cc)) @@ -244,7 +244,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - push!(pvols, PlacedVolume{T}(-1, transformation, daughter)) + push!(pvols, PlacedVolume(-1, transformation, daughter)) end end end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 41984e6..6f2e4cc 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -27,6 +27,7 @@ include("Tube.jl") include("Cone.jl") include("Polycone.jl") include("CutTube.jl") +using Unityper include("Boolean.jl") include("Trap.jl") include("Materials.jl") diff --git a/src/Navigators.jl b/src/Navigators.jl index b90ac38..ab1e1c3 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -158,10 +158,11 @@ function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir:: step = step_limit candidate = 0 #---Linear loop over the daughters------------------------------------------------- - for i in intersectedDaughters(nav, volume, point, dir) + @inbounds for i in intersectedDaughters(nav, volume, point, dir) pvol = volume.daughters[i] #---Assuming that it is not yet inside the daughter (otherwise it returns -1.) - dist = distanceToIn(pvol, point, dir) + xf = pvol.transformation + dist = distanceToIn(pvol.volume.shape, xf*point, xf*dir) dist < 0. && return (zero(T), pvol.idx ) if dist > 0. && dist != Inf && dist < step step = dist diff --git a/src/Volume.jl b/src/Volume.jl index 678bda7..00ad235 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -13,9 +13,7 @@ const Shape{T} = Union{NoShape{T}, Cone{T}, Polycone{T}, CutTube{T}, - BooleanUnion{T}, - BooleanIntersection{T}, - BooleanSubtraction{T}} where T<:AbstractFloat + AbstractBoolean{T}} where T<:AbstractFloat #---Volume------------------------------------------------------------------------- struct VolumeP{T<:AbstractFloat,PV} From af5f924e38074adb3dc21f85c0dc44b3ece5ed60 Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 16:34:33 -0400 Subject: [PATCH 15/16] Revert "Unityper does nothing" This reverts commit a90b0a1b9cd8ba8ee129365db1623f9c10726980. --- Project.toml | 1 - examples/XRay.jl | 2 +- src/Boolean.jl | 92 +++++++++++----------- src/Distance.jl | 23 ++++-- src/DistanceIn.jl | 192 ++++++++++++++++++++++----------------------- src/DistanceOut.jl | 118 ++++++++++++++-------------- src/GDML.jl | 4 +- src/Geom4hep.jl | 1 - src/Navigators.jl | 5 +- src/Volume.jl | 4 +- 10 files changed, 224 insertions(+), 218 deletions(-) diff --git a/Project.toml b/Project.toml index d1964ae..66ba0bb 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,6 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnionArrays = "d6dd79e4-993b-11e9-1366-0de1c9fe1122" -Unityper = "a7c27f48-0311-42f6-a7f8-2c11e75eb415" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/XRay.jl b/examples/XRay.jl index e72db49..7fd574f 100644 --- a/examples/XRay.jl +++ b/examples/XRay.jl @@ -8,7 +8,7 @@ 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 - world = vol + world = getWorld(vol) # Setup plane of points and results lower, upper = extent(world.shape) ix, iy, iz = idx[view, :] diff --git a/src/Boolean.jl b/src/Boolean.jl index 75f61f5..538903e 100644 --- a/src/Boolean.jl +++ b/src/Boolean.jl @@ -1,41 +1,31 @@ -@compactify begin - @abstract struct AbstractBoolean{T, SL, SR} <: AbstractShape{T} - left::SL = NoShape() - right::SR = NoShape() - transformation::Transformation3D{T} = Transformation3D{T}(0, 0, 0) - end - struct BooleanUnion{T, SL, SR} <: AbstractBoolean{T, SL, SR} end - struct BooleanSubtraction{T, SL, SR} <: AbstractBoolean{T, SL, SR} end - struct BooleanIntersection{T, SL, SR} <: AbstractBoolean{T, SL, SR} end +#---Boolean Types---------------------------------------------------------------------------------- +struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} + left::SL # the mother (or left) volume A in unplaced form + right::SR # (or right) volume B in placed form, acting on A with a boolean operation + transformation::Transformation3D{T} # placement of "right" with respect of "left" +end +struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} + left::SL # the mother (or left) volume A in unplaced form + right::SR # (or right) volume B in placed form, acting on A with a boolean operation + transformation::Transformation3D{T} # placement of "right" with respect of "left" +end +struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} + left::SL # the mother (or left) volume A in unplaced form + right::SR # (or right) volume B in placed form, acting on A with a boolean operation + transformation::Transformation3D{T} # placement of "right" with respect of "left" end + #---Constructor------------------------------------------------------------------------------------ function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanUnion{T,typeof(left),typeof(right)}(; left, right, transformation = place) + BooleanUnion{T,typeof(left),typeof(right)}(left, right, place) end function BooleanSubtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanSubtraction{T,typeof(left),typeof(right)}(; left, right, transformation = place) + BooleanSubtraction{T,typeof(left),typeof(right)}(left, right, place) end function BooleanIntersection(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat - BooleanIntersection{T,typeof(left),typeof(right)}(; left, right, transformation = place) + BooleanIntersection{T,typeof(left),typeof(right)}(left, right, place) end -#---Boolean Types---------------------------------------------------------------------------------- -# struct BooleanUnion{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} -# left::SL # the mother (or left) volume A in unplaced form -# right::SR # (or right) volume B in placed form, acting on A with a boolean operation -# transformation::Transformation3D{T} # placement of "right" with respect of "left" -# end -# struct BooleanSubtraction{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} -# left::SL # the mother (or left) volume A in unplaced form -# right::SR # (or right) volume B in placed form, acting on A with a boolean operation -# transformation::Transformation3D{T} # placement of "right" with respect of "left" -# end -# struct BooleanIntersection{T<:AbstractFloat, SL<:AbstractShape, SR<:AbstractShape} <: AbstractShape{T} -# left::SL # the mother (or left) volume A in unplaced form -# right::SR # (or right) volume B in placed form, acting on A with a boolean operation -# transformation::Transformation3D{T} # placement of "right" with respect of "left" -# end - #---Utilities--------------------------------------------------------------------------------------- #---Printing and Plotting--------------------------------------------------------------------------- @@ -73,27 +63,39 @@ function GeometryBasics.mesh(shape::BooleanIntersection{T, SL, SR}) where {T,SL, end #---Basic functions--------------------------------------------------------------------------------- -function extent(shape::AbstractBoolean) +function extent(shape::BooleanUnion{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} + (; left, right, transformation) = shape + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (min.(minLeft, minRight), max.(maxLeft, maxRight)) +end +function extent(shape::BooleanSubtraction{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} + extent(shape.left) +end +function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T}} where {T,SL,SR} (; left, right, transformation) = shape - @compactified shape::AbstractBoolean begin - BooleanUnion => begin - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (min.(minLeft, minRight), max.(maxLeft, maxRight)) - end - BooleanIntersection => begin - minLeft, maxLeft = extent(left) - minRight, maxRight = extent(right) .* Ref(transformation) - (max.(minLeft, minRight), min.(maxLeft, maxRight)) - end - BooleanSubtraction => extent(shape.left) - end + minLeft, maxLeft = extent(left) + minRight, maxRight = extent(right) .* Ref(transformation) + (max.(minLeft, minRight), min.(maxLeft, maxRight)) end -function safetyToOut(shape::AbstractBoolean) +function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToOut(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToOut(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end + +function safetyToIn(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} + return 0 +end +function safetyToIn(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 end -function safetyToIn(shape::AbstractBoolean) +function safetyToIn(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::T where {T,SL,SR} return 0 end diff --git a/src/Distance.jl b/src/Distance.jl index 37e11f0..146ec50 100644 --- a/src/Distance.jl +++ b/src/Distance.jl @@ -1,4 +1,4 @@ -function distanceToIn(shape, point::Point3{T}, dir)::T where T +function distanceToIn(shape, point, dir) if shape isa Trap distanceToIn_trap(shape, point, dir) elseif shape isa Trd @@ -13,11 +13,18 @@ function distanceToIn(shape, point::Point3{T}, dir)::T where T distanceToIn_polycone(shape, point, dir) elseif shape isa CutTube distanceToIn_cuttube(shape, point, dir) - elseif shape isa AbstractBoolean - distanceToIn_boolean(shape, point, dir) + elseif shape isa BooleanUnion + distanceToIn_booleanunion(shape, point, dir) + elseif shape isa BooleanSubtraction + distanceToIn_booleansubtraction(shape, point, dir) + elseif shape isa BooleanIntersection + distanceToIn_booleanintersection(shape, point, dir) + elseif shape isa PlacedVolume + xf = shape.transformation + distanceToIn_placedvolume(shape.volume.shape, xf*point, xf*dir) end end -function distanceToOut(shape, point::Point3{T}, dir)::T where T +function distanceToOut(shape, point, dir) if shape isa Trap distanceToOut_trap(shape, point, dir) elseif shape isa Trd @@ -32,7 +39,11 @@ function distanceToOut(shape, point::Point3{T}, dir)::T where T distanceToOut_polycone(shape, point, dir) elseif shape isa CutTube distanceToOut_cuttube(shape, point, dir) - elseif shape isa AbstractBoolean - distanceToOut_boolean(shape, point, dir) + elseif shape isa BooleanUnion + distanceToOut_booleanunion(shape, point, dir) + elseif shape isa BooleanSubtraction + distanceToOut_booleansubtraction(shape, point, dir) + elseif shape isa BooleanIntersection + distanceToOut_booleanintersection(shape, point, dir) end end diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index 0e652f6..4401214 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -1,104 +1,102 @@ +function distanceToIn_placedvolume(pvol, p::Point3{T}, d::Vector3{T})::T where T<:AbstractFloat + distanceToIn(pvol, p, d) +end ## Boolean -function distanceToIn_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} - @compactified shape::AbstractBoolean begin - BooleanUnion => begin - (; left, right, transformation) = shape - distA = distanceToIn(left, point, dir) - distB = distanceToIn(right, transformation * point, transformation * dir) - return min(distA, distB) +function distanceToIn_booleanunion(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_booleanintersection(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 - BooleanIntersection => begin - (; 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 + return dist +end +function distanceToIn_booleansubtraction(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 - BooleanSubtraction => begin - (; 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 - return dist + + # 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/DistanceOut.jl b/src/DistanceOut.jl index 6205b60..6ab14cf 100644 --- a/src/DistanceOut.jl +++ b/src/DistanceOut.jl @@ -1,74 +1,70 @@ #Boolean -function distanceToOut_boolean(shape::AbstractBoolean, point::Point3{T}, dir::Vector3{T})::T where {T} - @compactified shape::AbstractBoolean begin - BooleanUnion => begin - (; left, right, transformation) = shape - dist = T(0) - positionA = inside(left, point) - if positionA != kOutside # point inside A - while(true) - distA = distanceToOut(left, point, dir) +function distanceToOut_booleanunion(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 # 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 + if inside(right, lpoint) == kOutside 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) + break end - end - end - BooleanIntersection => begin - (; left, right, transformation) = shape - distA = distanceToOut(left, point, dir) - distB = distanceToOut(right, transformation * point, transformation * dir) - return min(distA, distB) - end - BooleanSubtraction => begin - (; left, right, transformation) = shape - distA = distanceToOut(left, point, dir) - distB = distanceToIn(right, transformation * point, transformation * dir) - return min(distA, distB) + end + return dist - kPushTolerance(T) + else + return T(-1) end - end + end +end +function distanceToOut_booleanintersection(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_booleansubtraction(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 #Plane diff --git a/src/GDML.jl b/src/GDML.jl index 0f36d69..73afee1 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -225,7 +225,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac volname = attrs["name"] shape = nothing material = nothing - pvols = [] + pvols = PlacedVolume{T}[] for cc in child_nodes(e) #--- loop over ref to solid, material and daughters if is_elementnode(cc) aa = attributes_dict(XMLElement(cc)) @@ -244,7 +244,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - push!(pvols, PlacedVolume(-1, transformation, daughter)) + push!(pvols, PlacedVolume{T}(-1, transformation, daughter)) end end end diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 6f2e4cc..41984e6 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -27,7 +27,6 @@ include("Tube.jl") include("Cone.jl") include("Polycone.jl") include("CutTube.jl") -using Unityper include("Boolean.jl") include("Trap.jl") include("Materials.jl") diff --git a/src/Navigators.jl b/src/Navigators.jl index ab1e1c3..b90ac38 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -158,11 +158,10 @@ function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir:: step = step_limit candidate = 0 #---Linear loop over the daughters------------------------------------------------- - @inbounds for i in intersectedDaughters(nav, volume, point, dir) + for i in intersectedDaughters(nav, volume, point, dir) pvol = volume.daughters[i] #---Assuming that it is not yet inside the daughter (otherwise it returns -1.) - xf = pvol.transformation - dist = distanceToIn(pvol.volume.shape, xf*point, xf*dir) + dist = distanceToIn(pvol, point, dir) dist < 0. && return (zero(T), pvol.idx ) if dist > 0. && dist != Inf && dist < step step = dist diff --git a/src/Volume.jl b/src/Volume.jl index 00ad235..678bda7 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -13,7 +13,9 @@ const Shape{T} = Union{NoShape{T}, Cone{T}, Polycone{T}, CutTube{T}, - AbstractBoolean{T}} where T<:AbstractFloat + BooleanUnion{T}, + BooleanIntersection{T}, + BooleanSubtraction{T}} where T<:AbstractFloat #---Volume------------------------------------------------------------------------- struct VolumeP{T<:AbstractFloat,PV} From 7d19ffc3aa2dbf5b49a55ba60ecac4d380986d2d Mon Sep 17 00:00:00 2001 From: Jerry Ling Date: Fri, 23 Sep 2022 17:21:40 -0400 Subject: [PATCH 16/16] remove some alias --- examples/XRay.jl | 2 +- src/DistanceIn.jl | 3 ++- src/GDML.jl | 14 +++++++------- src/Geom4hep.jl | 2 +- src/Navigators.jl | 36 ++++++++++++++++++------------------ src/Volume.jl | 31 ++++++++++++++++--------------- 6 files changed, 45 insertions(+), 43 deletions(-) diff --git a/examples/XRay.jl b/examples/XRay.jl index 7fd574f..b71cf62 100644 --- 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) diff --git a/src/DistanceIn.jl b/src/DistanceIn.jl index 4401214..0b52b38 100644 --- a/src/DistanceIn.jl +++ b/src/DistanceIn.jl @@ -59,7 +59,8 @@ function distanceToIn_booleanintersection(shape::BooleanIntersection{T, SL, SR}, end return dist end -function distanceToIn_booleansubtraction(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR} +function distanceToIn_booleansubtraction(shape::BooleanSubtraction{T, SL, SR}, + point::Point3{T}, dir::Vector3{T})::T where {T,SL,SR} (; left, right, transformation) = shape lpoint = transformation * point diff --git a/src/GDML.jl b/src/GDML.jl index 73afee1..829a56e 100644 --- a/src/GDML.jl +++ b/src/GDML.jl @@ -3,16 +3,16 @@ name = LightXML.name include("Units.jl") -struct GDMLDicts{T<:AbstractFloat} +struct GDMLDicts{T<:AbstractFloat, S} materials::Dict{String,AbstractMaterial{T}} solids::Dict{String,AbstractShape{T}} - volumes::Dict{String,Volume{T}} + volumes::Dict{String,VolumeP{T, S}} positions::Dict{String, Vector3{T}} rotations::Dict{String, RotXYZ{T}} - function GDMLDicts{T}() where T<:AbstractFloat + function GDMLDicts{T, S}() where {T<:AbstractFloat, S} new(Dict{String,AbstractMaterial{T}}(), Dict{String,AbstractShape{T}}(), - Dict{String,Volume{T}}(), + Dict{String,VolumeP{T, S}}(), Dict{String, Vector3{T}}(), Dict{String, RotXYZ{T}}() ) @@ -20,7 +20,7 @@ struct GDMLDicts{T<:AbstractFloat} end # process -function fillDefine!(dicts::GDMLDicts{T}, element::XMLElement) where T<:AbstractFloat +function fillDefine!(dicts::GDMLDicts{T, S}, element::XMLElement) where {T<:AbstractFloat, S} (; positions, rotations) = dicts for c in child_nodes(element) if is_elementnode(c) @@ -248,7 +248,7 @@ function fillVolumes!(dicts::GDMLDicts{T}, element::XMLElement) where T<:Abstrac end end end - volume = elemname == "volume" ? Volume{T}(volname, shape, material) : Assembly{T}(volname) + volume = elemname == "volume" ? VolumeP(volname, shape, material) : Assembly{T}(volname) for pvol in pvols placeDaughter!(volume, pvol.transformation, pvol.volume) end @@ -289,7 +289,7 @@ end # process the full file function processGDML(fname::String, ::Type{T}=Float64) where T - dicts = GDMLDicts{T}() + dicts = GDMLDicts{T, PlacedVolume{T}}() world = nothing xdoc = parse_file(fname) xroot = root(xdoc) diff --git a/src/Geom4hep.jl b/src/Geom4hep.jl index 41984e6..1366f01 100644 --- a/src/Geom4hep.jl +++ b/src/Geom4hep.jl @@ -7,7 +7,7 @@ export CutTube, Plane, Boolean, Trap, BooleanUnion export getindex, capacity, surface, extent, normal, distanceToOut, distanceToIn, inside, safetyToOut, safetyToIn export Material, Isotope, Element export Transformation3D, RotMatrix3, RotXYZ, one, isone, transform, hasrotation, hastranslation, inv, lmul! -export PlacedVolume, Volume, Assembly, placeDaughter!, draw, draw!, drawDistanceToIn, drawDistanceToOut, getWorld, children +export PlacedVolume, VolumeP, Assembly, placeDaughter!, draw, draw!, drawDistanceToIn, drawDistanceToOut, getWorld, children export AABB, area, volume, BVHParam, buildBVH, BVH, pvolindices, pushPvolIndices! export AbstractNavigator, TrivialNavigator, BVHNavigator, NavigatorState, computeStep!, locateGlobalPoint!, reset!, isInVolume, currentVolume, getClosestDaughter, containedDaughters export kTolerance diff --git a/src/Navigators.jl b/src/Navigators.jl index b90ac38..86b2e68 100644 --- a/src/Navigators.jl +++ b/src/Navigators.jl @@ -3,31 +3,31 @@ abstract type AbstractNavigator end abstract type AbstractNavigatorState end #---Define the various navigators------------------------------------------------------------------- -struct TrivialNavigator{T<:AbstractFloat} <: AbstractNavigator - world::Volume{T} +struct TrivialNavigator{T<:AbstractFloat, S} <: AbstractNavigator + world::VolumeP{T, S} end -function TrivialNavigator(world::Volume{T}) where T +function TrivialNavigator(world::VolumeP{T, PV}) where {T, PV} TrivialNavigator{T}(world) end #---BVHNavigator------------------------------------------------------------------------------------- -struct BVHNavigator{T<:AbstractFloat} <: AbstractNavigator - world::Volume{T} +struct BVHNavigator{T<:AbstractFloat, S} <: AbstractNavigator + world::VolumeP{T, S} param::BVHParam{T} bvhdict::Dict{UInt64,BVH{T}} # Used to cache the lists of indices to placed volumes (avoid useless allocations) pvolind::Vector{Int64} end -function BVHNavigator(world::Volume{T}; SAHfrac::T=T(8)) where T - nav = BVHNavigator{T}(world, BVHParam{T}(SAHfrac), Dict{UInt64,BVH{T}}(),Int64[]) +function BVHNavigator(world::VolumeP{T, S}; SAHfrac::T=T(8)) where {T, S} + nav = BVHNavigator{T, S}(world, BVHParam{T}(SAHfrac), Dict{UInt64,BVH{T}}(),Int64[]) #--- Create accelerationstructures--------------------------------- buildBVH(nav, world, Set{UInt64}()) sizehint!(nav.pvolind,1024) nav end -function buildBVH(nav::BVHNavigator{T}, vol::Volume{T}, set::Set{UInt64}) where T +function buildBVH(nav::BVHNavigator{T}, vol::VolumeP{T, S}, set::Set{UInt64}) where {T, S} id = objectid(vol) #---return immediatelly if BVH is already there---- id in set && return @@ -49,18 +49,18 @@ function Base.show(io::IO, nav::BVHNavigator{T}) where T end #---NavigatorState---------------------------------------------------------------------------------- -mutable struct NavigatorState{T<:AbstractFloat, NAV<:AbstractNavigator} <: AbstractNavigatorState +mutable struct NavigatorState{T<:AbstractFloat, S, NAV<:AbstractNavigator} <: AbstractNavigatorState navigator::NAV # Navigator to be used - topvol::Volume{T} # Typically the unplaced world - currvol::Volume{T} # the current volume + topvol::VolumeP{T,S} # Typically the unplaced world + currvol::VolumeP{T,S} # the current volume isinworld::Bool # inside world volume volstack::Vector{Int64} # keep the indexes of all daughters up to the current one tolocal::Vector{Transformation3D{T}} # Stack of transformations end #---Constructors------------------------------------------------------------------------------------ -function NavigatorState{T}(top::Volume{T}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where T - x = NavigatorState{T,typeof(nav)}(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) +function NavigatorState{T}(top::VolumeP{T, S}, nav::AbstractNavigator=TrivialNavigator{T}(top)) where {T, S} + x = NavigatorState(nav, top, top, false, Vector{Int64}(), Vector{Transformation3D{T}}()) sizehint!(x.volstack,16) sizehint!(x.tolocal,16) return x @@ -99,8 +99,8 @@ end end #---Basic loops implemented using acceleration structures------------------------------------------- -containedDaughters(::TrivialNavigator{T}, vol::Volume{T}, ::Point3{T}) where T = eachindex(vol.daughters) -function containedDaughters(nav::BVHNavigator{T}, vol::Volume{T}, point::Point3{T}) where T +containedDaughters(::TrivialNavigator{T}, vol::VolumeP, ::Point3{T}) where T = eachindex(vol.daughters) +function containedDaughters(nav::BVHNavigator{T}, vol::VolumeP, point::Point3{T}) where T id = objectid(vol) if haskey(nav.bvhdict, id) bvh = nav.bvhdict[id] @@ -111,8 +111,8 @@ function containedDaughters(nav::BVHNavigator{T}, vol::Volume{T}, point::Point3{ return nav.pvolind end -intersectedDaughters(::TrivialNavigator{T}, vol::Volume{T}, ::Point3{T}, ::Vector3{T}) where T = eachindex(vol.daughters) -function intersectedDaughters(nav::BVHNavigator{T}, vol::Volume{T}, point::Point3{T}, dir::Vector3{T}) where T +intersectedDaughters(::TrivialNavigator{T}, vol::VolumeP, ::Point3{T}, ::Vector3{T}) where T = eachindex(vol.daughters) +function intersectedDaughters(nav::BVHNavigator{T}, vol::VolumeP, point::Point3{T}, dir::Vector3{T}) where T id = objectid(vol) if haskey(nav.bvhdict, id) bvh = nav.bvhdict[id] @@ -154,7 +154,7 @@ end state.isinworld end -function getClosestDaughter(nav::NAV, volume::Volume{T}, point::Point3{T}, dir::Vector3{T}, step_limit::T ) where {T,NAV} +function getClosestDaughter(nav::NAV, volume::VolumeP, point::Point3{T}, dir::Vector3{T}, step_limit::T ) where {T,NAV} step = step_limit candidate = 0 #---Linear loop over the daughters------------------------------------------------- diff --git a/src/Volume.jl b/src/Volume.jl index 678bda7..2e058f1 100644 --- a/src/Volume.jl +++ b/src/Volume.jl @@ -33,24 +33,25 @@ struct PlacedVolume{T<:AbstractFloat} end #---Convenient Alias to simplify signatures--------------------------------------- -const Volume{T} = VolumeP{T,PlacedVolume{T}} where T<:AbstractFloat +# const Volume{T} = VolumeP{T,PlacedVolume{T}} where T<:AbstractFloat #---Constructor-------------------------------------------------------------------- -function Volume{T}(label::String, shape::Shape{T}, material::Material{T}) where T<:AbstractFloat - Volume{T}(label, shape, material, Vector{PlacedVolume{T}}()) # call the default constructor +function VolumeP(label::String, shape::Shape{T}, material::Material{T}) where T<:AbstractFloat + PV = PlacedVolume{T} + VolumeP{T, PV}(label, shape, material, Vector{PV}()) # call the default constructor end #---Utilities---------------------------------------------------------------------- -function Base.show(io::IO, vol::Volume{T}) where T - name = vol.label - print(io, "Volume{$T} name = $name") -end +# function Base.show(io::IO, vol::Volume{T}) where T +# name = vol.label +# print(io, "Volume{$T} name = $name") +# end -function AbstractTrees.children(vol::Volume{T}) where T - collect((d.volume for d in vol.daughters)) -end +# function AbstractTrees.children(vol::Volume{T}) where T +# collect((d.volume for d in vol.daughters)) +# end -function Base.getindex(vol::Volume{T}, indx...) where T +function Base.getindex(vol::VolumeP, indx...)::VolumeP v = vol for i in indx v = v.daughters[i].volume @@ -61,12 +62,12 @@ end function contains(pvol::PlacedVolume{T}, p::Point3{T})::Bool where T<:AbstractFloat inside(pvol.volume.shape, pvol.transformation * p) == kInside end -function contains(vol::Volume{T}, p::Point3{T})::Bool where T<:AbstractFloat +function contains(vol::VolumeP, p::Point3{T})::Bool where T<:AbstractFloat inside(vol.shape, p) == kInside end -function placeDaughter!(volume::Volume{T}, placement::Transformation3D{T}, subvol::Volume{T}) where T<:AbstractFloat +function placeDaughter!(volume::VolumeP, placement::Transformation3D{T}, subvol::VolumeP) where T<:AbstractFloat if subvol.shape isa NoShape for d in subvol.daughters push!(volume.daughters, PlacedVolume(length(volume.daughters)+1, d.transformation * placement, d.volume)) @@ -77,7 +78,7 @@ function placeDaughter!(volume::Volume{T}, placement::Transformation3D{T}, subvo volume end -function getWorld(vol::Volume{T}) where T<:AbstractFloat +function getWorld(vol::VolumeP{T, S}) where {T<:AbstractFloat, S} if vol.label == "world" return vol else @@ -85,7 +86,7 @@ function getWorld(vol::Volume{T}) where T<:AbstractFloat box = Box{T}((high - low)/2. .+ 0.1) # increase the bounding box tra = Transformation3D{T}(one(RotMatrix3{T}), -(high + low)/2.) mat = Material{T}("vacuum"; density=0) - world = Volume{T}("world", box, mat) + world = VolumeP("world", box, mat) placeDaughter!(world, tra, vol) return world end