[PYTHON] Using the 3D plot library MayaVi from Julia

[Caution] This article was published in December 2016. The code below is for Julia 0.6.

Julia advent calendar 2016 This is the article on the 6th day.

Let's try using a library called MayaVi from Julia, which has recently been provided with a conda package.

What is MayaVi

MayaVi is a package for 3D plots. The outline is as follows.

--A 3D plot package for Python developed by a scientific computing consultant Enthought. It is included in the integrated environment Enthough Canopy distributed by the company. It is also included in the Scientific and Scientific Calculation Python Package Python (x, y) for Windows. --The drawing engine is a library The Visualization Toolkit (VTK) written in C ++. vtk is suitable for visualization of large-scale and complex data. --MayaVi pronounces "Ma-ya-vee". In Sanskrit, it is equivalent to the English word "magical".

A library with a similar personality is ParaView (http://www.paraview.org). It's aimed at parallel decentralization and is probably more famous than MayaVi. I started using MayaVi first and have had a lot of trouble installing ParaView. The usability is similar, and MayaVi is recommended for personal projects.

Install mayavi's conda package

The handling of conda packages in Julia was explained in yesterday's article. (Add conda package from Julia) Use Conda.add () to add the mayavi package to the miniconda environment installed by Julia.

julia> using Conda

julia> Conda.add("mayavi")
Using Anaconda Cloud api site https://api.anaconda.org
Fetching package metadata .......
Solving package specifications: ..........

# All requested packages already installed.
# packages in environment at /Users/hs/.pyenv/versions/anaconda-2.4.0/envs/conda_jl:
...The following is omitted.

If you want to add mayavi to your anaconda environment, you can just hit the conda command from the shell (command line, terminal).

$ conda install mayavi -n conda_jl
Using Anaconda Cloud api site https://api.anaconda.org
Fetching package metadata .............
Solving package specifications: ..........

Package plan for installation in environment /Users/hs/.pyenv/versions/anaconda-2.4.0/envs/conda_jl:
...The following is omitted.

MayaVi functions gallery

First, try the example in MayaVi functions gallery.

julia> using PyCall

julia> @pyimport mayavi.mlab as mlab

julia> mlab.test_plot3d()
PyObject <mayavi.modules.glyph.Glyph object at 0x337fc95f0>

julia> mlab.show()

Source file above: plot3d () https://gist.github.com/d4578e1f4d22bbfa0f418f0caff239c7

Execution result (screenshot)

mayavi.mlab is a tool for simply plotting Python's numpy.array. If you know matplotlib, can you guess that it corresponds to matplotlib.pyplot?

When you start mlab.show (), the actual drawing will be done. At the top of the 3D plot, there is a menu bar with several icons, which allows you to rotate and scale the shape. The feature around here is called traits (no good translation found).

Now, the other eight examples in the Mayavi Functions gallery were also displayed successfully. Put the source file in gist. --Points (sphere) points3d (): https://gist.github.com/a3b0b95b58a42b7b010aea779536e6f3 --Image imshow (): https://gist.github.com/326f9b72a56f6187f0d62b7e9e3ca4ed --High surface surf (): https://gist.github.com/53d022bfcb9b268ea829f23f9e2aae4e --Contour line contour_surf (): https://gist.github.com/bc51e4baf2f09df352e8019e37f34bdf --Curved surface mesh (): https://gist.github.com/8fc119d17a4bd44cefdbb7439a3a028f --Bar chart barchart (): https://gist.github.com/89d6a5f3a6a321fd099d1b04e9cc7ea7 --Triangular mesh triangular_mesh (): https://gist.github.com/3c44621b47282921458803c2c15c1aaf --Isosurface contour3d (): https://gist.github.com/002372eecca7cc831764ea731acaf155

Execution result of test_mesh () (screenshot)

surface_from_irregular_data

Below are four example programs that provide tips and related tips for porting a Python Mayavi program to Julia. It's a program called surface_from_irregular_data.py.

--Python source: http://docs.enthought.com/mayavi/mayavi/auto/example_surface_from_irregular_data.html --Julia source below: https://gist.github.com/2cbe32c7b10a2b1c90657e06e9ba4237

function f(x, y)
  exp(-(x .^ 2 + y .^ 2))
end

srand(12345)
xs = 4.0 * (rand(500) - 0.5)
ys = 4.0 * (rand(500) - 0.5)
zs = f(xs, ys)

using PyCall
@pyimport mayavi.mlab as mlab
mlab.figure(1, fgcolor=(0, 0, 0), bgcolor=(1, 1, 1))

# Visualize the points
pts = mlab.points3d(xs, ys, zs, zs, scale_mode="none", scale_factor=0.2)

# Create and visualize the mesh
mesh = mlab.pipeline[:delaunay2d](pts)
surf = mlab.pipeline[:surface](mesh)

mlab.view(47, 57, 8.2, (0.1, 0.15, 0.14))
mlab.show()

Execution result (screenshot)

Extend an irregular 2D point (x, y) to a 3D point (x, y, z) using the function z = f (x, y). Draw them as points (spheres) (points3d). It also interpolates 3D points (delaunay2d) and draws its surface (surf). The process of adding some processing to the data is called a pipeline. Julia source is almost the same as Python source. I paid attention to the following points. --Python's mlab.pipeline.delaunay2d etc. cannot be called as it is, but it is calledmlab.pipeline [: delaunay2d]. --In function, corrected the exponentiation of each element to. ^ .

triangular_mesh

Next, try porting the program inside test_triangular_mesh to Julia.

--Python source: http://docs.enthought.com/mayavi/mayavi/auto/mlab_helper_functions.html?highlight=triangular%20mesh#triangular-mesh --Julia source below: https://gist.github.com/bc98e07f807ebc03f75837cb76117b79

# An example of a cone, ie a non-regular mesh defined by its triangles.
n = 8
t = linspace(-pi, pi, n)
xy = exp(im * t)
x = real(xy)
y = imag(xy)
z = zeros(n)

triangles = [ (0, i, i + 1) for i in 1:n-1 ]
unshift!(x,0.0)
unshift!(y,0.0)
unshift!(z,1.0)
t=collect(t)
unshift!(t,0.0)

using PyCall
@pyimport mayavi.mlab as mlab
mlab.triangular_mesh(x, y, z, triangles, scalars=PyObject(t))
mlab.show()

Execution result (screenshot) * The figure is reduced / rotated:

mlab.triangular_mesh is an instruction to draw (multiple) triangles. As arguments, give the coordinates of the point to be the vertex and the number of the vertex of the triangle. Functions imported with @pyimport pass array arguments as numpy.array to Python. If you want to pass it as a regular list, wrap it in PyObject (). Note that array subscripts start at 1 in Julia, compared to 0 in Python. The point numbers of the vertices indicated by triangles will be offset by one in Julia.

Here are some other Julia tips. --linspace (start, end, n) creates an arithmetic progression of n elements. To get the numbers. Use collect. --ʻIm is an imaginary unit. --zeros (n) creates an array of Float64withnelements. The value is0.0. --ʻUnshift! (v, e) adds the element ʻeto the beginning of the arrayv. Append at the end with push! (V, e) or ʻappend! (V, e). The array v is destroyed for any instruction. (It's ! at the end of the instruction)

spherical_harmonics

Now let's port ʻexample_spherical_harmonics.py` to Julia.

--Python source: http://docs.enthought.com/mayavi/mayavi/auto/example_spherical_harmonics.html --Julia source below: https://gist.github.com/60d97f20a19d820e299c253612728907

# phi, theta = np.mgrid[0:pi:101j, 0:2 * pi:101j]
phi   = [ u1 for u1 in linspace(0,pi,101), v1 in linspace(0,2*pi,101) ]
theta = [ v1 for u1 in linspace(0,pi,101), v1 in linspace(0,2*pi,101) ]

r = 0.3
x = r * sin(phi) .* cos(theta)
y = r * sin(phi) .* sin(theta)
z = r * cos(phi)

using PyCall
@pyimport mayavi.mlab as mlab
@pyimport scipy.special as spe

mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0,0,0), size=(400, 300))
mlab.clf()
# Represent spherical harmonics on the surface of the sphere
for n in 1:6-1, m in 0:n-1
  s = real( spe.sph_harm(m, n, theta, phi) )
  mlab.mesh(x - m, y - n, z, scalars=s, colormap="jet")
  s[s .< 0] *= 0.97
  s /= maximum(s)
  mlab.mesh(s .* x - m, s .* y - n, s .* z + 1.3, scalars=s, colormap="Spectral" )
end

mlab.view(90, 70, 6.2, (-1.3, -2.9, 0.25))
mlab.show()

Execution result (screenshot)

The calculation of the spherical harmonics is done by calling Python's scipy.special.sph_harm. Then draw a curved surface (mesh). You can draw well. Atomic orbital azimuth quantum numbers s, p, d ... that appear in physical chemistry. Some tips. --numpy.mgrid is a function that creates a direct product of 2D coordinates or more, but Julia does not have it. But given its implications, it's easy to implement with comprehension. --Julia's for statement can write multiple loops. What is written on the right is the inner loop. --Python's range (n) is equivalent to Julia's 0: n-1. Python's range (m, n) is equivalent to Julia's m: n-1.

simple structured grid

The last example is a little tricky. Let's port ʻexample_simple_structured_grid.py` to Julia.

--Python source: http://docs.enthought.com/mayavi/mayavi/auto/example_simple_structured_grid.html#example-simple-structured-grid --Julia Source: https://gist.github.com/d0c2c1c9a6fdb04c258f1961d559ee3b

# x, y, z = mgrid[1:6:11j, 0:4:13j, 0:3:6j]
x = [ x1 for x1 in linspace(1.0,6.0,11), y1 in linspace(0.0,4.0,13), z1 in linspace(0.0,3.0,6) ]
y = [ y1 for x1 in linspace(1.0,6.0,11), y1 in linspace(0.0,4.0,13), z1 in linspace(0.0,3.0,6) ]
z = [ z1 for x1 in linspace(1.0,6.0,11), y1 in linspace(0.0,4.0,13), z1 in linspace(0.0,3.0,6) ]
base=x[:,:,1] + y[:,:,1]
for i in 1:size(z)[3]
  z[:,:, i] = base[:,:] * 0.25 * (i-1)
end

pts=zeros(Float64, tuple(size(z)...,3))
pts[:,:,:,1] = x
pts[:,:,:,2] = y
pts[:,:,:,3] = z
scalars1 = x .* x + y .* y + z .* z

vectors1=zeros(Float64, tuple(size(z)...,3))
vectors1[:,:,:,1] = (4.0 - y * 2.0)
vectors1[:,:,:,2] = (x * 3.0 - 12.0)
vectors1[:,:,:,3] = sin(z * pi)

# pts = pts.transpose(2, 1, 0, 3).copy()
# pts= permutedims(pts, [3,2,1,4] )
# pts= reshape(pts, ( prod(size(pts)[1:3]), 3))

# vectors1= permutedims(vectors1, [3,2,1,4] )
# vectors1= reshape(vectors1, ( prod(size(vectors1)[1:3]), 3))

using PyCall
@pyimport tvtk.api as tvtk_api

# Create the dataset.vec
sg=tvtk_api.tvtk[:StructuredGrid](dimensions=size(x),points=pts)
sg[:point_data][:scalars] = vec(scalars1)
sg[:point_data][:scalars][:name] = "temperature"
sg[:point_data][:vectors] = vectors1
sg[:point_data][:vectors][:name] = "velocity"

@pyimport mayavi.mlab as mlab

d = mlab.pipeline[:add_dataset](sg)
gx = mlab.pipeline[:grid_plane](d)
gy = mlab.pipeline[:grid_plane](d)
gy[:grid_plane][:axis] = "y"
gz = mlab.pipeline[:grid_plane](d)
gz[:grid_plane][:axis] = "z"
iso = mlab.pipeline[:iso_surface](d)
iso[:contour][:maximum_contour] = 75.0
vec1 = mlab.pipeline[:vectors](d)
vec1[:glyph][:mask_input_points] = true
vec1[:glyph][:glyph][:scale_factor] = 1.5

mlab.show()

--Execution result (screenshot)

--Python execution result (screenshot)

Create a 3D grid (StructuredGrid). At each point, assign a scalar value (scalar1) and a vector value (vector1) (point_data). Pour this into the pipeline to draw the iso_surface and vectors. It also draws planes with x = 0, y = 0, z = 0 (grid_plane). When giving coordinate data to vtk, give an array so that x first, then y, and finally z move (column major). Since python-numpy is row-major, the original python source swaps the storage order. On the other hand, Julia is a column major, so you can leave it as it is. (Reference Row-major order and Column-major order) By the way, you can use Python's numpy.transpose instruction to transpose the axes of a multidimensional array. Julia's transpose only swaps rows and columns in a matrix (two-dimensional array), not for multidimensional arrays. Use permutedims to swap the axes of a multidimensional array in Julia. The following two are equivalent (in Julia, the axes also count from 1):

a = a.transpose(2, 1, 0, 3).copy()

a = permutedims(a, [3,2,1,4] )
#When rewriting an array
permutedims!(a, [3,2,1,4] )

Some tips. Tuple nesting does not expand. Add ... to expand. It is called splat construct.

julia> ((1,2),(3,4))
((1,2),(3,4))

julia> ((1,2)...,(3,4))
(1,2,(3,4))

julia> ((1,2),(3,4)...)
((1,2),3,4)

julia> ((1,2)...,(3,4)...)
(1,2,3,4)

For a multidimensional array ʻa, Python-numpy's ʻa.shape is Julia's size (a).

At the end

So far, I've introduced an example of using Mayavi with Julia in a hurry. In many cases, I looked at the Python source and showed that it could be rewritten almost mechanically. Now, you can also view Mayavi inside Jupyter. This will be introduced in another article. -> I wrote it. [Python, Julia] 3D display in Jupyter-Mayavi library