7

I'm currently processing strain data with python and use matplotlib (v. 1.5.1) to create various graphical outputs for finite strain ellipsoids.

Processing 1000s of ellipsoids parameters is quite fast (I'm reusing some sweet python code available here https://github.com/minillinim/ellipsoid/blob/master/ellipsoid.py), but the bottleneck in my workflow is related to the time it takes to plot the plethora of 3D objects in a 3d plot.

Below I've attached a small snippet of python code that computes and plot a bunch of random ellipsoids. While 'ellipNumber' is small it works like a charm. But, when it reaches 100 it takes much longer.. with 1000s I bet you won't have the patience to wait.

In 2D, I understand using a collection is the way to go to improve performance: How can I plot many thousands of circles quickly?

Assuming a collection is indeed the way to go, I looked around for an example and attempted to populate a Poly3DCollection with ellipsoid coordinates like they did here for polygons in 3D: Plotting 3D Polygons in python-matplotlib, but I had no luck with setting up the vertices based on the 2d x,y and z arrays.

Any suggestion/comment on how to improve the plotting performance of the ellipsoids would be very appreciated!

Cheers 

import numpy as np
from numpy import linalg
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib.colors as colors



fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')

# number of ellipsoids 
ellipNumber = 10

#set colour map so each ellipsoid as a unique colour
norm = colors.Normalize(vmin=0, vmax=ellipNumber)
cmap = cm.jet
m = cm.ScalarMappable(norm=norm, cmap=cmap)

#compute each and plot each ellipsoid iteratively
for indx in xrange(ellipNumber):
    # your ellispsoid and center in matrix form
    A = np.array([[np.random.random_sample(),0,0],
                  [0,np.random.random_sample(),0],
                  [0,0,np.random.random_sample()]])
    center = [indx*np.random.random_sample(),indx*np.random.random_sample(),indx*np.random.random_sample()]

    # find the rotation matrix and radii of the axes
    U, s, rotation = linalg.svd(A)
    radii = 1.0/np.sqrt(s) * 0.3 #reduce radii by factor 0.3 

    # calculate cartesian coordinates for the ellipsoid surface
    u = np.linspace(0.0, 2.0 * np.pi, 60)
    v = np.linspace(0.0, np.pi, 60)
    x = radii[0] * np.outer(np.cos(u), np.sin(v))
    y = radii[1] * np.outer(np.sin(u), np.sin(v))
    z = radii[2] * np.outer(np.ones_like(u), np.cos(v))

    for i in range(len(x)):
        for j in range(len(x)):
            [x[i,j],y[i,j],z[i,j]] = np.dot([x[i,j],y[i,j],z[i,j]], rotation) + center


    ax.plot_surface(x, y, z,  rstride=3, cstride=3,  color=m.to_rgba(indx), linewidth=0.1, alpha=1, shade=True)
plt.show()

3D plot with 10 random ellipsoids:

3D plot with 10 random ellipsoids

Community
  • 1
  • 1
Guillaume Duclaux
  • 71
  • 1
  • 4
  • Internally every `plot_surface` creates an `art3d.Poly3DCollection`. So you don't gain anything by doing manually, what is being done automatically. I also think that 100 objects is really at the edge of what one can grasp looking at a figure, so how can the need for 1000 objects be motivated? Could there be a better way of conveying the same information? – ImportanceOfBeingErnest Jan 31 '17 at 11:05
  • You have a very valid point, and I am also producing a range of maps and diagrams to illustrate strain patterns variation within my surface or volume of interest. In my field (geology) we normally measure finite strain in rocks using some markers, and represent the strain quantity with ellipsoids. So this is definitely one important way of illustrating the spatial variation of the deformation. – Guillaume Duclaux Jan 31 '17 at 11:38
  • I'm working with 3D models as big as 1024^3, and when studying a simple section of this volume I end-up with 1024^2 ellipsoids, many (let say 70%) are seldom strained and I can just threshold them out before plotting. I normally plot only one every 100ish ellipsoid, but I still end-up with about 3000 of them. The method works and I remain convinced it is valuable. Here, I was just hopping to improve the performance of the process. – Guillaume Duclaux Jan 31 '17 at 11:40
  • Having a closer look you may want to reduce the number of points being drawn. At the moment, you calculate 60 x 60 = 3600 points per ellipse. I would suggest to use only 9 points per dimension. This might require to set `stride` to 1 instead of 3 but you'd still gain a factor of 5 in speed. – ImportanceOfBeingErnest Jan 31 '17 at 14:28
  • 1
    Reducing the number of points per ellipsoid offers a massive improvement in the plotting time. With a stride of 1 the resulting object are smooth enough for my application. Thanks for the suggestions! – Guillaume Duclaux Feb 01 '17 at 08:49
  • A lighter solution would be to plot the eigenvectors scaled by their magnitude – Seph42 May 06 '20 at 20:39

0 Answers0