Does delete! or Observable @view work with Axis3?

[g2g]

As part of an essay (including a listing of the 3D PGA reference implementation ripga3d.jl in the appendix) proposing that Julia + Makie + projective geometric algebra (PGA) are an excellent combination for implementing intense geometry applications, I’m trying to port bivector.net’s 3D slicer application (e.g., for 3D printing) from JavaScript to Julia and Makie.

The moving horizontal slice looks good in the 2D animation but not in the 3D animation (because the previous horizontal slice in the 3D animation is not erased). I wasn’t able to get the Observable @view or the delete!() to erase the previous horizontal slice in Axis3. Any suggestions?

 # tslice.jl : test tetrahedron slice
#   to test Julia port of ganja.js 3D slicing example
#
# The test tetrahedron has three faces (the base face
# is empty), all edges of length 2.0, and a peak that
# is directly above (i.e., z offset) the origin. All the
# test tetrahedron's vertices and faces are defined,
# according to the wavefront 3D object file format,
# by seven lines of text (within the following 
# multiline comment): 
#=
v 1.0 -0.5773502692 0.0
v 0.0 1.15470053838 0.0
v -1.0 -0.5773502692 0.0
v 0.0 0.0 1.632993162
f 1 2 4
f 2 3 4
f 3 1 4
=#

using GLMakie, GLMakie.FileIO
using GeometryBasics
include("ripga3d.jl")

function tslice()
	# load and scale (by 13) faces of bunny object
	F = load("tetrahedron.obj")
	nFace = length(F)
	
	# define cut line segment buffer & references to it
	LS = fill(NaN32, 3*nFace, 3) 
	nLS = 0
	nPrevLS = 0
	s2d::Any = 0
	s3d::Any = 0

	# precalculate face vector and triangle side bivectors
	FS3 = Array{Float32,2}(undef,64,nFace)
	for iF = 1:nFace # for each triangle face in bunny object
		P1 = point(F[iF][1][1],F[iF][1][2],F[iF][1][3])
		P2 = point(F[iF][2][1],F[iF][2][2],F[iF][2][3])
		P3 = point(F[iF][3][1],F[iF][3][2],F[iF][3][3])
#		P12 = ga"P1∨P2"
#		P23 = ga"P2∨P3"
#		P31 = ga"P3∨P1"
#		FS3[ 1:16,iF] = ga"P1∨P2∨P3" # face
#		FS3[17:32,iF] = ga"FS3[1:16,iF]·P12" # side 1
#		FS3[33:48,iF] = ga"FS3[1:16,iF]·P23" # side 2
#		FS3[49:64,iF] = ga"FS3[1:16,iF]·P31" # side 3
		P12 = P1 & P2
		P23 = P2 & P3
		P31 = P3 & P1
		FS3[ 1:16,iF] = P1 & P2 & P3 # face
		FS3[17:32,iF] = FS3[1:16,iF] | P12 # side 1
		FS3[33:48,iF] = FS3[1:16,iF] | P23 # side 2
		FS3[49:64,iF] = FS3[1:16,iF] | P31 # side 3
	end

	# initialize figures
	fig = Figure(resolution = (1000, 500))
	ax3d = Axis3(fig[1,1],
		elevation = pi/16,
		azimuth = -3*pi/4,
		viewmode = :fit,
		aspect = (1,1,1))
	ax2d = Axis(fig[1,2],
		limits = (-1,1, -0.75,1.25),
		yticks = (-0.75:0.5:1.25),
		title = "cross section at horizontal slice plane",
		aspect = 1)
	m = mesh!(ax3d, F)
	
	# perform slicing
	zMax = 1.632993162
	nFrame = 360
	record(fig, "tslice.mp4", 1:nFrame) do frame
		ax3d.azimuth[] = -3pi/4 - 2pi*(frame-1)/nFrame
		zCut = zMax - abs(zMax - 2*zMax*(frame-1)/nFrame)
		cut = zCut*e0 - e3
		
		# erase old slice
		if nPrevLS > 0
			delete!(ax2d, s2d)
# NOTE: currently, delete! on ax3d causes the following error:
# ERROR: MethodError: no method matching delete!(::Axis3, ::Lines{Tuple{Vector{Point{3, Float32}}}}) 
#			delete!(ax3d, s3d)
			nPrevLS = 0
		end
		
		for iF = 1:nFace
			N = FS3[1:16,iF]
			side1 = FS3[17:32,iF]
			side2 = FS3[33:48,iF]
			side3 = FS3[49:64,iF]
#			lc = ga"N∧cut"
#			P1 = ga"lc∧side1" # line^plane->point
			lc = N ^ cut
			P1 = lc ^ side1 # line^plane->point
			nPoint = 0
			if P1[15] != 0
				nPoint += 1
				iLS = 3*nLS + nPoint
				LS[iLS,:] = [P1[14] P1[13] P1[12]] ./ P1[15]
			end
#			P2 = ga"lc∧side2" # line^plane->point
			P2 = lc ^ side2 # line^plane->point
			if P2[15] != 0
				nPoint += 1
				iLS = 3*nLS + nPoint
				LS[iLS,:] = [P2[14] P2[13] P2[12]] ./ P2[15]
			end
#			P3 = ga"lc∧side3" # line^plane->point
			P3 = lc ^ side3 # line^plane->point
			if P3[15] != 0
				nPoint += 1
				iLS = 3*nLS + nPoint
				LS[iLS,:] = [P3[14] P3[13] P3[12]] ./ P3[15]
			end
			if nPoint == 2
				nLS += 1
			elseif nPoint == 3
				iLS = 3*nLS + 3
				LS[iLS,:] .= NaN32
			end
		end # for each face

		# plot new slice
		if nLS > 0
			slice2d = @view LS[1:3*nLS,1:2]
			s2d = lines!(ax2d, slice2d, color=:black)
			slice3d = @view LS[1:3*nLS,:]
			s3d = lines!(ax3d, slice3d, color=:white)
			
			nPrevLS = nLS
			nLS = 0 # reset line segment buffer
		end
	end # for each video frame
end # tslice()

[jules]

Do not execute lines! in every frame, instead do it once at the beginning with an Observable that holds the first set of points (or an empty vector). Then compute your new points in each frame and update the observable. The steps are roughly like this:

points = Observable(Point2f[])
lines!(ax, points)
for frame in frames
    points[] = calculate_new_points(frame)
end

It doesn’t have to be a vector of Point2f but there might be an edge case if you start with an empty matrix, not sure. Maybe it works if you make it a 0 by 2 matrix, but I think the conversion pipeline has to recognize what kind of point the matrix stores, 2d or 3d.

If you want to keep using your LS array, you can wrap it in an Observable, mutate it like you’re already doing it, and then call notify on the observable when you’re done mutating. The inputs of the two lines can then be derived from this basis using the logic you have written out already. An example:

LS_obs = Observable(LS)
slice_2d = @lift @view $LS_obs[1:3*nLS, 1:2]

# or you could write equivalently
slice_2d = lift(LS_obs) do array
    @view array[1:3*nLS, 1:2]
end

lines!(ax2d, slice_2d)
for frame in frames
    mutate_somehow!(LS)
    notify(LS_obs) # now slice_2d will update with the new data and the lines plot will, too
end

[g2g]

@jules

An example:

Thanks for giving that example. (I would never have thought of adding a second layer of Observable.)

Removing the complexity of projective geometric algebra, here is an MWE.

# tmslice.jl : test tetrahedron slice
#   to test Makie 2d and 3d plotting
#
# The test tetrahedron has three faces (the base face
# is empty), all edges of length 2.0, and a peak that
# is directly above (i.e., z offset) the origin. All the
# test tetrahedron's vertices and faces are defined,
# according to the wavefront 3D object file format,
# by seven lines of text (within the following 
# multiline comment): 
#=
v 1.0 -0.5773502692 0.0
v 0.0 1.15470053838 0.0
v -1.0 -0.5773502692 0.0
v 0.0 0.0 1.632993162
f 1 2 4
f 2 3 4
f 3 1 4
=#

using GLMakie, GLMakie.FileIO

function tmslice()
	# load mesh of 3D object
	F = load("tetrahedron.obj")
	nFace = length(F)
	
	# initialize figures
	fig = Figure(resolution = (1000, 500))
	ax3d = Axis3(fig[1,1],
		elevation = pi/16,
		azimuth = -3*pi/4,
		viewmode = :fit,
		aspect = (1,1,1))
	ax2d = Axis(fig[1,2],
		limits = (-1,1, -0.75,1.25),
		yticks = (-0.75:0.5:1.25),
		title = "horizontal cross section",
		aspect = 1)
	m = mesh!(ax3d, F)
	
	# plot first slice
	nFrame = 360
	zMax = 1.632993162
	ZLR = LinRange(-zMax, zMax, nFrame+1)
	Z = zMax .- abs.(ZLR)
	rMax = 1.15470053838
	R = abs.(LinRange(-rMax, rMax, nFrame+1))
	THETA = LinRange(pi/2, 5*pi/2, 4)
	LS = fill(NaN32, 15, 3) # Line Segment buffer 
	LS[1:4,1] = R[1] .* cos.(THETA)
	LS[1:4,2] = R[1] .* sin.(THETA)
	LS[1:4,3] .= Z[1]
	LS_obs = Observable(LS)
	slice2d = @lift @view $LS_obs[1:5, 1:2]
	slice3d = @lift @view $LS_obs[1:5, :]
	lines!(ax2d, slice2d, color=:black)
	lines!(ax3d, slice3d, color=:black, linewidth=5)
	fig
	
	# record video of slicing
	record(fig, "tmslice.mp4", 1:nFrame) do iFrame
		ax3d.azimuth[] = -3pi/4 - 2pi*(iFrame-1)/nFrame
		LS[1:4,1] = R[iFrame] .* cos.(THETA)
		LS[1:4,2] = R[iFrame] .* sin.(THETA)
		LS[1:4,3] .= Z[iFrame]
		notify(LS_obs)
	end
end # tmslice()

I was concerned about plotting a variable number of line segments in the slices, but that works too!

# tm2slice.jl : test tetrahedron slice
#   to test Makie 2d and 3d plotting when the
#	number of line segments per slice varies
#
# The test tetrahedron has three faces (the base face
# is empty), all edges of length 2.0, and a peak that
# is directly above (i.e., z offset) the origin. All the
# test tetrahedron's vertices and faces are defined,
# according to the wavefront 3D object file format,
# by seven lines of text (within the following 
# multiline comment): 
#=
v 1.0 -0.5773502692 0.0
v 0.0 1.15470053838 0.0
v -1.0 -0.5773502692 0.0
v 0.0 0.0 1.632993162
f 1 2 4
f 2 3 4
f 3 1 4
=#

using GLMakie, GLMakie.FileIO

function tm2slice()
	# load mesh of 3D object
	F = load("tetrahedron.obj")
	nFace = length(F)
	
	# initialize figures
	fig = Figure(resolution = (1000, 500))
	ax3d = Axis3(fig[1,1],
		elevation = pi/16,
		azimuth = -3*pi/4,
		viewmode = :fit,
		aspect = (1,1,1))
	ax2d = Axis(fig[1,2],
		limits = (-1,1, -0.75,1.25),
		yticks = (-0.75:0.5:1.25),
		title = "horizontal cross section",
		aspect = 1)
	m = mesh!(ax3d, F)
	
	# plot first slice
	nFrame = 360
	zMax = 1.632993162
	ZLR = LinRange(-zMax, zMax, nFrame+1)
	Z = zMax .- abs.(ZLR)
	rMax = 1.15470053838
	R = abs.(LinRange(-rMax, rMax, nFrame+1))
	THETA = LinRange(pi/2, 5*pi/2, 4)
	LS = fill(NaN32, 15, 3) # Line Segment buffer 
	LS[1:4,1] = R[1] .* cos.(THETA)
	LS[1:4,2] = R[1] .* sin.(THETA)
	LS[1:4,3] .= Z[1]
	nTri = 1
	LS_obs = Observable(LS)
	slice2d = @lift @view $LS_obs[1:5*nTri, 1:2]
	slice3d = @lift @view $LS_obs[1:5*nTri, :]
	lines!(ax2d, slice2d, color=:black)
	lines!(ax3d, slice3d, color=:black, linewidth=5)
	fig
	
	# record video of slicing
	record(fig, "tm2slice.mp4", 1:nFrame) do iFrame
		ax3d.azimuth[] = -3pi/4 - 2pi*(iFrame-1)/nFrame
		
		iSet = div(iFrame, 8, RoundDown)
		nTri = mod(iSet, 3) + 1
		for iTri = 1:nTri
			i = 5*(iTri-1) + 1 # starting index of triangle
			s = 1 + (iTri-1)*0.10 # scale of triangle
			LS[i:i+3,1] = (R[iFrame]*s) .* cos.(THETA)
			LS[i:i+3,2] = (R[iFrame]*s) .* sin.(THETA)
			LS[i:i+3,3] .= Z[iFrame]
		end
		notify(LS_obs)
	end
end # tm2slice()