Skip to content

Commit

Permalink
recover again
Browse files Browse the repository at this point in the history
  • Loading branch information
@Moelf
Moelf committed Sep 22, 2022
1 parent eafe81e commit 99f4396
Show file tree
Hide file tree
Showing 9 changed files with 12 additions and 492 deletions.
159 changes: 0 additions & 159 deletions src/Cone.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,80 +263,6 @@ function distanceToConicalSurface(cone::Cone{T}, point::Point3{T}, dir::Vector3{
ok ? distance : T(Inf)
end

function distanceToOut(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

function isInsideR(cone::Cone{T}, hit::Point3{T}) where T<:AbstractFloat
(; rmin1, rmax1, rmin2, rmax2, Δϕ, ϕWedge,
Expand All @@ -346,91 +272,6 @@ function isInsideR(cone::Cone{T}, hit::Point3{T}) where T<:AbstractFloat
rminTol = rmin1 == rmin2 ? rmin1 + kTolerance(T) : innerOffset + innerSlope * hit[3] + kTolerance(T)
r2 >= rminTol * rminTol && r2 <= rmaxTol * rmaxTol
end


function distanceToIn(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
#---Drawing functions------------------------------------------------------------
function GeometryBasics.coordinates(cone::Cone{T}, facets=36) where {T<:AbstractFloat}
(; rmin1, rmax1, rmin2, rmax2, z, ϕ₀, Δϕ) = cone
Expand Down
49 changes: 0 additions & 49 deletions src/CutTube.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,52 +99,3 @@ function safetyToIn(ctube::CutTube{T}, point::Point3{T}) where T<:AbstractFloat

end

function distanceToOut(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

function distanceToIn(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

3 changes: 2 additions & 1 deletion src/Distance.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions src/DistanceIn.jl
View file
Original file line number Diff line number Diff line change
@@ -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}
Expand Down
79 changes: 0 additions & 79 deletions src/Polycone.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,24 +83,6 @@ 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}
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

function safetyToIn(pcone::Polycone{T,N}, point::Point3{T})::T where {T,N}
(; sections, zᵢ) = pcone
z = point[3]
Expand Down Expand Up @@ -164,64 +146,3 @@ function safetyToOut(pcone::Polycone{T,N}, point::Point3{T})::T where {T,N}
end
return minSafety
end

function distanceToOut(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
16 changes: 8 additions & 8 deletions src/Transformation3D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,14 @@ Transformation3D{T}(dx, dy, dz, rotx, roty, rotz) where T<:AbstractFloat = Trans
Transformation3D{T}(dx, dy, dz, rot::Rotation{3,T}) where T<:AbstractFloat = Transformation3D(RotMatrix3{T}(rot), Vector3{T}(dx, dy, dz))

# Transforms
@inline transform(t::Transformation3D{T}, p::Point3{T}) where T<:AbstractFloat = t.rotation * (p - t.translation)
@inline transform(t::Transformation3D{T}, d::Vector3{T}) where T<:AbstractFloat = t.rotation * d
@inline invtransform(t::Transformation3D{T}, p::Point3{T}) where T<:AbstractFloat = t.translation + (p' * t.rotation)'
@inline invtransform(t::Transformation3D{T}, d::Vector3{T}) where T<:AbstractFloat = (d' * t.rotation)'
@inline Base.:*(t::Transformation3D{T}, p::Point3{T}) where T<:AbstractFloat = transform(t,p)
@inline Base.:*(t::Transformation3D{T}, d::Vector3{T}) where T<:AbstractFloat = transform(t,d)
@inline Base.:*(p::Point3{T}, t::Transformation3D{T}) where T<:AbstractFloat = invtransform(t,p)
@inline Base.:*(d::Vector3{T}, t::Transformation3D{T}) where T<:AbstractFloat = invtransform(t,d)
transform(t::Transformation3D{T}, p::Point3{T}) where T<:AbstractFloat = t.rotation * (p - t.translation)
transform(t::Transformation3D{T}, d::Vector3{T}) where T<:AbstractFloat = t.rotation * d
invtransform(t::Transformation3D{T}, p::Point3{T}) where T<:AbstractFloat = t.translation + (p' * t.rotation)'
invtransform(t::Transformation3D{T}, d::Vector3{T}) where T<:AbstractFloat = (d' * t.rotation)'
Base.:*(t::Transformation3D{T}, p::Point3{T}) where T<:AbstractFloat = transform(t,p)
Base.:*(t::Transformation3D{T}, d::Vector3{T}) where T<:AbstractFloat = transform(t,d)
Base.:*(p::Point3{T}, t::Transformation3D{T}) where T<:AbstractFloat = invtransform(t,p)
Base.:*(d::Vector3{T}, t::Transformation3D{T}) where T<:AbstractFloat = invtransform(t,d)

function transform(v::Vector{Transformation3D{T}}, p::Point3{T}) where T<:AbstractFloat
for t in v
Expand Down
Loading

0 comments on commit 99f4396

Please sign in to comment.