用epicycles绘制图像：从轮廓到傅立叶级数
导言
在这个项目中，位图图案变成一个封闭的轨迹，然后变成旋转矢量的总和。 每个下一个向量的开始都在前一个向量的末尾。 链条的最后一端起到笔的作用，逐渐再现绘图。 向量的末端相对于它们的开始移动的圆圈称为epicycles。
该方法的基础是复杂的傅立叶级数。 任何足够规则的周期性轨迹都可以通过具有不同半径、初始角度和整数频率的谐波来描述。 有限的点集合是已知的，因此系数是通过离散傅立叶变换（dft）获得的。 旋转矢量的几何解释由[G.Sanderson，2019]解释（https://www.3blue1brown.com/lessons/fourier-series /);复变换描述于S.W.Smith,ch. 31。
主要的实际困难出现在转换之前：像素需要变成有序的行走，必须协调绘图的各个部分，并且必须选择运动参数。 DFT不知道如何自行查找线或确定其通过的顺序。 存储的谐波越多，所制备的轨迹的近似性越好，链的图像越丰富。
该方法适用于复杂数字和光谱分析的教育演示，标志和铭文的动画，以及扁平曲线的紧凑表示。 傅立叶系数在形状描述和识别中的应用具有独立的历史：参见[Zahn和Roskies，1972]（https://doi.org/10.1109/TC.1972.5008949 )。 在他们的工作中，使用了另一个曲线描述符;这里直接转换其点的坐标。
完成工作后，您将能够解释振幅与段长度，频率与转数，相位与初始角度的关系，并且还可以根据测量的误差选择谐波数量。
设置任务
**入口：**一个静态PNG，JPEG或BMP图像。 PNG是首选：白色背景上的黑线，深色背景上的白线，或透明背景上的对象。 彩色标志也被接受：不透明的图像被转换为灰度，并与透明度，您可以选择一个剪影使用阿尔法通道。 选择蒙版后不会保存源调色板。
它已经包含在包中 engee-poster.png 具有AnyMath徽标和名称的图像的示例。 对于新图像，通过AnyMath文件浏览器将文件上传到项目文件夹，并在掩码中指定名称。 路径字段选择上传的文件;在掩码中没有单独的窗口小部件用于发送文件。
**结果：**循环GIF动画的一个完整的爬行与选定的FPS和暂停完成的绘图。 在文件夹中 output 二进制掩码、阶段插图、最后一帧、CSV中的系数和实际错误报告也被保存。 主要参数汇总于表1。
表1-实验参数
|
团体 |
参数和含义 |
|
图像 |
路径；自动/暗/亮/阿尔法；自动或手动阈值；处理尺寸；边框或笔画中间 |
|
离散化;离散化 |
N=沿着路线长度均匀分布的512、1024、2048或4096个点 |
|
方法 |
精度模式:百分比从1到100;段数模式:K从1到300. 只考虑所选模式 |
|
颜色 |
背景;旋转段;绘制轮廓;背景上的初始轨迹;圆圈。 值由字符串#RRGGBB给出 |
|
动画制作 |
FPS从5到50;绘图持续时间从2到30秒;暂停从0到3秒;显示圆圈 |
为了节省合理的计算时间，文件大小（10MiB），解码像素数（1200万），处理掩码的侧面（800像素），遍历的复杂性（6000个顶点和40个部分），样本数（4096），段数（300）和GIF帧数（900）受到限 如果超过复杂度，则建议缩小图像;如果精度过高，则限制分段数量并显示实际达到的百分比。
重要的数学运算直接显示在单元格中。 像素边界遍历、骨架化、服务可视化和GIF记录在 epicycles_functions.jl. 参数采用原生[AnyMath code cell masks]设计（https://engee.com/helpcenter/stable/cn/feature/script-editor-masks.html )。 更改掩码后，执行它和所有后续单元格;对于完整的重新计算，请使用"全部运行"。
using Images, FileIO, FFTW, Plots, Statistics, Markdown
# @markdown启用课程功能并准备结果文件夹。
include(joinpath(@__DIR__, "epicycles_functions.jl"))
lesson_dir = @__DIR__
output_dir = joinpath(lesson_dir, "output")
mkpath(output_dir)
gr(format=:png)
nothing
输入图像处理
假设并获得二进制掩码
绘图应具有对比对象和相对均匀的背景。 要自动检测暗或亮的物体，假设图像的边缘主要属于背景。 如果对象接触大部分边缘，请手动选择模式。
在第一步骤中，形成二元掩模。：意思 true 对应于对象, false -背景。 如果有明显的透明度,自动模式使用一个阈值为0.5的alpha通道. 这对于所附的徽标很方便：其填充的颜色不会干扰轮廓提取。
亮度用于不透明图像。 通过Otsu方法找到自动阈值：选择直方图划分，最大化类间方差。 见[大津，1979]（https://doi.org/10.1109/TSMC.1979.4310076 )。 手动阈值从0设置为1。 包含少于10个像素的小的孤立组作为噪声被去除。
# @markdown**加载和选择对象**
путь_изображения = "engee-poster.png" # @param {type:"string",placeholder:"input_logo.png"}
выделение_объекта = "Авто" # @param["Auto"，"Dark"，"Light"，"Alpha"]
авто_порог = true # @param {type:"boolean"}
порог = 0.5 # @param {type:"slider",min:0.0,max:1.0,step:0.01}
размер_обработки = 608 # @param {type:"slider",min:128,max:800,step:16}
тип_контура = "Граница" # @param["边界"，"中风中间"]
nothing
在下一个单元格中，掩码模式转换为函数参数。 相对路径相对于项目文件夹进行解析。 仅当超出所选边时，图像才会缩小；小源不会放大。 然后以处理的图像的像素为单位测量坐标和误差。
mode = Dict("Авто"=>:auto, "Тёмный"=>:dark, "Светлый"=>:light, "Альфа"=>:alpha)[выделение_объекта]
input_path = isabspath(путь_изображения) ? путь_изображения : joinpath(lesson_dir, путь_изображения)
mask, image_meta = EpicycleLesson.load_mask(input_path;
foreground=mode, threshold=авто_порог ? nothing : порог,
max_dimension=размер_обработки)
image_meta
有序的轮廓和部分的连接
**"边框"**模式绕过与背景相邻的像素单元的边缘。 外部边界和孔的边界被保留。 粗圈描边可以有两个轮廓–这是遍历其边缘的预期结果。
"中划"模式**首先使用算法[Zhang–Suen,1984](https://doi.org/10.1145/357994.358023 ），然后跟踪骨架的像素链。 它是有用的手写铭文和线条画。 为提供的徽标选择边框旁路。 小于12像素的非常短的部分被排除在外。
Epicycle链的一端连续绘制一条路线。 因此，各个部分由桥梁连接。 我们从最长的轮廓开始，然后选择我们已经构建的路线的最近的一对顶点和剩余的部分之一。 我们插入一个过渡，完全绕过新的部分，并沿着相同的过渡返回。 打开的行程绕到两端并返回到入口点。 重复，直到所有部分都打开。
这是一种贪婪的连接方式：它给出了短的局部交叉，但不保证全局最短路线。 桥通过两次并保持目标轨迹的一部分。 它们将在完成的绘图中可见;裁剪它们意味着更改连接的任务。
using Plots, Markdown
paths, closed = тип_контура == "Граница" ? EpicycleLesson.boundary_paths(mask) : EpicycleLesson.centerline_paths(mask)
keep = [i for i in eachindex(paths) if length(paths[i]) >= 2 && EpicycleLesson.path_length(paths[i],closed[i]) >= 12]
paths, closed = paths[keep], closed[keep]
isempty(paths) && error("После выделения не осталось значимых контуров. Измените порог или размер обработки.")
center = complex(size(mask,2)/2, -size(mask,1)/2)
paths = [p .- center for p in paths]
route, bridges = EpicycleLesson.join_paths(paths,closed)
geometry_figure = EpicycleLesson.geometry_plot(mask,paths,closed,bridges)
savefig(geometry_figure,joinpath(output_dir,"01_contours.png")); display(geometry_figure)
Markdown.parse("*Рисунок 1 – Бинарная маска, выделенные контуры и соединительные мосты.*")
.png)
图1以黑色显示蒙版，边框以蓝色显示，桥点点缀。 对于附加的PNG，自动选择给出了12个轮廓和11个桥。 如果背景突出而不是可识别的轮廓，请在切换到傅立叶变换之前更改对象模式或阈值。
长度参数化和统一计数
像素顶点不能简单地以任何顺序转移到DFT。 首先，它们形成连续的路线，然后在路线上选择沿着行进长度的等距计数。 这消除了笔速度对源顶点密度的依赖性。
让 -封闭路线的连续顶点, -累积长度, -全长，以及 -计数的数量。 样品的长度和位置的定义在（1）中给出:
为 该点通过沿边（2）的线性插值找到:
在路线结束时，重复第一个顶点以描述最后一条边。 阵列中 没有重复的端点：周期参数位于半间隔上 . 它将在闭包可视化时出现。
# @markdown**目标路由计数数**
число_отсчётов = 2048 # @param [512,1024,2048,4096] {type:"raw"}
nothing
平面点写成复数 . 像素垂直坐标改变标志，使轴 它指向上方。 移动一半的宽度和高度中心的图像,以方便图形. DFT的恒定分量仍然由笔在路线上的平均位置确定。
电话号码 设置目标路由的表示的详细信息。 它没有指定GIF帧的数量。 对于一个小标志，从2048个点开始就足够了;增加 当路线中有小细节时，这是有意义的。
z = EpicycleLesson.resample_path(route, число_отсчётов)
N = length(z)
(отсчётов=N, участков=length(paths), мостов=length(bridges),
длина_маршрута_px=round(EpicycleLesson.path_length(route,false);digits=2))
离散傅立叶变换
复杂轨迹的分析与恢复
对于连续的周期性轨迹，复数傅立叶级数写为（3）。 这里的参数 -一个完整的旁路的比例，和频率 它是用每个周期的转数来衡量的:
我们的计算使用 离散点。 正向变换(4)计算系数，逆和(5)恢复节点中的值。:
指数符号和规范化必须相互匹配。 在朱莉娅 fft 计算总和而不除以 ，因此，单元格明确使用 fft(z) / N. 功能 ifft 它已经包含了一个分为 ;因此，对于这里选择的归一化，逆计算具有形式 ifft(C) * N. 定义由[AbstractFFTs documentation]检查（https://juliamath.github.io/AbstractFFTs ...jl/stable/api/);在Smith,ch. 31。
**FFT和DFT：**DFT是数学变换本身（4），FFT是计算它的快速方法。 直接求和需要订单 操作，FFT-顺序 . 这里使用FFTW包。
Julia数组从一个开始索引: C[1] 回应 . 其余的元素需要与符号频率进行比较。 用偶数 我们使用订单 . 负频率表示顺时针旋转，正频率表示逆时针旋转。 对于奈奎斯特频率，接受负符号，如 fftfreq.
信号 是复数的，因此，其系数一般不具有实信号的频谱的共轭对称性。 不可能只保存正频率或人为地将它们成对组合。
using FFTW, Plots, Markdown
C = fft(z) / N
frequencies = [q <= (N-1)÷2 ? q : q-N for q in 0:N-1]
offset = C[1]
spectrum_figure = EpicycleLesson.spectrum_plot(C,frequencies)
savefig(spectrum_figure,joinpath(output_dir,"02_spectrum.png")); display(spectrum_figure)
Markdown.parse("*Рисунок 2 – Спектр комплексной траектории и гармоники, упорядоченные по амплитуде.*")
.png)
图2显示了路由描述中涉及的频率。 常数分量 从图表中排除：它设置链的开始，不是旋转段。 其值等于平均值 ，其可能与图像的几何中心不同。
段数和测量精度
让 包含索引 存储的非恒定谐波。 截断轨迹由（6）定义:
我们按降序对谐波进行排序 . 这种选择节省了最多的能量，并最大限度地减少了RMS误差在初始离散节点在固定 . 与此同时，频率仍然是它们自己的：链中的顺序改变了中间环节的位置，但笔的数量和末端不依赖于它。
让我们介绍一个RMS错误 、路线规模 和精度指标 到(7):
例如，**99%**表示RMS误差不超过1% . 这是准备好的路线，包括桥梁的坐标的准确性的指标。 它不测量源中匹配像素的比例，也不限制最大局部误差。
谐波的正交性和Parseval相等给出（8）;有关原理的解释，请参阅Smith，ch。 10。 在这里，归一化已经被简化为复数系数（4）:
因此，对于所请求的百分比 找到最小值就足够了 ，其中光谱的尾部满足条件（9）:
"段数"模式直接设置 ;然后只计算百分比。 偏移量 它不包括在段数中。 如果给定百分比需要超过300个分段，则需要300个分段，并显示所达到的精度。 没有这样的限制的100％请求将需要高达 我会用机器精度重建节点。
# @markdown**谐波选择**
режим_точности = "По точности" # @param["按精度"，"按段数"]
точность_процентов = 99.5 # @param {type:"slider",min:1.0,max:100.0,step:0.1}
число_отрезков = 89 # @param {type:"slider",min:1,max:300,step:1}
nothing
下一个小区根据（9）调用一个短的选择过程，然后在节点中执行一个独立的反向计算。 该报告显示实际百分比、RMS和最大误差。 后者对于估计角度和交汇点很有用：即使有明显的局部偏差，平均误差也可能很小。
using FFTW, Statistics
selection = EpicycleLesson.select_harmonics(C; accuracy_percent=точность_процентов,
n_segments=режим_точности == "По числу отрезков" ? число_отрезков : nothing)
selected = selection.indices
C_kept = zeros(ComplexF64,N); C_kept[1] = offset; C_kept[selected] = C[selected]
node_errors = abs.(ifft(C_kept)*N .- z)
rms_error = sqrt(mean(abs2,node_errors))
scale = sqrt(mean(abs2,z .- offset))
actual_accuracy = clamp(100*(1-rms_error/scale),0,100)
model = (samples=z, offset=offset, coefficients=C[selected], frequencies=frequencies[selected],
accuracy=actual_accuracy, rms_error=rms_error, max_error=maximum(node_errors), capped=selection.capped)
(отрезков=length(selected), нужно_для_запроса=selection.needed, предел_сработал=selection.capped,
точность_процентов=actual_accuracy, RMS_px=rms_error, максимум_px=model.max_error)
检查恢复的路由
让我们比较图3中的目标路线和截断总和（6）。 为了显示，总和在8192点的更密集的网格上进行插值。 对于少量谐波，角度被平滑；异常值可能出现在转弯和桥梁附近。 提高精度会使所选路线更接近，但不会消除从图像中突出显示它的错误。
轮廓样本的数量和存储的谐波的数量解决了不同的问题：第一个描述了初始数据，第二个限制了它们的频谱表示的复杂性。
using Plots, Markdown
reconstructed = EpicycleLesson.reconstruct(model)
comparison_figure = EpicycleLesson.comparison_plot(route,reconstructed,model)
savefig(comparison_figure,joinpath(output_dir,"03_reconstruction.png")); display(comparison_figure)
Markdown.parse("*Рисунок 3 – Целевой маршрут с мостами и его восстановление выбранными гармониками.*")
.png)
绘制epicycles
旋转矢量和部分和链
欧拉公式 将每个谐波变成一个旋转矢量。 如果 -选定系数, -其频率，半径，初始相位和矢量由（10）确定:
链的开始是在 . 其余节点是部分和（11）:
圆 它有一个中心 和半径 . 中心本身与先前的链接一起移动;每个链接的长度保持不变。 cis(θ) 在Julia中，它计算一个复数 ，而 cumsum 建立部分总和。
对于段的完整绘图与 逆时针旋转三圈， -顺时针五个。 在绘图期间 秒物理角速度等于 弧度每秒。
颜色和单相查看
例如，使用六个十六进制数字设置五种颜色 101827. 符号 # 颜色之前添加一个代码。 "轮廓"颜色应用于笔迹;"背景上的轨迹"颜色应用于准备好的路线，包括桥梁。 将深色背景替换为浅色背景时，请选择对比鲜明的线条颜色。 隐藏圆圈会影响设计；所有选定的线段继续参与总和。
# @markdown**绘制颜色，十六进制代码**
цвет_фона = "f3f3f3" # @param {type:"string"}
цвет_отрезков = "87d1f1" # @param {type:"string"}
цвет_контура = "1b84b8" # @param {type:"string"}
цвет_траектории_на_фоне = "e5e5e5" # @param {type:"string"}
цвет_окружностей = "cccccc" # @param {type:"string"}
style = (background="#"*цвет_фона, segments="#"*цвет_отрезков, contour="#"*цвет_контура,
reference="#"*цвет_траектории_на_фоне, circles="#"*цвет_окружностей)
nothing
# @markdown**预览阶段和GIF设置**
фаза = 0.35 # @param {type:"slider",min:0.0,max:1.0,step:0.01}
показать_окружности = true # @param {type:"boolean"}
fps = 25 # @param {type:"slider",min:5,max:50,step:1}
длительность_секунд = 12.0 # @param {type:"slider",min:2.0,max:30.0,step:0.5}
пауза_секунд = 1.0 # @param {type:"slider",min:0.0,max:3.0,step:0.5}
nothing
在下面的单元格中，（10）和（11）是字面上实现的。 阵列 nodes 它包含链的开始和所有段的结束;最后一个元素设置笔的位置。 图4允许您在创建整个动画之前查看所选阶段。
using Plots, Markdown
vectors = model.coefficients .* cis.(2π .* model.frequencies .* фаза)
nodes = vcat(offset, offset .+ cumsum(vectors))
limits = EpicycleLesson.canvas_limits(model,route,reconstructed)
phase_figure = EpicycleLesson.frame_plot(model,route,reconstructed,фаза,style,limits;
show_circles=показать_окружности)
savefig(phase_figure,joinpath(output_dir,"04_epicycles.png")); display(phase_figure)
Markdown.parse("*Рисунок 4 – Цепочка эпициклов и след пера в выбранной фазе.*")
.png)
时间采样和GIF记录
让 -FPS, -所要求的绘图时间。 成形;成形 帧并包括根据(12)的周期的两端:
在第一帧 笔位于路线的开始，在最后 绘图关闭。 然后重复帧被添加到暂停。 Epicycles是为当前的一个绘制的 ，并且所述迹线取自直至该相位的级数的密集插值。 因此，轮廓在低FPS下不会变成粗糙的折线。
记录是通过 Plots.Animation, frame 和 gif 在[文档图]中描述（https://docs.juliaplots.org/stable/animations /)。 为帧选择图形后端GR。 对渲染性能的评估是必须的 每个圆上有固定数量的点；大 增加了链的时间和视觉密度。
GIF以0.01s的增量存储帧延迟。 因此，可以通过交替的相邻延迟来近似任意设定的FPS。 值25FPS每帧正好给出0.04秒，默认情况下被选中。 浏览器播放器还可以限制播放速度。
具有大频率模数的谐波快速旋转。 对于罕见的时间样本，由于频率的时间重叠，它们的运动可能看起来突然或具有可见的反向方向。 对于平滑度，增加绘图持续时间或FPS。 这是运动观察的一个属性;计算的笔位置遵循公式（6）。
使用标准参数，创建300个绘图帧和25个暂停帧：共325帧，13秒播放。 预览阶段不会影响GIF中的完整旁路。 下一个单元格将报表、系数和动画保存到文件夹中 output.
using Plots, Markdown
EpicycleLesson.save_results(model,mask,image_meta,output_dir)
gif_result = EpicycleLesson.render_gif(model,route,joinpath(output_dir,"epicycles.gif");
style, fps, duration=длительность_секунд, hold_seconds=пауза_секунд,
show_circles=показать_окружности)
display(gif_result.animation)
Markdown.parse("*Рисунок 5 – GIF: $(gif_result.frames) кадров; задано $(fps) FPS; файл output/epicycles.gif.*")
.gif)
结论
工作中构造了以下序列：位图→二进制掩码→有序轮廓→带桥的路线→均匀计数→DFT→谐波采样→epicycle链→GIF。 快速变换是在几何形状准备好之后使用的。 路线搜索及其参数化直接影响所得谱。
研究的基础是：平面坐标的复杂表示，欧拉公式，正向和反向DFT，幅度和相位，符号频率，常数分量，谐波正交性和Parseval相等，归一化RMS误差，长度参数化以及空间和时
**重复使用。**在脚本旁边上传一个新的PNG/JPEG/BMP，更改掩码中的路径并检查图1。 然后选择计数数和精度，检查图3，设置颜色和FPS并记录GIF。 服务档案 .jl 它应该留在脚本旁边。 对于没有训练文本的函数的多种使用，您可以连接模块 EpicycleLesson 从一个普通的朱莉娅剧本。
巩固的实验。
- 对于所附徽标，比较90％，95％和99％：记下段数，RMS和最大误差。 解释为什么将段数增加约100％变得昂贵。
- 修正段数并增加N.比较误差并检查目标轨迹本身改变了多少。
- 比较厚厚的手写字母上的"边框"和"笔画中间"。 说明区域边界与笔线的区别。
- 更改固定频谱的持续时间和FPS。 评估旋转的平滑度和已经恢复的足迹的形状。
**进一步开发：**根据路线的总长度选择桥梁;进入路段之间的笔升降模式;直接与SVG一起工作;将栅格边界平滑为DFT;比较幅度采样，同时保持低频;导出MP4以
周期清楚地显示了复杂的运动是如何建立在简单的谐波. 对于形状分析任务，相同的系数可以作为特征，而对于轨迹生成，它们可以作为坐标的方便的周期性表示。
资料来源一览表
文学与理论
- 史密斯S.W.科学家和工程师的数字信号处理指南。 加州技术出版社，1997年。 第31章：复傅里叶变换(PDF);[第10章：Parseval的关系](https://www.dspguide.com/ch10/7.htm ）;[第12章：FFT的工作原理]（https://www.dspguide.com/ch12/2.htm )。 全面的DFT，归一化，节能和快速算法。
- Zahn C.T.，Roskies R.Z.平面闭合曲线的傅立叶描述符。 计算机上的IEEE交易，C-21(3), 269-281, 1972. [DOI：10.1109/TC.1972.5008949]（https://doi.org/10.1109/TC.1972.5008949 );[文章全文(PDF)](https://www.cs.jhu.edu /~米莎/论文/Zahn72.pdf)。 谱描述符对平坦闭合曲线的应用。
- Otsu N.A Threshold Selection Method from Gray-Level Histograms。 IEEE关于系统、人和控制论的交易, 9(1), 62-66, 1979. [DOI：10.1109/TSMC.1979.4310076]（https://doi.org/10.1109/TSMC.1979.4310076).基于直方图的自动阈值选择。
- 张T.Y.，孙C.Y.一种快速并行算法细化数字模式。 ACM的通讯, 27(3), 236-239, 1984. [DOI：10.1145/357994.358023]（https://doi.org/10.1145/357994.358023）。二元笔划的骨架化。
- Sanderson G.（3Blue1Brown）。 但什么是傅立叶级数？ 从热流到画圈，2019年。 【作者图文并茂讲座】（https://www.3blue1brown.com/lessons/fourier-series /)。 傅立叶级数的几何解释。
实施文件
- AnyMath：[code cell masks]（https://engee.com/helpcenter/stable/cn/feature/script-editor-masks.html );文本标记。
- AbstractFFTs。jl:[fft、ifft和fftfreq的定义](https://juliamath.github.io/AbstractFFTs ...jl/stable/api/)。
- 情节。JL:创建动画。