holoviz / datashader

Quickly and accurately render even the largest data.
http://datashader.org
BSD 3-Clause "New" or "Revised" License
3.32k stars 365 forks source link

Trimesh: document associating data with simplices #1110

Open Huite opened 2 years ago

Huite commented 2 years ago

Is your feature request related to a problem? Please describe.

(I think this question applies to the Holoviews/Geoviews trimesh as well.)

Having the option to efficiently visualize large irregular meshes via datashader is quite useful. However, unless I'm mistaken, the trimesh currently only supports data on the nodes. This does the trick for e.g. finite element computation, but not for finite volume, where the data is associated with the simplex/face/cell (more similar to raster data).

(If my understanding is incorrect, I'd be happy to hear how to do it -- I think neither the docs/examples nor the source code show it in that case!)

Describe the solution you'd like

Ideally, I can just give some data that matches the length of the simplices, and get a triangle-by-triangle fill. For datashader, I suppose you could stick it in an additional column of the simplices dataframe.

With matplotlib, I can plot both data associated with the nodes via triplot or data on the simplices with PolyCollection (although the latter one seems a bit obscure).

Describe alternatives you've considered

Any alternative seems rather impractical (e.g. I could create a voronoi tesselation, then triangulate it to get data associated with the cell centroids to get something on the nodes).

Add any other context or screenshots about the feature request here.

I guess this relates somewhat to https://github.com/holoviz/holoviews/issues/3812.

It's relatively straightforward to triangulate a mesh of (convex) cells in a vectorized manner, storing the indices to grab the values, e.g.:

IntArray = np.ndarray
IntDtype = np.int64

def _triangulate(i: IntArray, j: IntArray, n_triangle_per_row: IntArray) -> IntArray:
    n_triangle = n_triangle_per_row.sum()
    n_face = len(i)
    index_first = np.argwhere(np.diff(i, prepend=-1) != 0)
    index_second = index_first + 1
    index_last = np.argwhere(np.diff(i, append=-1) != 0)

    first = np.full(n_face, False)
    first[index_first] = True
    second = np.full(n_face, True) & ~first
    second[index_last] = False
    third = np.full(n_face, True) & ~first
    third[index_second] = False

    triangles = np.empty((n_triangle, 3), IntDType)
    triangles[:, 0] = np.repeat(j[first], n_triangle_per_row)
    triangles[:, 1] = j[second]
    triangles[:, 2] = j[third]
    return triangles

def triangulate(face_node_connectivity: IntArray, fill_value: int):
    n_face, n_max = face_node_connectivity.shape

    if n_max == 3:
        triangles = face_node_connectivity.copy()
        return triangles, np.arange(n_face)

    valid = face_node_connectivity != fill_value
    n_per_row = valid.sum(axis=1)
    n_triangle_per_row = n_per_row - 2
    i = np.repeat(np.arange(n_face), n_per_row)
    j = face_node_connectivity.ravel()[valid.ravel()]
    triangles = _triangulate(i, j, n_triangle_per_row)

    triangle_face_connectivity = np.repeat(
        np.arange(n_face), repeats=n_triangle_per_row
    )
    return triangles, triangle_face_connectivity

This would allow plotting of any mesh (triangles, quads, or even voronoi cells).

It should be easier than the node based visualization: there's no need to interpolate values inside of the triangles, just sampling suffices.

ianthomas23 commented 2 years ago

Datashader already supports weights on the triangles as well as on the vertices. Here is an example showing both:

import datashader as ds
import pandas as pd

canvas = ds.Canvas(plot_width=200, plot_height=200)

# Vertex-based weights.
vertices = pd.DataFrame(
    [[0, 0, 1.0], [1, 0, 2.5], [0, 1, 1.5], [-1, 0, 0.5], [0, -1, 0.75]],
    columns=("x", "y", "vert_weight"))
triangles = pd.DataFrame(
    [[0, 1, 2], [0, 2, 3], [0, 3, 4], [0, 4, 1]],
    columns=["v0", "v1", "v2"])

agg = canvas.trimesh(vertices, triangles, agg=ds.sum("vert_weight"))
im = ds.transfer_functions.shade(agg, how="linear")
ds.utils.export_image(im, "vert_weight")

# Triangle-based weights.
vertices = pd.DataFrame(
    [[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1]],
    columns=("x", "y"))
triangles = pd.DataFrame(
    [[0, 1, 2, 0.5], [0, 2, 3, 1.0], [0, 3, 4, 1.5], [0, 4, 1, 2.0]],
    columns=["v0", "v1", "v2", "tri_weight"])

agg = canvas.trimesh(vertices, triangles, agg=ds.sum("tri_weight"))
im = ds.transfer_functions.shade(agg, how="linear")
ds.utils.export_image(im, "tri_weight")

and the images produced: vert_weight tri_weight

Weights are taken from the vertices DataFrame if there are more than 2 columns, otherwise they are taken from the triangles DataFrame.

ianthomas23 commented 2 years ago

Certainly the docs could be more explicit about this.

jbednar commented 2 years ago

A PR adding the above example to the Datashader docs would be very welcome, @Huite !

Huite commented 2 years ago

Great! -- I'll put the PR on my (somewhat crowded) to do list!