Message Boards Message Boards

Coin Flip: Simulate Collisions of 3D Rigid Bodies

Tim Laska, A Priori, Ltd
Posted 4 years ago

Rendered Coin Flip

I have used the free program Blender v2.79b to simulate the handling of 100s of complex shapes through a geometrically complex industrial machine with many moving parts including vibrating elements. So, it should be able to handle a "coin flip". I believe Blender still uses the Bullet Physics Engine as its solver. I should warn you that collision simulation can get difficult and that there are a lot of tricks of the trade you must learn to be accurate and fast in a general case.

Blender has a python interface and it can be run as a background task (Bullet also has a python interface, but I am not familiar with its operation). Since Mathematica can create text files with StringTemplate and can execute system commands, we should be able to create a python script to drive a Blender simulation.

Python Script Generation

Blender has a fairly well-documented API and there many resources one can find online to generate python script. 

import bpy
from math import pi

for o in bpy.data.objects:
    if o.type == 'MESH' or o.type == 'EMPTY':
        o.select = True
    else:
        o.select = False

# Delete all objects in the scene
bpy.ops.object.delete()

# Add the floor
bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "BOX"
boxObj.name = "Ground"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "CYLINDER"
bpy.context.object.rigid_body.friction = 0.25
bpy.context.object.rigid_body.restitution = 0.75
boxObj.name = "Coin"
# Set reference to the coin
coin = bpy.data.objects["Coin"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Advance two frames and add translational and rotational motion
sce.frame_set(3)
# Translate up a little
coin.location.z = 3.45
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = 1
coin.rotation_euler.y = 0.1
coin.rotation_euler.z = 0.1
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
            'point_cache': bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print("          X                  Y                  Z                  RX                  RY                  RZ")
print(tr.x, tr.y, tr.z, eu.x, eu.y, eu.z)

if eu.x > pi / 2.:
    print("Coin is heads")
else:
    print("Coin is tails")

The script above does is the following:

  1. Selects and deletes all extraneous objects in the scene.
  2. Adds a passive cube (note that cubes are better than planes for collisions) for the ground.
  3. Adds an active cylindrical coin.
  • Sets friction to 0.25
  • Sets restitution to 0.75
  1. Adds keyframe animate to set initial velocity and rotation.
  2. Prints out and tests if coin has rotated more than $\frac{\pi}{2}$ radians to determine heads/tails.

Driving from Mathematica

We can create a parametric model in Mathematica by replacing hard coded parameters with template variables using `` delimiters as in the createCoinFlip function.

createCoinFlip[z_, rx_, ry_, rz_, friction_, restitution_] := 
 StringTemplate["import bpy
from math import pi

for o in bpy.data.objects:
    if o.type == 'MESH' or o.type == 'EMPTY':
        o.select = True
    else:
        o.select = False

# Delete all objects in the scene
bpy.ops.object.delete()

# Add the floor
bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"BOX\"
boxObj.name = \"Ground\"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, \
location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"CYLINDER\"
bpy.context.object.rigid_body.friction = `friction`
bpy.context.object.rigid_body.restitution = `restitution`
boxObj.name = \"Coin\"
# Set reference to the coin
coin = bpy.data.objects[\"Coin\"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Advance two frames and add translational and rotational motion
sce.frame_set(3)
# Translate up a little
coin.location.z = `z`
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = `rx`
coin.rotation_euler.y = `ry`
coin.rotation_euler.z = `rz`
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
            'point_cache': \
bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print(\"\
          X                  Y                  Z                  RX                  R\
Y                  RZ\")
print(tr.x, tr.y, tr.z, eu.x, eu.y, eu.z)

if eu.x > pi / 2.:
    print(\"Coin flip result is heads\")
else:
    print(\"Coin flip result is tails\")
"][<|"z" -> z, "rx" -> rx, "ry" -> ry, "rz" -> rz, 
"friction" -> friction, "restitution" -> restitution|>]

Blender will send a lot of information to standard out. We can parse this output with Find to extract a line of interest. Putting it all together, the following will create a python script, run Blender in the background, and parse the output.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlip[3.95, 1, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --background --python coinflip.py >>" <> outputfile];
stext = OpenRead[outputfile];
Find[stext, "Coin"]
Close[stext];
DeleteFile[outputfile]
(* Coin is tails *)

You can visualize the results of the simulation by removing the "--background" and repeating the step above.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlip[3.45, 1, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --python coinflip.py >>" <> outputfile];
stext = OpenRead[outputfile];
Find[stext, "Coin"]
Close[stext];
DeleteFile[outputfile]

If you left click anywhere on the screen and hit the play button, then you should see the following:

Blender Coin Flip

Rendering in Blender

You can take advantage of Blender's photorealistic rendering capability to create a nice animation if desired.

Rendered Coin Flip

Additional Post Processing in Mathematica

Blender is geared more towards the artist whereas Mathematica is geared more towards the physicist. We can find synergy when we combine the strengths of both tools.

What follows is a simple example of how to conduct additional post processing on a Blender simulation in Mathematica.

First, let's modify the python generation script to give the positions and orientations of the coin at each frame (we will insert a string "PosRot" to identify the appropriate lines).

createCoinFlipTransform[z_, rx_, ry_, rz_, friction_, restitution_] :=
  StringTemplate["import bpy
from math import pi

for o in bpy.data.objects:
    if o.type == 'MESH' or o.type == 'EMPTY':
        o.select = True
    else:
        o.select = False

# Delete all objects in the scene
bpy.ops.object.delete()

# Add the floor
bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"BOX\"
boxObj.name = \"Ground\"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, \
location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
cylObj = bpy.context.active_object
cylObj.rigid_body.collision_shape = \"CYLINDER\"
bpy.context.object.rigid_body.friction = `friction`
bpy.context.object.rigid_body.restitution = `restitution`
cylObj.name = \"Coin\"
# Set reference to the coin
coin = bpy.data.objects[\"Coin\"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Advance two frames and add translational and rotational motion
sce.frame_set(3)
# Translate up a little
coin.location.z = `z`
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = `rx`
coin.rotation_euler.y = `ry`
coin.rotation_euler.z = `rz`
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
            'point_cache': \
bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()

for i in range(250):
    sce.frame_set(i)
    tr = coin.matrix_world.translation
    eu = coin.matrix_world.to_euler()
    print(\"PosRot\",tr.x, tr.y, tr.z, eu.x , eu.y , eu.z )
"][<|"z" -> z, "rx" -> rx, "ry" -> ry, "rz" -> rz, 
"friction" -> friction, "restitution" -> restitution|>]

We can extract the positions and orientations of the simulation with the following code.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlipTransform[4, -Pi 0.75, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --background --python coinflip.py >>" <> outputfile];
stext = OpenRead[outputfile];
data = ToExpression@StringSplit[#] & /@ FindList[stext, "PosRot"];
{tx, ty, tz, rx, ry, rz} = Transpose@data[[All, {2, 3, 4, 5, 6, 7}]];
Close[stext];
DeleteFile[outputfile]

We can define a cuboid and cylinder that have the same dimensions as the Blender simulation and we can create a transformation function with the following code.

box = {Cuboid[{-5, -5, -0.5}, {5, 5, 0.5}]};
cyl = {Cylinder[{{0, 0, -0.05}, {0, 0, 0.05}}, 1], 
   AbsolutePointSize[10], 
   Opacity[1], {Black, Point[{0, 0, 0}]}, {Red, 
    Point[{1, 0, 0}]}, {Green, Point[{0, 1, 0}]}, {Blue, 
    Point[{0, 0, 1}]}};
m = IdentityMatrix[4];
m[[1 ;; 3, 1 ;; 3]] = EulerMatrix[{a, b, c}, {1, 2, 3}];
m[[1 ;; 3, -1]] = {x, y, z};
transform[a_, b_, c_, x_, y_, z_] = TransformationFunction[m];

Now, we can combine plots of position and orientation (or other quantities like angular momentum) into a Manipulate[] function.

Manipulate[
 Column[{Row[{ListPlot[{tx[[1 ;; i]], ty[[1 ;; i]], tz[[1 ;; i]]}, 
      Filling -> Axis, ImageSize -> {200, 200}, PlotRange -> All, 
      PlotLegends -> {"tx", "ty", "tz"}],
     ListPlot[{rx[[1 ;; i]], ry[[1 ;; i]], rz[[1 ;; i]]}, 
      Filling -> Axis, ImageSize -> {200, 200}, PlotRange -> All, 
      PlotLegends -> {"rx", "ry", "rz"}]}], 
   Graphics3D[{{Opacity[0.75], Red, box}, 
     GeometricTransformation[{Opacity[.85], Yellow, cyl}, 
      transform[rx[[i]], ry[[i]], rz[[i]], tx[[i]], ty[[i]], 
       tz[[i]]]]}, SphericalRegion -> True,  Boxed -> False, 
    ImageSize -> {400, 400}]}], {i, 1, 250, 1}]

Physics Graphs and 3D Animation

The original version of this post can be found HERE.
POSTED BY: Tim Laska
1 Reply
Sort By:

enter image description here - Congratulations! This post is now featured in our Staff Pick column as distinguished by a badge on your profile of a Featured Contributor! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract