Drawing parametric surfaces in Plots.jl

General Usage
1

Does anyone know how to make a 3D plot of surfaces which are parametric? I mean e.g. sphere or Möbius strip. Looking at help and my own experiments, function surface draws on rectangles seem to not provide option of giving x and y as 2D matrices. However some old discussions, e.g.

https://stackoverflow.com/questions/44608878/plots-jl-turn-off-axis-and-grid-lines

suggest that this option was available some time ago. Now surface function allows to provide x,y,z as lists, but forcing them to make correct mesh seems to be quite unpredictable.

I’ve done this in Makie, but I am curious if I missed some option in Plots.

Julia v1.3.1, Plots v0.29.8

3 Likes
2

I have this plot recipe which was working last I checked, but only with the two backends. I don’t know how the rendering would work with more complicated surfaces, but cones and spheres work fine with pyplot.

"""
    plot_parametric_surface(F, [xlim=(-5,5)], [ylim=(-5,5)], [nx=50], [ny=50])
Plot the parameterized surface described by `F:R^2 -> R^3`. This works with `pyplot` and `plotly`, but not `gr`.
A surface `(x,y,f(x,y))` can be directly plotted with `surface` as:

surface(xs, ys, f)

where `xs` and `ys` provide a grid to plot over.
For parametrically described surfaces, the above doesn't work. This function provides an interface:

F(u,v) = [u*cos(v), u*sin(v), u] # a cone
plot_parametric_surface(F, xlim=(0, 1), ylim=(0,2pi))

The values of `xlim` are used as the range to plot over; `nx` specifies the number of points to use. Similarly for `ylim`.
"""
@userplot Plot_Parametric_Surface
@recipe function __(r::Plot_Parametric_Surface; nx=50, ny=50)

    F = first(r.args)

    xlim = get(plotattributes,:xlims, (-5,5))
    ylim = get(plotattributes,:ylims, (-5,5))



    us = range(xlim[1], xlim[2], length=nx)
    vs = range(ylim[1], ylim[2], length=ny)
    ws = unzip(F.(us, vs'))

    seriestype := :surface
    x := ws[1]
    y := ws[2]
    z := ws[3]
    xlims := extrema(ws[1])
    ylims := extrema(ws[2])
    ws

end

unzip(vs::Vector) = Tuple([[vs[i][j] for i in eachindex(vs)] for j in eachindex(vs[1])])
unzip(v,vs...) = unzip([v, vs...])
1 Like
3

This code seems similar to my (failed) approach. However, for me it does not work at all. " Couldn’t process recipe args" and the rest of the error is gibberish. I do not have experience with user recipes, so I am not sure what is wrong here.

4

Sorry, I should have unpacked the recipe. Something like this should work, using the unzip function defined previously:

F(u,v)  =  [cos(u)*sin(v),  sin(u)*sin(v), cos(v)]
vs =  us = range(0,2pi, length=100)
surface(unzip(F.(us, vs')...))
5

It doesn’t work properly. It has the same problem as my attempts - the mesh is off. This is not visible with this many points, but change to lenght = 20 and you will see this is not a sphere, this is some sort of ball shaped zig zag.

6

Yeah, that’s right. Looks terrible with n=5. Sorry. I don’t have thoughts on fix for that, as I had only tested this for simpler surfaces.

7

Thought maybe the plot_parametric_surface function from SymPy might help here, but you can try this and see it also does quite poorly:

using SymPy, PyPlot
@vars u v
SymPy.plot_parametric_surface((sin(u)*sin(v),cos(u)*sin(v), cos(v)), (u,0,2pi), (v,0,2pi)) # ok
SymPy.plot_parametric_surface((sin(u)*sin(v),cos(u)*sin(v), cos(v)), (u,0,2pi), (v,0,2pi), nb_of_points_u=10, nb_of_points_v=10) # not  even a hint of a  sphere...
8

@bachrathyd has some nice 3D plots of implicit functions in https://github.com/bachrathyd/MDBM.jl

IntervalConstraintPropagation.jl can so this by approximating the surface with boxes, which you can view with Makie.jl.

9 
using Plots
plotly() # or pyplot() - gr() does not work for me

# spherical: (radius r, inclination θ, azimuth φ)
X(r,theta,phi) = r * sin(theta) * sin(phi)
Y(r,theta,phi) = r * sin(theta) * cos(phi)
Z(r,theta,phi) = r * cos(theta)

thetas = range(0, stop=pi,   length=50)
phis   = range(0, stop=pi/2, length=50)

xs = [X(1, theta, phi) for theta in thetas, phi in phis] 
ys = [Y(1, theta, phi) for theta in thetas, phi in phis]
zs = [Z(1, theta, phi) for theta in thetas, phi in phis]

surface(xs, ys, zs)

surface

3 Likes
10

For me it doesn’t work in GR as well, in old versions it didn’t throw and error, but the result was wrong, now it throws an error (if I remember correctly, I don’t have Julia 1.3 now). So I guess we should just use new versions of Plots with plotly for parametric surfaces.

11

You can also use Makie, which does this quite nicely. Here is a different, but similar example:

using AbstractPlotting, GLMakie
r1, r2 = 2, 1/2
r(u,v) = ((r1 + r2*cos(v))*cos(u), (r1 + r2*cos(v))*sin(u), r2*sin(v)) # torus
us = vs = range(0, 2pi, length=25)
xs, ys, zs = [[pt[i] for pt in r.(us, vs')] for i in 1:3]
surface(xs, ys, zs)
2 Likes
12

I know, to quote my first post

:slight_smile:

13

Thank you for this example! One thing I don’t understand is why do we need to transpose vs? And do you have any good resource to learn about how to make these plots?

1 Like
14

The transpose is used to create a grid of values as opposed to a vector when Julia’s broadcasting is used. This example might be instructive:

julia> ys = xs = -1:1
-1:1

julia> tuple.(xs, ys)
3-element Array{Tuple{Int64,Int64},1}:
 (-1, -1)
 (0, 0)
 (1, 1)

julia> tuple.(xs, ys')
3×3 Array{Tuple{Int64,Int64},2}:
 (-1, -1)  (-1, 0)  (-1, 1)
 (0, -1)   (0, 0)   (0, 1)
 (1, -1)   (1, 0)   (1, 1)

I’m not sure of a good resource beyond the documentation of the plotting package for these plots specifically. Other plotting packages outside of Julia may create matrices for the u and v using a call to meshgrid, but this isn’t necessary for Makie or Plots, where the u and v values may be a vector and the f(u,v) values a matrix.

3 Likes