[ANN] Meshes.jl - Computational Geometry in Julia

[juliohm]

Meshes.jl v0.16 is out with new features: Taubin smoothing, general shape intersection with the GJK algorithm and the ability to store Tables.jl of data for any rank of the mesh.

• rank 0 (vertices)
• rank 1 (segments)
• rank 2 (polygons)
• rank 3 (polyhedrons)

To illustrate the Taubin smoothing method, consider the following mesh:

It can be smoothed with:

smesh = smooth(mesh, TaubinSmoothing(100, weights=:cotangent))

and this is the result:

Unlike refinement methods illustrated in previous posts, smoothing methods do not change the number of vertices of the mesh, they just move the coordinates to remove high-frequencies.

[heliosdrm]

@juliohm

If we could write small packages for specific file formats such as OBJ and OFF similar to PlyIO.jl, then these packages would give the user full flexibility to manipulate the contents of the file. At the end MeshIO.jl could just depend on these packages to do the basic read with a single load command.

There is something that I’m missing. If I understand correctly, you are suggesting to create packages that read OBJ, OFF, etc. files and give Meshes.Mesh objects, which Meshes.jl can handle, is that right? But on the other hand, there is MeshIO.jl, which already provides functions to read those files, but their output are GeometryBasics.AbstractMesh objects.

Now, how shall this duality be managed? Since Meshes.jl is a fork of GeometryBasics.jl – and looking into its documentation, it seems that there is a lot of great work invested on it, is Meshes.jl meant to be merged back on GeometryBasics.jl in the future? Or maybe it will replace it, and then MeshIO.jl should change its output to Meshes.Mesh?

[juliohm]

@heliosdrm

If I understand correctly, you are suggesting to create packages that read OBJ, OFF, etc. files and give Meshes.Mesh objects, which Meshes.jl can handle, is that right?

No, I am suggesting to write packages that are agnostic to geometry and mesh types like PlyIO.jl does. They return a simple dictionary structure that one can use to query coordinates and connectivities from the file. After these building blocks are ready, it shouldn’t be difficult to write another layer of software in MeshIO.jl to consume these dictionary-like structures and return instances of specific mesh types.

@heliosdrm

But on the other hand, there is MeshIO.jl, which already provides functions to read those files, but their output are GeometryBasics.AbstractMesh objects.

You will notice for example that reading binary PLY files with MeshIO.jl is currently a bit problematic. PlyIO.jl can read binary files without problems. That is why I think these efforts could be split into smaller packages where the maintainers know the format in depth. MeshIO.jl could just be a wrapper package that puts together all these loaders into a common API. But I can’t speak for the whole community, this is just my opinion.

@heliosdrm

Now, how shall this duality be managed? Since Meshes.jl is a fork of GeometryBasics.jl – and looking into its documentation, it seems that there is a lot of great work invested on it, is Meshes.jl meant to be merged back on GeometryBasics.jl in the future?

Meshes.jl is a completely different codebase, and if you open the source code you will see the different design choices. The fork existed historically because I wanted to make sure that Meshes.jl features are a super-set of GeometryBasics.jl features in case the visualization team behind Makie.jl needed a fresh approach to achieve their own goals. It is very reasonable to start from a common ground to make sure that nothing is missing in terms of visualization requirements.

@heliosdrm

Or maybe it will replace it, and then MeshIO.jl should change its output to Meshes.Mesh ?

Meshes.jl is a project on its own, and I have to work on it as part of my research. If it is gonna replace any other project in the long-term that is not up to me. I certainly see value in a coherent picture where advanced meshing algorithms are developed in Meshes.jl and where advanced visualizations are developed in Makie.jl in parallel. In the meantime, people can help with both projects in multiple ways.

[heliosdrm]

I think I’m starting to understand. I was confused by the reference in the README of MeshIO.jl to Meshes.jl, which points to the package announced in this post. That made me think that the type of mesh objects that MeshIO deals with, should be interoperable with the types of this new package, but I was failing to find how.

However, it seems that when that line of the README was written (6 years ago), it referred to a different package, possibly the deprecated package that was renamed to OldMeshes.jl.

So, the current Meshes.jl is related to those other projects, and partly based on them, but departs from them in a fashion that in this moment makes it incompatible with MeshIO.jl.

[jzr]

@paulmelis

I’m not familiar with Zulip, but it seems like a Slack clone, so does it really facilitate discussions as well as Discourse?

Zulip is not a Slack clone. Instead,

Zulip is an open source alternative to Slack which has

• Better Markdown support
• Julia Syntax highlighting
• Unlimited message history for us and other open source communities
• Beautiful built in LaTeX support for displaying math
• An approach to conversation threads that may be initially a bit different from how Slack handles it, but I think much better once you use it a bit
• Really responsive and helpful developers. If you have a problem or complaint with Zulip you can open an issue on their Git repo and expect a quick response! This really contrasted with my experience with Slack at various times, especially when Slack started screwing with text formatting.

Come join the Julia Zulip chat!

Hey everyone, just wanted to advertise the Julia Zulip chat server. I think it’s great but there’s a chicken and egg problem where it’s hard to convince people to use it because there’s currently few active users. Zulip is an open source alternative to Slack which has Better Markdown support Julia Syntax highlighting Unlimited message history for us and other open source communities Beautiful built in LaTeX support for displaying math An approach to conversation threads that may be initially …

(But there is also a time-limited “ephemeral” channel for some informal chats, to emulate the slack experience.) Discourse is nice for a forum feel, but zulip is nice for a more lightweight feel while retaining the powerful features of a discussion system. (Just for an illustrative example, many more messages fit on a Zulip screen than on a Discourse screen.)

I really like Zulip, and a bunch of our most distinguished community members have moved their contributions from Slack to Zulip because they value the accessibility, threading, searchability, formatting, open-source community, data ownership, and so on. I hope that trend continues.

[juliohm]

The master branch of MeshViz.jl already supports visualization of n-gon meshes with Makie.jl:

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)

Here is a large example with 435545 vertices and 871306 triangles:

mesh = readply("dragon.ply")

viz(mesh)

We can refine the mesh and obtain 2613810 vertices and 2613918 quadrangles in 20s:

refined = refine(mesh, CatmullClark())

and visualize in 5s:

viz(refined)

Tomorrow I will work on adding an option to show the edges of the triangles, quadrangles, etc., which should be straightforward given the data structures in Meshes.jl.

[juliohm]

The solution was easier to write than I thought:

The same visualization I showed before in MeshLab now with Makie

[jzr]

I’m curious why is he upside down.

[juliohm]

I don’t know, I just loaded the PLY file and plotted. Maybe the coordinates are flipped in the file.

[sdanisch]

Yeah half of the formats have z up and the other half has y up or something like that ^^

[paulmelis]

@juliohm

I don’t know, I just loaded the PLY file and plotted. Maybe the coordinates are flipped in the file.

The axis labels are also somewhat confusing, as +Y is upwards in the image, but the labels are upside-down?

[juliohm]

This issue is unrelated to the recipes above. Perhaps you can ask for support in the Makie forum or open an issue on GitHub.

[juliohm]

A super exciting update. We can now visualize any type of mesh or collection of geometries from Meshes.jl using a consistent and unified set of attributes in MeshViz.jl. Let me illustrate the features that are already available in the latest stable releases of Meshes.jl + MeshViz.jl so that you understand how general this approach is.

MeshViz.jl provides a single plot recipe called viz, which currently has the following attributes:

julia> using MeshViz

help> viz
  viz(object)

  Visualize Meshes.jl object with various options:

    •  elementcolor - color of the elements (e.g. triangles)

    •  facetcolor - color of the facets (e.g. edges)

    •  vertexcolor - color of the vertices (i.e. points)

    •  showvertices - tells whether or not to show the vertices

    •  showfacets - tells whether or not to show the facets

    •  variable - informs which variable to visualize

    •  decimation - decimation tolerance for polygons

Most often we will need to use elementcolor to color the elements of the mesh or collections of geometries. Let’s show some examples from the MeshViz.jl README:

using Meshes, MeshViz
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("dragon.ply")

viz(mesh,
  showfacets=false,
  elementcolor=1:nelements(mesh),
  colormap=:Spectral
)

grid = CartesianGrid(10,10)

viz(grid)

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

grid = CartesianGrid(10,10,10)

viz(grid)

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

We can also plot cloropleth plots at scale. This is possible because we use sophisticated Meshes.jl algorithms to simplify (or decimate) the polygonal chains before calling Makie.

To give an example, we can automatically download a very detailed map of Brazil, and plot it:

using GeoTables

# Brazil states as Meshes.jl polygons
BRA = GeoTables.gadm("BRA", children=true)

viz(BRA.geometry)

The colors are set like before using the elementcolor attribute:

viz(BRA.geometry, elementcolor=1:length(BRA.geometry))

The good news is that we can visualize any place in the world basically because we are using GADM.jl within GeoTables.jl to convert shapefiles and geopackages to Meshes.jl geometries automatically for the user. Notice how we can get the boundaries of cities inside of Rio de Janeiro:

RIO = GeoTables.gadm("BRA", "Rio de Janeiro", children=true)

viz(RIO.geometry, decimation=0.001)

The decimation tolerance controls the simplification step. The smaller the decimation the more we preserve the original boundaries of the shapes. The plan is to set this automatically in the future based on the zoom level of the visualization.

GeoStats.jl + Meshes.jl

Now the part I am really excited about We can do geostatistical estimation, simulation, learning on these domains using GeoStats.jl and produce very sophisticated geospatial workflows:

# find row of table with Amazon geometry
ind = findfirst(==("Amazonas"), BRA.NAME_1)

# simplify geometry (thousands of vertices) before meshing
shape = simplify(BRA.geometry[ind], DouglasPeucker(0.05))

mesh = discretize(shape, FIST())

viz(mesh)

Let’s for example simulate a Gaussian field on this mesh:

using GeoStats

# ask for one realization of Gaussian variable z
problem = SimulationProblem(mesh, :z => Float64, 1)

# simple LU solver
solver = LUGS(:z => (variogram=SphericalVariogram(range=15.,),))

# simulate Gaussian field on the mesh!
sol = solve(problem, solver)

viz(sol[1], variable=:z,
  colormap = :inferno,
  showfacets=false
)

Notice that there are visualization artifacts due to the way Makie.jl interpolates the colors, and of course the mesh is also too coarse for any realistic interpolation. We are working on fixing the interpolation issue in the visualization pipeline and also working on more discretization methods that guarantee semi-regular meshes

Just for fun, let’s also create Gaussian data over a Cartesian grid (a.k.a. image data, raster data):

grid = CartesianGrid((-55.,-20.), (-45.,-10.), dims=(200,200))

problem2 = SimulationProblem(grid, :z => Float64, 1)

# for grids we can use the FFT solver
solver2 = FFTGS(:z => (variogram=SphericalVariogram(range=5.,),))

sol2 = solve(problem2, solver2)

Let’s put everything together to see how things look:

viz(BRA.geometry)

viz!(sol[1], variable=:z,
  colormap = :inferno,
  showfacets=false
)

viz!(sol2[1], variable=:z,
  colormap = :inferno,
  showfacets=false
)

If you work with geospatial data like myself, you will see the value of this unified interface. We are trying to remove all the pain people have with GIS software, and are getting very close from a very unique experience. People from other languages will wish they had something so general

This is the tip of the iceberg. I will be presenting some cool features in my JuliaCon talk next month.

[Humphrey_Lee]

@juliohm

We support basic plotting features via Plots.jl recipes, but the plan is to migrate these recipes to Makie.jl when its recipe system becomes available. Help is appreciated, especially if you are familiar with the AbstractPlotting.jl internals. In the long-run, it would be ideal to see more synergy between the Makie.jl stack and Meshes.jl (cc: @sdanisch @jkrumbiegel).

Any update on the “migration” to Makie (GSOC)? Specifically, how to make the Geostats.jl example work with Makie.jl?

[juliohm]

@Humphrey_Lee please watch our JuliaCon2021 talk:

Geostatistical Learning | Júlio Hoffimann | JuliaCon 2021

and check the accompanying notebook:

GitHub - juliohm/JuliaCon2021: Pluto presentation for JuliaCon2021 talk

We are hosting Makie recipes in MeshViz.jl temporarily while the Makie recipe system is being refactored with lightweight dependencies.

[Humphrey_Lee]

Thanks. Will check them out.

[juliohm]

Meshes.jl v0.18 is out with new algorithms, tons of improvements and bug fixes:

https://juliageometry.github.io/Meshes.jl/stable

A few breaking changes that deserve extra attention:

• Coordinates of points constructed with Integer literals in expressions such as Point(1, 2) are converted to Float64. Most algorithms in geometric processing assume a continuous R^n and passing discrete Integer coordinates just leads to InexactError everywhere. This will improve the experience of first-time users of the language who are not familiar with its rich type system.

• Intersection algorithms between two geometries g1 and g2 now allow a first argument function intersecttype(f, g1, g2) that is applied to the return value to force type stability. This is Andrey’s trick to handle heterogeneous return types. Although we’ve implemented this feature in this release, we are still considering alternative designs for type-stable intersections in future releases. This issue is being discussed here.

If you would like to get involved, check the good first issue label in our issue tracker.

[PetrKryslUCSD]

BTW at least a few CAD softwares (perhaps not all) use integers to store geometry. It is easier to build exact predicates that way.

[juliohm]

We may need to consider custom types for coordinates that are not Integer but can be used for exact representation. Or maybe jump directly to interval arithmetic which sounds more promising in general.

[BambOoxX]

Hey, I’ve seen above that Meshes was a fork from GeometryBasics. However it seems the metadata part of GeometryBasics was not kept. What is the correct practice to add metadata to Meshes meshes ?