Add Polycone

Project.toml

@@ -63,7 +63,7 @@ LinearAlgebra = "<0.0.1, 1"
 OrderedCollections = "1.1"
 ParallelProcessingTools = "0.4"
 Parameters = "0.12"
-PolygonOps = "0.1"
+PolygonOps = "0.1.1"
 Polynomials = "2, 3, 4"
 ProgressMeter = "1.5"
 RadiationDetectorSignals = "0.3.5"

docs/src/man/csg.md

@@ -73,6 +73,22 @@ cone = CSG.Geometry(T, cfn)
 plot(cone)
 ````
 
+#### Polycone
+
+
+````@example primitives
+cfn = joinpath(path_to_example_primitives_config_files, "Polycone.yaml")
+print(open(f -> read(f, String), cfn))
+````
+
+Load the primitive from the configuration file via `CSG.Geometry`
+
+````@example primitives
+polycone = CSG.Geometry(T, cfn)
+plot(polycone)
+````
+
+
 ### Ellipsoid
 #### Sphere
 

examples/example_config_files/ivc_splitted_config/detector_geometry.yaml

@@ -1,26 +1,7 @@
 geometry:
-  difference:
-    - tube:
-        r: 35
-        h: 80
-        origin:
-          z: 40
-    - cone:
-        r:
-          bottom:
-            from: 35
-            to: 36
-          top:
-            from: 23.71
-            to: 36
-        h: 64
-        origin:
-          z: 52
-    - tube:
-        r: 5
-        h: 80
-        origin:
-          z: 65
+  polycone:
+    r: [0, 35, 35, 24.42, 5, 5, 0]
+    z: [0, 0, 20, 80, 80, 25, 25]
 contacts:
 - material: HPGe
   id: 1

examples/example_primitive_files/Polycone.yaml

@@ -0,0 +1,3 @@
+polycone:
+  r: [0, 35, 35, 25, 5, 5, 0, 0]  
+  z: [0, 0, 20, 80, 80, 25, 25, 0]

ext/Geant4/io_gdml.jl

@@ -1,7 +1,7 @@
 using SolidStateDetectors
 using SolidStateDetectors: Cylindrical, Cartesian, World, Simulation, SSDFloat
 using SolidStateDetectors.ConstructiveSolidGeometry: CSGUnion, CSGDifference, CSGIntersection, 
-    Box, Cone, Ellipsoid, Torus, RegularPrism, show_CSG_tree, AbstractVolumePrimitive, AbstractGeometry,
+    Box, Cone, Ellipsoid, Torus, RegularPrism, Polycone, AbstractVolumePrimitive, AbstractGeometry,
     AbstractConstructiveGeometry, CylindricalPoint, origin, rotation
 using IntervalSets
 using OrderedCollections: OrderedDict
@@ -171,6 +171,14 @@ function has_volume(e::RegularPrism, v::Bool = false)
     return true
 end
 
+function has_volume(p::Polycone, v::Bool = false)
+    if isapprox(PolygonOps.area(tuple.(p.r, p.z)), 0)
+        v && @warn "Polycone: The points passed result in a zero-volume Polycone"
+        return false
+    end
+    return true
+end
+
 
 # Check whether Boolean solid has any volume at all
 function has_volume(e::AbstractConstructiveGeometry, v::Bool = false)
@@ -374,6 +382,28 @@ function parse_geometry(e::RegularPrism{T, <:Any, N}, x_solids::XMLElement, x_de
     end
 end
 
+function parse_geometry(p::Polycone, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing
+    if has_volume(p, v)
+        y = new_child(x_solids, "genericPolycone")
+
+        for i in eachindex(p.r)
+            rzpoint = new_child(y, "rzpoint")
+            set_attributes(rzpoint, OrderedDict(
+                "r" => p.r[i],
+                "z" => p.z[i]
+            ))
+        end
+
+        set_attributes(y, OrderedDict(
+            "name" => pf * string(id),
+            "startphi" => 0,
+            "deltaphi" => 360,
+            "lunit" => SolidStateDetectors.internal_length_unit,
+            "aunit" => "deg"
+        ))
+    end
+end
+
 
 function parse_geometry(e::World{T, 3, Cylindrical}, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing where {T}
     w = new_child(x_solids, "tube")

ext/SolidStateDetectorsGeant4Ext.jl

@@ -8,6 +8,7 @@ using Geant4.SystemOfUnits
 using LightXML
 using LinearAlgebra
 using Parameters
+using PolygonOps
 using ProgressMeter
 using RadiationDetectorSignals
 using StaticArrays

src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl

@@ -50,9 +50,9 @@ module ConstructiveSolidGeometry
 
     abstract type AbstractConstructiveGeometry{T} <: AbstractGeometry{T} end
 
-    _csg_convert_args(eltype::Type{T}, r::Real) where T = convert(T, r) 
-    _csg_convert_args(eltype::Type{T}, r::Tuple) where T = broadcast(x -> _csg_convert_args(T, x), r)
-    _csg_convert_args(eltype::Type{T}, r::Nothing) where T = nothing
+    _csg_convert_args(::Type{T}, r::Real) where T = convert(T, r) 
+    _csg_convert_args(::Type{T}, r::Tuple) where T = broadcast(x -> _csg_convert_args(T, x), r)
+    _csg_convert_args(::Type{T}, r::Nothing) where T = nothing
 
     _csg_get_promoted_eltype(::Type{T}) where {T <: AbstractArray} = eltype(T)
     _csg_get_promoted_eltype(::Type{T}) where {T <: Real} = T

src/ConstructiveSolidGeometry/IO.jl

@@ -10,6 +10,7 @@ const CSG_dict = Dict{String, Any}(
     "sphere" => Ellipsoid,
     "box" => Box,
     "torus" => Torus,
+    "polycone" => Polycone,
     "TriangularPrism" => TriangularPrism,
     "QuadranglePrism" => QuadranglePrism,
     "PentagonalPrism" => PentagonalPrism,
@@ -59,7 +60,7 @@ end
     T(to_internal_units(float(x)))
 end
 @inline _parse_value(::Type{T}, s::String, unit::Unitful.Units) where {T} = _parse_value(T, uparse(s), unit)
-@inline _parse_value(::Type{T}, a::AbstractVector, unit::Unitful.Units) where {T} = _parse_value.(T, a, unit)
+@inline _parse_value(::Type{T}, a::Union{<:AbstractVector, <:Tuple}, unit::Unitful.Units) where {T} = _parse_value.(T, a, unit)
 
 # parses dictionary entries of type {"from": ..., "to": ... } to a Tuple of the interval boundaries
 function _parse_interval_from_to(::Type{T}, dict::AbstractDict, unit::Unitful.Units)::Tuple{T,T} where {T}

src/ConstructiveSolidGeometry/VolumePrimitives/Polycone.jl

@@ -0,0 +1,184 @@
+"""
+    struct Polyone{T,CO,N,TP} <: AbstractVolumePrimitive{T, CO}
+
+Volume primitive describing a [Polycone](@ref), similar to the G4Polycone defined in Geant4.
+
+## Parametric types
+* `T`: Precision type.
+* `CO`: Describes whether the surface belongs to the primitive. 
+    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
+    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
+* `N`: Integer describing the number of corners.
+* `TP`: Type of the angular range `φ`.
+    * `TP == Nothing`: Full 2π Cone.
+
+## Fields
+* `r::NTuple{N,T}`: `r`-coordinates of the corners of the polycone.
+* `z::NTuple{N,T}`: `z`-coordinates of the corners of the polycone.
+* `φ::TP`: Range in polar angle `φ` over which the `Cone` extends (in radians).
+* `origin::CartesianPoint{T}`: The position of the center of the `Cone` at the middle height.
+* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Cone` around its `origin`.
+
+## Definition in Configuration File
+
+A `Polycone` is defined in the configuration file as part of the `geometry` field 
+of an object through the field `polycone`.
+
+Example definitions of a polycone looks like this:
+```yaml
+polycone:
+  r: [0, 35, 35, 24.42, 5, 5, 0, 0]
+  z: [0, 0, 20, 80, 80, 25, 25, 0]
+  origin:
+    z: 1.0   # => origin = [0.0, 0.0, 1.0]
+```
+This is a polycone describing the semiconductor of the example inverted coaxial detector.
+
+See also [Constructive Solid Geometry (CSG)](@ref).
+"""
+struct Polycone{T,CO,N,TP<:Nothing} <: AbstractVolumePrimitive{T,CO}
+    r::NTuple{N,T}
+    z::NTuple{N,T}
+    φ::TP
+
+    origin::CartesianPoint{T}
+    rotation::SMatrix{3,3,T,9}
+
+    function Polycone{T,CO}(r::Union{<:AbstractArray, <:Tuple}, z::Union{<:AbstractArray, <:Tuple}, φ, origin, rotation) where {T, CO}
+        nr = length(r)
+        nz = length(z)
+        nr != nz && throw(ConfigFileError("In PolyCone: r and z must have the same length."))
+        # convert to Tuple and determine the length N
+        _r, _z = if first(r) != last(r) || first(z) != last(z)
+            nr += 1
+            nz += 1
+            (r..., first(r)), (z..., first(z))
+        else
+            (r...,), (z...,)
+        end
+        # sort the points counter-clockwise in the r-z-plane
+        if PolygonOps.area(tuple.(_r,_z)) < 0
+            _r, _z = reverse(_r), reverse!(_z)
+        end
+        new{T,CO,nr,typeof(φ)}(_r, _z, φ, origin, rotation)
+    end
+end
+
+#Type conversion happens here
+function Polycone{T}(CO, r, z, φ, origin, rotation) where {T}
+    _r = _csg_convert_args.(T, r)
+    _z = _csg_convert_args.(T, z)
+    _φ = _csg_convert_args(T, φ)
+    Polycone{T,CO}(_r, _z, _φ, origin, rotation)
+end
+
+#Type promotion happens here
+function Polycone(CO, r::TR, z::TZ, φ::TP, origin::PT, rotation::ROT) where {TR, TZ, TP, PT, ROT}
+    eltypes = _csg_get_promoted_eltype.((TR, TZ, TP, PT, ROT))
+    T = float(promote_type(eltypes...))
+    Polycone{T}(CO, r, z, φ, origin, rotation)
+end
+
+function Polycone(::Type{CO}=ClosedPrimitive;
+    r = (0,1,1,0),
+    z = (0,0,1,1),
+    φ = nothing,
+    origin = zero(CartesianPoint{Int}), 
+    rotation = one(SMatrix{3, 3, Int, 9})
+) where {CO}
+    Polycone(CO, r, z, φ, origin, rotation)
+end
+
+function Polycone{T}(::Type{CO}=ClosedPrimitive;
+    r = (0.,1.,1.,0.),
+    z = (0.,0.,1.,1.),
+    φ = nothing,
+    origin = zero(CartesianPoint{Float64}), 
+    rotation = one(SMatrix{3, 3, Float64, 9})
+) where {T, CO}
+    Polycone{T}(CO, r, z, φ, origin, rotation)
+end
+
+Polycone{T,CO,N,TP}( c::Polycone{T,CO,N,TP}; COT = CO,
+            origin::CartesianPoint{T} = c.origin,
+            rotation::SMatrix{3,3,T,9} = c.rotation) where {T, CO<:Union{ClosedPrimitive, OpenPrimitive}, N, TP} =
+    Polycone{T,COT}(c.r, c.z, c.φ, origin, rotation)
+
+
+####################################################################
+####################################################################
+
+
+function _inpolygon(pt::Tuple{T,T}, polygon::NTuple{N,Tuple{T,T}}; csgtol::T = csg_default_tol(T)) where {N,T <: Real}
+    # use the normal PolygonOps.inpolygon method to check if the point is inside
+    PolygonOps.inpolygon(pt, polygon, in = true, on = true, out = false) && return true
+    # if it is flagged as outside, check if the distance to any of the edges is smaller than csgtol
+    rp,zp = pt
+    @inbounds for i in Base.OneTo(N-1)
+        r1,z1 = polygon[i]
+	r2,z2 = polygon[i+1]
+        r1 == r2 && z1 == z2 && continue
+        s = clamp(((zp - z1) * (z2 - z1) + (rp - r1) * (r2 - r1)) / ((z2 - z1)^2 + (r2 - r1)^2), zero(T), one(T))
+        if hypot(rp - r1 - s * (r2 - r1), zp - z1 - s * (z2 - z1)) <= csgtol
+            return true
+        end
+    end
+    # if not, then it is really outside
+    return false
+end
+
+function _in(pt::CartesianPoint{T}, c::Polycone{<:Any, ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
+    return (isnothing(c.φ) || _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)) && 
+        _inpolygon((hypot(pt.x, pt.y), pt.z), tuple.(c.r, c.z); csgtol) 
+end
+
+function _in(pt::CartesianPoint{T}, c::Polycone{<:Any, OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
+    return (isnothing(c.φ) || _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)) && 
+        PolygonOps.inpolygon((hypot(pt.x, pt.y), pt.z), tuple.(c.r, c.z), in = true, on = iszero(pt.x) && iszero(pt.y), out = false)
+end
+
+function surfaces(c::Polycone{T,<:Any,N,Nothing}) where {T,N}
+    s = []
+    @inbounds for i in Base.OneTo(N-1)
+        r1::T = c.r[i]
+        r2::T = c.r[i+1]
+        ## skip the surfaces on the z-axis
+        iszero(r1) && iszero(r2) && continue
+        z1::T = c.z[i]
+        z2::T = c.z[i+1]
+        origin = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), (z1+z2)/2), c) 
+        vol = if z1 == z2
+            EllipticalSurface{T}(r = r1 <= r2 ? (r1,r2) : (r2,r1), origin = origin, rotation = r1 <= r2 ? c.rotation : -c.rotation * RotZ{T}(π))
+        else
+            FullConeMantle{T, z1 < z2 ? (:inwards) : (:outwards)}(z1 < z2 ? (r1,r2) : (r2,r1), c.φ, abs(z1-z2)/2, origin, c.rotation)
+        end
+        push!(s, vol)
+    end
+    (s...,)
+end
+
+
+function Geometry(::Type{T}, ::Type{Polycone}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
+    length_unit = input_units.length
+    origin = get_origin(T, dict, length_unit)
+    rotation = get_rotation(T, dict, input_units.angle)
+    haskey(dict, "r") || throw(ConfigFileError("Polycone needs entry `r`."))
+    haskey(dict, "z") || throw(ConfigFileError("Polycone needs entry `z`."))
+    r = _parse_value(T, dict["r"], length_unit)
+    z = _parse_value(T, dict["z"], length_unit)
+    polycone = Polycone{T, ClosedPrimitive}(r, z, nothing, origin, rotation)
+    transform(polycone, transformations)
+end
+
+
+function Dictionary(c::Polycone{T})::OrderedDict{String, Any} where {T}
+    dict = OrderedDict{String, Any}()
+    @assert isnothing(c.φ) "Polycone needs `φ` field to be nothing."
+    dict["r"] = c.r
+    dict["z"] = c.z
+    if c.origin != zero(CartesianVector{T}) dict["origin"] = c.origin end
+    if c.rotation != one(SMatrix{3,3,T,9})  dict["rotation"] = Dictionary(c.rotation) end 
+    OrderedDict{String, Any}("polycone" => dict)
+end
+
+extremum(c::Polycone{T}) where {T} = hypot(max(abs.(extrema(c.r))...), max(abs.(extrema(c.z))...))

src/ConstructiveSolidGeometry/VolumePrimitives/VolumePrimitives.jl

@@ -1,3 +1,9 @@
+# All volume primitives need to have the following methods defined: 
+# * _in(::CartesianPoint, ::VolumePrimitive)  (in the local coordinate system)
+# * surfaces(::VolumePrimitive) --> all surfaces of the volume primitive as Tuple
+# * Geometry(...) --> a method to read the geometry from a config file
+# * Dictionary(...) --> a method to convert the primitive to a config 
+
 ClosedPrimitive(p::VP) where {VP <:AbstractVolumePrimitive} = VP(p, COT = ClosedPrimitive)
 OpenPrimitive(p::VP) where {VP <: AbstractVolumePrimitive} = VP(p, COT = OpenPrimitive)
 switchClosedOpen(p::AbstractVolumePrimitive{<:Any,ClosedPrimitive}) = OpenPrimitive(p)
@@ -19,4 +25,5 @@ include("Box.jl")
 include("Cone.jl")
 include("Ellipsoid.jl")
 include("RegularPrism.jl")
+include("Polycone.jl")
 include("Torus.jl")

test/ConstructiveSolidGeometry/CSG_IO.jl

@@ -30,6 +30,11 @@ example_primitive_dir = joinpath(@__DIR__, "../../examples/example_primitive_fil
         @test ellipsoid_full_sphere isa CSG.FullSphere{T}
     end
 
+    @testset "Polycone" begin
+        polycone = Geometry(T, joinpath(example_primitive_dir, "Polycone.yaml"))
+        @test polycone isa CSG.Polycone{T}
+    end
+
     @testset "RegularPrism" begin
         hexagon = Geometry(T, joinpath(example_primitive_dir, "RegularPrism_hexagon.yaml"))
         @test hexagon isa CSG.HexagonalPrism{T}

test/ConstructiveSolidGeometry/CSG_primitives.jl

@@ -164,6 +164,22 @@ no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(Carte
         rot_torus = @inferred CSG.Torus{Float64}(φ=π/2,θ=nothing,rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0) * SMatrix{3}(cos(π/4),sin(π/4),0,-sin(π/4),cos(π/4),0,0,0,1))
         @test tuple_torus ==rot_torus
     end    
+    @testset "Polycone" begin
+        polycone1 = CSG.Polycone(CSG.ClosedPrimitive,r=Float32[0,2,4,0,0],z=Float32[0,1,2,3,0],origin=zero(CartesianPoint{Float16}),rotation=one(SMatrix{3, 3, Float16, 9}))
+        polycone2 = CSG.Polycone{Float32}(r = [0,2,4,0], z = [0,1,2,3]) # omit last entry
+        @test polycone1 == polycone2
+
+        dict = Dict("polycone" => Dict(
+            "r" => (0,2,4,0,0),
+            "z" => (0,1,2,3,0),
+        ))
+        polycone3 = Geometry(T,dict,default_units,no_translations)
+        @test dict == Dictionary(polycone3)
+
+        polycone_open = CSG.Polycone(CSG.OpenPrimitive,r=Float32[0,2,4,0,0],z=Float32[0,1,2,3,0],origin=zero(CartesianPoint{Float16}),rotation=one(SMatrix{3, 3, Float16, 9}))
+        @test  in(CartesianPoint{Float32}(4,0,2), polycone1)
+        @test !in(CartesianPoint{Float32}(4,0,2), polycone_open)
+    end
     @testset "Ellipsoid" begin
         ellip1 = @inferred CSG.Ellipsoid(CSG.ClosedPrimitive,r=1f0,origin = zero(CartesianPoint{Float16}),rotation = one(SMatrix{3, 3, Float16, 9}))
         ellip2 = @inferred CSG.Ellipsoid{Float32}(r=1.0)