Spiro2SVG：使用贝塞尔曲线将曲线图转换为 SVG

• 下载源代码 - 75.9 KB
• 下载演示 - 41.7 KB

基于开源项目 Spiro2SVG。本文将讨论该项目的一部分，仅关注标准的回转仪形状。

引言

有几十种免费程序可用于生成回转仪形状。您可以在 Code Project、GitHub、SourceForge、Android Apps 或 Launchpad 上找到它们。大多数这些程序有两个共同点：它们通常只允许您以位图形式保存图像，并且它们通常使用一系列许多（也许是数百个）短直线来渲染图像，有时也称为折线。这是旧式数控机床的遗留问题，当时的硬件只能渲染直线，而样条只能通过 法国曲线 产生。Spiro2SVG项目旨在超越这一点，进入一种本质上是矢量图形的表示，仅使用样条。该程序会将通常以XML格式保存的回转仪数据文件转换为矢量图形格式（SVG），并使用三次贝塞尔曲线。这会导致文件大小更小，图像始终平滑，即使放大也是如此；并且是一种非常适合网站的文件格式。它还避免了一个非常棘手的问题，即需要多少条短直线才能产生高质量的图像？

使用程序

当前版本将转换由程序 SpiroGraph_1.0.2.1 生成的*.spiro文件。此程序附带一个示例文件库以及一个在线 图像库。这些图像遵循 方程

x = (a + b) cos(t) + c cos((1 + a/b)t)	(1)
y = (a + b) sin(t) + c sin((1 + a/b)t)	(2)

这些形状称为外摆线（b > 0）或内摆线（b < 0）。本系列的后续文章将介绍其他类型的形状的转换。

可以通过双击 Spiro2SVG.jar 文件来启动程序。它会要求您浏览并定位一个 .spiro 文件，其中一些文件包含在 .\Samples 文件夹中。这将显示图像的配置。对于每个对象，都会显示以下字段：

public static final String[][] spiroNames = new String[][] 
            {{"StatorRadius"}, {"RotorRadius"}, {"NumRotations"}, {"AnglesPerCycle"},
             {"RotorSlide"}, {"OriginX"}, {"OriginY"}, {"InitialAngle"},
             {"PenDistance"}, {"Lock"}, {"CurvePenWidth"}, {"Zoom"},
             {"Argb"}, {"FillArgb"}, {"FillMode"}, {"Edit Drawing Style"}};

对于每个对象，“编辑绘图样式”的底部字段可以单独编辑为“点”、“线”或“贝塞尔”。使用“线”应该会生成与原始 SpiroGraph_1.0.2.1 程序相同的图像。使用“贝塞尔”将基于较少数量的贝塞尔拟合点生成平滑插值。在某些情况下，使用“线”渲染方法可能实际上是有利的，如果图像的原始作者为了产生有些锯齿状的外观而故意使用了少量线（AnglesPerCycle）。否则，可以使用显示中的默认设置。默认设置基于这样的假设：任何每瓣线数少于 20 的对象，其意图是用“线”绘制成锯齿状；否则，将使用“贝塞尔”曲线平滑绘制。 （一瓣是回转仪的最小唯一部分，通常每个回转仪的革命中有 a/b 个瓣。）一旦选择了渲染方法，请单击“确定”以显示保存 SVG 文件的浏览对话框。该程序还可以从 CLI 静默运行，命令为 'java -jar Spiro2SVG.jar -?' 以查看可用选项。如果您想使用批处理文件转换整个目录，或者如果您希望捕获可能生成的任何错误消息，这可能会很有用。在这种情况下，可以使用选项 (-p, -l, -b) 将渲染方法强制对所有对象都相同，或者可以使用上述默认方法。

程序计算顺序

1. main.main(String[] args) - 这将启动程序，可以使用 CLI 或 JFileChooser 来指定输入数据文件。浏览 Samples 文件夹以查找合适的文件。
2. SpiroParse.parse_spiro_file(String fname, String m_style) - 这将读取数据文件并填充数组 main.rowNames 和 main.rowData，这些数组包含每个对象的字段名称和字段值。
3. SpiroConfig.showDialog(JDialog parent, String fname) - 如果程序不是静默运行的，此对话框允许选择绘图样式。此对话框以相当通用的方式编写，以允许使用来自不同类型数据文件的不同字段名称集。
4. SpiroWrite.write_svg_file(String fname, boolean m_isspiro) - 这会启动写入序列，通过写入 <svg> 标头和一些单独的 <path> 属性（如填充、描边和描边宽度）来开始写入。这些将构成 SVG 样式属性。详细的路径信息将在下面计算。
5. SpiroWrite.convertspiroParms(int index, FileWriter out) - 对于每个回转仪对象，将传入的数据字段转换为与公式 (1) 一致的格式。这是为了允许未来可能使用不同的数据源。
6. SpiroWrite.getPath(FileWriter out, PathIterator pit) - 这是一个通用的转换器，用于将 PathIterator（它操作 Java Shape 对象）转换为文本格式，该格式将指定 SVG <path> 元素中的“d”参数。在 Java 中，路径信息的流程如下：
   - 在最底层，我们有 Rectangle2D、Line2D 或 CubicCurve2D 对象，具体取决于选择的绘图样式。这些对象是在例程 SpiroWrite.getspiroShape(int index, double a, double b, double c) 中创建的。
   - 这些被连接成一个 Path2D 对象以形成复杂的形状。
   - 然后我们使用 AffineTransform.createTransformedShape(Shape pSrc) 操作最终路径，该操作将对每个对象执行整体定位、缩放和可能的旋转。
   - 这会产生一个 Shape 对象，然后我们使用 Shape.getPathIterator() 来检查它，以生成有关路径每个段的详细信息。
   - PathIterator 信息然后馈送到 SVG <path> 的“d”参数中。
7. SpiroWrite.getspiroShape(int index, double a, double b, double c) - 这是 Java 路径信息的来源。它将返回 Line2D、Rectangle2D 或 CubicCurve2D 段的序列。在前两种情况下，处理在此点完成。但是，对于贝塞尔绘图样式，有趣的部分才刚刚开始。贝塞尔拟合过程是分两步进行的。首先，我们选择一组有趣的 t 值，用于公式 (1) 来指定每个段的起点和终点，使用 SpiroCalc.get_t_values(t_values, a, b, c)。然后，我们使用 SpiroCalc.getBezier(t_old, t_new) 将贝塞尔曲线拟合到每个段。这将在下面描述。

在内摆线（b < 0）上查找 t 值

首先是一些假设：我们假设 a > 0（定子半径）和 b < 0 用于内摆线，因此转子在定子内部滚动。我们还假设 c > 0，尽管实际上这并不重要。更改 c 的符号等效于将整个对象旋转 πb/a，这不会以任何方式改变形状。为了拟合内摆线，只需要分析第一个瓣，该瓣的范围是 t ∈ (0, 2πb/a)。后续的瓣将只是第一个瓣的旋转克隆。此外，每个瓣都有反射对称性，因此如果 t 对应于曲线的特定特征，那么 2πb/a-t 也将对应于瓣的反射版本的相同特征。我们希望选择 t 值来对应容易识别的特征，例如最大半径、最小半径、曲率为零的拐点以及运动暂时静止的尖点。此外，我们需要施加一个非平凡的约束，即最大弧角不能大于 90 度。通过弧角，我们指的是段的起点处的切向量与终点处的切向量之间的角度差。如果该角度大于 90 度，则无法获得可靠的贝塞尔拟合。

接下来是一些定义：我们将斜率定义为 dy/dx = ẏ/ẋ，其中 ẋ 表示 dx/dt。但是，如果该点是静止的，意味着 x 和 y 是 t 的二次函数而不是线性函数，那么我们将 dy/dx 定义为 ÿ/ẍ。曲率 定义 为

κ = (ẋÿ - ẏẍ) / (ẋ2 + ẏ2)3/2

(3)

在尝试计算上面描述的拟合点之前，我们需要首先对内摆线形状进行分类，以确定这些特征是否存在。分类取决于 c 的大小。有三个特殊值 c

• c = 0，这将产生一个圆。
• c = b/(1 + a/b)，这是产生拐点的最小 c 值。
• c = -b，这将产生一个内摆线，它有一个尖点（静止点）。

对于 c > b/(1 + a/b)，每个瓣将有两个拐点。
对于 c > -b，拐点消失，被一个额外的环（凸形）取代，这将需要创建两个新的拟合点，这两个点位于距离最大半径点 90 度的地方，因此它们发生在新的环的最大宽度点。对于内摆线（第一个瓣），这是斜率为零的点。

有了这个分类，我们可以按如下方式计算拟合点：

1. 最大半径将始终出现在 t = 0
2. 最小半径将始终出现在 t = πb/a
3. 尖点（如果存在）将出现在 t = 0
4. 拐点（如果存在）将出现在
   cos(at/b) = (1 + c²(1 + a/b)/b²) b / c / (2 + a/b)
5. 额外的环的最大宽度点，如果存在，将是方程的解
   cos(t) = -(c/b) cos((1+a/b)*t)
   该方程使用牛顿-拉夫森方法在 calc_cos_t(double fAmp, double fFreq) 中求解。

下图显示了五个典型的内摆线，参数值为 a = 210, b = -70, c = {10, 30, 50, 70, 90}。这五种形状说明了以下情况：无拐点、两个拐点、一个尖点和一个额外的环。在每种情况下，第一个瓣的拟合点显示为红色矩形。

在 Eptrochoid（b > 0）上查找 t 值

在 Eptrochoid 上查找合适的 t 值的过程比内摆线稍复杂，如下所示。对于内摆线，如果我们考虑 a/b = N 的情况，其中 N 是整数，如果我们计算切线在整个曲线回到起点时所做的完整旋转次数，我们会发现切向量执行了 N-1 次旋转。唯一瓣的数量是 N。因此，如果我们实现一个每瓣产生四个拟合点的过程，那么我们可以保证每段的弧角永远不会超过 90 度，这是必需的。这是假设点选择得当，这很容易做到。然而，对于 Eptrochoid，切向量的旋转次数是 N+1，所以每瓣四个拟合点将不足够，尤其是在 N 较小的情况下。

除了之前为内摆线定义的拟合点外，我们现在考虑一些新的候选点：

6. 如果 c > b，则每个唯一瓣中将有两个环，一个大的（凸形）通过最大半径，一个小的（凹形）通过最小半径。尝试在这两个环的最大宽度点进行拟合是有吸引力的，这将比以前多一个拟合点。然而，对于大的环，存在一个问题，即最大宽度点有时会消失。很难事先预测何时会发生这种情况。我们只知道，随着 c 的减小，大环的最大宽度点将与拐点合并并消失，留下一个无解的方程。由于我们无法预测何时会发生这种情况，因此我们需要一种替代方法。替代方法是首先找到小（凹形）环的最大宽度点，然后选择一个距离该点 90 度的另一个拟合点，使其切线与最小半径点的切线平行。这保证存在（事实上，即使小环不存在，它也存在），并且满足每段弧角不超过 90 度的要求。此解决方案将使用 calc_sin_t(double fAmp, double fFreq) 找到。
7. 实验发现，如果一个端点是尖点或者甚至接近尖点（c ≈ b），并且总弧角接近 90 度，那么解决方案可能会表现不佳。这通常表现为负的贝塞尔臂长（见下文）。为防止这种情况，我们可能在最接近尖点的两个点之间插入一个插值 t 值。一个额外的拟合点可以大大提高这种情况下的精度。
8. 在另一种情况，即不存在尖点，但我们处于拐点或甚至接近拐点（c ≈ b/(1 + a/b)）时，已经发现解决方案可能同样表现不佳，因此我们也在此处插入一个额外的插值 t 值。
9. 在 a < 2b（瓣数小于 2）的情况下，有必要在大（凸形）环的最大宽度点插入一个新的 t 值，以避免弧角大于 90 度。

下图显示了四个典型的 Eptrochoid，参数值为 a = 150, b = 50, c = {10, 30, 50, 70}。这四种形状说明了以下情况：无拐点、两个拐点、一个尖点和一个额外的环。在每种情况下，第一个瓣的拟合点显示为红色矩形。

将三次贝塞尔曲线拟合到回转仪图形

三次贝塞尔曲线在 t ∈ (0,1) 的范围内由以下方程参数化定义：

x = x0 (1-t)3 + 3 x1 t(1-t)2 + 3 x2 t2(1-t) + x3 t3	(4)
y = y0 (1-t)3 + 3 y1 t(1-t)2 + 3 y2 t2(1-t) + y3 t3	(5)

该曲线有八个未知数：其中四个将使用起点和终点位置指定。另外两个将通过匹配从 (x₀, y₀) 到 (x₁, y₁) 以及从 (x₃, y₃) 到 (x₂, y₂) 的贝塞尔控制臂的斜率来指定。最后，公式 (3) 中的起点和终点曲率将用于指定贝塞尔控制臂的长度。拟合方程的实际求解将在 main.calcBezier(Point2D.Double[][] ptSpiro, double t1, double t2, double max_v) 中执行。但是，在此之前，有必要在 SpiroCalc.getBezier 中进行一些预处理，以确定两个端点的曲率和切线方向。（这两个计算都取决于端点是否静止，因此有两种不同的公式可以完成此操作。）切线角度很重要，原因如下：如果我们使用原始曲线方向求解方程，我们将不得不处理大量特殊情况，具体取决于切线角度在开始或结束时是否垂直，以及开始或结束时的曲率是否为零，或者是否为尖点等。为了尝试减少特殊情况的数量，我们将重新定向所有曲线段，使其起始和结束切线永远不会垂直，通过使用介于两个方向之间的方向。这是因为曲率与方向无关。解决方案完成后，我们将在使用之前将解决方案转换回原始方向。接下来，我们定义数组

Point2D.Double[][] ptSpiro = new Point2D.Double[3][2];

这包含了关于每个单独回转仪段的所有必要信息。第一个索引 (0-2) 指定这是回转仪的位置、相对于 t 的一阶导数还是二阶导数，第二个索引 (0-1) 指定起始点或终点。这就是我们将尝试通过调整贝塞尔曲线来匹配的信息。ptSpiro 数据被发送到 main.calcBezier。该例程将使用例程 getrotX() 和 getrotY() 将数据转换为新的旋转框架，然后再求解方程。

在拟合贝塞尔曲线时，唯一 nontrivial 的部分是拟合曲率的过程。在此之前，公式 (3) 中 κ 的原始定义将稍微简化。表达式中的分母可以用以下替代形式表示：

(ẋ² + ẏ²)^3/2 = ẋ³ (1 + α²)^3/2

其中 α 是斜率 dy/dx。但是，我们将约束贝塞尔斜率等于回转仪斜率，因此涉及 α 的项对于这两个对象都是通用的，可以提取出来。这意味着我们将匹配一个简化的曲率表达式，形式为：

Cu = (ÿ - α ẍ) / ẋ²

其中 Cu 是一个常数，将使用回转仪参数进行评估，右侧是基于贝塞尔参数的函数。将其显式写出，我们得到以下两个约束：

Cu₀ = (2/3) ((y₀ - 2y₁ + y₂) - α₀ (x₀ - 2x₁ + x₂)) / (x₁ - x₀)²
Cu₁ = (2/3) ((y₁ - 2y₂ + y₃) - α₁ (x₁ - 2x₂ + x₃)) / (x₃ - x₂)²

其中 Cu_i 指的是起点或终点处的简化的回转仪曲率，α_i 是相应的斜率：

α₀ = (y₁ - y₀) / (x₁ - x₀)
α₁ = (y₃ - y₂) / (x₃ - x₂)

现在我们引入未知贝塞尔臂长，使用符号：

δ₀ = x₁ - x₀
δ₁ = x₃ - x₂

这将产生最终方程组：

(3/2) Cu0 δ02 =   (y3 - y0) - α0 (x3 - x0) - δ1 (α1 - α0)	(6)
(3/2) Cu1 δ12 = -(y3 - y0) + α1 (x3 - x0) - δ0 (α1 - α0)	(7)

在这些方程中，只有臂长 δ₀ 和 δ₁ 是未知的。现在我们可以识别一些特殊情况：

• α₀ = α₁ : 端点斜率平行。在这种情况下，方程是解耦的，可以独立求解。
• Cu_i = 0 : 如果端点曲率为零，则该方程对于 δ 是线性的，解法也很直接。
• 一般情况：在这种情况下，我们将公式 (6) 对于 δ₁ 求解为 δ₀ 的函数，并将其代入公式 (7) 以获得 δ₀ 的四次方程。这将在 solve_quartic() 中解析求解。在四个解中，我们选择符号正确的那个，因为我们期望两个贝塞尔臂都朝内，相互指向。

在接受曲线拟合结果之前，还有一个最终问题需要处理。如果端点的曲率为零，或者端点是静止的，并且弧角接近 90 度，则可能发生贝塞尔臂长略为负值，指向错误的方向。如果发生这种情况，臂长将被任意设置为零，这意味着另一个端点的曲率将不完全正确。然而，这种情况相当罕见，而且比贝塞尔臂方向反转可能引起的明显缺陷要好。最后，作为一种诊断工具，如果您想在视觉上定性地确认解决方案，可以使用节点工具在 Inkscape 中轻松显示贝塞尔控制臂。

结论

总而言之，我们得到了类似这样的结果，这是一个基于 SpiroGraph_1.0.2.1 的示例文件 *00_Animation_Simple.spiro* 的 SVG 文件。