Skip to content

Commit 9de89eb

Browse files
MoelfMoelf
committed
recover original performance
1 parent 2f2a458 commit 9de89eb

File tree

11 files changed

+237
-213
lines changed

11 files changed

+237
-213
lines changed

‎src/Boolean.jl‎

Lines changed: 0 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -80,69 +80,6 @@ function extent(shape::BooleanIntersection{T, SL, SR})::Tuple{Point3{T},Point3{T
8080
end
8181

8282

83-
function inside(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
84-
(; left, right, transformation) = shape
85-
lpoint = transformation * point
86-
87-
positionA = inside(left, point)
88-
positionA == kInside && return kInside
89-
90-
positionB = inside(right, lpoint)
91-
positionB == kInside && return kInside
92-
93-
if positionA == kSurface && positionB == kSurface
94-
normalA = normal(left, point)
95-
normalB = normal(right, lpoint) * transformation
96-
if dot(normalA, normalB) < 0
97-
return kInside # touching solids -)(-
98-
else
99-
return kSurface # overlapping solids =))
100-
end
101-
elseif positionA == kSurface || positionB == kSurface
102-
return kSurface
103-
else
104-
return kOutside
105-
end
106-
end
107-
108-
function inside(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
109-
(; left, right, transformation) = shape
110-
lpoint = transformation * point
111-
112-
positionA = inside(left, point)
113-
positionA == kOutside && return kOutside
114-
115-
positionB = inside(right, lpoint)
116-
if positionA == kInside && positionB == kInside
117-
return kInside
118-
elseif positionA == kInside && positionB == kSurface ||
119-
positionA == kSurface && positionB == kInside ||
120-
positionA == kSurface && positionB == kSurface
121-
return kSurface
122-
else
123-
return kOutside
124-
end
125-
end
126-
127-
function inside(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
128-
(; left, right, transformation) = shape
129-
lpoint = transformation * point
130-
131-
positionA = inside(left, point)
132-
positionA == kOutside && return kOutside
133-
134-
positionB = inside(right, lpoint)
135-
136-
if positionA == kInside && positionB == kOutside
137-
return kInside;
138-
elseif positionA == kInside && positionB == kSurface ||
139-
positionB == kOutside && positionA == kSurface
140-
return kSurface
141-
else
142-
return kOutside
143-
end
144-
end
145-
14683
function safetyToOut(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::T where {T,SL,SR}
14784
return 0
14885
end

‎src/Box.jl‎

Lines changed: 0 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -28,19 +28,6 @@ function normal(box::Box{T}, point::Point3{T}) where T
2828
end
2929
end
3030

31-
function inside(box::Box{T}, point::Point3{T})::Int64 where T<:AbstractFloat
32-
dist = -Inf
33-
for i in 1:3
34-
d = abs(point[i]) - box.fDimensions[i]
35-
if d > dist
36-
dist = d
37-
end
38-
end
39-
abs(dist) <= kTolerance(T)/2 ? kSurface : dist < 0.0 ? kInside : kOutside
40-
#dist = maximum(abs.(point) - box.fDimensions)
41-
#isapprox(dist, 0.0, atol = kTolerance(T)/2) ? kSurface : dist < 0.0 ? kInside : kOutside
42-
end
43-
4431
function safetyToOut(box::Box{T}, point::Point3{T}) where T<:AbstractFloat
4532
minimum(box.fDimensions - abs.(point))
4633
end

‎src/Cone.jl‎

Lines changed: 0 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -93,38 +93,6 @@ function extent(cone::Cone{T})::Tuple{Point3{T},Point3{T}} where T<:AbstractFloa
9393
extent(Tube{T}(min(rmin1,rmin2), max(rmax1,rmax2),z, ϕ₀, Δϕ))
9494
end
9595

96-
function inside(cone::Cone{T}, point::Point3{T})::Int64 where T<:AbstractFloat
97-
x, y, z = point
98-
(; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge,
99-
outerSlope, outerOffset, innerSlope, innerOffset) = cone
100-
101-
# Check Z
102-
outside = abs(z) > cone.z + kTolerance(T)/2
103-
outside && return kOutside
104-
cinside = abs(z) < cone.z - kTolerance(T)/2
105-
106-
# Check on RMax
107-
r2 = x * x + y * y
108-
rmax = rmax1 == rmax2 ? rmax1 : outerOffset + outerSlope * z
109-
outside |= r2 > rmax * rmax + kTolerance(T) * rmax
110-
outside && return kOutside
111-
cinside &= r2 < rmax * rmax - kTolerance(T) * rmax
112-
113-
# Check on RMin
114-
if rmin1 > 0. || rmin2 > 0.
115-
rmin = rmin1 == rmin2 ? rmin1 : innerOffset + innerSlope * z
116-
outside |= r2 <= rmin * rmin - kTolerance(T) * rmin
117-
outside && return kOutside
118-
cinside &= r2 > rmin * rmin + kTolerance(T) * rmin
119-
end
120-
121-
# Check on Phi
122-
if Δϕ < 2π
123-
outside |= isOutside(ϕWedge, x, y)
124-
cinside &= isInside(ϕWedge, x, y)
125-
end
126-
return outside ? kOutside : cinside ? kInside : kSurface
127-
end
12896

12997
function safetyToIn(cone::Cone{T}, point::Point3{T}) where T<:AbstractFloat
13098
x,y,z = point

‎src/CutTube.jl‎

Lines changed: 0 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -95,23 +95,6 @@ function extent(ctube::CutTube{T})::Tuple{Point3{T},Point3{T}} where T<:Abstract
9595
return (Point3{T}(aMin), Point3{T}(aMax))
9696
end
9797

98-
function inside(ctube::CutTube{T}, point::Point3{T})::Int64 where T<:AbstractFloat
99-
bot, top = ctube.planes
100-
tub = ctube.tube
101-
x, y, z = point
102-
103-
# Check the cut planes first
104-
a = safety(bot, point)
105-
b = safety(top, point)
106-
outside = a > kTolerance(T)/2 || b > kTolerance(T)/2
107-
outside && return kOutside
108-
cinside = a < -kTolerance(T)/2 && b < -kTolerance(T)/2
109-
inplanes = cinside ? kInside : kSurface
110-
# Check the tube
111-
intube = inside(tub, point)
112-
return inplanes == kSurface && intube != kOutside ? inplanes : intube
113-
end
114-
11598
function safetyToIn(ctube::CutTube{T}, point::Point3{T}) where T<:AbstractFloat
11699

117100
end

‎src/Geom4hep.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ include("GDML.jl")
3737
include("./DistanceIn.jl")
3838
include("./DistanceOut.jl")
3939
include("./Distance.jl")
40+
include("./Inside.jl")
4041
include("Benchmark.jl")
4142
#include("CuGeom.jl") # PackageCompiler fails?
4243

‎src/Inside.jl‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
function inside(shape, point)
2+
@nospecialize shape
3+
if shape isa Trap
4+
inside_trap(shape, point)
5+
elseif shape isa Trd
6+
inside_trd(shape, point)
7+
elseif shape isa Cone
8+
inside_cone(shape, point)
9+
elseif shape isa Box
10+
inside_box(shape, point)
11+
elseif shape isa Tube
12+
inside_tube(shape, point)
13+
elseif shape isa Polycone
14+
inside_polycone(shape, point)
15+
elseif shape isa CutTube
16+
inside_cuttube(shape, point)
17+
elseif shape isa BooleanUnion
18+
inside_booleanunion(shape, point)
19+
elseif shape isa BooleanSubtraction
20+
inside_booleansubtraction(shape, point)
21+
elseif shape isa BooleanIntersection
22+
inside_booleanintersection(shape, point)
23+
end
24+
end
25+
function inside_booleanunion(shape::BooleanUnion{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
26+
(; left, right, transformation) = shape
27+
lpoint = transformation * point
28+
29+
positionA = inside(left, point)
30+
positionA == kInside && return kInside
31+
32+
positionB = inside(right, lpoint)
33+
positionB == kInside && return kInside
34+
35+
if positionA == kSurface && positionB == kSurface
36+
normalA = normal(left, point)
37+
normalB = normal(right, lpoint) * transformation
38+
if dot(normalA, normalB) < 0
39+
return kInside # touching solids -)(-
40+
else
41+
return kSurface # overlapping solids =))
42+
end
43+
elseif positionA == kSurface || positionB == kSurface
44+
return kSurface
45+
else
46+
return kOutside
47+
end
48+
end
49+
50+
function inside_booleanintersection(shape::BooleanIntersection{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
51+
(; left, right, transformation) = shape
52+
lpoint = transformation * point
53+
54+
positionA = inside(left, point)
55+
positionA == kOutside && return kOutside
56+
57+
positionB = inside(right, lpoint)
58+
if positionA == kInside && positionB == kInside
59+
return kInside
60+
elseif positionA == kInside && positionB == kSurface ||
61+
positionA == kSurface && positionB == kInside ||
62+
positionA == kSurface && positionB == kSurface
63+
return kSurface
64+
else
65+
return kOutside
66+
end
67+
end
68+
69+
function inside_booleansubtraction(shape::BooleanSubtraction{T, SL, SR}, point::Point3{T})::Int64 where {T,SL,SR}
70+
(; left, right, transformation) = shape
71+
lpoint = transformation * point
72+
73+
positionA = inside(left, point)
74+
positionA == kOutside && return kOutside
75+
76+
positionB = inside(right, lpoint)
77+
78+
if positionA == kInside && positionB == kOutside
79+
return kInside;
80+
elseif positionA == kInside && positionB == kSurface ||
81+
positionB == kOutside && positionA == kSurface
82+
return kSurface
83+
else
84+
return kOutside
85+
end
86+
end
87+
88+
89+
function inside_trap(trap::Trap{T}, point::Point3{T}) where T<:AbstractFloat
90+
z = point[3]
91+
planes = trap.planes
92+
# inside z?
93+
outside = abs(z) > trap.z + kTolerance(T)/2
94+
cinside = abs(z) < trap.z - kTolerance(T)/2
95+
96+
for i in 1:4
97+
s = safety(planes[i], point) # positive if the point is on the outside halfspace, negative otherwise.
98+
outside |= s > kTolerance(T)/2
99+
cinside &= s < -kTolerance(T)/2
100+
end
101+
cinside ? kInside : outside ? kOutside : kSurface
102+
end
103+
function inside_cone(cone::Cone{T}, point::Point3{T}) where T<:AbstractFloat
104+
x, y, z = point
105+
(; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge,
106+
outerSlope, outerOffset, innerSlope, innerOffset) = cone
107+
108+
# Check Z
109+
outside = abs(z) > cone.z + kTolerance(T)/2
110+
outside && return kOutside
111+
cinside = abs(z) < cone.z - kTolerance(T)/2
112+
113+
# Check on RMax
114+
r2 = x * x + y * y
115+
rmax = rmax1 == rmax2 ? rmax1 : outerOffset + outerSlope * z
116+
outside |= r2 > rmax * rmax + kTolerance(T) * rmax
117+
outside && return kOutside
118+
cinside &= r2 < rmax * rmax - kTolerance(T) * rmax
119+
120+
# Check on RMin
121+
if rmin1 > 0. || rmin2 > 0.
122+
rmin = rmin1 == rmin2 ? rmin1 : innerOffset + innerSlope * z
123+
outside |= r2 <= rmin * rmin - kTolerance(T) * rmin
124+
outside && return kOutside
125+
cinside &= r2 > rmin * rmin + kTolerance(T) * rmin
126+
end
127+
128+
# Check on Phi
129+
if Δϕ < 2π
130+
outside |= isOutside(ϕWedge, x, y)
131+
cinside &= isInside(ϕWedge, x, y)
132+
end
133+
return outside ? kOutside : cinside ? kInside : kSurface
134+
end
135+
function inside_polycone(pcone::Polycone{T,N}, point::Point3{T}) where {T,N}
136+
x, y, z = point
137+
indexLow = getSectionIndex(pcone, z - kTolerance(T))
138+
indexHigh = getSectionIndex(pcone, z + kTolerance(T))
139+
indexLow < 0 && indexHigh < 0 && return kOutside
140+
indexLow < 0 && indexHigh == 1 && return inside(pcone.sections[1], point-Vector3{T}(0,0,pcone.zᵢ[1]))
141+
indexHigh < 0 && indexLow == N && return inside(pcone.sections[N], point-Vector3{T}(0,0,pcone.zᵢ[N]))
142+
if indexLow == indexHigh
143+
return inside(pcone.sections[indexLow], point-Vector3{T}(0,0,pcone.zᵢ[indexLow]))
144+
else
145+
insideLow = inside(pcone.sections[indexLow], point-Vector3{T}(0,0,pcone.zᵢ[indexLow]))
146+
insideHigh = inside(pcone.sections[indexHigh], point-Vector3{T}(0,0,pcone.zᵢ[indexHigh]))
147+
insideLow == kSurface && insideHigh == kSurface && return kInside
148+
insideLow == kOutside && insideHigh == kOutside && return kOutside
149+
return kSurface
150+
end
151+
end
152+
function inside_box(box::Box{T}, point::Point3{T}) where T<:AbstractFloat
153+
dist = -Inf
154+
for i in 1:3
155+
d = abs(point[i]) - box.fDimensions[i]
156+
if d > dist
157+
dist = d
158+
end
159+
end
160+
abs(dist) <= kTolerance(T)/2 ? kSurface : dist < 0.0 ? kInside : kOutside
161+
#dist = maximum(abs.(point) - box.fDimensions)
162+
#isapprox(dist, 0.0, atol = kTolerance(T)/2) ? kSurface : dist < 0.0 ? kInside : kOutside
163+
end
164+
function inside_trd(trd::Trd{T}, point::Point3{T}) where T<:AbstractFloat
165+
x, y, z = point
166+
167+
# inside z?
168+
outside = abs(z) > (trd.z + kTolerance(T)/2)
169+
inside = abs(z) < (trd.z - kTolerance(T)/2)
170+
171+
# inside x?
172+
c = cross(abs(x)-trd.x1, z+trd.z, trd.x2-trd.x1, 2*trd.z)
173+
outside |= (c < -kTolerance(T)/2)
174+
inside &= (c > kTolerance(T)/2)
175+
176+
# inside y?
177+
c = cross(abs(y)-trd.y1, z+trd.z, trd.y2-trd.y1, 2*trd.z)
178+
outside |= (c < -kTolerance(T)/2)
179+
inside &= (c > kTolerance(T)/2)
180+
181+
# return
182+
inside ? kInside : outside ? kOutside : kSurface
183+
end
184+
function inside_cuttube(ctube::CutTube{T}, point::Point3{T}) where T<:AbstractFloat
185+
bot, top = ctube.planes
186+
tub = ctube.tube
187+
x, y, z = point
188+
189+
# Check the cut planes first
190+
a = safety(bot, point)
191+
b = safety(top, point)
192+
outside = a > kTolerance(T)/2 || b > kTolerance(T)/2
193+
outside && return kOutside
194+
cinside = a < -kTolerance(T)/2 && b < -kTolerance(T)/2
195+
inplanes = cinside ? kInside : kSurface
196+
# Check the tube
197+
intube = inside(tub, point)
198+
return inplanes == kSurface && intube != kOutside ? inplanes : intube
199+
end
200+
function inside_tube(tub::Tube{T}, point::Point3{T}) where T<:AbstractFloat
201+
x, y, z = point
202+
203+
# Check Z
204+
outside = abs(z) > tub.z + kTolerance(T)/2
205+
outside && return kOutside
206+
cinside = abs(z) < tub.z - kTolerance(T)/2
207+
208+
# Check on RMax
209+
r2 = x * x + y * y
210+
outside |= r2 > tub.rmax2 + kTolerance(T) * tub.rmax
211+
outside && return kOutside
212+
cinside &= r2 < tub.rmax2 - kTolerance(T) * tub.rmax
213+
214+
# Check on RMin
215+
if tub.rmin > 0.
216+
outside |= r2 <= tub.rmin2 - kTolerance(T) * tub.rmin
217+
outside && return kOutside
218+
cinside &= r2 > tub.rmin2 + kTolerance(T) * tub.rmin
219+
end
220+
221+
# Check on Phi
222+
if tub.Δϕ < 2π
223+
outside |= isOutside(tub.ϕWedge, x, y)
224+
cinside &= isInside(tub.ϕWedge, x, y)
225+
end
226+
227+
return outside ? kOutside : cinside ? kInside : kSurface
228+
end

0 commit comments

Comments
 (0)