Mesh3D Only Plotting Half

Nkrause · January 3, 2021, 6:13am

Hello, I’m trying to plot a 3d surface. I’ve used a marching cubes algorithm to generate a tuple of points. When I plot them in scatter() I get the following image, which is what I want and expect.

scatter([a for (a,_,_) in points],[b for (_,b,_) in points], [c for (_,_,c) in points] )

However when I try and use Mesh3d to generate the mesh surface I obtain the following.

mesh3d([a for (a,_,_) in points], [b for (_,b,_) in points], [c for (_,_,c) in points])

And only half of the points appear to be plotted. Is there a reason for this? I’m using Plots with a PlotlyJS backend. The entire setup is running in a pluto notebook, and I can append the code if someone wants to see the whole thing. The MWE is immediately below, along with a .csv version of the points tuple variable beneath it that one should be able to copy and paste!

using Plots
	Pkg.add("PlotlyJS")
Pkg.add("QuadGK")
	using QuadGK
	plotly()
Pkg.add("Meshing")
	using Meshing

function unknotintegralfunction(a,b,c,t)
	1/(sqrt(((a.-cos.(t))^2).+((b.-sin.(t))^2).+c^2)) #using .math means element wise
end

function unknotpotential(a,b,c,A) 
	ans, err = quadgk(t->unknotintegralfunction(a,b,c,t),0,(2*pi),order=quadorder) 
	# we'll use the default errors to start
	A[floor(Int,(10*a+21)),floor(Int,(10*b+21)),floor(Int,(10*c+21))]=ans
	return [a,b,c,ans,err]
end

inflateunknot(f,N,p,A)=[f(x,y,z,A) for x in -1*N:p:N, y in -1*N:p:N, z in -1*N:p:N if z ≠ 0 || (x*x +y*y)≠ 1 ]
A = zeros(41,41,41)
ans = inflateunknot(unknotpotential,2,0.1,A)
c=A[20,20,20]+0.5
points,surfaces = isosurface(A,MarchingCubes(iso=c), origin = Vec(1,1,1), widths = Vec(4.0,4.0,4.0)) 

scatter([a for (a,_,_) in points],[b for (_,b,_) in points], [c for (_,_,c) in points] )
mesh3d([a for (a,_,_) in points], [b for (_,b,_) in points], [c for (_,_,c) in points])

dpsanders · January 3, 2021, 6:36am

Please provide a MWE (minimal working example that people can copy and paste to reproduce the problem, including all data and using calls). See Please read: make it easier to help you. In particular you haven’t said which of the many plotting libraries you are using.

Nkrause · January 3, 2021, 6:42am

My apologies, I’ve edited the original question to include the code generating the data to be plotted and the plots themselves along with details about the plotting libraries.

dpsanders · January 3, 2021, 6:51am

Thanks! You might want to reduce that down to normal Julia code, rather than a Pluto notebook, and reduce it down to the bare minimum that reproduces the problem

Nkrause · January 3, 2021, 6:58am

Ah, will do! Was worried there might be something in the pluto notebook causing the issue itself. I’ll update tomorrow with a stripped down example.

Nkrause · January 3, 2021, 6:53pm

I’ve edited the question to include an MWE, I followed the guide you linked to to output a copy-pastable version of my points variable. However it’s too many points to just copy paste and include in the example, not too sure how I should go about providing the variable values.

dpsanders · January 3, 2021, 7:23pm

I don’t see the points variable. It’s fine / better to generate that using code, not copy paste the data.

Nkrause · January 3, 2021, 8:14pm

Hello again! I appreciate your patience, I’ve updated the code with what I think is enough to generate the dataset. It works on my machine so hopefully it reproduces on yours!

rafael.guerra · January 3, 2021, 9:46pm

The code provided would not run in my machine “as is”. The edited code below runs. Need to set alphahull in mesh3d.

using Plots, PlotlyJS, QuadGK, Meshing, StaticArrays
plotly()

function unknotintegralfunction(a,b,c,t)
	1/(sqrt( (a - cos(t))^2 + (b - sin(t))^2 + c^2 ))
end

function unknotpotential(a,b,c,A) 
	ans0, err = quadgk(t->unknotintegralfunction(a,b,c,t),0,2π,order=quadorder) 
	A[floor(Int,(10*a+21)),floor(Int,(10*b+21)),floor(Int,(10*c+21))] = ans0
	return [a,b,c,ans0,err]
end

inflateunknot(f,N,p,A) = [f(x,y,z,A) for x in -N:p:N, y in -N:p:N, z in -N:p:N if z ≠ 0 || (x*x +y*y)≠ 1 ]

quadorder = 20
A = zeros(41,41,41)
ans = inflateunknot(unknotpotential,2,0.1,A)
c = A[20,20,20] + 0.5
points,surfaces = isosurface(A,MarchingCubes(iso=c), origin = SVector(1,1,1), widths = SVector(4.0,4.0,4.0))

scatter([a for (a,_,_) in points],[b for (_,b,_) in points], [c for (_,_,c) in points] )

PlotlyJS.plot(
    PlotlyJS.mesh3d(
        x = [a[1] for a in points],
        y = [a[2] for a in points],
        z = [a[3] for a in points],
        alphahull = 5
    )
)

Nkrause · January 3, 2021, 9:54pm

Thank you Rafael, I appreciate the time and effort. I’ll take a look at your code and see where I made some mistakes in my own. With regards to the mesh, is it possible that I’m feeding too many points to Mesh3d? Would the algorithm have an easier time with a sparser dataset?

dpsanders · January 3, 2021, 11:01pm

It looks to me like the points are not being generated in the right order for the meshing algorithm to parse.

What happens if you replace plotly by plotlyjs? IIUC, the former is remote and the latter is local.

Nkrause · January 4, 2021, 12:56am

switching to PlotlyJS creates a number of errors that take me to the following webpage.

It discusses Jupyter notebooks but makes no mention of Pluto. Installing the dependencies as below unfortunately does not resolve the error.

begin
	Pkg.add("WebIO")
	using WebIO
end

Is the PlotlyJS backend not supported for Pluto? If so I might stick with a 3d scatter plot for now rather than worrying about trying to get mesh3d working in a notebook.

dpsanders · January 4, 2021, 5:37am

You could try with WGLMakie.

rafael.guerra · January 4, 2021, 7:01am

The default alphahull parameter was not appropriate. Please check my edited post above.

Nkrause · January 4, 2021, 4:40pm

Thank you Rafael!

empet · January 6, 2021, 3:35pm

@Nkrause WebIO is not enabled in Pluto: https://github.com/fonsp/Pluto.jl/issues/299