Difference between revisions of "Animated 3D plotting with Blender"

From Penguin Development
Jump to navigationJump to search
(Created page with "[https://www.blender.org/ Blender] is an extremely versatile 3D creation and animation suite. Since it is fully scriptable in [https://www.python.org/ Python], Blender may be...")
(No difference)

Revision as of 07:08, 9 November 2016

Blender is an extremely versatile 3D creation and animation suite. Since it is fully scriptable in Python, Blender may be used to generate animated 3-dimensional plots of data or mathematical functions. Below is an example of one way to generate such a plot.


Important information

The script below has been tested with Blender 2.78a for Linux. It should generally be compatible with other versions and operating systems. The script works with the pristine Blender startup file; heavily modified startup files may break certain assumptions.

To run the script, simply call blender --python blenderplot-ani.py

To create the plot in a different base file, instead use blender basefile.blend --python blenderplot-ani.py. The order of arguments matters here.

Note that the animation data will not be saved with your .blend file! This means you must run the Python script on a pristine "background" .blend file each time you wish to examine the geometry or create a render. If the .blend file is saved, only the first frame will be captured.

The script

Code for blenderplot-ani.py follows.

  1 #!/bin/true
  2 # vim: se fo=tcroq tw=78 :
  3 # Simple animated 3D plot example using Blender ( https://www.blender.org/ ).
  4 # Given the function f(k, x, t)=exp(ikx-iωt), plots Re(f) against x and k,
  5 # with colour given by Im(f) and t being the time.
  6 #
  7 # For pedagogical purposes, this just computes f at each frame (twice: once
  8 # for the vertex positions and once for their colours). This is horribly
  9 # inefficient; it would be much better to generate a 3D array for f(x, k, t)
 10 # once and slice this array for each frame -- however, this is left as an
 11 # exercise to the reader.
 12 #
 13 # To run, call
 14 #   blender --python blenderplot-ani.py
 15 
 16 import os.path
 17 
 18 import numpy as np
 19 
 20 import bpy
 21 
 22 ### Begin user settings
 23 omega = 1
 24 font = '/usr/share/fonts/cm-unicode/cmunti.ttf' # Must be a unicode font!
 25 ### End user settings
 26 
 27 ### Begin generic Blender rendering code
 28 # Global object counter.
 29 obj_ind = 10000
 30 
 31 plot_id = None
 32 
 33 line_material = bpy.data.materials.new('line')
 34 line_material.diffuse_color = (0, 0, 0)
 35 line_material.diffuse_shader = 'LAMBERT'
 36 line_material.specular_color = (0, 0, 0)
 37 line_material.specular_shader = 'COOKTORR'
 38 line_material.use_shadows = False
 39 line_material.use_cast_shadows = False
 40 line_material.use_raytrace = True
 41 line_material.ambient = 0
 42 
 43 text_material = bpy.data.materials.new('text')
 44 text_material.diffuse_color = (.15, .05, .035)
 45 text_material.diffuse_shader = 'OREN_NAYAR'
 46 text_material.diffuse_intensity = .9
 47 text_material.roughness = 2
 48 text_material.specular_color = (.6, .2, .1)
 49 text_material.specular_shader = 'PHONG'
 50 text_material.specular_hardness = 80
 51 text_material.specular_intensity = .85
 52 text_material.use_shadows = True
 53 text_material.use_cast_shadows = False
 54 text_material.use_raytrace = True
 55 text_material.raytrace_mirror.use = True
 56 text_material.mirror_color = (.7, .3, .15)
 57 text_material.raytrace_mirror.reflect_factor = .3
 58 text_material.emit = 0
 59 text_material.ambient = 0
 60 
 61 plot_material = bpy.data.materials.new('plot')
 62 plot_material.specular_color = (.5, .5, .5)
 63 plot_material.specular_shader = 'COOKTORR'
 64 plot_material.specular_intensity = .2
 65 plot_material.use_shadows = True
 66 plot_material.use_transparent_shadows = True
 67 plot_material.use_raytrace = True
 68 plot_material.use_transparency = True
 69 plot_material.transparency_method = 'RAYTRACE'
 70 plot_material.alpha = .95
 71 plot_material.specular_alpha = 1
 72 plot_material.raytrace_transparency.depth = 5
 73 plot_material.use_vertex_color_paint = True
 74 
 75 text_font = bpy.data.fonts.load(os.path.abspath(os.path.expanduser(font)))
 76 
 77 def heatmap(heat):
 78     """Heat map: given a "heat" between 0 and 1, return a tuple of RGB
 79     values."""
 80     r = np.max((2*heat-1., 0))
 81     b = np.max((1.-2*heat, 0))
 82     return (r, 1.-r-b, b)
 83 
 84 def zheat(z, zmin, zmax, **kwargs):
 85     """Colour a vertex based on its z-height compared to the minimum and
 86     maximum z-values that occur."""
 87     return heatmap((z-zmin)/(zmax-zmin if zmin != zmax else 1))
 88 
 89 def plot_function(x, y, func, auto_axes = True, xmarks=None,
 90         ymarks = None, zmarks = None, labels = None, thickness = 0.025,
 91         text_rot = None, colourfunc=zheat, zmin = None, zmax = None):
 92     """Plot the function (lambda) func of x and y. The resulting surface is
 93     smooth-shaded. Vertices may be coloured according to colourfunc: this is a
 94     function that accepts the following parameters and returns an RGB-tuple:
 95         x, y, z, xmin, xmax, ymin, ymax, zmin, zmax"""
 96     global obj_ind, plot_id
 97 
 98     ids = {
 99         'axes': [],
100         'axis_labels': [],
101         'xmarks': [],
102         'ymarks': [],
103         'zmarks': [],
104         'xlabels': [],
105         'ylabels': [],
106         'zlabels': []
107     }
108 
109     if text_rot is None:
110         text_rot = np.array((0, 0, 0))
111 
112     if plot_id is None:
113         obj_id = 'plot_{}'.format(obj_ind)
114         obj_ind += 1
115         # Generate all vertices in the plot at z = 0
116         verts = [(i, j, 0) for i in x for j in y]
117         faces = []
118         count = 0
119         # Build faces from the vertices
120         for i in range(len(y)*(len(x)-1)):
121             if count < len(y)-1:
122                 faces.append((i, i+1, i+len(y)+1, i+len(y)))
123                 count += 1
124             else:
125                 count = 0
126 
127         # Create a mesh and an object at the origin
128         mesh = bpy.data.meshes.new(obj_id)
129         obj = bpy.data.objects.new(obj_id, mesh)
130         obj.location = (0, 0, 0)
131         bpy.context.scene.objects.link(obj)
132         mesh.from_pydata(verts, [], faces)
133         mesh.update(calc_edges=True)
134 
135         # Create a new vertex colour map
136         colours = obj.data.vertex_colors.new()
137 
138         # Set material
139         obj.data.materials.append(plot_material)
140 
141         # Smooth-shade polygons
142         for pol in obj.data.polygons:
143             pol.use_smooth = True
144     else:
145         obj = bpy.data.objects[plot_id]
146         colours = obj.data.vertex_colors.active
147 
148     verts = obj.data.vertices
149 
150     # Move vertices to their correct position
151     for v in verts:
152         v.co.z = func(v.co.x, v.co.y)
153     obj.data.update(calc_edges=True)
154 
155     sv = sorted([(v.co.x, v.co.y, v.co.z) for v in verts], key=lambda q: q[2])
156 
157     # Colour vertices
158     for pol in obj.data.polygons:
159         for idx in pol.loop_indices:
160             co = obj.data.vertices[obj.data.loops[idx].vertex_index].co
161             colours.data[idx].color = colourfunc(x=co.x, y=co.y, z=co.z,
162                 xmin=np.min(x), xmax=np.max(x),
163                 ymin=np.min(y), ymax=np.max(y),
164                 zmin=sv[0][2], zmax=sv[-1][2])
165 
166     if auto_axes and plot_id is None:
167         # Axes
168         ids['axes'].append(add_line(np.array((min(x), min(y), 0)),
169             np.array((max(x), min(y), 0)), thickness, False))
170         ids['axes'].append(add_line(np.array((max(x), min(y), 0)),
171             np.array((max(x), max(y), 0)), thickness, False))
172         ids['axes'].append(add_line(
173             np.array((min(x), min(y), sv[0][2] if zmin is None else zmin)),
174             np.array((min(x), min(y), sv[-1][2] if zmax is None else zmax)),
175             thickness, False))
176         # Axis marks
177         if xmarks is not None:
178             for pos, label in xmarks:
179                 p = np.array((pos, min(y), 0))
180                 ids['xmarks'].append(add_line(p,
181                     p-np.array((0, 1.5*thickness, 0)), thickness, False))
182                 if label is not None and len(label) > 0:
183                     ids['xlabels'].append(add_text(
184                         p-np.array((0, 7*thickness, 0)), label,
185                         thickness, text_rot))
186         if ymarks is not None:
187             for pos, label in ymarks:
188                 p = np.array((max(x), pos, 0))
189                 ids['ymarks'].append(add_line(p,
190                     p+np.array((1.5*thickness, 0, 0)), thickness, False))
191                 if label is not None and len(label) > 0:
192                     ids[ 'ylabels'].append(add_text(
193                         p+np.array((7*thickness, 0, 0)), label,
194                         thickness, text_rot))
195         if zmarks is not None:
196             for pos, label in zmarks:
197                 p = np.array((min(x), min(y), pos))
198                 ids['zmarks'].append(add_line(p,
199                     p-np.array((1.5*thickness/np.sqrt(2),
200                         1.5*thickness/np.sqrt(2), 0)), thickness, False))
201                 if label is not None and len(label) > 0:
202                     ids['zlabels'].append(add_text(
203                         p-np.array((7*thickness, 7*thickness, 0))/np.sqrt(2),
204                         label, thickness, text_rot))
205 
206         # Axis labels
207         if labels is not None:
208             ids['axis_labels'].append(add_text(
209                 np.array((max(x)+8*thickness, min(y), 0)),
210                 labels[0], 2*thickness, text_rot))
211             ids['axis_labels'].append(add_text(
212                 np.array((max(x), max(y)+8*thickness, 0)),
213                 labels[1], 2*thickness, text_rot))
214             ids['axis_labels'].append(add_text(
215                 np.array((min(x), min(y),
216                     (sv[-1][2] if zmax is None else zmax) + 8*thickness)),
217                 labels[2], 2*thickness, text_rot))
218 
219     if plot_id is None:
220         plot_id = obj_id
221 
222     ids['plot'] = plot_id
223 
224     return ids
225 
226 def add_text(r, text, size=0.025, rotation=None):
227     """Add text at the position r. Size is a relative parameter; use
228     trial-and-error here."""
229     global obj_ind
230     obj_id = 'text_{}'.format(obj_ind)
231     obj_ind += 1
232     rot = [np.pi/2, 0, 0] if rotation is None else rotation.tolist()
233     bpy.ops.object.text_add(location=r.tolist(), rotation=rot)
234     obj = bpy.context.active_object
235     obj.name = obj_id
236     obj.data.name = obj_id
237     obj.data.body = text
238 
239     # Set the font
240     obj.data.font = text_font
241 
242     obj.data.offset_x = -2*size
243     obj.data.offset_y = -2*size
244     obj.data.shear = 0.0
245     obj.data.size = 8*size
246     obj.data.space_character = 1
247     obj.data.space_word = 4*size
248     obj.data.extrude = size/3
249 
250     obj.data.materials.append(text_material)
251 
252     return obj_id
253 
254 def add_line(r1, r2, w=0.01, rel_w=True):
255     """Add a "line" (cylinder) between the points r1 and r2. The width is
256     either w (if rel_w is False) or w*|r2-r1| (if rel_w is True)."""
257     global obj_ind
258     obj_id = 'line_{}'.format(obj_ind)
259     obj_ind += 1
260     rc = (r2+r1)/2 # Centroid
261     rr = r2-rc # Position of r2 relative to centroid
262     r = np.sqrt(np.sum(rr**2))
263     theta = np.arccos(rr[2]/r)
264     phi = np.arctan2(rr[1], rr[0])
265     bpy.ops.mesh.primitive_cylinder_add(vertices=16,
266             radius=.5*w*(r if rel_w else 1), depth=2*r,
267             location=rc.tolist(), rotation=(0, theta, phi))
268     obj = bpy.context.active_object
269     obj.name = obj_id
270     for pol in obj.data.polygons:
271         pol.use_smooth = True
272     obj.data.materials.append(line_material)
273 
274     return obj_id
275 ### End generic Blender rendering code
276 
277 def frame_change(scene):
278     """Update the plot for a given frame."""
279     frame = min(max(scene.frame_current, 0), n_frames - 1)
280     plot_function(x, k, funcgen(t[frame]), colourfunc=colgen(t[frame]))
281     # Update the t-indicator
282     bpy.data.objects[ttext_id].data.body = 't = {: >4.3f}'.format(t[frame])
283 
284 if __name__ == '__main__':
285     # Set up x, y and t data
286     n_frames = 51
287     nx = 101
288     nk = 101
289     xscale = 1/2
290     kscale = 1/3
291     zscale = 2
292     x = np.linspace(0, 10, nx)*xscale
293     k = np.linspace(0, 4*np.pi, nk)*kscale
294     t = np.linspace(0, 10, n_frames)
295 
296     # Function to plot
297     func = lambda x, k, t: np.exp(1j*(k*x-omega*t))*zscale
298     # Generator for plottable function f(x, k) at time t
299     funcgen = lambda t: lambda x, k: np.real(func(x, k, t))
300     # Colour generator
301     colgen = lambda t: lambda x, y, **kwargs: \
302             heatmap((np.imag(func(x, y, t))+1)/2)
303 
304     # Absolute z range for axes
305     azmin = -1*zscale
306     azmax = 1*zscale
307 
308     # Axis marks
309     xmarks = [(i*xscale, str(i)) for i in range(1, 11)]
310     kmarks = [(j*np.pi/2*kscale, '{}π/2'.format(j if j != 1 else '') \
311         if j % 2 == 1 else '{}π'.format(j // 2 if j != 2 else '')) \
312         for j in range(1, 9)]
313     zmarks = [(z/10*zscale, str(z/10)) for z in range(-10, 12, 2)]
314 
315     # Hide the 吸牛 splash screen
316     bpy.context.user_preferences.view.show_splash = False
317 
318     # Remove existing meshes
319     for item in bpy.context.scene.objects:
320         if item.type == 'MESH':
321             bpy.context.scene.objects.unlink(item)
322     for item in bpy.data.objects:
323         if item.type == 'MESH':
324             bpy.data.objects.remove(item)
325     for item in bpy.data.meshes:
326         bpy.data.meshes.remove(item)
327 
328     # Set the camera position
329     bpy.data.objects['Camera'].location = (11, -6, 5.5)
330     bpy.data.objects['Camera'].rotation_euler = (1.1, 0, 0.8)
331 
332     # Let all texts face the camera
333     rot = np.array(bpy.data.objects['Camera'].rotation_euler)
334 
335     # Initial t=0 plot; this also sets up the axes.
336     print(plot_function(x, k, funcgen(0),
337         True, labels=('x', 'k', 'z'), text_rot=rot, xmarks=xmarks,
338         ymarks=kmarks, zmarks=zmarks, zmin=azmin, zmax=azmax,
339         colourfunc=colgen(0)))
340 
341     # Text label indicating current time
342     ttext_id = add_text(np.array((5.6, -1.2, 0)), 't = {: >4.2f}'.format(0),
343         size=0.05, rotation=rot)
344 
345     # Set min/max/current frame
346     bpy.data.scenes['Scene'].frame_start = 0
347     bpy.data.scenes['Scene'].frame_end = n_frames - 1
348     bpy.data.scenes['Scene'].frame_current = 0
349 
350     # Add frame change handler. This is what makes the animation happen!
351     bpy.app.handlers.frame_change_pre.append(frame_change)
352 
353     # Add some environment lighting
354     wld = bpy.data.worlds['World']
355     wld.light_settings.use_environment_light = True
356     wld.light_settings.environment_energy = .5
357 
358     # Add a white backdrop plane
359     plane_material = bpy.data.materials.new('backdrop')
360     plane_material.diffuse_color = (1, 1, 1)
361     plane_material.use_shadeless = True
362     bpy.ops.mesh.primitive_plane_add(location=(0, 0, -5))
363     bpy.context.active_object.scale = (50, 50, 0)
364     bpy.context.active_object.data.materials.append(plane_material)

Output video

File:Animated 3D plot.ogv