一、问题引人
工作上遇到一个算法问题,需要从一张图片的中心点开始,生成一条螺旋线并沿着螺旋线进行特征提取。参考一些资料,我选择了一条比较常用的螺旋线方程,如下所示:
上式中x和y表示坐标点,t表示时间,可见x和y都是关于t的函数。当我对t按照步长为1来递增时,效果如下:
t以步长为1进行增长
在生成螺旋线时,发现相邻两个点之间的距离越来越大,而我需要寻找的是螺旋线上等间距分布的点集。
二、数学建模
为了解决上述的问题,要设计了一个算法,在螺旋线方程的基础上引入一个限定条件,求出t的“合理值”使得相邻两个点之间距离相等,算法思路如下:
步骤一:给定t一个初始值,比如t等于4,带入螺旋线方程求出p(x,y);
步骤二:以p为圆心,r为半径生成一个小圆,求出圆与螺旋线之间的交点;
步骤三:选择距离上一个圆心最近的交点作为新的圆心,再次生成一个圆,圆的半径等于r,继续求圆与螺旋线之间的交点,以此类推,直到生成的点集个数达到一定阈值。
如图所示,红色点表示圆心,蓝色代表小圆,为了更好的示意,我将圆缩小了,实际上每个红色的点恰好分布在前一个圆的边上,所有红色的点集合就是要求的螺旋线上等间距分布的点。假设圆心的坐标用(a,b)表示,半径为r,每次求圆与螺旋线的交点都是求一个方程组的解。
乍一看觉得这个方程很好解,只要把x和y带入最后一个式子,只剩一个自变量,然而经过漫长的求解最终幡然醒悟,这个方程没有解析解(即无法用解析式可以表达的解),只能用数值逼近的方法求近似解。世界上有很多方程都没有解析解,但是解确实存在,这点确实也让人感到好奇和感叹,并非所有的问题都能依靠公式来表达,反而复杂的问题却可以用简单的二分法来找答案。
三、二分法求解方程
将x、y带入最后一个式子,化简得到关于t的方程:
f(t) = t * t - 2 * t * (a * cos(t) + b * sin(t)) + a * a + b * b - r * r;
导函数为
f’(t) = 2 * (t - (a + b * t) * cos(t) - (b - a * t) * sin(t));
观察函数图像和导函数图像,发现
函数的性质
1.存在多个根
2.随着x增大,y震荡递增且振幅越来越大
3.极值点均匀分布,周期为π
算法步骤
步骤一:从t = 0开始,以固定步长的形式搜索第一个极值点t1和第二个极值点t2(导函数为0作为判断条件);
步骤二:T = t2 - t1;
步骤三:第三个极值点t3 = t2 + T,第四个极值点t4 = t3 + T,以此类推,求出前k个极值点,直到连续3个极值点大于都大于0为止;
步骤四:将前k个极值点位置分成k - 1个区间,[t1,t2)、[t2,t3) …
步骤五:对每个区间应用二分查找法来找根。
二分法参考链接:https://blog.youkuaiyun.com/u013232740/article/details/43342165/
四、代码
/*************************************************************//**
* @file Bisection.cpp
* @brief this class provides interfaces for bisection algorithm
* @author Qingquan Zhou (zhouqingquan118@163.com)
* @date 10/20/2022 (M/D/Y)
*
***************************************************************/
#include <iostream>
#include <vector>
#include <cmath>
#include "Bisection.hpp"
double Bisection::function(double x, const double a, const double b