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

Trimesh rasterization seems slightly off #806

Closed
ronaldsipkema opened this issue Oct 9, 2019 · 8 comments · Fixed by #1092
Closed

Trimesh rasterization seems slightly off #806

ronaldsipkema opened this issue Oct 9, 2019 · 8 comments · Fixed by #1092
Labels
Milestone

Comments

@ronaldsipkema
Copy link

ronaldsipkema commented Oct 9, 2019

Description

While playing around with Canvas.trimesh and rasterizing simple geomtery it seems that the rule described in the documentation ("a pixel (bin) is treated as belonging to a given triangle if its center falls either inside that triangle or along its top or left edge") is not adhered to, i.e., some pixels that do not belong to a triangle get a value, and the other way around.

rasterized_circle

How to reproduce

import pandas as pd
import shapely.geometry
from scipy.spatial import Delaunay
import datashader as ds
import matplotlib.pyplot as plt


def triangulate(vertices, x="x", y="y"):
    triang = Delaunay(vertices[[x,y]].values)
    return pd.DataFrame(triang.simplices, columns=['v0', 'v1', 'v2'])


shape = shapely.geometry.Point(0, 0).buffer(1)
vertices = pd.DataFrame(shape.exterior.coords, columns=('x', 'y'))
vertices['z'] = 1.0

simplices = triangulate(vertices)

canvas = ds.Canvas(plot_width=10, plot_height=10,
                x_range=(-1, 1), y_range=(-1, 1))

canvas.trimesh(vertices, simplices, agg=ds.any()).plot(alpha=0.5)
plt.plot(vertices.x, vertices.y)

@jbednar
Copy link
Member

jbednar commented Oct 23, 2019

Thanks for the report! I can see that the array pattern shown does not match the Shapely circle, but as the triangles are not visible I can't see how the triangles relate to the array pattern, so I can't see where the discrepancy arises.

@jbednar
Copy link
Member

jbednar commented Oct 23, 2019

Oh, and I wasn't able to run the code above, due to a Pandas issue:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-d9e37701cbad> in <module>
     12 
     13 shape = shapely.geometry.Point(0, 0).buffer(1)
---> 14 vertices = pd.DataFrame(shape.exterior.coords, columns=('x', 'y'))
     15 vertices['z'] = 1.0
     16 

~/miniconda3/envs/pyviz/lib/python3.7/site-packages/pandas/core/frame.py in __init__(self, data, index, columns, dtype, copy)
    420                                          dtype=values.dtype, copy=False)
    421             else:
--> 422                 raise ValueError('DataFrame constructor not properly called!')
    423 
    424         NDFrame.__init__(self, mgr, fastpath=True)

ValueError: DataFrame constructor not properly called!

@ronaldsipkema
Copy link
Author

I can see that the array pattern shown does not match the Shapely circle, but as the triangles are not visible I can't see how the triangles relate to the array pattern, so I can't see where the discrepancy arises.

If you change the last line of the script to plt.triplot(vertices.x, vertices.y, simplices, alpha=0.5) you'll get the following image:

with_triangles

@ronaldsipkema
Copy link
Author

Oh, and I wasn't able to run the code above, due to a Pandas issue:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-d9e37701cbad> in <module>
     12 
     13 shape = shapely.geometry.Point(0, 0).buffer(1)
---> 14 vertices = pd.DataFrame(shape.exterior.coords, columns=('x', 'y'))
     15 vertices['z'] = 1.0
     16 

~/miniconda3/envs/pyviz/lib/python3.7/site-packages/pandas/core/frame.py in __init__(self, data, index, columns, dtype, copy)
    420                                          dtype=values.dtype, copy=False)
    421             else:
--> 422                 raise ValueError('DataFrame constructor not properly called!')
    423 
    424         NDFrame.__init__(self, mgr, fastpath=True)

ValueError: DataFrame constructor not properly called!

I just tested in a newly created conda environment (with pandas 0.25.2, shapely 1.6.4) and everything works fine. Which versions are you using? I'm afraid I can't directly come up with why this would stop working without knowing the pandas versions and / or what the value of shape.exterior.coords looks like on your system.

@jbednar
Copy link
Member

jbednar commented Oct 28, 2019

I was able to run it after updating Pandas, and I've simplified it to what looks like a genuine off-by-a-half-pixel error to me:

import pandas as pd, datashader as ds, holoviews as hv
hv.extension("bokeh")

coords    = [(1.0, 0.0, 1.0), (0, -1.0, 1.0), (-0.6, -0.7, 1.0), (1.0, 0.0, 1.0)]
vertices  = pd.DataFrame(coords, columns=('x', 'y', 'z'))
simplices = pd.DataFrame(dict(v0=[1],v1=[0],v2=[2]))
canvas    = ds.Canvas(plot_width=10, plot_height=10, x_range=(-1, 1), y_range=(-1, 1))
agg       = canvas.trimesh(vertices, simplices, agg=ds.any())

opts=hv.opts.Image(cmap="greys", show_grid=True, xticks=10, yticks=10, alpha=0.5)
hv.Image(agg).opts(opts) * hv.Path(coords)

image

I.e., the whole image seems to be rendered up and to the left of where it should be, and in particular the pixel on the bottom row between -0.2 and 0 should be on. Seems like that implies either that the various test_trimesh_* functions are either inadequate or contain incorrect test data.

@ianthomas23
Copy link
Member

Upon investigation, the code is doing what is intended but not what the user expects. There is not a constant half a pixel offset. The location of each of the (x, y) vertices is snapped to the middle of the nearest pixel before the triangles are rendered, so each vertex is offset by between minus half a pixel and plus half a pixel in each of the x and y directions. This is best illustrated through an example based on @jbednar's example above:

import datashader as ds
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
import pandas as pd

coords         = np.asarray([(1.0, 0.0, 1.0), (0.0, -1.0, 1.0), (-0.6, -0.7, 1.0), (1.0, 0.0, 1.0)])
snapped_coords = np.asarray([(0.9, 0.1),      (0.1, -0.9),      (-0.7, -0.7),      (0.9, 0.1     )])

vertices = pd.DataFrame(coords, columns=('x', 'y', 'z'))
triang = mtri.Triangulation(vertices["x"], vertices["y"])
tris = pd.DataFrame(triang.triangles, columns=["v0", "v1", "v2"])

canvas = ds.Canvas(plot_width=10, plot_height=10, x_range=(-1, 1), y_range=(-1, 1))
agg = canvas.trimesh(vertices, tris, agg=ds.count())

fig, ax = plt.subplots(1, 1)
agg.plot(ax=ax, cmap="Greys", alpha=0.3)

# Grid
x, y = np.meshgrid(np.linspace(canvas.x_range[0], canvas.x_range[1], 1 + canvas.plot_width),
                   np.linspace(canvas.y_range[0], canvas.y_range[1], 1 + canvas.plot_height))
ax.plot(x, y, x.T, y.T, c='k', alpha=0.1)

ax.plot(coords[:, 0], coords[:, 1], label="unsnapped")
ax.plot(snapped_coords[:, 0], snapped_coords[:, 1], 'o-', label="snapped")
ax.legend()
plt.show()

The output produced is:
trimesh
The user-specified triangle is drawn in blue and the snapped triangle in orange. The rendering is correct for the orange (snapped) triangle as the pixels whose centres are inside the orange triangle are rendered in grey. But the user is expecting those pixels whose centres are inside the blue (unsnapped) triangle to be rendered, and is disappointed with the results.

To fix this it is necessary to change the rendering code to use unsnapped rather than snapped vertices. This will need some investigation into whether it is a minor alteration or a major rewrite.

@ianthomas23
Copy link
Member

The problems in the OP's circle plot are a side-effect of this. The vertices on the right of the circle are all snapped to the middle of the last x-pixel, and when aggregating we don't include pixels whose centres are exactly on the right hand triangle edge so the right-hand column is empty. The same applies to the bottom row.

This is most easily illustrated using a trimesh that is a square entirely covering the canvas:

import datashader as ds
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

coords = np.asarray([(0.0, 0.0, 1.0), (1.0, 0.0, 1.0), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0)])
vertices = pd.DataFrame(coords, columns=('x', 'y', 'z'))
triangles = [[0, 1, 2], [0, 2, 3]]
tris = pd.DataFrame(triangles, columns=["v0", "v1", "v2"])

canvas = ds.Canvas(plot_width=4, plot_height=4, x_range=(0, 1), y_range=(0, 1))
agg = canvas.trimesh(vertices, tris, agg=ds.count())

fig, ax = plt.subplots()
agg.plot(ax=ax, cmap="Greys", alpha=0.3)

# Grid
x, y = np.meshgrid(np.linspace(canvas.x_range[0], canvas.x_range[1], 1 + canvas.plot_width),
                   np.linspace(canvas.y_range[0], canvas.y_range[1], 1 + canvas.plot_height))
ax.plot(x, y, x.T, y.T, c='k', alpha=0.1)

plt.show()

trimesh

@jbednar
Copy link
Member

jbednar commented Jun 7, 2022

Ouch! Good sleuthing. Elsewhere in Datashader we enforce a special case including the bottom and right edge of the canvas because of how surprising it is to people that that final pixel is not included (even though it's a fully consistent definition of gridded coordinates). Here, I do think that unsnapped pixels are what people would most want, so if it's feasible to avoid snapping them, I think we should. Even without snapping, getting full coverage for that last example would probably take a special case, and I do think it's a useful special case.

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

Successfully merging a pull request may close this issue.

4 participants