Spiro2SVG II: Roulettes, Lissajous Figures, and Farris Wheels to SVG

• Download demo - 44.9 KB
• Download source - 79.5 KB

Based on the Open Source project Spiro2SVG. This article is a follow-up to the article Spiro2SVG - Convert Spirographs to SVG using Bezier Curves. We will develop more general methods of converting roulettes into SVG, including Lissajous Figures and Farris Wheels.

Introduction

Previously, a method was described for converting standard spirograph shapes (hypotrochoids and epitrochoids) into SVG using a Bezier curve-fit method. That method was fairly specialized, not suitable for other types of roulette shapes. The current article will discuss the conversion of two new types of shapes, using more general methods. The first type comes from the program SpiroJ.jar. This program comes complete with a library of sample files, most of which are included here in the Samples folder. This program also uses two rotors, as before, but the amplitude and frequency in the x-direction are allowed to be different from the amplitude and frequency in the y-direction. This allows for the creation of a wide variety of shapes, most of which do not have any rotational symmetry, except perhaps a single mirror symmetry in either the x- or y-direction. A subset of this type of shape is the Lissajous Figure, which mimics the behavior of an oscilloscope. The second type of new shape is the Farris Wheel, described in an article by Frank Farris. This uses three rotors instead of two, and the rotors are once again "round", with the x-and y-components being the same for both amplitide and frequency; but with arbitrary phase shifts in each rotor. These shapes can be generated by the program GeometryDemonstration.exe which is part of the article "Spirograph Shapes: WPF Bezier Shapes from Math Formulae" by Ken Johnson.

The main thrust of this article is to develop more general methods of handling the above shapes, with an emphasis on saving the data in an SVG file rather than displaying it; and on doing so using cubic Bezier curves, so that the user does not need to make any decisions as to how many fit points are needed or what quality is desired. We will also describe methods to preserve whatever rotational symmetry may exist in the original figure, even if that symmetry is not known ahead of time.

Using the Program

The program usage is essentially the same as the previous paper: "Convert Spirographs to SVG". Clicking on the Spiro2SVG.jar file will bring up a File Open dialog. Browse to the Samples folder and open an .xml file, not a .spiro file as before. The .xml files come in two types. The first type is the format used by SpiroJ, for two rotors. Here is the complete file for the image meduza.svg which is shown above.

<?xml version="1.0" encoding="UTF-8"?>
<spiro>
   <editor type="advanced" />
   <shape linewidth="1" linecolor="#000000" fillcolor="#DDDDFF" />
   <generator steps="600" />
   <operator type="rotor">
      <radius x="50.0" y="50.0" />
      <frequency x="3.0" y="3.0" />
   </operator>
   <operator type="rotor">
      <radius x="30.0" y="30.0" />
      <frequency x="11.0" y="8.0" />
   </operator>
</spiro>

Note that the frequency is different in the x- and y-directions, and that the file contains the parameter "generator steps," which is the number of straight lines to use. We will not use this parameter if we choose the Bezier rendering method.

The second type of .xml file contains the format needed for the Farris Wheel. Unfortunately, I have not found any program that actually saves these data files for future use. Therefore, I have arbitrarily created a format which mimics the SpiroJ format. Here is the xml data file for the image Farris_1_7_-17.svg which is shown above:

<?xml version="1.0" encoding="UTF-8"?>
<spiro>
   <editor type="Farris" />
   <shape linewidth="1.7" linecolor="#00006A" fillcolor="#B8E0CC" />
   <generator steps="600" />
   <operator type="rotor">
      <radius r="120" />
      <frequency w="1" />
      <phase phi="0" />
   </operator>
   <operator type="rotor">
      <radius r="60" />
      <frequency w="7" />
      <phase phi="0" />
   </operator>
   <operator type="rotor">
      <radius r="40" />
      <frequency w="-17" />
      <phase phi="90" />
   </operator>
</spiro>

This data file can be used as a template for creating data files for other types of Farris Wheels, with different frequencies and amplitudes. The relevant input data can be copied from the GeometryDemonstration program once you have displayed an interesting shape using that program. Note that this data file has three rotors instead of two, and no distinction between x and y. We will also be assuming throughout that the frequencies are always integers, so that the objects will become closed shapes within a finite number of revolutions; and we will specify the phase "phi" in units of degrees.

After loading the input file, you will have the choice of rendering method, as before, and the choice of output svg file name. In this article, we will discuss only the Bezier rendering method, since the other two are straightforward and proceed the same as in the first article.

Program Calculation Sequence

The calculation sequence is as described in the first article, except that it sometimes splits and then merges again as different file formats are encountered. The flow of information is shown in this figure:

If a .spiro file was chosen, it will be processed as described in the first article. If an .xml file was chosen, the routine SpiroParse.parse_spiroJ_file() will decide what kind it is. If it is of type SpiroJ (two rotors), the following field names will be used in the SpiroConfig dialog:

public static final String[][] SpiroJNames = new String[][]
           {{"Radius_x1"}, {"Radius_y1"}, {"Frequency_x1"}, {"Frequency_y1"},
            {"Radius_x2"}, {"Radius_y2"}, {"Frequency_x2"}, {"Frequency_y2"},
            {"generator_steps"}, {"line_width"}, {"line_color"}, {"fill_color"},
            {"Edit Drawing Style"}};

If it is a Farris Wheel (three rotors), the field names will be:

public static final String[][] FarrisNames = new String[][]
           {{"Radius_1"}, {"Frequency_1"}, {"Phase_1"}, {"Radius_2"},
            {"Frequency_2"}, {"Phase_2"}, {"Radius_3"}, {"Frequency_3"}, {"Phase_3"},
            {"generator_steps"}, {"line_width"}, {"line_color"}, {"fill_color"},
            {"Edit Drawing Style"}};

For the SpiroJ shapes, the routine SpiroWrite.convertspiroJParms() will make a final determination of the type of the shape. It is possible for a SpiroJ shape to be identical to a standard spirograph object, if the frequencies and amplitudes in the x- and y-directions are the same. In this case, the object will be re-directed into the standard spirograph processing routine, getspiroShape, since it is specialized to deal with this shape. Otherwise it will be processed by the routine getspiroJShape which is more generic in nature. For the Farris Wheels there is no decision to be made at this point, since it is clear that they must be processed separately, in the routine getFarrisShape.

We now assume that the Bezier rendering method has been chosen. The process of fitting a Bezier curve to a known spirograph segment is the same as in the first article, but the choice of appropriate t values is quite different, and needs to be described.

Choosing t values for a SpiroJ Shape

The SpiroJ object obeys the equations:

x = r_x1 cos(w_x1t) + r_x2 cos(w_x2t)
y = r_y1 sin(w_y1t) + r_y2 sin(w_y2t)

In general this object will not have any rotational symmetry, and there is no simple way to sub-classify these shapes, other than the two special cases mentioned above, which are the standard spirograph and the Lissajous figure. Therefore, we will adopt a very simple-minded generic approach in this section. The primary requirement is that the arc angle between start- and end-point can never be allowed to be greater than 90 degrees. To achieve this, we fit all points which have vertical slope, as well as horizontal slope, and all inflection points. The vertical slope constraint leads to this equation for t:

The horizontal slope constraint leads to this equation for t:

r_y1 w_y1 cos(w_y1t) + r_y2 w_y2 cos(w_y2t) = 0

The inflection point is determined as follows: from Eq.(3) in the first article, we have an expression for curvature. Setting this to zero yields the constraint: ẋÿ - ẏẍ = 0. This leads to the following equation for t:

This equation is different than the previous ones because it contains a product of two trigonometric functions. However, we can simplify it, so that it does not contain any products, using the identities:

We now need a generic solver that can set the sum of an arbitrary number of sine and cosine waves to zero. The previous article used the routines SpiroCalc.calc_cos_t() and SpiroCalc.calc_sin_t() to solve similar equations, but they are too specialized for this problem. In order to solve these equations we need a general method of representing the equation data: we use the array

rotors = new double[][] {{amplitude, frequency, phase}, ...};

where the first subscript refers to the index number of the rotor. For example, Eq.(2) will yield eight rotors coming from eight unique frequencies. The second subscript refers to the three parameters: amplitude, frequency and phase; where phase will be specified in radians and will be used primarily to distinguish between sine and cosine waves. With this notation, Eq.(1) is represented by:

rotors = new double[][] {{rx1*wx1, wx1, -Math.PI/2},
                         {rx2*wx2, wx2, -Math.PI/2}};

while Eq.(2), after simplification using Eq.(3), is represented by:

rotors = new double[][] {{rx1*ry1*wx1*wy1*(wx1 + wy1)/2, wx1 - wy1, 0},
                         {rx1*ry1*wx1*wy1*(wx1 - wy1)/2, wx1 + wy1, 0},
                         {rx1*ry2*wx1*wy2*(wx1 + wy2)/2, wx1 - wy2, 0},
                         {rx1*ry2*wx1*wy2*(wx1 - wy2)/2, wx1 + wy2, 0},
                         {rx2*ry1*wx2*wy1*(wx2 + wy1)/2, wx2 - wy1, 0},
                         {rx2*ry1*wx2*wy1*(wx2 - wy1)/2, wx2 + wy1, 0},
                         {rx2*ry2*wx2*wy2*(wx2 + wy2)/2, wx2 - wy2, 0},
                         {rx2*ry2*wx2*wy2*(wx2 - wy2)/2, wx2 + wy2, 0}};

This data is sent to SpiroJCalc.solve_cos_t(double r[][], double t0) where the roots will be found using the Newton-Raphson method. The main problem encountered here is the choice of t0, the initial estimate of the root. There will be many roots to find, and there is always the possibility of converging to the wrong root and accidentally missing the one that was hoped for. We avoid this by using a certain amount of overkill: for every single rotor (cosine wave), and for every single possible zero of this individual rotor within the range of interest, we perform a separate call to the solver. This should bring us reasonably close to the actual root, if it exists, assuming that the other rotors are not varying too quickly. The code for calling the solver is:

for (i = 0; i < rotors.length; i++)
   if ((Math.abs(rotors[i][0]) > TOL) && (Math.abs(rotors[i][1]) > TOL))
      for (int j = 0; j < Math.round(Math.abs(2*rotors[i][1])); j++)
         N = main.insert_t_value(N, N, t, 
             solve_cos_t(rotors, Math.PI*(j + 0.5)/Math.abs(rotors[i][1])));

There is no theoretical guarantee that this method will find all the roots, but so far it has not yet failed. What it does do, however, is produce a large number of duplicate, redundant, roots. These are detected and removed in the routine SpiroJCalc.sort_t_values(double[] t, int n).

Given a set of t values to be used as fit points, we can produce the Bezier curve fit the same as before, using the SpiroJCalc.getBezier() routine. A typical SVG result is shown below using the SpiroJ file holub.xml. The diagram also shows the fit points in red. Note that there is a fundamental difference in strategy here, compared to the previous article: previously, when using standard spirograph shapes, we knew a priori what the symmetry was, meaning how many unique lobes there were, so we were able to fit just the first lobe and then replicate the result as often as needed. For the current shapes we have no knowledge of what the symmetry may or may not be, so we have to explicitly calculate all the fit points for the entire object. This is not a problem for the general SpiroJ shape which normally has no symmetry anyways, but it will potentially be a problem for the Farris Wheel discussed below, which actually does have rotational symmetry, except that the symmetry is not explicitly known or used by the fitting method. For that case, we will need to develop a method that will automatically respect the object's symmetry even if we do not know what the degree of the symmetry is. But first, we need to make a side-trip to investigate a specialized feature of the Lissajous shape.

Curvature at Lissajous Stationary Points

The method developed above for SpiroJ shapes works well for the general shape that has no symmetry. However, it fails in the special case of a Lissajous figure that reverses its' direction of motion. This yields a stationary point that must be treated specially. An example of this is given below. This shows two superimposed figures: the first one has w_x1 = 1, w_y1 = 3, which leads to a closed figure where the motion never stops. The second one has w_x1 = 2, w_y1 = 3, which leads to a figure where the motion stops and reverses itself, so that it looks like a single curve, but is actually two paths superimposed on each other, with stationary points at each end. These stationary points are different from the ones encountered in the normal spirograph. In the normal spirograph the stationary points were actually cusps where the curvature had a singularity as |c| became equal to |b|. The spirograph curvature at this point could be either ±∞ depending on the relative magnitude of c versus b. For the Lissajous stationary points, on the other hand, the curvature is well-behaved but must be evaluated analytically, with some care, since both numerator and denominator are approaching zero. Without loss of generality, we can represent the Lissajous figure using the simplified notation:

x = r_x cos(w_xt)
y = r_y sin(w_yt)

The stationary points will occur when ẋ = 0 and ẏ = 0 simultaneously. This requires that:

w_xt = Nπ
w_yt = (M + ½)π

where M and N are arbitrary integers. Eliminating t between these equations, we get a new constraint on the frequencies:

It will not always be possible to solve these equations for an arbitrary w_x and w_y. For example, Eq.(4) never has a solution if w_x is odd. However, if a solution exists, then we wish to perform a Taylor series expansion about that specific t value in order to capture the leading term in the expressions for ẋ, ẍ, etc., so that we can calculate the curvature. We perform the Taylor expansion in terms of the variable Δt = t - Nπ/w_x = t - (M + ½)π/w_y. This yields

ẋ ≈ -(-1)^N r_x w_x (w_x Δt - (w_x Δt)³/6)
ẏ ≈ -(-1)^M r_y w_y (w_y Δt - (w_y Δt)³/6)

Substituting this and the corresponding ẍ and ÿ into the expression for curvature, Eq.(3) of the first article, we get

κ_stationary = (-1)^N+M r_x r_y w_x² w_y² (w_x² - w_y²) / 3 / (r_x² w_x⁴ + r_y² w_y⁴)^3/2

where we have kept only the leading terms, which were of order Δt³ in both numerator and denominator, and which therefore cancel perfectly as Δt aproaches zero. We use this result in SpiroJCalc.getBezier() to calculate curvature at stationary points. This yields the result shown below as a blue curve, for w_x1 = 2, w_y1 = 3.

Finding Extrema of Farris Wheel Curvature

The Farris Wheel presents a new challenge because it normally will possess some rotational symmetry, but it is clear that our current fitting method will not display that symmetry. Our method currently fits all vertical and horizontal slopes. However, if the Farris wheel has six-fold symmetry, and if we fit one lobe using vertical and horizontal slopes, and if we then rotate those fitted points by sixty degrees, then it is clear that they will not correspond to horizontal and vertical slopes on the new lobe, so the new lobe will not be a perfect clone of the original one, rotated by 60 degrees. This is the sense in which we wish to respect the symmetry of the object, but the problem is that the exact nature of the symmetry may not be immediately obvious, and the shape is sufficiently complex that it is not clear how to make explicit use of the symmetry. Similarly, for the standard spirograph shape it was easy to sub-classify the shape and develop special techniques for each sub-category; for the Farris Wheel this is probably not feasible. So we need to use shape measurements that will automatically preserve the symmetry if it exists, without imposing it explicitly. One such measure is the property curvature, which is invariant with respect to orientation. Therefore any measure derived from curvature will also be invariant. This leads us to consider using points of maximum or minimum curvature as fit points.

We begin by searching for four types of fit points:

• points of maximum absolute value of curvature
• points of minimum absolute value of curvature
• inflection points (zero curvature)
• points whose tangent is 90 degrees perpendicular to a point of maximum curvature

As a side note, we find that this set of criteria is essentially identical to what was originally done for the hypotrochoid, because it turns out that, for a standard spirograph, a point of maximum radius also happens to be a point of maximum curvature as well, and similarly for minimum radius (curvature). Therefore the current set of criteria is a natural generalization of the method previously used for standard spirographs.

To calculate the inflection points, we use the same method as for the SpiroJ shape, and also use Eq.(3), to get the set of rotors:

double Sum23 = r1*r1*w1*w1*w1 + r2*r2*w2*w2*w2 + r3*r3*w3*w3*w3;
double A12   = r1*r2*w1*w2;
double A23   = r2*r3*w2*w3;
double A31   = r3*r1*w3*w1;
rotors = new double[][] {{Sum23, 0, 0},
                         {A12*(w1 + w2), w1 - w2, phi1 - phi2},
                         {A23*(w2 + w3), w2 - w3, phi2 - phi3},
                         {A31*(w3 + w1), w3 - w1, phi3 - phi1}};

Next, the points of maximum/minimum curvature are obtained by setting dκ/dt = 0, which yields the equation:

The individual terms in this expression are given by:

ẋÿ - ẏẍ = r₁² w₁³ + r₂² w₂³ + r₃² w₃³
           + r₁ r₂ w₁ w₂ (w₁ + w₂) cos(θ₁₂)
           + r₂ r₃ w₂ w₃ (w₂ + w₃) cos(θ₂₃)
           + r₃ r₁ w₃ w₁ (w₃ + w₁) cos(θ₃₁)

ẋ² + ẏ² = r₁² w₁² + r₂² w₂² + r₃² w₃²
            + 2 r₁ r₂ w₁ w₂ cos(θ₁₂)
            + 2 r₂ r₃ w₂ w₃ cos(θ₂₃)
            + 2 r₃ r₁ w₃ w₁ cos(θ₃₁)

where θ_ij = (w_i - w_j)t + phi_i - phi_j. Taking first derivatives of these terms and substituting into Eq.(5) we obtain an equation for t, which can be cast into the form of a standard set of twelve rotors using Eq.(3). The rotors are a bit too complex to be shown here, but they are defined explicitly in FarrisCalc.get_t_values() and solved in SpiroJCalc.solve_cos_t().

Finally, the points that are rotated 90 degrees relative to the point of maximum curvature are obtained by scanning through the previously obtained fit points and searching for an angle change greater than 90 degrees. If it occurs, then the tangent angle of the point of highest curvature is calculated and the search begins for a fit point that is perpendicular to it on either side. As an initial estimate we use a t value that is the average of the t value at the point of highest curvature and its nearest neighbor on either side. So far, this initial estimate has always led to convergence on the desired result. A typical Farris Wheel is shown below, using frequencies of (-2, 5, 19) which yield a 7-fold rotational symmetry. The fit points are shown in red. Note that it shows the desired symmetry, in that the fit points for one lobe can be rotated by (1/7)^th of a revolution to superimpose directly on the fit points for the next lobe.

Stationary Points in the Farris Wheel

The Farris Wheel can also show cusps which need to be treated specially because the curvature is not defined at these points. These cusps occur when the motion is temporarily stationary. They are analogous to the cusps that occur in the normal spirograph when |c| = |b| (i.e. when a hypotrochoid turns into a hypocycloid). However, for a Farris Wheel, these stationary points can occur over a range of different radii, not just a single radius, because we have one more degree of freedom in a Farris Wheel. The best way to describe these cusps is to calculate the behavior of the velocity of a point on the Farris Wheel. The velocity vector of a Farris Wheel point is itself a Farris Wheel as well, with new parameters. If the position of the original Farris Wheel point was given by a set of rotors of the type:

rotors = new double[][] {{r1, w1, phi1}, ...};

then the velocity vector of the same point is given by the vector sum of rotors of the type:

rotors = new double[][] {{r1*w1, w1, phi1 + Math.PI/2}, ...};

Using this new velocity representation, we can now see that a stationary point in the original Farris Wheel, Eq.(6), will occur whenever the three velocity vectors, Eq.(7), add up to zero. This is not always theoretically possible. For example, if we arbitrarily rename the three rotors in such a way that |r₁ w₁| > |r₂ w₂| > |r₃ w₃|, then it will be theoretically possible for the velocity to be zero only if

|r₁ w₁| < |r₂ w₂| + |r₃ w₃|.

This defines maximum and minimum values that a wheel radius can have, assuming we hold the other two radii fixed. If this inequality is satisfied, so that a stationary point is theoretically possible, then we can calculate the specific conditions under which it actually will occur. The specific conditions will show up as constraints on the phase shifts phi_i; for example:

r₁² w₁² = r₂² w₂² + r₃² w₃² + 2 r₂ r₃ w₂ w₃ cos(phi₂ - phi₃)

where analogous equations will also apply to (phi₁ - phi₂) and (phi₃ - phi₁). All three of these constraints can be viewed as expressions of the Cosine Law which applies to the vector sum of three vectors that add up to zero. Armed with these equations we can now build a cusp detector for the Farris Wheel to determine a priori when a cusp exists. This is done in the routine FarrisCalc.calc_t_cusp(). It is necessary to predict these cusps up front, because otherwise we would likely get an infinite or undefined curvature at this point. If a cusp is detected, then we will arbitrarily set the Bezier arm length to zero at this point, the same as was done for the standard spirograph.

Under normal circumstances, the presence of cusps of this type is quite rare, but it is interesting to see what they look like, by creating them deliberately. The figure below shows an animation based on a series of 32 images which each contain cusps. They were produced by taking the Farris Wheel from the top figure, with frequencies (1, 7, -17), and varying the smallest wheel radius to produce a sequence of different figures with cusps, from the smallest possible radius to the largest radius, with a number of examples in between. This was done just as a test, to confirm that the Bezier curve calculation remains well-behaved in these cases.

Future Work

The methods used for the Farris Wheel are more generally applicable than the methods used for either the SpiroJ shapes or the standard spirograph shapes. It seems quite likely that it would be possible to standardize the code to a very high degree by applying the curvature maximization/minimization code to all the shapes described here, to produce a single unified approach to all three cases. However, that would involve a large amount of rewriting of code that is already working well, so I think I will leave that as an exercise for an interested reader.