Animated 3D plotting with Blender

From Penguin Development
Jump to: navigation, search

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 / Animated_3D_plot.ogv