Hexahedron volume plot

[raman_kumar]

How to volume plot a Hexahedron for colors given at each vertices?

using Meshes, CoordRefSystems ,GLMakie, LinearAlgebra, Unitful

r =  Float32[1.6, 1.754088, 1.9230156, 2.1082115, 2.311243, 2.5338273, 2.7778475, 3.0453684, 3.3386526, 3.6601818, 4.012676, 4.3991165, 4.8227735, 5.287231, 5.7964177, 6.354642, 6.9666257, 7.637547, 8.373081]
θ =  Float32[0.0, 0.09817477, 0.19634955, 0.2945243, 0.3926991]
ϕ =  Float32[0.0, 0.3926991, 0.7853982, 1.1780972, 1.5707964, 1.9634954, 2.3561945, 2.7488935, 3.1415927, 3.5342917, 3.9269907, 4.3196898, 4.712389, 5.105088, 5.497787, 5.8904862, 6.2831855]

g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ)
hxdrn = element(g,1)
spherpts = coords.(vertices(hxdrn))
cartpts =  convert.(Cartesian, spherpts)
viz(hxdrn, alpha=0.3)

[raman_kumar]

@ffreyer What stops us to volume plot a Hexahedron for colors given at each vertices? Is there some limitation in OpenGL or in Makie.jl itself? Can it be solved with Vulkan?

[ffreyer]

Ignoring GPU rendering, I think the main thing you need is a function that tells you the color which should go at a given position. For a uniform 3D cartesian grid like a 3D array provides that’s relatively easy. You just find the 8 grid points around the sample position and nest linear interpolations. A tetrahedron is probably also some combination of linear interpolations. Beyond that, no clue. If you want to try working it out, here’s roughly what :absorptionrgba would look like for a single Rect3:

using Colors
using GeometryBasics

function get_interpolated_color(sample_pos, rect_positions, rect_colors)
    # This relies on the order of vertices
    mini = first(rect_positions)
    maxi = last(rect_positions)

    if sample_pos in Rect3f(mini, maxi .- mini)
        color = mapreduce(+, rect_positions, rect_colors) do vertex_pos, color
            interp = abs.(sample_pos - vertex_pos) ./ (maxi - mini)
            return color * prod(interp)
        end

        return color
    else
        return Vec4f(0,0,0,0)
    end
end

rect_positions = coordinates(Rect3f(-0.5, -0.5, -0.5, 1, 1, 1))
rect_colors = [Vec4f(0.5 + 0.5x, 0.5 + 0.5y, 0.5 + 0.5z, 0.1) for (x, y, z) in rect_positions]

width = 400
height = 300
img = [RGBA(0,0,0,0) for x in 1:width, y in 1:height]

R = Mat3f(0.7886752, -0.21132484, -0.57735026, -0.21132484, 0.7886752, -0.57735026, 0.57735026, 0.57735026, 0.5773503)
steps = 100
step_weight = 1.0
for (i, x) in enumerate(range(-1, 1, width))
    @info i
    for (j, y) in enumerate(range(-1, 1, height))
        accumulated_color = Vec3f(0)
        remaining_intensity = 1.0
        for z in range(-1, 1, steps)
            pos = R * Point3f(x, y, z)
            sample_rgba = get_interpolated_color(pos, rect_positions, rect_colors)

            sample_weight = step_weight * sample_rgba[4]
            remaining_intensity *= (1 - sample_weight)
            sample_color = sample_weight * remaining_intensity * Vec3f(sample_rgba)
            accumulated_color += sample_color

            remaining_intensity < 0.01 && break
        end
        img[i, j] = RGBA(clamp.(accumulated_color, 0, 1)..., 1 - remaining_intensity)
    end
end

using GLMakie
image(img)

You’d have to swap out rect_positions and rect_colors for the vertex coordinates and colors of you mesh and write a get_interpolated_color() function that returns the appropriate colors for your mesh. If you have that though, you could also just resample your data into 3D array and run the ray tracing part on the GPU with Makie volume plots.

When it comes to doing this on the GPU, you also have the problem that this doesn’t fit in the normal rendering pipeline. Regardless of which Graphics API you use, the normal rendering pipeline will have a vertex shader which acts on one vertex each, then the API will combine that data into faces and figure out which pixels are covered by them and then run a fragment shader to draw those pixels.
Makie essentially draws the bounding box of the volume plot and does ray tracing inside it for the fragment shader. To do the ray tracing you’ll need all volume data. If it’s on a regular cartesian grid you can simply put it in a 3D texture. That’ll also handle the interpolation for you. For anything else you’ll need to do the interpolation yourself and pack and unpack it yourself. (This could be a bit simpler with compute shaders but Makie doesn’t support them and you’d still need to handle data passthrough and interpolation.)

What stops us from doing this is mainly that it’s a lot of effort for something very specific.

[raman_kumar]

@ffreyer

I think the main thing you need is a function that tells you the color which should go at a given position.

I have found a way to get the color value for each vertex position of a mesh element. I took average of color values of all elements adjacent to a vertex as color of that vertex.

@raman_kumar

vertex_color = Float64[]
for i in 1:nelements(g)
    el= element(g,i)
    bd = Boundary{3,0}(topology(g))(i)	#Ntuple of vertex ids of i-th element
    for id in bd
		adels= Coboundary{0,3}(topology(g))(id)	#Ntuple of element ids adjacent to this vertex
		clct = collect(Element_color[(elem2cart(topology(g),idx))...] for idx in adels)	#position of i-th element in grid	
		push!(vertex_color, sum(clct)/length(clct))
	end
	@inline viz!(el, color =vertex_color, alpha=0.2)
	#@show clr
end

Now i have color value for each vertex in mesh element.

• Can it be done that i subdivide/tesselate these Hexahedron elements into small cuboids and interpolate these colors to cuboids vertices and get volume plot of these cuboids and then stitch together ? My mesh has these two types of many Hexahedrons.
  
  
  for element(g,1)
  
  
  for element(g,2)

[raman_kumar]

@juliohm

• How to find cuboid that can occupy the maximum or large volume inside the elements(Hexahedrons) of discretized dscel?
• I have taken vertex color of elements of RectilinearGrid g as average of color of elements surrounding that vertex.
• Color at the vertices of these cuboids can be found by interpolating the vertex color clr of that element of RectilinearGrid g.
• After that i will volume!() plot these cuboids and have full stitched plot.

using Meshes, CoordRefSystems ,GLMakie

r =  Float32[1.6, 1.754088, 1.9230156, 2.1082115, 2.311243, 2.5338273, 2.7778475, 3.0453684, 3.3386526, 3.6601818, 4.012676, 4.3991165, 4.8227735, 5.287231, 5.7964177, 6.354642, 6.9666257, 7.637547, 8.373081]
θ =  Float32[0.0, 0.09817477, 0.19634955, 0.2945243, 0.3926991]
ϕ =  Float32[0.0, 0.3926991, 0.7853982, 1.1780972, 1.5707964, 1.9634954, 2.3561945, 2.7488935, 3.1415927, 3.5342917, 3.9269907, 4.3196898, 4.712389, 5.105088, 5.497787, 5.8904862, 6.2831855]

g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ)

clr = [-129.9932861328125, -123.21377563476562, -140.03717041015625, -145.09832763671875, -113.99639892578125, -110.56198120117188, -125.68373107910156, -126.87786865234375]
el = element(g,1)
dscel = discretize(el, RegularDiscretization(10,10))
viz(dscel, alpha=0.5)

discretised dscel

[raman_kumar]

1.In this code:

• I am trying to find a boundingbox for an element of mesh and create a RegularGrid of dims = (10, 10, 10) on that box.
• I will then keep only those cuboids of this RegularGrid that are inside the element.
• I already have color values of element’s vertices. Then I will do interpolation of color values on these cuboid`s vertices and get volume plot.

using Meshes, CoordRefSystems ,GLMakie, LinearAlgebra, Unitful

r =  Float32[1.6, 1.754088, 1.9230156, 2.1082115, 2.311243, 2.5338273, 2.7778475, 3.0453684, 3.3386526, 3.6601818, 4.012676, 4.3991165, 4.8227735, 5.287231, 5.7964177, 6.354642, 6.9666257, 7.637547, 8.373081]
θ =  Float32[0.0, 0.09817477, 0.19634955, 0.2945243, 0.3926991]
ϕ =  Float32[0.0, 0.3926991, 0.7853982, 1.1780972, 1.5707964, 1.9634954, 2.3561945, 2.7488935, 3.1415927, 3.5342917, 3.9269907, 4.3196898, 4.712389, 5.105088, 5.497787, 5.8904862, 6.2831855]

g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ)
hxdrn = element(g,1)
spherpts = coords.(vertices(hxdrn))
cartpts =  convert.(Cartesian, spherpts)
points = [Meshes.Point(Float64.(ustrip.((c.x, c.y, c.z)))...) for c in cartpts]
sel =Hexahedron(tuple(points...))
bb = Meshes.boundingbox(sel)
dbb = RegularGrid(bb.min, bb.max, dims = (10, 10, 10))
cuboids = []
for cl in nelements(dbb)
	c = vertex(element(dbb,cl),1)
	if c ∈ sel     # keep cuboid if inside the element
		push!(cuboids, cl)
	end
end

viz!(dbb, alpha=0.4, showsegments = true, pointsize =20, pointcolor =:Green)
viz!(sel)

But this code c∈sel is giving error.

ERROR: MethodError: no method matching ntuple(::Meshes.var"#761#762"{Meshes.Point{…}}, ::Quantity{Float64, 𝐋, Unitful.FreeUnits{…}})
The function ntuple exists, but no method is defined for this combination of argument types..

---

2.Another problem is that

src/predicates/issubset.jl

using Meshes
h = Hexahedron((0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1))
h2 = Hexahedron((0,0,0),(2,0,0),(2,2,0),(0,2,0),(0,0,2),(2,0,2),(2,2,2),(0,2,2))
h ⊆ h2

gives following error.

ERROR: MethodError: no method matching isconvex(::Hexahedron{CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}, 𝔼{3}})
The function isconvex exists, but no method is defined for this combination of argument types.

[raman_kumar]

using GLMakie, Meshes
l =Meshes.Point(0,0,1.6)
h =Meshes.Point(0.172,0.065,1.75)
dbb = RegularGrid(l, h, dims = (10, 10, 10))
viz(dbb, color=1:nelements(dbb), alpha=0.4, showsegments=true)

work fine while

viz(dbb, color=1:nvertices(dbb), alpha=0.4, showsegments=true)

gives error

ERROR: not implemented
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35