Pushing the recipe system with advanced mesh plots

[juliohm]

According to the documentation, we could simply add a method to convert_arguments to get a plot for a custom type:

using Meshes, AbstractPlotting

import AbstractPlotting: convert_arguments

# convert a collection of points to a vector of SVector
convert_arguments(pset::PointSet) = ([coordinates(pset[i]) for i in nelements(pset)],)

When I try to plot however, I get an error:

using Makie

pset = PointSet(rand(Meshes.Point2, 100))

plot(pset)

ERROR: PlotMethodError: no plot method for arguments (::PointSet{2, Float64}). To support these arguments, define
  plot!(::Combined{Any, S} where S<:Tuple{PointSet{2, Float64}})
Available methods are:
  plot!(plot::Combined{Any, var"#s199"} where var"#s199"<:Tuple{AbstractVector{var"#s199"} where var"#s199"<:Complex}) in AbstractPlotting at /home/juliohm/.julia/packages/AbstractPlotting/gBTAw/src/basic_recipes/convenience_functions.jl:1

Stacktrace:
 [1] _plot!(p::Combined{Any, Tuple{PointSet{2, Float64}}})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/interfaces.jl:630
 [2] plot!(p::Combined{Any, Tuple{PointSet{2, Float64}}})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/interfaces.jl:625
 [3] plot!(scene::Scene, P::Type{Combined{Any, Tuple{PointSet{2, Float64}}}}, attributes::Attributes, input::Tuple{Observable{PointSet{2, Float64}}}, args::Observable{Tuple{PointSet{2, Float64}}})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/interfaces.jl:711
 [4] plot!(scene::Scene, P::Type{Any}, attributes::Attributes, args::PointSet{2, Float64}; kw_attributes::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:show_axis,), Tuple{Bool}}})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/interfaces.jl:622
 [5] plot(P::Type{Any}, args::PointSet{2, Float64}; axis::NamedTuple{(), Tuple{}}, figure::NamedTuple{(), Tuple{}}, kw_attributes::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/figureplotting.jl:21
 [6] plot
   @ ~/.julia/packages/AbstractPlotting/gBTAw/src/figureplotting.jl:18 [inlined]
 [7] plot(args::PointSet{2, Float64}; attributes::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/recipes.jl:12
 [8] plot(args::PointSet{2, Float64})
   @ AbstractPlotting ~/.julia/packages/AbstractPlotting/gBTAw/src/recipes.jl:12
 [9] top-level scope
   @ REPL[4]:1

What am I missing?

[SimonDanisch]

The docs seem outdated, since the first argument isn't actually optional.
This would be the correct definition:

function convert_arguments(p::AbstractPlotting.PointBased, pset::PointSet)
    vecs = [coordinates(pset[i]) for i in 1:nelements(pset)]
    return convert_arguments(p, vecs) # recursive convert, since we need to convert to points
end

[juliohm]

I still get the same error here. Maybe the version?

I am on AbstractPlotting.jl v0.15.26 and Makie.jl v0.12.

[juliohm]

Also am I correct to say that all the recipe-related code lives in AbstractPlotting.jl nowadays? I would like to contribute to it to make sure that everything works smoothly with custom types, and that we can exchange data at a fundamental layer with built-in Julia Arrays.

[SimonDanisch]

Ah try scatter(pset) ... plot(...) hasn't gotten as much love from me since I like to be explicit about the desired plot type.

[juliohm]

So what would be the path to contributing to the plot recipe system? Any guidelines?

[juliohm]

If you have a set of issues opened that are related to the recipe system, I can try to address some of them.

[SimonDanisch]

So, as you may have noticed, recipes have been a bit neglected.
I myself want to slowly get back to figure out what the biggest pain points are right now.
In theory, most of the basics should still work as described in the docs, but may have bitrotten a little.
So how about you try to get your Meshes.jl recipes going and document problems & fix docs along the way?
As this issue nicely demonstrates, the error messages and the integration with the plot command are also in dire need of improvement ;)
I think, e.g. to fix plot(pset) one needs to overload plottype, to call convert_arguments before giving up. Although that's a bit finicky, if we want to avoid calling convert_arguments a second time in the plot function.

Feel free to ping me on slack or github to discuss whatever problems you're running into.

[juliohm]

That sounds like a good plan. I will put up a repository with just the recipes and then we can evolve it alongside the recipe system itself and its documentation.

[juliohm]

I have created a playground package as planned in https://github.com/JuliaGeometry/MeshViz.jl

It would be great if we could fix the recipe system to work out of the box with a single convert_arguments definition. I will try to read the AbstractPlotting.jl source to understand the pipeline from types to visualizations. Is #409 up-to-date?

[juliohm]

With Makie.jl v0.13.4 the recipes defined in MeshViz.jl are working as expected. Following the plan, I am trying to write recipes for Meshes.jl types to push the recipe system to its limits.

This is the current version of the module:

module MeshViz

using Meshes

import Makie
import Makie: convert_arguments
import Makie: plottype

# ---------
# PointSet
# ---------

plottype(::PointSet) = Makie.Scatter

convert_arguments(P::Type{<:Makie.Scatter}, pset::PointSet) =
  convert_arguments(P, coordinates.(pset))

# -----------
# SimpleMesh
# -----------

plottype(::SimpleMesh) = Makie.Mesh

function convert_arguments(P::Type{<:Makie.Mesh}, mesh::SimpleMesh)
  # retrieve geometry + topology
  verts = vertices(mesh)
  topo  = topology(mesh)
  elems = elements(topo)

  # convert to Julia arrays
  coords = reduce(hcat, coordinates(v) for v in verts)'
  connec = reduce(hcat, collect(indices(e)) for e in elems)'

  convert_arguments(P, coords, connec)
end

end

So plotting a SimpleMesh is straightforward:

julia> using Meshes, MeshViz

julia> using GLMakie

julia> mesh = SimpleMesh(Point2[(0,0),(1,0),(1,1),(0,1)], connect.([(1,2,3),(3,4,1)]))
2 SimpleMesh{2,Float64}
  4 vertices
    └─Point(0.0, 0.0)
    └─Point(1.0, 0.0)
    └─Point(1.0, 1.0)
    └─Point(0.0, 1.0)
  2 elements
    └─Triangle(1, 2, 3)
    └─Triangle(3, 4, 1)

julia> plot(mesh)

But I can't find an attribute to show the wireframe (i.e. edges of the triangles). Is it currently possible?

To give some context, here is how you can get access to the edges of the mesh without repetition:

julia> # convert to half-edge data structure for topological relations
julia> t = convert(HalfEdgeTopology, topology(mesh))
HalfEdgeTopology(HalfEdge[HalfEdge(1, 1), HalfEdge(2, nothing), HalfEdge(3, 1), HalfEdge(1, 2), HalfEdge(4, 2), HalfEdge(1, nothing), HalfEdge(3, 2), HalfEdge(4, nothing), HalfEdge(2, 1), HalfEdge(3, nothing)], Dict(2 => 4, 1 => 1), Dict(4 => 5, 2 => 9, 3 => 3, 1 => 1), Dict((3, 2) => 5, (1, 2) => 1, (3, 1) => 2, (1, 3) => 2, (4, 1) => 3, (2, 1) => 1, (1, 4) => 3, (3, 4) => 4, (4, 3) => 4, (2, 3) => 5…))

julia> # boundary relation that maps edges to vertices
julia> ∂ = Boundary{1,0}(t)
Boundary{1, 0, HalfEdgeTopology}(HalfEdgeTopology(HalfEdge[HalfEdge(1, 1), HalfEdge(2, nothing), HalfEdge(3, 1), HalfEdge(1, 2), HalfEdge(4, 2), HalfEdge(1, nothing), HalfEdge(3, 2), HalfEdge(4, nothing), HalfEdge(2, 1), HalfEdge(3, nothing)], Dict(2 => 4, 1 => 1), Dict(4 => 5, 2 => 9, 3 => 3, 1 => 1), Dict((3, 2) => 5, (1, 2) => 1, (3, 1) => 2, (1, 3) => 2, (4, 1) => 3, (2, 1) => 1, (1, 4) => 3, (3, 4) => 4, (4, 3) => 4, (2, 3) => 5…)))

julia> # loop over edges and print incident vertices
julia> for i in 1:nfacets(t)
         @show ∂(i)
       end
∂(i) = [1, 2]
∂(i) = [3, 1]
∂(i) = [4, 1]
∂(i) = [3, 4]
∂(i) = [2, 3]

With the pairs of indices above I can use the global vector of vertices to create a lines plot. Can you help explain how this can be done with the current version of the recipe system?

[SimonDanisch]

You should be able to plot linesegments(view(vertices, indices)) efficiently.
You can overload plot!(meshplot::Mesh{Tuple{args...}}) to make mesh(your_mesh; your_kw...) work. The your_kw will end up in meshplot.
See the poly recipe for Mesh:
https://github.com/JuliaPlots/Makie.jl/blob/master/src/basic_recipes/poly.jl#L52

[juliohm]

You mean defining a full recipe instead of a type recipe, right? So that I can call multiple commands instead of just converting the input type to coordinates and connectivities.

[juliohm]

I managed to write full type recipes and they are working nicely:

https://github.com/JuliaGeometry/MeshViz.jl/blob/aa6516ad941cb25e37e7f14510d6aff409a4a231/src/MeshViz.jl#L1-L66

I have a few questions regarding recipe attributes:

1. Where in the source code I can find an up-to-date list of attributes with short descriptions? I noticed that https://makie.juliaplots.org/dev/generated/plot-attributes.html#List-of-attributes contains a partial list, but attributes such as overdraw are not documented yet.
2. When we create a new recipe, we pass a list of attributes to it. What exactly is this list? Is it a default list of values that users can change in the call site? Or is it a fixed set of values that users can't change? I noticed that we can grab the attribute value from the theme of the current scene by calling Makie.theme(scene, :color) for example. I wonder what happens if I omit all attributes in the recipe definition. Will their values reproduce the default theme or the current theme?

[juliohm]

3. What is the meaning of the patch prefix in :color vs. :patchcolor, :strokewidth vs. :patchstrokewidth?

[jkrumbiegel]

It's a separate attribute because often patch-like plots (poly, density, etc) should have different colors than lines or markers. Often this is a difference in alpha only.

[juliohm]

Got it, thanks. I am able to plot triangle meshes, it looks great:

using Meshes, MeshViz # master
import GLMakie

using PlyIO

function readply(fname)
  ply = load_ply(fname)
  x = ply["vertex"]["x"]
  y = ply["vertex"]["y"]
  z = ply["vertex"]["z"]
  points = Meshes.Point.(x, y, z)
  connec = [connect(Tuple(c.+1)) for c in ply["face"]["vertex_indices"]]
  SimpleMesh(points, connec)
end

mesh = readply("beethoven.ply")

viz(mesh)

mesh = readply("iphi.ply")

viz(mesh)

Now I will try to add an option for showing the edges of the triangles. What is the recommended method to add new attributes in a plot recipe? For example, suppose I want to add a keyword option showfacets=true, does it suffice to add it to the Attributes list? I will give it a try, but let me know if there is a preferred method.

[juliohm]

Another quick question. Is it correct to say that OpenGL and related libraries know how to plot triangles efficiently but not necessarily quadrangles, pentagons, etc? I need to split the n-gons into triangles before calling the Makie.mesh! function, correct?

[SimonDanisch]

Another quick question. Is it correct to say that OpenGL and related libraries know how to plot triangles efficiently but not necessarily quadrangles, pentagons, etc?

Yeah, it will all be triangles ;) It handles all the ngons from GeometryBasics as it knows how to convert them to triangles, so I guess you'd need to do the same for Meshes.jl.

[juliohm]

What about the wireframe? Is it a similar situation where OpenGL provides an efficient routine for plotting the edges given a list of triangles? Or in this case we really need to provide a list of segments?

[SimonDanisch]

Or in this case we really need to provide a list of segments?

Well, you need to draw lines, so no triangles ;) But yes, that means you need to generate the list of linesegments!

[juliohm]

Awesome, I will work on it, it should be straightforward with the half-edge structure. Will report the results soon.

[juliohm]

Everything is working now in the master branch with full support for n-gon meshes including the display of edges:

Tomorrow I will work on adding support for plotting with colors on faces.

[juliohm]

Recipes with color working nicely too:

using Meshes, MeshViz
import GLMakie

mesh = readply("dragon.ply")

viz(mesh, elementcolor=1:nelements(mesh))

I tried playing with the interpolate option but it didn't change the result. Is there a method to disable interpolation completely so that we can see sharp transitions between triangles?

[SimonDanisch]

Not really. The problem is, that we're coloring the vertices, which get interpolated by default.
What we can do though is, to create a mesh with triangles that don't share any vertices and then give each vertex the same color for the same triangle.
I looked into it some time ago, and I didn't find any more efficient solution for cell shading with OpenGL at least...

[juliohm]

Got it, thanks. Before researching this direction, I would like to ask if it would be ok to refactor some of the existing recipes like Mesh, Poly, Wireframe, and Linesegments so that they all accept the same shape of coordinates and connectivities. For example, it would be nice if all recipes accepted a list of vertices as a Vector{<:Tuple} and a list of connectivities as a Vector{<:Tuple} as well. That way we could eliminate the reshaping that is currently necessary when we combine multiple recipes in the same plot. Would a PR in that direction be welcome?

[juliohm]

I've finished writing most of the recipes I need. You can find examples here: https://discourse.julialang.org/t/ann-announcing-meshes-jl/53973/33

I wonder if there is a simple method to find out the coordinates of the corners of the axis, as well as the pixel size. The ultimate goal is to be able to filter geometries inside a given zoom level, and after filtering is done, do a simplification of geometries based on the pixel size so that very big geometries in terms of vertices are simplified to smaller geometries when seen from far away. This is what advanced GIS software do to avoid plotting thousands of thousands of vertices of polygons every time.

[SimonDanisch]

Very nice :)

Yes, there are lots of improvements we can introduce, if we really want to enable high performance display of large meshes.
The challenge will be, to do most of the work on the GPU, which I think will be necessary to really get state of the art performance.
But maybe we can already start figuring out the design of the culling pipeline while running the algorithms on the CPU.

[ffreyer]

Can't you get a size estimate quite easily via boundingbox?

[juliohm]

Can't you get a size estimate quite easily via boundingbox?

I need the limits of the axis in data coordinates. I could then use these limits to filter the geometries that are visible given the current zoom level. Makes sense?

[SimonDanisch]

I guess for 3D we'll need the camera Frustum..
For 2D we can use something like this:
#973 (comment)

[juliohm]

Nice, I will take a look at the 2D case when I find some time. I need to switch tasks to handle other priorities now, but will come back at some point. I fully agree that some of the algorithms in Meshes.jl could be dispatched on the GPU when one is available for state of the art speed 🚀

[juliohm]

Closing as solved. The recipes are working, there is a lot of code repetition, but I think things will get better in the future with the recipe system.