Two animation scripts - and why Movimentum is still necessary

Here is a small tool chain to create animations for mechanical explanations. And here is an explanation of the compiler that creates ImageMagick (or other) calls from animation scripts.

Is there still a need to complete the Movimentum project (which is getting more and more complex)?

Well, there is! Let me explain this with two very simple mechanisms,

• the "Scotch yoke mechanism" and
• the simple crank (mechanism)

which also serve as examples of animation scripts in my "tiny animation language." Here are diagrams of both:

Scotch-yoke mechanism

Crank mechanism

Getting these mechanisms to turn their wheels will also be a short tutorial about how to write animation scripts.

Scotch-yoke is easy

Let us start is the animation for Scotch-yoke. Here is the result of what we'll develop:

The "model" for this gear is easy. We overlay three pieces:

• the fixed blue block (file Block.PNG),
• the turning red wheel (file Wheel.PNG), and
• the orange slider (file Slider.PNG).

The blue block is always at the same place:

  Block.PNG -geometry +0+0 -composite

The wheel is simply rotated by some angle !A (whose actual values will be set later). The center of the rotation must be found by inspecting the image with some graphics viewer that can show coordinates:

  ( Wheel.PNG -distort SRT "619,265,!A" ) -geometry +0+0 -composite

Finally, we have to move the slider by some x distance. This distance must vary like a sine wave. My little language can create sine waves only for vector variables, therefore we define such a variable %S, but use only its x component:

  Slider.PNG -geometry +%S:x+0 -composite

Overlaying these parts yields the following call to convert—the result is written to f_*.PNG:

$ ... convert -size 824x528 xc:white \
  Block.PNG -geometry +0+0 -composite \
  ( Wheel.PNG -distort SRT "619,265,!A" ) -geometry +0+0 -composite \
  Slider.PNG -geometry +%S:x+0 -composite \
  f_!T.PNG

When I created this animation the first time, I used 825x529 as size. The libx264 codec, however, complained that width and height were not even, so I changed them to 824x528—and it worked!

What about the variables' values?

Increasing !A is easy. We choose a rotation of 5 degrees per frame, starting at 0:

  !A 0 5 360

Moving %S is a little awkward. Here is one(!) way of doing it: %S starts at x coordinate 0 (because this is the initial shift of the slider), and has the same y as the rotation center of the wheel. Moreoever, we want to rotate it about some point which we call %C. Therefore, we have

  %S 0 0 %C 5

Now, we need to place %C. As the %S anchor of the slider is at its left end (the red circle), %C must be the rotation center of that end. Therefore, its coordinates are not (619,0)—the center of the wheel—, but (R,0), where R is the pixel radius of the wheel. Measuring this in some graphics program, I found R=118, and so the coordinates for %C are:

  %C 118 0 0 0

How many frames should we produce? 360 divided by 5 is 72, so 72 frames make one revolution. Two rotations would be fine, so we create 150 frames. Thus, our animation script is:

  !T 10000 1 .

  !A 0 5 360
  %C 118 0 0 0
  %S 0 0 %C 5
  $ 150 convert -size 824x528 xc:white \
    Block.PNG -geometry +0+0 -composite \
    ( Wheel.PNG -distort SRT "619,265,!A" ) -geometry +0+0 -composite \
    Slider.PNG -geometry +%S:x+0 -composite \
    e\f_!T.png

... and this is exactly what produces the animation shown above!

Actually, I like to start animations without any movement; and then accelerate the mechanism. Here is a script that does this:

!T 20000 1 .

# We start with 12 frames with increments of 0 degrees, i.e., nothing moves:
!A 0 0 360
%C 118 0 0 0
%S 0 0 %C 0
$ 12 convert -size 824x528 xc:white \
  Block.PNG -geometry +0+0 -composite \
  ( Wheel.PNG -distort SRT "619,265,!A" ) -geometry +0+0 -composite \
  Slider.PNG -geometry +%S:x+0 -composite \
  e\f_!T.png

# Next, we set the increment to 1 degrees for 4 frames ...
!A . 1 .
%S . . %C 1
$ 4

# ... then to 2 ...
!A . 2 .
%S . . %C 2
$ 4

# ... then to 3 ...
!A . 3 .
%S . . %C 3
$ 4

# ... then to 4:
!A . 4 .
%S . . %C 4
$ 4

# Then, we run with increments of 5 degrees for 150 frames.
!A . 5 .
%S . . %C 5
$ 150

# For deceleration, we reduce the increments to 3 for 4 frames ...
!A . 3 .
%S . . %C 3
$ 4

# ... then to 1 ...
!A . 1 .
%S . . %C 1
$ 4

# Finally, we end with 12 frames without movement:
!A . 0 .
%S . . %C 0
$ 12

And here is the result:

So far so good. But what about the crank?

The crank is difficult

We place the block and the wheel as before:

  Block.PNG -geometry +0+0 -composite
  ( Wheel.PNG -distort SRT "619,265,!A" ) -geometry +0+0 -composite

The driving rod must turn at its right end; and moreover, it must change its angle. In my script language, we need something like

  ( Rod.PNG -distort SRT "619,265,!R" ) -geometry +%S:x+%S:y -composite

where !R is the rod's angle. We can reuse %S from the Scotch-yoke mechanism ... I won't explain why, just try it out. Leaving !R at zero for the moment, we get a crank with "flying rod":

  !T 30000 1 .

  !R 0 0 360
  !A 0 5 360
  %C 118 0 0 0
  %S 0 0 %C 5
  $ 50 convert -size 824x528 xc:white \
    Block.PNG -geometry +0+0 -composite \
    ( Wheel.PNG -distort SRT "619,265,!A" ) -geometry +0+0 -composite \
    ( Rod.PNG -distort SRT "501,265,!R" ) -geometry +%S:x+%S:y -composite \
    e\f_!T.png

But what can we do about the rod's angle? Well—nothing. The "piston motion equations" for the crank mechanism are not simple linear or rotational formulas. Here is the equation for x in the "angle domain":

I could add functionality for this equation (and maybe I will, if I create the third animation with a crank). But of course there will be another linkage which still cannot be computed; and if I add that, there will be one more, and so on. Movimentum, on the other hand, as a general equation solver should solve almost all equations for mechanical linkages, without any need to patch the animation compiler.

What can we do right now? We can interpolate. Let us divide the circle into eight segments, each of which has 72/8 = 9 frames. At the end of the first segment, the right end of the rod will be 0.707r above the center line. The angle in the last frame must be asin(0.707 · r / d), where r is the radius of the crank and d the length of the rod. In our case, we find an angle of asin(0.707 · 118 / 447) = 10.75°. If we distribute this evenly over the 9 frames, we find an angle of 1.19° per frame. So we try the following:

!T 40000 1 .

!R 0 -1.19 360
!A 0 5 360
%C 118 0 0 0
%S 0 0 %C 5
$ 9

# stop at end
!R . 0 .
!A . 0 .
%S . . %C 0
$ 10

That looks almost acceptable. Ok, the left end goes a little bit up and down—but we'll leave it at that.

At the end of the next octant, the angle is asin(118 / 447) = 15.3°. We have to distribute 15.3° – 10.75° over 9 frames, which gives us 0.51° per frame for that octant:

!T 50000 1 .

!R 0 -1.19 360
!A 0 5 360
%C 118 0 0 0
%S 0 0 %C 5
$ 9

!R . -0.51 .
$ 9

# stop at end
!R . 0 .
!A . 0 .
%S . . %C 0
$ 10

For the next octants, we can copy the results which gives us the following animation:

!T 60000 1 .

!R 0 -1.19 360
!A 0 5 360
%C 118 0 0 0
%S 0 0 %C 5
$ 9

!R . -0.51 .
$ 9

!R . 0.51 .
$ 9

!R . 1.19 .
$ 18

!R . 0.51 .
$ 9

!R . -0.51 .
$ 9

!R . -1.19 .
$ 9

# stop at end
!R . 0 .
!A . 0 .
%S . . %C 0
$ 10

Mhm—not as good as needed. Of course, subdividing the circle in more segments is possible—but becomes more and more cumbersome. So, maybe a new version of an animation language would be helpful if cranks appear at more and more places.

I hope that this short introduction to my tool chain for simple mechaincal animations has been interesting. Of course, there is much more capable software around (ImageMagick itself can do some animation, then there's Flash programming, modern 3D CAD tools and what-have-you), but for me, it was fun to develop this tooling with a 100-line compiler, some existing command-line tools, and a simple CAD software at the front end!

A tiny language for mechanical animations - the AWK script

The AWK compiler for my little animation language is quite simple. Right now, it has 104 lines, and I do not intend to expand it. You can download the 3kB file here: www.haraldmmueller.de/eb/animate.awk.

1. Overall structure

The overall structure of the script is simple: There is a single block for each of the 5 possible commands.

/^#/      { ...
          }
/^!/      { ...
          }
/^%/      { ...
          }
/^~/      { ...
          }
/^\$.../  { ...
          }

The state of the script is contained in the following eleven variables:

1. rem: the prefix for shell comments, to be set in the gawk call.
2. valuS[v]: value of scalar variable v.
3. stepS[v]: increment of scalar variable v.
4. moduS[v]: modulus of scalar variable v.
5. frmtS[v]: print format of scalar variable v (necessary because modulus controls leading zeros).
6. valuX[v]: x value of vector variable v.
7. valuY[v]: x value of vector variable v.
8. stepXOrM[v]: either a number → increment for x value of vector variable v; or a string starting with % → anchor for second format of vector increment.
9. stepYOrA[v]: increment for y value of vector variable v (if stepXOrM[v] is a number); or angle increment (if stepXOrM[v] starts with %).
10. currcmd: The previous command. 
11. out: The flag remembering the ~ state.

The following sections explain the five code pieces.

2. Explicit Comments

This is easy: We print rem, the line number, and the text; and then continue with the next line ("next"):

/^#/  { print rem " (" NR ":" $0 ")"
        next
      }

3.Scalar variable assignment

This straightforward code consists of four if statements:

• The first one checks that there are exactly three parameters present (a rudimentary attempt at error handling—more of this is needed ...). The print statements uses gawk's extension > "/dev/stderr" which works also under Windows.
• The other three store the respective values in valuS, stepS, and moduS, unless a dot was specified. frmtS is set when the value or the modulus is set.0 + ... is the AWK idiom for converting a string to a number:

 
/^!/  { if (NF != 4) {
          print NR ": ERROR - exactly 3 parameters expected in " $0 > "/dev/stderr"
          next
        }
        if ($2 != ".") {
          valuS[$1] = 0 + $2
          frmtS[$1] = "%0d"
        }
        if ($3 != ".") {
          stepS[$1] = 0 + $3
        }
        if ($4 != ".") {
          moduS[$1] = 0 + $4
          frmtS[$1] = "%0" int(log(moduS[$1]) / log(10) + 1) "d"
        }
        print rem " v=" valuS[$1] ...more debugging output...
        next
      }

4. Vector variable assignment

Also vector variable assignments just store the values:

/^%/  { if (NF != 5) {
          print NR ": ERROR - exactly 4 parameters expected in " $0 > "/dev/stderr"
          next
        }
        if ($2 != ".") {
          valuX[$1] = 0 + $2
        }
        if ($3 != ".") {
          valuY[$1] = 0 + $3
        }
        if ($4 != ".") {
          stepXOrM[$1] = $4
        }
        if ($5 != ".") {
          stepYOrA[$1] = 0 + $5
        }
        print rem " v=" valuS[$1] ":" valuY[$1] " ...more debugging output...
        next
      }

5. Handling ~

This is simple:

/^~/    { out = $2
          next
        }

6. Command handling

Obviously, this is the interesting part. A $ line must contain a subsequent number of steps. Therefore, the script checks not only for a $ at the beginning of the line, but also for a subsequent number. The steps are remembered, then we collect the command if it is present. At last, the command is emitted multiple times. Here is the overall structure:

/^\$[ \t]+[0-9]+[ \t]+/    {
        steps = 0 + $2               
        ...remove leading $ and $2
        if ($0 != "") {
          ...collect lines with trailing \ into currcmd
        }
        for (i = 0; i < steps; i++) {
          cmd = currcmd
          ...interpolate scalar variables in cmd
          ...interpolate vector variables in cmd
          print cmd
          ...increment scalar variables
          ...increment vector variables
        }
      }      

Removing the leading $ and count is easy in AWK:

        sub(/^\$[ \t]+[0-9]+[ \t]+/, "")

If there is a command after the count, we collect it in variable currcmd. Handling trailing \ for command concatenation requires five more lines:

        if ($0 != "") {
          currcmd = $0
          while (/[\\]$/) {
            sub(/[\\]$/, " ")     
            getline
            currcmd = currcmd $0
          }
        }

Inside the emitting loop, we handle variables, print the output line, and finally increment the variables.

Replacing scalars has a few special points:

• First, all values are rounded to the nearest integer. The reason is that the images work with pixels.
• If modulus is set, we take the variable value modulo that modulus
• And finally, we format the value using the current format.

          for (v in valuS) {
            vv = int(valuS[v]+0.5)
            if (moduS[v] > 0) vv %= moduS[v]
            vv = sprintf(frmtS[v], vv);
            gsub(v, vv, cmd)
          }         

Replacing vectors is actually easier than scalars, as there are no formats and modulus for vectors. Also here, we round to the nearest integer:
         
          for (v in valuX) {
            vx = int(valuX[v]+0.5);
            vy = int(valuY[v]+0.5);
            gsub(v ":x", vx, cmd)
            gsub(v ":y", vy, cmd)
          }

After (possibly) printing the command, we have to increment the values. For scalars, this is a "one-liner":

          for (v in valuS) {
            valuS[v] += stepS[v]
          }

For vectors, we have two variants. One is easy—we just add the increments. The 0 + is needed here because stepXOrM is not coerced to a number when reading because it might be a %name:

          for (v in valuX) {
            if (stepXOrM[v] ~ /^%/) {
              ...angle increment - see below...
            } else {
              valuX[v] += 0 + stepXOrM[v]
              valuY[v] += stepYOrA[v]
            }
          }

Finally, we need some simple trigonometry for the angle increment. Here is the algorithm:

• First, we need the coordinates of the rotation center (mx and my).
• Then, we compute the delta vector between the current location and the rotation center (dx and dy).
• We convert this to angular form (r and phi) and immediately add the angle increment, after converting it from degrees to radians.
• Finally, we compute the new location by transforming and adding the new delta vector to the rotation center:

              mx = valuX[stepXOrM[v]]
              my = valuY[stepXOrM[v]]
              dx = valuX[v] - mx
              dy = valuY[v] - my
              r = sqrt(dx*dx + dy*dy)
              phi = atan2(dy, dx) + (stepYOrA[v] / 180 * 3.1416)
              valuX[v] = mx + r * cos(phi)
              valuY[v] = my + r * sin(phi)

And this is it!

The next posting will contain two examples of simple animations and, additionally, some arguments on why cranks are difficult.

A tiny language for mechanical animations using ImageMagick

My Movimentum project is currently fast asleep—but in a discussion about railroad barriers, I felt that I needed to explain a certain mechanical linkage with animations. Christmas season is a good time to experiment with new ideas (at least new to me), so I created three "mechanical explanations" for an old Austrian gear that drove (and drives) railroad barriers. You can see the results in that posting (in German). If you are interested in what it explains, just ask here—then I'll give you a quick overview about what this is all about. However, this posting here is really about the tool chain I used to create these animations.

Here is the second animation I made:

Here is an overview of my process. Below, I will share my experiences for each step:

1. I create the parts of the animation with a CAD tool (I use an older version of Caddy++).
2. The result is exported ("printed") into a single PNG file with PDFCreator.
3. The PNG file is split into multiple PNGs, each of which gets background transparency (I use an old version of Paint Shop Pro here).
4. Next, the images need "location calibration": In step 3., they do not come out exactly overlapping. The small corrections to move them to the right position are done here, using ImageMagick.
5. Next, I write the animation script. For this, I defined a very simple language that is interpreted by a 100-line AWK program (I use Gnu's Gawk for this). The description of this language is the major part of this posting.
6. The result of this is a batch file that calls ImageMagick's convert program to assemble the images into frames.
7. At last, ffmpeg is used to create an animation file from the frames.

For the record, I work with Windows. However, I think most of this also applies e.g. to Linux.

1. Creating the parts

Most probably, the parts have to be aligned in the animation. For this, it is very useful to add cross hairs as anchor points at critical locations, e.g., at the center of rotating items.

Moreover, to place the anchors at a predictable point, I put each part into a box of equal size, where the anchors are at the same position relative to the left upper box corner. In addition, I create an empty box that contains only the cross hairs of an anchor.

Here is a part of such a drawing. The lower left box is the empty one.

2. Creating PNG files

Instead of snipping the parts images from the CAD program or a PDF rendering, PDFCreator can be installed as a virtual printer that directly creates PNG files (despite its name). This avoids color and other distortions. Currently, I create a single PNG file with all the parts—this is faster than creating a single PNG file per part by printing from the CAD program, because I have to adjust each parts file anyway in the next step.

Already here, I select the correct scaling of the PNG—at least with Caddy++, this is possible in the print dialog. This avoids any later scalings which will destroy the cross hairs of the anchor points.

PDFCreator creates the PNG in portrait orientation even if you "print" in landscape—I do not know whether this is a fault of PDFCreator or Caddy++ or whatever else. In Windows, the created PNG is immediately opened in the Photo Viewer, so it is only a single click to rotate it to the intended orientation.

3. Splitting the PNG image and adding global transparency

Using Paint Shop Pro, I select the content of each box in the result and paste each part into a new PNG file. Unfortunately, box selection uses screen coordinates in Paint Shop Pro. Therefore, in the new files, the pasted parts may not lie at perfectly the same relative spot—I consistently got offsets of ±1 for both coordinates. In step 4 I explain how I correct this—but maybe there is a better way ... Because of this correction, I save the parts in files called ...Approx.PNG; the corrected files will have names without Approx.

While I have the new part file open, I make the background transparent, because this is what is needed for stacking the parts in the animation. Also the empty box gets its own "part file" with transparent background.

4. "Location calibration"

As explained above, the anchor points of the parts files will not be exactly at the same relative position to the upper left corner. This is a nuisance, as ImageMagick uses the upper left corner as a reference point for its computations, whereas we want to use our anchors. I repair this problem manually as follows:

• With ImageMagick, I overlay the the empty box over each part.
• In the resulting image, I check the distance of the part anchors vs. the empty box anchors.
• Again using ImageMagick, I correct the location of the part.

Here are the concrete steps for this on the command line, using ImageMagick's convert:

• Overlay and view with convert Part1Approx.PNG Empty.PNG -geometry +0+0 -composite temp.png
• Viewing temp.png reveals that the anchors can e.g. be superimposed precisely with +1+-2 instead of +0+0. This is a trial-and-error process. Nicely enough, Windows's Photo Viewer will refresh the image immediately, so this is a very quick process.
• Knowing the size of the empty box (e.g. use ImageMagick's identify—here I use 565x612 as example), invert the correction from +1+-2 to +-1+2, because we now want to shift the part's image, not the empty box. Then create the perfect part file via convert -size 565x612 xc:none Part1Approx.PNG -geometry +-1+2 -composite Part1.PNG.

5. The animation script

Now we come to the central piece: The animation script. The basic idea is that we create a separate ImageMagick convert command for each frame. Each such command will place the parts at specific places. But of course, instead of writing single commands for thousands of frames, I write patterns that contain variables that are expanded many times. The script is then run through an Gnu AWK program that emits the convert commands that compute the frames. Here is the call I use:

  gawk -v rem=REM -f animate.awk script.txt > script.bat

I will explain the "script compiler" animate.awk in a separate posting.

Update: During documenting the small compiler, I found I could make the language a little bit terser: The = step definition has been merged into the $ and + commands.

5.1. Command patterns and variables - overview

Here is such a command pattern that overlays a Part1 at a fixed place and a Part2 after a rotation and a shift (you can look up convert's options in ImageMagick's documentation):

... convert -size 900x800 xc:white \
    Part1.PNG -geometry +50+600 -composite \
    ( Part2.PNG -distort SRT "333,329,-!A2" ) \
      -geometry +%2:x+%2:y -composite \
  out\f_!T.png

The variables in this pattern are the strings starting with ! or %:

• !... is a scalar variable, Very often this is used for angles (like !A2 above), but sometimes also for numbering of files (like the !T at the end of the command).
• %... is a vector variable. Its components are accessed with %...:x and %...:y. In the command above, a vector variable %2 is used to shift the rotated Part2.

5.2. Script segments

A complete script consists of script segments where each segment describes a homogeneous movement. Each segment consists of two parts:

• A (possibly empty) sequence of variable assignments.
• A command pattern (convert patterns usually become very long, therefore they can be broken over multiple lines with trailing backslashes, as shown above).

Here is an example of such a segment:

  !T 10000 1 .
  $ 4 convert -size 900x800 xc:white \
      Text1.PNG -geometry +50+600 -composite out\f_!T.png

• !T 10000 1 . is an assignment to the variable !T: Its initial value is 10000, and it is incremented by 1 for each emitted convert. The last parameter (which is a single dot here) will be explained later.
• $ 4 convert ... means that the convert command is to be emitted 4 times, with variables substituted in each command.

This script segment will emit four lines where the variable !T is replaced with 10000, 10001, 10002 and 10003 in turn:

  convert -size .... out\f_10000.png
  convert -size .... out\f_10001.png
  convert -size .... out\f_10002.png
  convert -size .... out\f_10003.png  

5.3. The details of scalar variables

The assignment to a scalar variable has the following form:

  !name  initial  increment  modulus

name can consist of any characters you like that are not special in regexps (reason: In the AWK script, regexp replacement is used to interpolate the variables).

The three parameters can either be numbers (decimals are ok) or single dots (.). The dot means that the value remains as it was before. This is a central concept that helps to write short scripts. E.g., here is how you move something linearly, then let it stop:

  # 15 frames that reduce !X from 200 to some value.
  !X 200 -10 . 
  $ 15 convert ...

  # 10 frames without movement
  !X . 0 .
  $ 10 convert ...

In the second assignment, we left !X at whatever value it had after the first segment completed by specifying initial as . (dot). However, we set the increment to 0, so that this variable does not change its value in the 10 frames of this segment.

The modulus is sometimes useful for angle variables: The replacement value for a pattern is the variable's value modulo the modulus. I needed this because I precomputed rotated images of a part into 360 files called something like Part1_000.PNG, Part1_001.PNG,..., Part1_359.PNG. Each file contained the part rotated by so-and-so many degrees. For rotating the part, I could now write something like

  !A 350 2 360
  $ 20 convert ... Part_!A.PNG -geometry +0+0 -composite ...

Starting at 350 (degrees), the variable !A will be incremented by 2 in each of 20 frames. After 5 frames, !A will be 360 and continue with 362, 364 etc. up to 390. Without 360, the calls using Part_360.PNG etc. will fail, as these files do not exist. However, with modules 360, the emitted commands will contain Part_000.PNG, Part_002.PNG etc. The modulus, if defined, also determines the number of digits in the output number: Because modulus is 360 here, we get three-digit numbers.

5.4. The details of vector variables

A vector variable consists of two components, x and y. There are two possible assignments to vector variables. The first form does what you expect:

  %name  initialX  initialY  incrementX  incrementY

The variable is set and/or incremented as specified. Instead of a number, each argument can again be a single dot, which keeps the current value.

The other form is necessary for rotational movements:

  %name  initialX  initialY  %other  degrees

This means that the change to %name is a rotation around %other by the angle specified as the last argument. %other cannot be a dot here (the % is the marker that this form is intended), but the other three values can be numbers or dots. If  %other itself changes, the results are somewhat unpredictable, as the order of variable computations is not defined. If someone has the idea to try this—well, it's on your own risk ...

5.5. The details of commands

There are two possible forms for commands. One is

  $ count command-text

This creates count instances of command-text. The variables are replaced in each instance, and the result is written to the output. As explained above, long commands can be wrapped with \ at the lines' ends (the \ must actually be the last character on the line—no trailing white space!).

text can be anything your operating systems understands. I have, up to now, only called ImageMagick's convert here. However, one could e.g. write intermediary files and use them later in the frame sequence, e.g. to improve the performance of frame generation (thousands of calls to convert are not that fast ...).

It makes much sense to write the output of commands—especially the created frame images—to a specific output directory, to keep things cleanly separated.

The other command form omits the command-text:

  $ count

This means that the previous command is to be repeated in this segment. This is a very important concept: It allows one to specify the "model" in a single convert command that places all the parts at various places; and then have this assembly move simply by maniulating variables. One could even argue that a good animation script sets up a model once at the beginning; and then only modifies variables.

In addition to the $ lines, it is possible to switch off the output with

  ~ -

and switch it on at a specific point with

  ~ +

This is very helpful for debugging: It creates a script that will start "later in the animation", so that the script will only compute certain frames that are to be fine-tuned. This is preferable to (and much quicker than) editing the script each time it has been written before running it.

5.6. Comments

There are two possible forms of comments. One is

  # text

The text of such a comment is copied into the output, with a leading comment marker that must be passed to the compiler. This is helpful for debugging, because you can line up the script with the output.

The other form is everything else which is not understood by animation.awk, i.e., every line that does not start with !, %, ~, $, or #. Such lines are simply ignored. I use this feature to embed the script into a text that explains the animations, which will later become a web page. However, one could call this an unfortunate idea (even though it is along the lines of Donald Knuth's Literate Programming).

6. Running the batch file

This is as easy as it gets. You just run it.

7. Creating the animation

At last, ffmpeg is used to create an animation file from the frames. Here is the command line I use:

ffmpeg -f image2 -r 12 -i f_1%04d.png 
             -vcodec libx264 -pix_fmt yuv420p -y out.mp4

Of course, you can use other formats than mp4 for your animation (especially if you do not like proprietory formats). For mp4, the -pix_fmt yuv420p is important. The default pixel format, whose name I forgot, is not understood by the majority of programs or web sites out there.

If you want to create multiple animations from a single script (which happens if you embed your script code into running text, as I do), it makes sense to start the frame files with new numbers for each script, using !T 10000 . ., !T 20000 . ., !T 30000 . . etc. Then, you can run ffmpeg with file patterns f_1%04d, f_2%04d, f_3%04d etc.

In my next posting, I explain the small compiler (which you can download here).
In a subsequent posting, I will give a few animation examples and some background on why cranks are difficult.

Movimentum - Polynomial folding creates our second clip

We happily replace the old ConstantFoldingVisitor with the new PolynomialFoldingVisitor everywhere and expect that everything still works—or rather, works even better. After all, the new visitor does all of what the old one did (constants are polynomials!), and more. But alas, a handful of tests fails. For example, one test tells us that

Expected: <{EqualsZeroConstraint}0 = y+3*-y+y+6>
But was:  <{EqualsZeroConstraint}0 = ((-y+6))>

However, this was to be expected: Where we formerly got an expression tree, we now get a cleanly folded polynomial. Therefore, we change this test, together with the few others where the behavior has changed.

But the crucial test is our rotating bar. Here is, once more, its movimentum script (I have replaced the former x variable with len, because this value should be equal to the length of the bar):

    .config(8, 600, 400);
    WH : .bar C = [0,0] P = [-60,80];
    // length of bar: 100 - Pythagoras!

    @0
        WH.C     = [200,200];
        WH.P     = WH.C + [len,0].r(360*.t + 22.5);
    @2;

As before, the solver fails already for the first frame. But this time, we have much more manageable constraints:

{EqualsZeroConstraint}0 = ((-0.923879532511287*len))+-sqrt((-0.146446609406726*len^2+10000))
{EqualsZeroConstraint}0 = ((-0.38268343236509*len))+-sqrt((-0.853553390593274*len^2+10000))
{EqualsZeroConstraint}0 = -sqrt((-len^2+10000))

Which rules should we define to solve these constraints? As always, there's more than one way to skin the cat—not that I would actually consider that!—, but I opt for very "atomic" rules. Here is my suggestion:

• 0 = –E can be rewritten to 0 = E;
• 0 = sqrt(E) can also be rewritten to 0 = E;
• Finally, 0 = P[V] should solve the equation for V. Right now, we will only consider linear and quadratic equations here, because they suffice to solve our current rotating bar problem.

Here is the code for the first rule:

    var e = new TypeMatchTemplate<IAbstractExpr>();
    new RuleAction<ScalarConstraintMatcher>("0=-E",
        new EqualsZeroConstraintTemplate(-e).GetMatchDelegate(),
        matcher => matcher != null,
        (currNode, matcher, matchedConstraint) =>
            new SolverNode(currNode.Constraints
                                    .Except(matchedConstraint)
                                    .Union(new[] {
                                      new EqualsZeroConstraint(matcher & e),
                                    }),
                            currNode.ClosedVariables,
                            currNode));

(Compared to the rewrites in the PolynomialFoldingVisitor, this code seems long-winded. I'll leave it like this right now, but maybe I'll extract a few methods when the code gets on my nerves).

The second rule requires almost the same code, only the template and the actual match different:

    ...
    var t = new UnaryExpressionTemplate(new PositiveSquareroot(), e);
    new RuleAction<ScalarConstraintMatcher>("0=sqrt(E)",
        new EqualsZeroConstraintTemplate(t).GetMatchDelegate(),
        ...

Finally, we must find zeros for polynomials—a textbook job. Here is the rule—the crucial missing part is the method FindZeros, which must return the collection of all value where the poynomial intersects the x axis:

    var p = new TypeMatchTemplate<IPolynomial>();
    new RuleAction<ScalarConstraintMatcher>("0=P[V]",
        new EqualsZeroConstraintTemplate(p).GetMatchDelegate(),
        matcher => matcher != null,
        (currNode, matcher, matchedConstraint) =>
            FindZeros(matcher & p).Select(zero =>
                currNode.CloseVariable(
                    (matcher & p).Var,
                    Polynomial.CreateConstant(zero),
                    matchedConstraint)
                 )
            );

Here is the full code for FindZeros:

    private static IEnumerable<double> FindZeros(IPolynomial polynomial) {
        switch (polynomial.Degree) {
            case 0:
                throw new ArgumentException("Cannot find zeros for constant");
            case 1: {
                // 0 = ax + b --> x = -b/a                
                var a = polynomial.Coefficient(1);
                var b = polynomial.Coefficient(0);
                return new[] { -b / a };
            }
            case 2: {
                // 0 = ax² + bx + c --> with D = b² - 4ac,
                // we get x1,2 = (-b+-sqrt D)/2a.
                // Attention: These formulas are not good
                // from a numerical standpoint! - see Wikipedia.
                double a = polynomial.Coefficient(2);
                double b = polynomial.Coefficient(1);
                double c = polynomial.Coefficient(0);
                double d = b * b - 4 * a * c;
                double sqrtD = d.Near(0) ? 0 : Math.Sqrt(d));
                return new[] { (-b + sqrtD) / (2 * a),
                               (-b - sqrtD) / (2 * a) };
            }
            default:
                throw new NotImplementedException( "Cannot find zeros for " + 
                                "polynomial of degree " + polynomial.Degree);
        }
    }

And this time, all the constraints are solved! We have another clip:




But what's that? The bar does not rotate! Rather, it flaps around helplessly only on the left side. What is wrong? To find this out, I wrote a new test that asserts each endpoint coordinate against a hand-computed value. Here is its core:

foreach (var f in frames) {
    IDictionary<string, IDictionary<string, ConstVector>>
        anchorLocations = f.SolveConstraints(10000, ref previousState);

    double alpha = f.FrameNo / 8.0 * 2 * Math.PI - Math.PI / 8;

    double expectedX = 200 + 100 * Math.Cos(alpha);
    double expectedY = 200 + 100 * Math.Sin(alpha);
    ConstVector whpLocation = anchorLocations["WH"]["P"];
    Assert.AreEqual(expectedX, whpLocation.X, 1e-2);
    Assert.AreEqual(expectedY, whpLocation.Y, 1e-2);
}

Here is its output before it fails:

++++ solution node ++++
SolverNode 1021 < 0<=C-1020 < 0<=C-1019 < 0<=C-1018 < 0<=C-1017 < 0=C-1015 < 0=P[V]...
  len:=-100
  WH.C.X:=200
  WH.C.Y:=200
  WH.C.Z:=0
  WH.P.X:=107,612046748871
  WH.P.Y:=161,731656763491
  WH.P.Z:=0
  ZERO:=0

  Expected: 292.38795325112869d +/- 0.01d
  But was:  107.61204674887132d

Looking at the variables in the last node, we see that len is minus 100—i.e., our solver found a solution where the bar has a negative length! That can only mean that the bar is oriented into the opposite direction of what we expect! We did not consider this when writing the script ... but the solution is easy: We add a constraint that restricts the length to positive values:

    ...;
    len       > 0; 
    ...;

And now the test that checks the coordinates happily runs green!

When we add this new constraint also to the clip-producing test, we finally get our second working clip out of Movimentum:



We are happy!

However, there is a usability issue lurking behind the problem with the negative length: How is the script writer ever going to find such missing constraints? It seems one needs a script debugging concept on the user level if it should be possible to create even moderately complex animation scripts. Mhm—one more item for the todo list.

Movimentum - Polynomial folding 3

5. Plus

Using "ASCII art," here is a table that shows how addition must work. Again, the goal is to merge polynomial subexpressions of a single variable at the left, and collect everything else in a single "rest" expression on the right. Be aware that E and F might also be polynomials, if their variables are different from V!

    //        | Q[V]+F              Q[V]            F
    // -------+------------------------------------------------
    // P[V]+E | {P[V]+Q[V]}+(E+F)   {P[V]+Q[V]}+E   P[V]+(E+F)
    //        |
    // P[V]   | {P[V]+Q[V]}+F       {P[V]+Q[V]}     P[V]+E
    //        |
    // E      | Q[V]+(E+F)          Q[V]+E          E+F

Why don't I consider entries with unary minus, or even square roots?

• Square roots, of course, cannot be merged together by addition, so handling them separately will not create smaller, more compact expressions.
• For unary minus, the argument is a little trickier: First, there cannot be a –P[V] subexpression anywhere—it would have been rewritten to {–P[V]}, i.e., a simple polynomial with inverted coefficients, on a lower level. The E's and F's could be expressions with a leading unary minus, but moving their minus upwards will not reduce the complexity of the expression in general. In specific cases, such a rewrite might get rid of expressions—e.g., if E happens to be –F or the like. But we do not optimize such "rare occurrences" here.

Polynomial folding is a heuristics. Like constant folding, it tries to standardize expressions quickly. However, if it oversees certain simplifications, this is not a big problem, as long as we do not rely on simplification guarantees.

Coding the table for the plus operator is easy—for example, here are the three entries for the first line of the table:

    new StandardExpressionRewrite("(P+E)+(Q+F)", _plusRewrites,
        (p + e) + (q + f),
        HaveSameVar(p, q),
        m => ((m & p) + (m & q)) + ((m & e) + (m & f)));
    //                PO         AE         AE

    new StandardExpressionRewrite("(P+E)+Q", _plusRewrites,
        (p + e) + q,
        HaveSameVar(p, q),
        m => ((m & p) + (m & q)) + (m & e));

    new StandardExpressionRewrite("(P+E)+F", _plusRewrites,
        (p + e) + f,
        m => (m & p) + ((m & e) + (m & f)));
    }

Below the first rewrite rule, I have marked the + operators with their "type":

• PO is a polynomial operator, i.e., it will merge its two operands into a single polynomial.
• AE, on the other hand, is an AbstractExpr operator, i.e., it will create an expression tree (a BinaryExpression).

The actual visiting code is similar to the square and unary minus operators:

    public IAbstractExpr Visit(Plus op, IAbstractExpr lhs,
                               IAbstractExpr rhs, Ignore p) {
        return Rewrite(new BinaryExpression(lhs, op, rhs),
                       "Cannot handle " + op, _plusRewrites);
    }

6. Times

My first attempt for multiplication was the same as for addition: Create a table that contains all combinations of "interesting expression structures," and explicitly define the result. From the square operator, we learn that we should handle square roots separately if they occur on both sides of the operator. In that first attempt, I identified 10 subexpression forms that I needed to distinguish: P+√R, P+–√R, P+√E, P+–√E, P+E, P+–E, P, C, –E, and general E. This gives us a table of 100 entries (or even more, because we have to distinguish (P+√R)·(Q+√R) from (P+√R)·(Q+√S), where all of P, Q, R, S are polynomials of the same variable). This is too much work—specifying, coding, and last, but not least, testing.

I thought a little bit and came up with the idea of stage-wise processing—it is what we do manually to "simplify" such expressions. For example, if we have (x + √4x) · (x - √x), we will

• first multiply the two terms which gives us a sum of four terms;
• then simplify each term;
• and finally the whole sum.

The result from the first step can be a sum of up to four terms, i.e., an expression tree of depth at most 3. Various factors in the sums might need simplification—e.g. in our example, where √4x · √x will have to be merged to √4x², which means that we have to run the tree through folding again at least three levels down from the current root. Instead of introducing a special mechanism for this, we run the visitor over the complete result—this will uselessly visit expressions below depth 3, but we will have to optimize away useless visits anyway somewhen in the future, so we will ignore this overhead right now.

The three steps above are only a rough plan. Actually, we need to consider a little bit more what "simplifications" we have to apply to the result of the first step. Here is my suggestion for a three-step algorithm.

In the first version of this posting, I bragged "similarly to all previous algorithms, I do not give any proof why this accomplishes what we want; it is "more or less obvious," ...". Well, in addition to being obvious, it was also wrong. The tables for Step III claimed, among others, that P√Q could be transformed to √(P²Q). However, this cannot be true: The latter is never negative (as we define sqrt to be the positive square root), but the former can of course be negative—for example, if the polynomial P is x – 1, and x happens to be zero. Here is a better table; but still no proof—I risk keeping this a mere software project, not a scientific undertaking, although I'd like to do all these proofs. Maybe I find time for them somewhen later ...

    // Step I: Distribute.
    //
    //          Q+F                 F
    //
    //  P+E     {PQ}+Q*E+P*F+E*F    P*F+E*F
    //
    //  E       Q*E+E*F             =
    //

    // Step II: Lift unary minus to top.
    //
    //          -F       F
    //
    //  -E      E*F      -(E*F)
    //
    //  E       -(E*F)   =
    //

    // Step III: Pull up square roots as far as possible.
    //
    //          √Q              √F              Q       D               F
    //
    //  √P      P=Q:  P         √(P*F)          =       D<0: -√{PD²}    =
    //          else: √{PQ}                             else: √{PD²}
    //
    //  √E      √(Q*E)          E=F:  E         =       D<0: -√({D²}*E) =
    //                          else: √(E*F)            else: √({D²}*E)
    //
    //  P       =               =               {P*Q}   {P*D}           =
    //                                                X
    //  C       C<0: -√({C²Q}   C<0: -√({C²}*F) {C*Q}   {C*D}           =  
    //          else: √({C²Q}   else: √({C²}*F)
    //
    //  E       =             =                 =       =               =
    //

An equals sign, in the tables above, simple means that no rewrite is done for the indicated product. Together, these three steps require 3 + 3 + 15 = 21 rules—a lot less than the 100 or more needed in my first attempt!
The four rewritings around the X can be handled by the same rule, which simply matches the product of two polynomials—also constants are polynomials!

Here is the code for Step I: The rules are almost like earlier ones—but in addition, there is a call to RecursivelyFold on the result:

    new StandardExpressionRewrite("(P+E)*(Q+F)", _timesRewrites,
        (p + e) * (q + f),
        HaveSameVar(p, q),
        m => RecursivelyFold((m & p) * (m & q)
                            + (m & q) * (m & e)
                            + (m & p) * (m & f)
                            + (m & e) * (m & f)));
    new StandardExpressionRewrite("(P+E)*F", _timesRewrites,
        (p + e) * f,
        m => RecursivelyFold((m & p) * (m & f)
                            + (m & e) * (m & f)));
    new StandardExpressionRewrite("E*(Q+F)", _timesRewrites,
        e * (q + f),
        m => RecursivelyFold((m & q) * (m & e)
                            + (m & e) * (m & f)));

RecursivelyFold is a single line method:

    private static readonly PolynomialFoldingVisitor
        _recursiveFolder = new PolynomialFoldingVisitor();

    private static IAbstractExpr RecursivelyFold(AbstractExpr expr) {
        return expr.Accept(_recursiveFolder);
    }

Here is an excerpt from Step III:

    new StandardExpressionRewrite("√P*√P", _timesRewrites,
        sqrtP * sqrtP,
        m => RecursivelyFold(m & p));
    new StandardExpressionRewrite("√P*√Q", _timesRewrites,
        sqrtP * sqrtQ,
        HaveSameVar(p, q),
        m => PosSqrtAndRecursivelyFold((m & p) * (m & q)));
    new StandardExpressionRewrite("√P*√F", _timesRewrites,
        sqrtP * sqrtF,
        m => PosSqrtAndRecursivelyFold((m & p) * (m & f)));
    // √P*Q --> =
    new StandardExpressionRewrite("√P*D,D<0", _timesRewrites,
        sqrtP * d, m => (m & d).Value < 0,
        m => -PosSqrtAndRecursivelyFold(
                    (m & d).Value * (m & d).Value * (m & p)
        ));

PosSqrtAndRecursivelyFold wraps the parameter in a PositiveSquareroot after recursively folding it:

    private static IAbstractExpr PosSqrtAndRecursivelyFold(AbstractExpr expr) {
        return new UnaryExpression(RecursivelyFold(expr), new PositiveSquareroot());
    }

If you look closely, you will notice that I put all the rules for multiplication into a single list, called _timesRewrites. This means that later stages, i.e. the RecursiveFold calls, will run over the rules for earlier stages again. This is actually useless and could be optimized away, but I guess that that would require even more complex algorithms than the ones I invent here. So I leave it as it is.

Of course, we also need an implementation of the Visit method: But it is similar to the one for the + operator.

Together with more than one hundred unit tests, this code rewrites binary expressions with a Times operator to our intended P+E format!

6. Divide

Division, finally, is very simple:

• If the division is by a constant, we replace it with a multiplication of the reciprocal value.
• Otherwise, we leave the expression alone (we do not try polynomial division or the like!):

    new StandardExpressionRewrite("E/C", _divideRewrites,
        e / c,
        m => RecursivelyFold(
              (m & e) * Polynomial.CreateConstant(1 / (m & c).Value)
            )
    );
    new StandardExpressionRewrite("E", _divideRewrites,
        e,
        m => m & e);

The corresponding Visit method is, of course, like the ones for plus and times.

And this concludes the visitor for rewriting polynomials.