[ANN] Announcing Meshes.jl

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.

1 Like

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.

(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.

4 Likes

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.

3 Likes

The solution was easier to write than I thought:

The same visualization I showed before in MeshLab now with Makie :julia:

8 Likes

I’m curious why is he upside down.

1 Like

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

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

1 Like

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

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

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 :heart_eyes: 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 :slight_smile:

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 :julia:

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

47 Likes

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

@Humphrey_Lee please watch our JuliaCon2021 talk:

and check the accompanying notebook:

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

2 Likes

Thanks. Will check them out.

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.

4 Likes

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

1 Like

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.

3 Likes

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 ?

Check the Geospatial Data Science with Julia book:

https://juliaearth.github.io/geospatial-data-science-with-julia

Use a GeoTable.

2 Likes

Thanks for the tip ! I just went through the book, and I feel like I’m missing some things.

  1. It feels like everything is based on the data values at geometry points, but no data can be linked directly with the elements, is that so ?
  2. How do you dispatch plotting recipes based on a given georef. In my current practice, I have a dedicated Mesh-like type that stores all the info and then I use a custom recipe based on various Makie basic recipes to do the job. What would be the Meshes to achieve that ? Dispatch on the column names of the geotable ?

The geometry column can contain any geometry, not just points.

The recipes are dispatched on Domain types which are collections of geometries.

Please let us know if something is not clear. Better ask on Zulip to avoid spamming everyone in this thread.

1 Like