Skip to content
35 changes: 20 additions & 15 deletions src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
Original file line number Diff line number Diff line change
Expand Up @@ -125,10 +125,8 @@ end
@inline Base.adjoint(pt::CartesianPoint) = CartesianPoint(adjoint(pt.x), adjoint(pt.y), adjoint(pt.z))

# Barycentric combination
function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}; dims = :) where T
# Equivalent to A .- Ref(cartesian_zero), but allocation-free:
B = reinterpret(CartesianVector{T}, A) # ToDo: Avoid this memory allocation.
return cartesian_zero + mean(B; dims = dims)
function Statistics.mean(A::AbstractArray{<:CartesianPoint}; dims = :)
barycenter(A; dims)
end

function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}, weights::StatsBase.AbstractWeights; dims = :) where T
Expand All @@ -137,6 +135,12 @@ function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}, weights::StatsBa
return cartesian_zero + mean(B, weights; dims = dims)
end

barycenter(X; kwargs...) = _barycenter_impl(eltype(X), X; kwargs...)
function _barycenter_impl(::Type{T}, points; kwargs...) where {T<:AbstractCoordinatePoint}
mean_vector = mean(Base.Fix2(-, cartesian_zero), CartesianPoint.(points); kwargs...)
return convert(T, cartesian_zero + mean_vector)::T
end

# Base.broadcastable(p::CartesianPoint) = (p.x, p.y, p.z)


Expand Down Expand Up @@ -275,6 +279,14 @@ function Base.convert(::Type{CylindricalPoint{T}}, pt::CylindricalPoint{U}) wher
return CylindricalPoint{T}(convert(T, pt.r), convert(T, pt.φ), convert(T, pt.z))
end

function Base.convert(::Type{CylindricalPoint{T}}, pt::AbstractCoordinatePoint)::CylindricalPoint{T} where {T}
return CylindricalPoint(pt)
end

function Base.convert(::Type{CartesianPoint{T}}, pt::AbstractCoordinatePoint)::CartesianPoint{T} where {T}
return CartesianPoint(pt)
end

AbstractCoordinatePoint{T, Cylindrical}(r::Real, φ::Real, z::Real) where T = CylindricalPoint{T}(r, φ, z)

Base.:(==)(a::CylindricalPoint, b::CylindricalPoint) = a.r == b.r && a.φ == b.φ && a.z == b.z
Expand All @@ -298,19 +310,12 @@ Base.iszero(pt::CylindricalPoint) = iszero(pt.r) && iszero(pt.z)
@inline Base.adjoint(pt::CylindricalPoint) = CylindricalPoint(adjoint(pt.r), adjoint(pt.φ), adjoint(pt.z))

# Barycentric combination
function Statistics.mean(A::AbstractArray{<:CylindricalPoint}; dims = :)
CylindricalPoint(mean(CartesianPoint, A; dims = dims))
function Statistics.mean(A::AbstractArray{<:CylindricalPoint}; dims = :)
B = CartesianPoint.(A)
CylindricalPoint(mean(B; dims = dims))
end

function Statistics.mean(A::AbstractArray{<:CylindricalPoint}, weights::StatsBase.AbstractWeights; dims = :)
B = CartesianPoint.(A)
CylindricalPoint(mean(B, weights; dims = dims))
end


# function _Δφ(φ1::T, φ2::T)::T where {T}
# δφ = mod(φ2 - φ1, T(2π))
# min(δφ, T(2π) - δφ)
# end
#
# _φNear(φ::Real, φMin::T, φMax::T) where {T} = _Δφ(T(φ),φMin) ≤ _Δφ(T(φ),φMax) ? φMin : φMax
end
4 changes: 2 additions & 2 deletions src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ function _sample_excluding_border(t::Triangle{T}, spacing::T)::Vector{CartesianP
p2 = t.points[2] + (push/sφ)*(-u/nu + w/nw)
if (p2-p1) ⋅ u > 0
su = norm(p2-p1)/nu
[(su*a*u + su*b*v) + p1
[p1 + (su*a*u + su*b*v)
for a in range(0, stop = 1, length = max(2, 1 + Int(ceil(su*nu*sθ/spacing))))
for b in range(0, stop = 1 - a, length = max(2, 1 + Int(ceil(su*nv*(1 - a)/spacing))))]
else
Expand All @@ -55,7 +55,7 @@ function sample(p::Polygon{N,T}, spacing::T)::Vector{CartesianPoint{T}} where {N
end

function extremum(p::Polygon{N,T})::T where {N,T}
c = sum(p.points)/N
c = barycenter(p.points)
m = maximum([norm(point - c) for point in p.points])
end

Expand Down
6 changes: 3 additions & 3 deletions src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl
Original file line number Diff line number Diff line change
Expand Up @@ -145,12 +145,12 @@ function vertices(b::Box{T}) where {T}
return _transform_into_global_coordinate_system.((
CartesianPoint{T}(-b.hX, -b.hY, -b.hZ),
CartesianPoint{T}(+b.hX, -b.hY, -b.hZ),
CartesianPoint{T}(-b.hX, +b.hY, -b.hZ),
CartesianPoint{T}(+b.hX, +b.hY, -b.hZ),
CartesianPoint{T}(-b.hX, +b.hY, -b.hZ),
CartesianPoint{T}(-b.hX, -b.hY, +b.hZ),
CartesianPoint{T}(+b.hX, -b.hY, +b.hZ),
CartesianPoint{T}(-b.hX, +b.hY, +b.hZ),
CartesianPoint{T}(+b.hX, +b.hY, +b.hZ)
CartesianPoint{T}(+b.hX, +b.hY, +b.hZ),
CartesianPoint{T}(-b.hX, +b.hY, +b.hZ)
), Ref(b))
end

Expand Down
9 changes: 9 additions & 0 deletions src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,15 @@ end

@recipe function f(csg::AbstractConstructiveGeometry{T}, primitive::AbstractPrimitive{T}, axis::Symbol, slice::T, st::Symbol; n_samples = 40, CSG_scale = missing, projection = :none, full_det = false, linewidth = :auto) where {T}
if st == :samplesurface
label --> l
seriesalpha --> 0.2
seriescolor --> 1
markerstrokewidth --> 0
markersize --> 3
if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
aspect_ratio --> 1.0
end
seriestype := :scatter
CSG_scale = ismissing(CSG_scale) ? get_scale(csg) : CSG_scale
spacing = T(CSG_scale/n_samples)
pts_nounits = filter(p -> in(p,csg), sample(primitive, spacing))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,14 @@ end
mesh(s, n_arc, n_vert_lines)
elseif st == :samplesurface
label --> l
seriesalpha --> 0.4
seriesalpha --> 0.2
seriescolor --> 1
markerstrokewidth --> 0
markersize --> 3
if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
aspect_ratio --> 1.0
end
seriestype := :scatter
sample(s, extremum(s)/n_samples)
elseif st == :slice
@warn ":slice is not supported for AbstractSurfacePrimitive. Use :csg, :wireframe, :mesh3d, or :samplesurface."
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,15 @@ end

@recipe function f(p::AbstractPrimitive{T}, axis::Symbol, slice::T, st::Symbol; n_samples = 40, CSG_scale = missing, projection = :none, full_det = false, linewidth = :auto) where {T}
if st == :samplesurface
label --> l
seriesalpha --> 0.2
seriescolor --> 1
markerstrokewidth --> 0
markersize --> 3
if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
aspect_ratio --> 1.0
end
seriestype := :scatter
CSG_scale = ismissing(CSG_scale) ? extremum(p) : CSG_scale
spacing = T(CSG_scale/n_samples)
sample(p, spacing)
Expand Down
20 changes: 19 additions & 1 deletion test/ConstructiveSolidGeometry/CSG_coordinates.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,11 @@
using Test

using SolidStateDetectors.ConstructiveSolidGeometry: CartesianPoint, CartesianVector,
CartesianZero, cartesian_zero, CylindricalPoint, LocalAffineFrame, global_frame, frame_transformation
CartesianZero, cartesian_zero, CylindricalPoint, LocalAffineFrame, global_frame, frame_transformation, barycenter
using StaticArrays: Size, SVector, SMatrix
using InverseFunctions: inverse
using Statistics: mean

import Unitful

@testset "points_and_vectors" begin
Expand Down Expand Up @@ -44,6 +46,14 @@ import Unitful


@test @inferred(CartesianZero{Float32}() * Unitful.u"mm") === CartesianZero{typeof(zero(Float32) * Unitful.u"mm")}()

A = [CartesianPoint{Float32}(x,0,0) for x in -2:2]
@test isapprox(mean(A), CartesianPoint{Float32}(0,0,0))
@test isapprox(barycenter(A), CartesianPoint{Float32}(0,0,0))

S = SVector{length(A)}(A)
@test_broken isapprox(mean(S), CartesianPoint{Float32}(0,0,0))
@test isapprox(barycenter(S), CartesianPoint{Float32}(0,0,0))
# ToDo: Add more tests!
end

Expand All @@ -68,5 +78,13 @@ import Unitful
@test @inferred(CylindricalPoint(a[1], a[2], a[3])) === a
@test @inferred(CylindricalPoint(a[1], a[2], a[3])) == a
@test @inferred(CylindricalPoint(a[1], a[2], a[3])) ≈ a

A = [CylindricalPoint{Float32}(x,0,0) for x in -2:2]
@test isapprox(mean(A), CylindricalPoint{Float32}(0,0,0))
@test isapprox(barycenter(A), CylindricalPoint{Float32}(0,0,0))

S = SVector{length(A)}(A)
@test_broken isapprox(mean(S), CylindricalPoint{Float32}(0,0,0))
@test isapprox(barycenter(S), CylindricalPoint{Float32}(0,0,0))
end
end
Loading