Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mesh Plotting #218

Merged
merged 15 commits into from
Oct 14, 2021
Merged

Mesh Plotting #218

merged 15 commits into from
Oct 14, 2021

Conversation

hervasa2
Copy link
Collaborator

Added mesh plotting for ConeMantle. Supported in pyplot() and gr(). Note that the only way to do this in gr() is using st = :mesh3d which is "experimental". Previously st = :surface was able to properly color the face of a triangle in gr(), now it needs 5 or more points to work, and no longer respects the boundaries set by the points. Therefore :mesh3d must be used. For consistency when surface() or equivalently plot(; st = :surface) is called, pyplot() uses :surface and gr() uses :mesh3d to produce a visually equivalent result. The experimental :mesh3d gives the warning (directly from GR, not from SSD) GR: mesh3d is experimental (no face colors) for each element in the series. As stated, there are no face colors yet, and therefore it will look like a dense wireframe for now.

To achieve the desired plot mesh(cm::ConeMantle{T}) is called. This will return a Mesh{T}. A plot recipe for Mesh is provided. Its fields are 3 matrices, which pyplot() :surface takes as input.
With, gr() this does not work. So we turn to :mesh3d. However, :mesh3d takes arrays as inputs: a list of points just like :path3d. The keyword connections is used to tell how these points assemble into triangles (whose faces will be filled - someday). The function polymesh is provided to turn a Mesh{T} into a PolyMesh{T,4} whose fields are arrays of 4 points. Mesh{T} naturally deconstructs into quadrangels instead of triangles and therefore PolyMesh{T,4} is used. In the future PolyMesh{T,3} will be useful for other meshes which deconstruct more naturally into triangles.

3 examples follow

T = Float64
cylindermantle = CSG.ConeMantle{T,Tuple{T,T},Nothing,:outwards}((1,1),nothing,1,CartesianPoint{T}(0,0,0),RotX(0))
conemantle_tranformed = CSG.ConeMantle{T,Tuple{T,T},Nothing,:outwards}((1,2),nothing,1,CartesianPoint{T}(1,0,0),RotY(π/2))
conemantle_partial = CSG.ConeMantle{T,Tuple{T,T},Tuple{T,T},:outwards}((3,2),(0,3π/2),1,CartesianPoint{T}(0,0,0),RotX(0))

pyplot(); surface(cylindermantle)
Screen Shot 2021-09-16 at 5 59 47 PM

gr(); surface(cylindermantle)
Screen Shot 2021-09-16 at 6 01 11 PM

pyplot(); surface(conemantle_tranformed)
Screen Shot 2021-09-16 at 6 00 30 PM

gr(); surface(conemantle_tranformed)
Screen Shot 2021-09-16 at 6 01 34 PM

pyplot(); surface(conemantle_partial)
Screen Shot 2021-09-16 at 6 08 08 PM
gr(); surface(conemantle_partial)
Screen Shot 2021-09-16 at 6 01 57 PM

Comment on lines 17 to 21
struct PolyMesh{T,N}
x::Vector{SVector{N,T}}
y::Vector{SVector{N,T}}
z::Vector{SVector{N,T}}
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just to clarify: PolyMesh is an auxiliary struct for plotting using mesh3d in GR?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For now this it is its only purpose. However, I will need PolyMesh{T,3} for when I plot surfaces for which I cannot write a analytical mesh. These arbitrary meshes will be generated by the JuliaGeometry/VoronoiDelaunay.jl package and then stored in a PolyMesh{T,3} (which contains vectors of triangle vertices. PolyMesh{T,4} contains vectors of quadrangle vertices). Note that pyplot() understands PolyMesh{T} aswell. So in theory we do not need to check for the backend and can always pass a PolyMesh. But we save time by not doing this in pyplot(). Another reason to use Mesh{T} is that pyplot() looks really nice with this structure.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This all is very annoying... I don't mean you implementation but this whole Plots.jl, GR/PyPlot stuff for 3D plots...

Meshes.jl actually seems to make really nice progress and provides a mesh struct we can use.

using Plots
using Meshes

T = Float64
p0 = Point3([0,0,0])
px = Point3([1,0,0])
py = Point3([0,1,0])
pz = Point3([0,0,1])

# mesh of three triangles, just as an example
mesh = SimpleMesh([p0, px, py, pz], [connect((1, 2, 3)), connect((1, 2, 4))])

plot(mesh)

This produces a wireframe in GR. With PyPlot it throws an error....
I would favor to use the mesh from Meshes.jl and, as they also provide already plot recipes for Plots.jl
(even though they plot usually with MeshViz.jl), talk to them regarding your workaround for plotting nice surfaces via PyPlot to their package. They suggest to contact them via https://julialang.zulipchat.com/#narrow/stream/275558-meshes.2Ejl before making issues on GitHub. Maybe we could also improve/contribute to Plots.jl directly.
And we might talk to Josef Heinen regarding face support from GR.jl.

I know this will increase the time until everything is through, but I think this is the better way and will reduce the amount of code in our package as we would only need to implement the connections for the vertices of our volume primitives to create the meshes.

end
end
if !haskey(plotattributes, :show_normal) || plotattributes[:show_normal]
if !haskey(plotattributes, :show_normal) || plotattributes[:show_normal]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is !haskey intended, i.e. default to plot it if not specified?
I would personally prefer haskey and only plot it if explicitly specified.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line is not my code, but I agree

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, we should set the default to false.

@fhagemann
Copy link
Collaborator

Will this PR introduce mesh plotting only for ConeMantle or are you planning to add mesh methods for all other primitives in this PR as well?

@hervasa2
Copy link
Collaborator Author

Will this PR introduce mesh plotting only for ConeMantle or are you planning to add mesh methods for all other primitives in this PR as well?

This PR is meant to add mesh methods for all surface primitives. However, I started it with only one, to see if you liked the way its implemented before I proceed

@hervasa2
Copy link
Collaborator Author

  • Meshing methods for all instances of all surface primitives. This includes 2D meshing for the arbitrary Polygon{N,T}. Support for gr() and pyplot()

  • in(p::Polygon{4,T}) is now generalized for to in(p::Polygon{N,T}).

The current status of plotting now is:

  • plot(s::SurfacePrimitive) relies on method lines(s). lines(s) is not available for all instances of a given SurfacePrimitive. eg: no method matching lines(::SolidStateDetectors.ConstructiveSolidGeometry.TorusMantle{Float64, Nothing, Tuple{Float64, Float64}, :outwards})

  • surface(s::SurfacePrimitive) relies on methods mesh(s) or trimesh(s). It should be complete for all instances of all SurfacePrimitives. Pending review.

  • scatter(s::SurfacePrimitive) will rely on method sample(s). Previously implemented as SSD_style = :samplesurface. Not started.

Some pyplot() examples follow. All of them have also been tested to work in gr().

Screen Shot 2021-09-19 at 3 06 45 PM
Screen Shot 2021-09-19 at 3 28 05 PM
Screen Shot 2021-09-19 at 3 02 21 PM
Screen Shot 2021-09-19 at 3 02 35 PM
Screen Shot 2021-09-19 at 3 03 07 PM

Copy link
Collaborator

@lmh91 lmh91 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we really should talk to the Plots.jl, GR.jl & PyPlot.jl people and make this right.
Even if this needs some time.
But there seems plenty of other people wanting this feature. So it would be a really nice contribution to the Julia Ecosystem. And it would reduce the amount of code we have to maintain.

Project.toml Outdated
@@ -30,6 +30,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
VoronoiDelaunay = "72f80fcb-8c52-57d9-aff0-40c1a3526986"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not see the reason why we would need this dependency.

It also does not seem to be really maintained. Last commit in January and the last release is almost one year old. Also the last issues did not received any comments.

src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl Outdated Show resolved Hide resolved
Comment on lines 17 to 21
struct PolyMesh{T,N}
x::Vector{SVector{N,T}}
y::Vector{SVector{N,T}}
z::Vector{SVector{N,T}}
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This all is very annoying... I don't mean you implementation but this whole Plots.jl, GR/PyPlot stuff for 3D plots...

Meshes.jl actually seems to make really nice progress and provides a mesh struct we can use.

using Plots
using Meshes

T = Float64
p0 = Point3([0,0,0])
px = Point3([1,0,0])
py = Point3([0,1,0])
pz = Point3([0,0,1])

# mesh of three triangles, just as an example
mesh = SimpleMesh([p0, px, py, pz], [connect((1, 2, 3)), connect((1, 2, 4))])

plot(mesh)

This produces a wireframe in GR. With PyPlot it throws an error....
I would favor to use the mesh from Meshes.jl and, as they also provide already plot recipes for Plots.jl
(even though they plot usually with MeshViz.jl), talk to them regarding your workaround for plotting nice surfaces via PyPlot to their package. They suggest to contact them via https://julialang.zulipchat.com/#narrow/stream/275558-meshes.2Ejl before making issues on GitHub. Maybe we could also improve/contribute to Plots.jl directly.
And we might talk to Josef Heinen regarding face support from GR.jl.

I know this will increase the time until everything is through, but I think this is the better way and will reduce the amount of code in our package as we would only need to implement the connections for the vertices of our volume primitives to create the meshes.

end
end
if !haskey(plotattributes, :show_normal) || plotattributes[:show_normal]
if !haskey(plotattributes, :show_normal) || plotattributes[:show_normal]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, we should set the default to false.

@lmh91
Copy link
Collaborator

lmh91 commented Sep 22, 2021

I added mesh plotting via the seriestype :mesh3d into Plots.jl: JuliaPlots/Plots.jl#3835

So, once Plots.jl makes a new release we can remove the special case for the PyPlot backend and do not have to write backend dependent code.
Locally, I just commented out those cases and was able to reproduce your plot pyplot(); surface(cylindermantle) with
pyplot(); surface(cylindermantle, linewidth = 0).

I am in contact with Josef Heinen (GR.jl) to get support of filled faces also in GR.
I will add this into Plots.jl as well as soon as it is in GR.jl.

@oschulz
Copy link
Member

oschulz commented Sep 23, 2021

Awesome!

@lmh91
Copy link
Collaborator

lmh91 commented Oct 7, 2021

Okay, so now it is (will be) also supported by GR :)

JuliaPlots/Plots.jl#3868

@lmh91
Copy link
Collaborator

lmh91 commented Oct 8, 2021

So (once the PR in Plots.jl is through) we only need one struct for a mesh:

I would suggest this:

struct mesh{T}
    x::Array{T}
    y::Array{T}
    z::Array{T}
    connections::Vector{Vector{T}} # We probably could also use `VectorOfVectors` from `ArraysOfArrays.jl`
end

so we can use a recipe like

@recipe function f(m::Mesh)
     seriestype := :mesh3d
     connections =: m.connections
     m.x, m.y, m.z
end

Regarding the structure of connections (They can be given in different ways, see):
It is a Vector{Vector{T}}. Each element (::Vector{T}) describes one Polygon. So we are not limited to triangles:

T = Float64
p0 = T[0,0,0]
p1 = T[1,0,0]
p2 = T[0,1,0]
p3 = T[1,1,0]
p4 = T[1/2,1/2,1]
pts = [p0, p1, p2, p3, p4]
x, y, z = broadcast(i -> getindex.(pts, i), (1,2,3))
connections = [ # Pyramide
     [1, 2, 4, 3], # Quadrangle 
     [1, 2, 5], # Triangle
     [2, 4, 5], # Triangle
     [4, 3, 5], # Triangle
     [3, 1, 5]  # Triangle
]

@hervasa2
Copy link
Collaborator Author

Backend independent meshing implementation using development branch of Plots (https://github.com/lmh91/Plots.jl/tree/gr_mesh3d) by @lmh91.

  • Now there is only one mesh struct used by both pyplot and gr.
  • Ugly rotations of vertical surfaces is no longer needed.
  • gr now provides face colors
  • The connections field of Mesh{T} contains the indices of the points which will form a polygon in the mesh. All meshes use quadrilaterals except Polygon which uses itself.
conemantle_tranformed = CSG.ConeMantle{T,Tuple{T,T},Nothing,:outwards}((1,2),nothing,1,CartesianPoint{T}(1,0,0),RotY/2))
es = CSG.EllipticalSurface{T,Tuple{T,T},Nothing}((1,2),nothing, zero(CartesianPoint{T}), RotX/6))

n = 6
polygon = CSG.Polygon{n,T}([RotX/4)*Pt(cos(i), sin(i),0) for i in 0:2π/n:2π-2π/n])
 
em = CSG.EllipsoidMantle{T,T,Tuple{T,T},Tuple{T,T},:ouwards}(1,(0,3π/2),(0/2), zero(CartesianPoint{T}), RotX(0))
 
tm = CSG.TorusMantle{T,Tuple{T,T},Tuple{T,T},:outwards}(2,1,(0,3π/2),(0,3π/2),CartesianPoint{T}(1,1,1), RotX/4))

Screen Shot 2021-10-11 at 2 46 11 PM
Screen Shot 2021-10-11 at 2 46 29 PM

Screen Shot 2021-10-11 at 2 46 55 PM
Screen Shot 2021-10-11 at 2 47 05 PM

Screen Shot 2021-10-11 at 3 36 29 PM
Screen Shot 2021-10-11 at 3 36 40 PM

Screen Shot 2021-10-11 at 6 24 50 PM
Screen Shot 2021-10-11 at 6 24 37 PM

Screen Shot 2021-10-11 at 6 31 32 PM
Screen Shot 2021-10-11 at 6 24 02 PM

Copy link
Collaborator

@lmh91 lmh91 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice! Just some minor things, see comments.

src/ConstructiveSolidGeometry/plotting/Meshing.jl Outdated Show resolved Hide resolved
src/ConstructiveSolidGeometry/plotting/Meshing.jl Outdated Show resolved Hide resolved
@lmh91
Copy link
Collaborator

lmh91 commented Oct 13, 2021

The PR in Plots.jl it through and the release v1.22.6 is on its way.
It should be available soon.

@fhagemann
Copy link
Collaborator

I tried this out but it didn't work (my Plots was not up to date).
We should somehow adapt the requirements such that the needed version of Plots is loaded.
Any ideas how to do this?

@lmh91
Copy link
Collaborator

lmh91 commented Oct 14, 2021

I tried this out but it didn't work (my Plots was not up to date). We should somehow adapt the requirements such that the needed version of Plots is loaded. Any ideas how to do this?

I asked in the Julia Plotting channel on Slack about this. Maybe there is a way.

@lmh91 lmh91 merged commit 418da02 into JuliaPhysics:v0.7 Oct 14, 2021
@fhagemann
Copy link
Collaborator

How do I plot the detector using meshes?

Comment on lines 15 to 25
function _get_n_points_in_arc_φ(p::AbstractSurfacePrimitive, n::Int64)::Int64
φMin, φMax = get_φ_limits(p)
f = (φMax - φMin)/(2π)
Int(ceil(n*f))
end

function _get_n_points_in_arc_θ(p::AbstractSurfacePrimitive, n::Int64)::Int64
θMin, θMax = get_θ_limits(p)
f = (θMax - θMin)/(2π)
Int(ceil(n*f))
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am getting an error that these functions are not defined when trying to plot a detector.
I guess they need to go before all the includes in this file.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems to be on your side. I was just able to plot a detector locally.
Also, the docs build went through where a detector is plotted.
I used:

using Plots, SolidStateDetectors
T = Float32 
sim = Simulation{T}(SSD_examples[:InvertedCoax])
plot(sim.detector)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants