新手 HTML Java

Spiro2SVG II：轮子线、李萨如图和摩天轮到SVG

Alvin Penner

5.00/5 (3投票s)

2017年1月3日

CPOL

19分钟阅读

9701

341

一个文件转换器，用于使用贝塞尔曲线将通用的轮子线数据文件转换为SVG图形。

基于开源项目Spiro2SVG。本文是文章Spiro2SVG - 使用贝塞尔曲线将棘轮图转换为SVG的续篇。我们将开发更通用的方法，将轮子线转换为SVG，包括李萨如图和摩天轮。

引言

之前，我们介绍了一种使用贝塞尔曲线拟合方法将标准的棘轮图形状（内摆线和外摆线）转换为SVG的方法。该方法相当专业化，不适用于其他类型的轮子线形状。本文将使用更通用的方法讨论两种新类型形状的转换。第一种类型来自程序SpiroJ.jar。该程序附带了一个样本文件库，其中大部分文件包含在“Samples”文件夹中。该程序也使用了两个转子，与之前一样，但x方向的幅度和频率可以与y方向的幅度和频率不同。这允许创建各种形状，其中大多数形状没有旋转对称性，可能只有x或y方向的单一镜像对称性。这类形状的一个子集是李萨如图，它模仿了示波器的行为。第二种新形状是摩天轮，由Frank Farris在文章中描述。它使用了三个转子而不是两个，并且转子再次是“圆的”，x和y分量的幅度和频率相同；但每个转子都有任意的相位偏移。这些形状可以通过Ken Johnson的《棘轮图形状：来自数学公式的WPF贝塞尔形状》一文的程序GeometryDemonstration.exe生成。

本文的主要目的是开发处理上述形状的更通用方法，重点是将数据保存到SVG文件中而不是显示它；并使用三次贝塞尔曲线进行保存，这样用户就不需要决定需要多少拟合点或期望的质量。我们还将描述保留原始图形中可能存在的任何旋转对称性的方法，即使这种对称性是事先未知的。

使用程序

程序用法与上一篇文章基本相同：《转换棘轮图到SVG》。单击Spiro2SVG.jar文件将弹出文件打开对话框。浏览到Samples文件夹并打开一个.xml文件，而不是像以前那样打开.spiro文件。.xml文件有两种类型。第一种类型是SpiroJ使用的格式，用于两个转子。这是上面所示的meduza.svg图像的完整文件。

<?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>

请注意，x和y方向的频率不同，并且文件包含“生成器步数”参数，即要使用的直线数量。如果我们选择贝塞尔渲染方法，则不会使用此参数。

第二种类型的.xml文件包含摩天轮所需的格式。不幸的是，我还没有找到任何实际保存这些数据文件以备将来使用的程序。因此，我任意创建了一种模仿SpiroJ格式的格式。这是上面所示的Farris_1_7_-17.svg图像的xml数据文件。

<?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>

此数据文件可用作创建其他类型摩天轮数据文件的模板，具有不同的频率和幅度。一旦您使用GeometryDemonstration程序显示了有趣的形状，就可以从中复制相关的输入数据。请注意，此数据文件有三个转子而不是两个，并且x和y之间没有区别。我们还将一直假设频率始终是整数，这样对象将在有限的旋转次数内成为闭合形状；并且我们将以度为单位指定相位“phi”。

加载输入文件后，您将可以选择渲染方法，与之前一样，以及选择输出svg文件名。在本文中，我们将只讨论贝塞尔渲染方法，因为其他两种方法很简单，并且与第一篇文章中的Proceed方式相同。

程序计算顺序

计算顺序与第一篇文章中的描述相同，只是它有时会拆分然后再次合并，因为遇到了不同的文件格式。信息流显示在此图中。

如果选择了.spiro文件，它将按第一篇文章中的描述进行处理。如果选择了.xml文件，则例程SpiroParse.parse_spiroJ_file()将决定其类型。如果它是SpiroJ类型（两个转子），则SpiroConfig对话框将使用以下字段名称。

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"}};

如果它是摩天轮（三个转子），字段名称将是。

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"}};

对于SpiroJ形状，例程SpiroWrite.convertspiroJParms()将最终确定形状的类型。SpiroJ形状可能与标准棘轮图对象相同，如果x和y方向的频率和幅度相同。在这种情况下，该对象将被重定向到标准的棘轮图处理例程getspiroShape，因为它专门用于处理该形状。否则，它将由例程getspiroJShape处理，该例程更具通用性。对于摩天轮，此时无需做出任何决定，因为它清楚它们必须分开处理，在例程getFarrisShape中。

我们现在假设已选择贝塞尔渲染方法。将贝塞尔曲线拟合到已知棘轮图段的过程与第一篇文章相同，但适当的t值的选择大不相同，需要进行描述。

为SpiroJ形状选择t值

SpiroJ对象遵循以下方程

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

通常，该对象将没有旋转对称性，并且没有简单的方法对其进行子分类，除了上面提到的两个特例，即标准棘轮图和李萨如图。因此，我们将在本节中采用一种非常简单通用的方法。主要要求是，起点和终点之间的弧角绝不能大于90度。为了实现这一点，我们拟合所有具有垂直斜率的点，以及水平斜率和所有拐点。垂直斜率约束导致t的方程如下：

r_x1 w_x1 sin(w_x1t) + r_x2 w_x2 sin(w_x2t) = 0

(1)

水平斜率约束导致t的方程如下：

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

拐点确定如下：从第一篇文章的公式(3)中，我们得到了曲率的表达式。将其设为零会产生约束：ẋÿ - ẏẍ = 0。这导致t的方程如下：

r_x1 r_y1 w_x1 w_y1 (w_y1 sin(w_x1t) sin(w_y1t) + w_x1 cos(w_x1t) cos(w_y1t)) +
r_x1 r_y2 w_x1 w_y2 (w_y2 sin(w_x1t) sin(w_y2t) + w_x1 cos(w_x1t) cos(w_y2t)) +
r_x2 r_y1 w_x2 w_y1 (w_y1 sin(w_x2t) sin(w_y1t) + w_x2 cos(w_x2t) cos(w_y1t)) +
r_x2 r_y2 w_x2 w_y2 (w_y2 sin(w_x2t) sin(w_y2t) + w_x2 cos(w_x2t) cos(w_y2t)) = 0	(2)

这个方程与之前的方程不同，因为它包含三角函数乘积。但是，我们可以使用以下恒等式对其进行简化，使其不包含任何乘积：

sin(A) sin(B) = (cos(A - B) - cos(A + B))/2	(3)
cos(A) cos(B) = (cos(A - B) + cos(A + B))/2

现在我们需要一个通用的求解器，可以将任意数量的正弦和余弦波的总和设为零。上一篇文章使用了例程SpiroCalc.calc_cos_t()和SpiroCalc.calc_sin_t()来求解类似方程，但它们对于这个问题来说过于专业化。为了求解这些方程，我们需要一种通用的方法来表示方程数据：我们使用数组

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

其中第一个下标表示转子的索引号。例如，方程(2)将产生八个转子，来自八个唯一的频率。第二个下标表示三个参数：幅值、频率和相位；其中相位将以弧度表示，并将主要用于区分正弦波和余弦波。使用这种表示法，方程(1)表示为

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

而方程(2)，在使用公式(3)简化后，表示为

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}};

这些数据被发送到SpiroJCalc.solve_cos_t(double r[][], double t0)，其中将使用牛顿-拉弗森方法找到根。这里遇到的主要问题是t0的选择，即根的初始估计。会有很多根需要找到，并且总是存在收敛到错误根并偶然错过所需根的可能性。我们通过使用一定的“过度杀伤”来避免这种情况：对于每一个单独的转子（余弦波），以及这个单独转子在感兴趣范围内的每一个可能的零点，我们都会单独调用求解器。这应该会让我们非常接近实际根，如果它存在的话，假设其他转子变化不是太快。调用求解器的代码是：

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])));

没有理论保证这种方法会找到所有根，但到目前为止还没有失败。然而，它确实会产生大量的重复、冗余的根。这些根在例程SpiroJCalc.sort_t_values(double[] t, int n)中被检测和删除。

给定一组要用作拟合点的t值，我们可以像以前一样，使用SpiroJCalc.getBezier()例程来生成贝塞尔曲线拟合。下面显示了使用SpiroJ文件holub.xml的一个典型SVG结果。图也用红色显示了拟合点。请注意，这里的策略与上一篇文章存在根本性差异：以前，在使用标准的棘轮图形状时，我们*先验地*知道对称性是什么，这意味着有多少个唯一的瓣，所以我们能够只拟合第一个瓣，然后根据需要复制结果。对于当前的形状，我们不知道对称性可能是什么，所以我们必须显式计算整个对象的拟合点。这对于通常没有对称性的通用SpiroJ形状来说不是问题，但对于下面讨论的摩天轮来说，这可能是一个问题，摩天轮实际上确实具有旋转对称性，但拟合方法并没有明确知道或使用这种对称性。对于这种情况，我们将需要开发一种方法，即使我们不知道对称性的程度，也能自动尊重对象的对称性。但首先，我们需要进行一次短途旅行，研究李萨如图的一个特殊特征。

李萨如图静止点的曲率

上面为SpiroJ形状开发的方法适用于没有对称性的通用形状。然而，它在李萨如图反转其运动方向的特殊情况下会失败。这会产生一个必须特殊处理的静止点。下面给出了一个例子。这显示了两个重叠的图形：第一个图形的w_x1 = 1, w_y1 = 3，导致一个闭合图形，运动从未停止。第二个图形的w_x1 = 2, w_y1 = 3，导致一个图形，运动停止并反向，因此它看起来像一条曲线，但实际上是两条路径重叠，在两端都有静止点。这些静止点与正常棘轮图中遇到的静止点不同。在正常棘轮图中，静止点实际上是尖点，当|c|等于|b|时，曲率会有一个奇点。对于李萨如图静止点，曲率是良好的，但必须小心地进行解析评估，因为分子和分母都趋于零。在不失一般性的情况下，我们可以使用简化的表示法来表示李萨如图

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

当ẋ = 0且ẏ = 0同时发生时，将出现静止点。这要求

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

其中M和N是任意整数。在这些方程之间消去t，我们得到一个关于频率的新约束

2Nw_y = (2M + 1)w_x

(4)

对于任意的w_x和w_y，不一定总能解出这些方程。例如，如果w_x是奇数，公式(4)将永远没有解。但是，如果存在解，那么我们希望在那个特定的t值周围进行泰勒级数展开，以捕获ẋ, ẍ等表达式中的首项，这样我们就可以计算曲率。我们进行泰勒展开，变量为Δt = t - Nπ/w_x = t - (M + ½)π/w_y。这导致

ẋ ≈ -(-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)

将这个和对应的ẍ和ÿ代入曲率的表达式，即第一篇文章的公式(3)，我们得到

κ_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

其中我们只保留了首项，它们在分子和分母中都是Δt³阶，因此在Δt趋近于零时完美抵消。我们在SpiroJCalc.getBezier()中使用这个结果来计算静止点的曲率。这产生了下面用蓝色曲线显示的结果，对于w_x1 = 2, w_y1 = 3。

寻找摩天轮曲率的极值

摩天轮带来了新的挑战，因为它通常具有某种旋转对称性，但显然我们当前的拟合方法不会显示这种对称性。我们的方法当前拟合所有垂直和水平斜率。然而，如果摩天轮具有六重对称性，并且如果我们使用垂直和水平斜率来拟合一个瓣，然后将这些拟合点旋转六十度，那么很明显它们不会对应于新瓣上的水平和垂直斜率，所以新瓣不会是原始瓣旋转60度的完美克隆。这就是我们希望尊重对象对称性的意义，但问题在于对称性的确切性质可能不明显，并且形状足够复杂，以至于不清楚如何明确利用对称性。同样，对于标准棘轮图形状，很容易将形状进行子分类并为每个子类别开发特殊技术；对于摩天轮，这可能不可行。所以我们需要使用形状测量方法，这些方法可以在不明确强加对称性的情况下自动保留对称性（如果存在）。一种这样的测量是曲率，它相对于方向是不变的。因此，从曲率派生的任何度量也将是不变的。这使我们考虑使用最大或最小曲率点作为拟合点。

我们开始搜索四种类型的拟合点：

曲率绝对值最大的点
曲率绝对值最小的点
拐点（曲率为零）
切线方向与曲率最大点成90度垂直的点

顺便说一句，我们发现这组标准与最初为内摆线所做的基本相同，因为事实证明，对于标准棘轮图，最大半径点碰巧也是最大曲率点，最小半径（曲率）点也是如此。因此，当前的标准是先前用于标准棘轮图方法的自然推广。

为了计算拐点，我们使用与SpiroJ形状相同的方法，并使用公式(3)得到转子集

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}};

接下来，通过将dκ/dt设为0来获得最大/最小曲率点，得到方程

(ẋ² + ẏ²) d(ẋÿ - ẏẍ)/dt = (3/2) (ẋÿ - ẏẍ) d(ẋ² + ẏ²)/dt

(5)

这个表达式中的各个项由下式给出：

ẋÿ - ẏẍ = 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(θ₃₁)

其中 θ_ij = (w_i - w_j)t + phi_i - phi_j。将这些项进行一阶导数并代入公式(5)得到一个关于t的方程，该方程可以使用公式(3)转化为一组标准的十二个转子。这些转子太复杂，无法在此处显示，但它们在FarrisCalc.get_t_values()中明确定义，并在SpiroJCalc.solve_cos_t()中求解。

最后，相对于曲率最大点旋转90度的点是通过扫描先前获得的拟合点并寻找角度变化大于90度的点来获得的。如果发生这种情况，则计算最高曲率点的切线角度，并开始在该侧或另一侧寻找与它垂直的拟合点。作为初始估计，我们使用一个t值，该t值是最高曲率点处t值与其任一侧最近邻的t值的平均值。到目前为止，这个初始估计一直导致收敛到所需结果。下面显示了一个典型的摩天轮，使用频率（-2, 5, 19），这产生了7重旋转对称性。拟合点用红色显示。请注意，它显示了所需的对称性，即一个瓣的拟合点可以旋转（1/7）圈来直接叠加到下一个瓣的拟合点上。

摩天轮中的静止点

摩天轮也可以显示尖点，这些尖点需要特殊处理，因为在这些点处曲率是未定义的。这些尖点发生在运动暂时静止时。它们类似于正常棘轮图中|c| = |b|时出现的尖点（即当内摆线变成内摆线时）。但是，对于摩天轮，这些静止点可能出现在一系列不同的半径上，而不仅仅是单个半径，因为摩天轮有一个额外的自由度。描述这些尖点的最好方法是计算摩天轮上点的速度行为。摩天轮点的速度向量本身也是一个摩天轮，具有新的参数。如果原始摩天轮点的位置由以下类型的转子集给出

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

(6)

那么同一点的速度向量由以下类型的转子向量和给出

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

(7)

使用这个新的速度表示法，我们现在可以看到，原始摩天轮（方程(6)）中的静止点发生在所有三个速度向量（方程(7)）相加为零时。这在理论上并不总是可能的。例如，如果我们任意重命名三个转子，使得|r₁ w₁| > |r₂ w₂| > |r₃ w₃|，那么速度可能为零仅当

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

这定义了一个车轮半径可能存在的最大值和最小值，假设我们保持其他两个半径不变。如果满足此不等式，使得理论上可能存在静止点，那么我们可以计算它实际发生的特定条件。特定条件将显示为对相位偏移phi_i的约束；例如

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

其中类似的方程也适用于(phi₁ - phi₂)和(phi₃ - phi₁)。这三个约束都可以看作是余弦定律的应用，它适用于相加为零的三个向量的向量和。有了这些方程，我们现在就可以构建一个摩天轮的尖点检测器，以便*先验地*确定何时存在尖点。这在例程FarrisCalc.calc_t_cusp()中完成。有必要提前预测这些尖点，否则我们可能会在该点处得到无限大或未定义的曲率。如果检测到尖点，我们将任意将该点的贝塞尔臂长设为零，这与标准棘轮图的做法相同。

在正常情况下，这种类型的尖点的存在相当罕见，但通过故意创建它们，看看它们是什么样子很有趣。下图显示了一个动画，基于一系列32张图像，每张图像都包含尖点。它们是通过取顶部图中的摩天轮，频率为（1, 7, -17），并改变最小的车轮半径来生成一系列具有尖点的不同图形，从可能的最小半径到最大半径，中间有几个例子。这只是作为一个测试，以确认贝塞尔曲线计算在这些情况下仍然是良好的。

未来工作

摩天轮使用的方法比SpiroJ形状或标准棘轮图形状使用的方法更具通用性。通过将曲率最大化/最小化代码应用于这里描述的所有形状，可以高度标准化代码，以产生一种统一的方法处理所有三种情况。然而，这需要对已经工作良好的代码进行大量重写，所以我认为我会将它留给感兴趣的读者作为练习。