Fix mesh/surface with NaN points

[asinghvi17]

Description

• Skip mesh cells which have at least 1 point which has some NaN component.
• Allow NaNs to be passed on from surface to mesh in CairoMakie.

Fixes MakieOrg/GeoMakie.jl#133

Type of change

Delete options that do not apply:

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)

Checklist

• Added an entry in NEWS.md (for new features and breaking changes)
• Added or changed relevant sections in the documentation
• Added unit tests for new algorithms, conversion methods, etc.
• Added reference image tests for new plotting functions, recipes, visual options, etc.

[MakieBot]

Compile Times benchmark

Note, that these numbers may fluctuate on the CI servers, so take them with a grain of salt. All benchmark results are based on the mean time and negative percent mean faster than the base branch. Note, that GLMakie + WGLMakie run on an emulated GPU, so the runtime benchmark is much slower. Results are from running:

using_time = @ctime using Backend
# Compile time
create_time = @ctime fig = scatter(1:4; color=1:4, colormap=:turbo, markersize=20, visible=true)
display_time = @ctime Makie.colorbuffer(display(fig))
# Runtime
create_time = @benchmark fig = scatter(1:4; color=1:4, colormap=:turbo, markersize=20, visible=true)
display_time = @benchmark Makie.colorbuffer(display(fig))

	using	create	display	create	display
GLMakie	12.23s (12.17, 12.34) 0.06+-	1.08s (1.07, 1.09) 0.01+-	812.51ms (793.46, 876.82) 29.01+-	11.25ms (11.12, 11.47) 0.11+-	90.98ms (90.49, 91.35) 0.33+-
master	12.18s (12.05, 12.35) 0.09+-	1.09s (1.08, 1.10) 0.01+-	809.86ms (797.70, 839.66) 14.19+-	11.24ms (11.16, 11.31) 0.05+-	91.05ms (90.25, 92.15) 0.65+-
evaluation	+0.40%, 0.05s invariant (0.63d, 0.27p, 0.08std)	-0.21%, -0.0s invariant (-0.29d, 0.60p, 0.01std)	+0.33%, 2.65ms invariant (0.12d, 0.83p, 21.60std)	+0.12%, 0.01ms invariant (0.15d, 0.78p, 0.08std)	-0.08%, -0.07ms invariant (-0.14d, 0.80p, 0.49std)
CairoMakie	12.37s (12.15, 12.99) 0.28+-	1.40s (1.37, 1.42) 0.02+-	289.72ms (281.54, 298.91) 5.50+-	12.90ms (12.76, 13.12) 0.11+-	7.07ms (6.99, 7.15) 0.06+-
master	12.39s (12.20, 12.65) 0.16+-	1.40s (1.38, 1.42) 0.02+-	286.45ms (282.60, 291.38) 3.35+-	12.44ms (12.28, 12.68) 0.16+-	7.08ms (6.99, 7.17) 0.08+-
evaluation	-0.17%, -0.02s invariant (-0.09d, 0.87p, 0.22std)	+0.12%, 0.0s invariant (0.09d, 0.87p, 0.02std)	+1.13%, 3.27ms invariant (0.72d, 0.21p, 4.43std)	+3.50%, 0.45ms slower X (3.26d, 0.00p, 0.14std)	-0.04%, -0.0ms invariant (-0.04d, 0.95p, 0.07std)
WGLMakie	12.28s (12.16, 12.46) 0.11+-	1.36s (1.34, 1.41) 0.02+-	11.84s (11.58, 11.95) 0.12+-	14.78ms (13.41, 19.20) 1.98+-	1.20s (1.12, 1.51) 0.14+-
master	12.30s (12.24, 12.40) 0.07+-	1.38s (1.33, 1.43) 0.03+-	11.93s (11.84, 12.11) 0.10+-	14.54ms (13.53, 15.48) 0.73+-	1.16s (1.13, 1.18) 0.02+-
evaluation	-0.21%, -0.03s invariant (-0.28d, 0.61p, 0.09std)	-1.08%, -0.01s invariant (-0.56d, 0.32p, 0.03std)	-0.70%, -0.08s invariant (-0.75d, 0.19p, 0.11std)	+1.61%, 0.24ms invariant (0.16d, 0.77p, 1.36std)	+3.52%, 0.04s invariant (0.43d, 0.45p, 0.08std)

[ffreyer]

This generates some black triangles around the NaN values, probably because normals at the edge become NaN:

It also happens in GLMakie to a lesser degree. I think there NaN points are ignored when generating normals.

Makie.jl/GLMakie/assets/shader/util.vert

Lines 264 to 288 in db55264

	vec3 normal_from_points(
	vec3 s0, vec3 s1, vec3 s2, vec3 s3, vec3 s4,
	vec2 uv, vec2 off1, vec2 off2, vec2 off3, vec2 off4
	){
	vec3 result = vec3(0);
	if(isinbounds(off1) && isinbounds(off2))
	{
	result += cross(s2-s0, s1-s0);
	}
	if(isinbounds(off2) && isinbounds(off3))
	{
	result += cross(s3-s0, s2-s0);
	}
	if(isinbounds(off3) && isinbounds(off4))
	{
	result += cross(s4-s0, s3-s0);
	}
	if(isinbounds(off4) && isinbounds(off1))
	{
	result += cross(s1-s0, s4-s0);
	}
	// normal should be zero, but needs to be here, because the dead-code
	// elimanation of GLSL is overly enthusiastic
	return normalize(result);
	}

.

[asinghvi17]

Good point - will probably reimplement the normal finding in CairoMakie then, for ease of use.

[SimonDanisch]

Would be nice to put that into makie, since Wglmakie also calculates normals on the cpu

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[asinghvi17]

CI is green now, we should be good to go!

[ffreyer]

This produces a NaN vector with one NaN input because centroid cmbines all 3:

julia> Makie.nan_aware_orthogonal_vector(Point3f(0), Point3f(1,0,0), Point3f(NaN))
3-element Vec{3, Float32} with indices SOneTo(3):
 NaN
 NaN
 NaN

I'm not sure why this doesn't affect things down the line, but you can probably just do

any(isnan, (v1, v2, v3)) && return Vec3f(0)
return cross(v2 - v1, v3 - v1)

instead?

[SimonDanisch]

any(isnan, (v1, v2, v3)) this is in the drawing code, so I guess that's why it isn't a problem

[ffreyer]

Ah, no, it is affecting things down the line.

If we have two faces

a --- b
| \   |
|   \ |
c --- d

and only b is NaN, orthogonal_vector will produce NaN as the face normal of (a, b, d). This is getting added to the normal at vertex a, b and d via

for i in 1:length(F)
    fi = face[i]
    normals_result[fi] = normals_result[fi] + n
end

When we iterate the faces, (a, b, d) but (a, c, d) is not. It seems like the shading is doing something differently though, so it's not as obvious.

---

With a plain color you can see the effect

zs = rand(10, 10)
ns = copy(zs)
ns[4, 3:6] .= NaN
f, a, p = surface(1..10, 1..10, ns, colormap = [:lightblue, :lightblue])
scatter!(a, [Point3f(i, j, zs[i, j]) for i in 1:10 for j in 1:10], color = isnan.(ns'[:]), depth_shift = -1f-3)
wireframe!(a, CairoMakie.surface2mesh(to_value.(p.converted)...), depth_shift = -1f-3, color = :black)
f

If you check the number of NaN normals you get 15.

m = CairoMakie.surface2mesh(to_value.(p.converted)...)
isnan.(m.normals) |> sum

If you mark them you get this (not showing position with NaN values in z)

The bottom right corner is probably not in a triangle with a NaN position, so the normal is clean. The two on left - no idea...

With

function nan_aware_orthogonal_vector(v1, v2, v3)
    (isnan(v1) || isnan(v2) || isnan(v3)) && return Vec3f(0)
    Vec3f(cross(v2 - v1, v3 - v1))
end

you get

in CairoMakie (with NaN normals still marked, new zs/ns)

Looking at this in more detail also brings up two more questions:

• Why is GLMakie extending the surface past the vertex positions around NaN?
• Should CairoMakie be iterating quad faces to avoid getting triangles in corners? (I think we could just filter faces containing a NaN coordinate for this in surface2mesh)

[asinghvi17]

I had tried quad faces for surface earlier, but file sizes ballooned so it turned out to be counterproductive.

[ffreyer]

Huh, good to know. The changes I pushed just remove some faces so it shouldn't be a problem

[ffreyer]

This can look quite a bit different from GLMakie in extreme cases, e.g.

surface([0,1,1], [0, 1], [0 0; 1 0; 0 1], colormap = [:lightblue, :lightblue])

(CairoMakie | GLMakie)

but I'm not sure whether CairoMakie or GLMakie should be seen as the problem here.

In terms of NaN values this produces nicer results than GLMakie (no darkening) and it also happens to clean up the little artifact you get when rendering a sphere as a surface plot:

[ffreyer]

I pushed a fix for the NaN normals and removed the triangles by filtering. I now get

with

zs = rand(10, 10)
ns = copy(zs)
ns[4, 3:6] .= NaN
f, a, p = surface(1..10, 1..10, ns, colormap = [:lightblue, :lightblue])
m = CairoMakie.surface2mesh(to_value.(p.converted)...)
scatter!(a, m.position, color = isnan.(m.normals), depth_shift = -1f-3)
wireframe!(a, m, depth_shift = -1f-3, color = :black)
f

I think with that I'm happy. Now GLMakie needs to catch up.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[asinghvi17]

Thanks! Yeah, in surfaces we really should skip any cell with a NaN. Do you mind if I add your example as a reference test @ffreyer?

[ffreyer]

I think this may just need a refimg update

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[SimonDanisch]

Damn, this introduces a regression for heatmaps & surface:


The recorded version also seems wrong, but I diffusely remember giving up debugging the CI GPU, since it seemed to be an obscure glsl issue:

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[SimonDanisch]

WGLMakie flakiness should be fixed by #2661

[SimonDanisch]

Hm, it's this change that makes the test fail.. It looks a bit like the colors get moved out of the surface (left: heatmap, right: surface):
PR:

With the old vec2(index2D) / (vec2(dims)-1.0):

Considering, that the old behavior maps more directly to the displayed matrix, why was this changed?

[ffreyer]

With vec2(index2D) / (vec2(dims)-1.0) you get values between 0 and 1, so the first index will sample the left edge of a pixel in a texture and the last index will sample the right edge of a pixel.

With vec2(index2D) / (vec2(dims)) you'd always sample the left edge of a pixel. With the 0.5 shift that becomes the center of a pixel.

I think I made this change to fix the discrepancy between GLMakies generated mesh and the wireframe representation of CairoMakies mesh: (right)

[ffreyer]

I think is also an expected/desired change.

The color data we have here is:

 1.0  NaN        NaN        5.0
 2.5    2.66667    2.83333  3.0

In surface these should map directly to vertices and we should be interpolating between them. Assuming quad faces, this should generate all NaN's.

With vertices marked as + and the positions where the texture is sampled without interpolation marked with pink dots, this should be the before:

And this the after:

[SimonDanisch]

Hm, but we shouldn't render the edges pretty much outside the visible field?
I think, if we do this, we'll need to move the interpolation center to the middle of the quad!

[SimonDanisch]

Maybe to get this branch merged, revert things for now, so that we dont have a regression?

[ffreyer]

My understanding of surface plots is that it's a mesh generating method where vertex positions follow from xs and ys and the colors/values passed should exactly map to those positions. That's not the case on master. That's a bug imo.

Moving the color to the center of a quad doesn't work because we have N colors for N vertices. We could move the full vertex to quad center but then we would need to generate an N+1 grid to keep quad sizes regular.

Maybe it would be useful to compare to some other plotting libraries? To me it should just be like my after example.

[ffreyer]

If we look at this with automatically chosen color values, i.e. z values, I would expect the first color in the colormap to directly correspond to the lowest z value, the last to the highest and anything in between to be smoothly interpolated. That's not the case on master, but is the case in the pr.

This example makes the difference very obvious:

scene = Scene()
surface!(scene, [1 2 3; 1 2 3; 1 2 3], colormap = [:black, :white, :black], shading = false)
cam3d!(scene, projectiontype = Makie.Orthographic)
center!(scene)
axis3d!(scene)
lines!(scene, Point3f[(3, 3, 1), (3, 3, 2), (3, 3, 3)], color = [:black, :white, :black], linewidth = 50)
update_cam!(scene, Vec3f(5.5, 5, 2), Vec3f(1.5, 1.5, 2))
scene

The line to the left reaches full black at the line start/end, i.e. at minimum and maximum z. The surface on the other hand at z values significantly larger than the minimum and significantly lower than the maximum.

On another note, why is surface offset like heatmap? As a continuous thing I would expect it to use axes(zs, 1) and axes(zs, 2), not some larger range like heatmap.

[ffreyer]

Maybe it would be useful to compare to some other plotting libraries? To me it should just be like my after example.

MatPlotLibs surf seems to be discrete so it's not comparable.

Plots with GR backend

Plots with Plotly backend

Plots with PythonPlot backend

Plots with Gaston backend

All of these are

Plots.surface([1 1 1; 2 2 2; 3 3 3], colormap = cgrad([:black, :white, :black], 3))

---

This PR

Out of those, none do our 0.0:1.5:3.0 x and y range, half do flat face coloring, half do smooth coloring and all match min/max z to the first/last color in the colormap like in this pr.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[github-actions]

Missing reference images

Found 1 new images without existing references.
Upload new reference images before merging this PR.

[ffreyer]

Docs failure should be unrelated. The code it highlights works fine locally

[ffreyer]

Closing this since it's in the beta branch