平面匹配间的变换约束

1 约定

这里约定两个同纬度向量矩阵乘之后为1\(\times\)1大小时为一个数字。

2 空间中两个平面匹配建立关系
在空间中有平面A(\(\mathbf{n}_1,d_1\))与平面B(\(\mathbf{n}_2,d_2\)),这里\(\mathbf{n}_1\)为平面A的单位法向量,\(\mathbf{n}_2\)为平面B的单位法向量。

现在已知,平面A经过(\(R,\mathbf{t}\))变换后与平面B完全重合。对于一个平面的方向向量,在经过变换后必然满足如下(1)与(2)之一的关系(可以通过点的变换推导出来)  
(1)$$
\mathbf{n}_2=R\mathbf{n}_1
$$
或者  
(2)$$
\mathbf{n}_2=-R\mathbf{n}_1$$

现在平面A上取一个点\(\mathbf{p}\),则根据A平面方程有  
(3)$$
\mathbf{n}_1^T\mathbf{p}+d_1=0
$$
同样地,\(\mathbf{p}\)经过变换后,必在平面B上,即  
(4)$$
\mathbf{n}_2^T(R\mathbf{p}+\mathbf{t})+d_2=0
$$
将(1)(3)代入(4)有  
(5)$$
\mathbf{n}_2^T\mathbf{t}=d_1-d_2
$$
将(2)(3)代入(4)有  
(6)$$
\mathbf{n}_2^T\mathbf{t}=-d_1-d_2
$$
因此,对于匹配的两个平面,其相对测量有如下关系

2.1 法向严格一致
如果要求匹配的平面,法向严格一致,则使用如下约束  
(7)$$
\begin{aligned}
\mathbf{n}_2&=R\mathbf{n}_1
\\
\mathbf{n}_2^T\mathbf{t}&=d_1-d_2
\end{aligned}$$

2.2 法向严格相反
如果要求匹配的平面,法向严格相反,则使用如下约束  
(8)$$
\begin{aligned}
\mathbf{n}_2&=-R\mathbf{n}_1
\\
\mathbf{n}_2^T\mathbf{t}&=-d_1-d_2
\end{aligned}$$

2.3 任意平面匹配
有时,我们构建了平面,提供了平面法向参数,但是我们无法确定两个平面法向的一致性(对于同一个平面,不同使用者构建的平面参数的法向量可能完全相反),这时我们只要求平面完全重合(即,在一个平面上任意取一点,经过变换后的点位,一定落在另一个平面上),这时依据此来建立约束。

这里建立约束主要使用下面(15)或(19)或(23)的约束方程来建立。

这里展开公式(4)有  
(9)$$
\begin{aligned}
0&=\mathbf{n}_2^T(R\mathbf{p}+\mathbf{t})+d_2
\\
&=(R^T\mathbf{n}_2)^T\mathbf{p}+\mathbf{n}_2^T\mathbf{t}+d_2
\end{aligned}
$$
这里设置如下中间变量  
(10)$$
\begin{aligned}
\mathbf{h}&=R^T\mathbf{n}_2
\\
g&=\mathbf{n}_2^T\mathbf{t}+d_2
\\
\mathbf{p}&= \begin{bmatrix}  
x
\\
y
\\
z
\end{bmatrix}  \\
\mathbf{n}_1&= \begin{bmatrix}  
n_{1,0}
\\
n_{1,1}
\\
n_{1,2}
\end{bmatrix} 
\end{aligned}
$$
将公式(10)带入(3)(9)有  
(11)$$
\begin{cases}
n_{1,0}x+n_{1,1}y+n_{1,2}z+d_1=0
\\
h_0x+h_1y+h_2z+g=0
\end{cases}$$

2.3.1 消去x
如果\(n_{1,0}\)为向量\(\mathbf{n}_1\)中的绝对值最大的元素,则在(11)中的约束变为  
(12)$$
\begin{cases}
x+a_xy+b_xz+c_x=0
\\
h_0x+h_1y+h_2z+g=0
\end{cases}
$$
上面公式中\(a_x\)、\(b_x\)、\(c_x\)定义如下  
(13)$$
\begin{aligned}
a_x &=\dfrac{n_{1,1}}{n_{1,0}}
\\
b_x &=\dfrac{n_{1,2}}{n_{1,0}}
\\
c_x &=\dfrac{d_1}{n_{1,0}}
\end{aligned}
$$
在(12)中消去\(x\)有  
(14)$$
(h_1-h_0a_x)y+(h_2-h_0b_x)z+(g-h_0c_x)=0
$$
由于\(y\)、\(z\)取任意值,均满足(14)算式,因此可得到如下约束方程  
(15)$$
\begin{cases}
h_1=h_0a_x
\\
h_2=h_0b_x
\\
g=h_0c_x
\end{cases}$$

2.3.2 消去y
如果\(n_{1,1}\)为向量\(\mathbf{n}_1\)中的绝对值最大的元素,则在(11)中的约束变为  
(16)$$
\begin{cases}
a_yx+y+b_yz+c_y=0
\\
h_0x+h_1y+h_2z+g=0
\end{cases}
$$
上面公式中\(a_y\)、\(b_y\)、\(c_y\)定义如下  
(17)$$
\begin{aligned}
a_y &=\dfrac{n_{1,0}}{n_{1,1}}
\\
b_y &=\dfrac{n_{1,2}}{n_{1,1}}
\\
c_y &=\dfrac{d_1}{n_{1,1}}
\end{aligned}
$$
在(12)中消去\(y\)有  
(18)$$
(h_0-h_1a_y)x+(h_2-h_1b_y)z+(g-h_1c_y)=0
$$
由于\(x\)、\(z\)取任意值,均满足(18)算式,因此可得到如下约束方程  
(19)$$
\begin{cases}
h_0=h_1a_y
\\
h_2=h_1b_y
\\
g=h_1c_y
\end{cases}$$

2.3.3 消去z
如果\(n_{1,2}\)为向量\(\mathbf{n}_1\)中的绝对值最大的元素,则在(11)中的约束变为  
(20)$$
\begin{cases}
a_zx+y+b_zy+c_z=0
\\
h_0x+h_1y+h_2z+g=0
\end{cases}
$$
上面公式中\(a_z\)、\(b_z\)、\(c_z\)定义如下  
(21)
$$
\begin{aligned}
a_z &=\dfrac{n_{1,0}}{n_{1,2}}
\\
b_z &=\dfrac{n_{1,1}}{n_{1,2}}
\\
c_z &=\dfrac{d_1}{n_{1,2}}
\end{aligned}
$$
在(12)中消去\(z\)有  
(22)
$$
(h_0-h_2a_z)x+(h_1-h_2b_z)y+(g-h_2c_z)=0
$$
由于\(x\)、\(y\)取任意值,均满足(22)算式,因此可得到如下约束方程  
(23)$$
\begin{cases}
h_0=h_2a_z
\\
h_1=h_2b_z
\\
g=h_2c_z
\end{cases}$$

PDF页面合并提取小工具

1. 背景

有时我们需要提取PDF文件指定页面,有时我们需要将一些PDF文件进行合并。一个比较常见的场景是,某天我们需要将一个扫描文件(图片形式), 将它嵌入或者追加到一个已有PDF文件里,这时就可以使用这个小工具进行处理。

 

2. 解决

除了商业软件可解决背景里提到的问题外,一些免费的在线PDF工具也提供了类似功能,如果想使用一个免费且又不想文件外流的工具,那就可以使用本页面提供的这个小工具。这个小工具只提供PDF文件的分割以及合并。如果会一点编程知识,使用Python或者其它语言库也能快速解决。这个工具没有多大技术含量,就是调用了PDFSharp提供的几个接口封装而成。

 

3. 下载

PDF合并提取小工具

C++线性规划库

1 背景

学习工作中的很多问题归结为线性规划问题。因此,对于线性规划的求解显得很有必要。

现在已经有很多开源的库或者商业软件可以求解这类问题,由于自己根据单纯形法封装了线性规划求解函数在自己的软件MathSword,为了方便C++使用,特地封装成C++可直接调用的类。

 

2 下载

封装的类在lineProgSolver.h头文件,可点击LineProg下载,对于代码里使用到的库MathLib.dll, 可到这里下载对于版本的MathSword, 解压后从解压目录Calculator下获取。

 

3 例子

以Beale问题为例

求如下目标函数的最大值

f = 0.75 * x1 – 20 * x2 + 0.5 * x3 – 6 * x4

其约束条件如下

0.25 * x1 – 8 * x2 – x3 + 9 * x4 ≤ 0

0.5 * x1 – 12 * x2 – 0.5 * x3 + 3 * x4 ≤0

x3 ≤ 1

x1>=0,x2>=0,x3>=0,x4>=0

使用如下代码求解

int main()

{

    using namespace LineProg;

    LineProgSolver solver("D:\\SoftWare\\MathSword\\Calculator\\MathLib.dll");

    SolverVar x1, x2, x3, x4;

    //1. 如果多次调用,每次调用前clear一下

    solver.clear();

    //2. 添加约束: 0.25 * x1 – 8 * x2 – x3 + 9 * x4 ≤ 0

    solver.addConstraint(0.25 * x1 – 8 * x2 – x3 + 9 * x4, EquationType_E::LessEqual_E, 0);

    //3. 添加约束:0.5 * x1 – 12 * x2 – 0.5 * x3 + 3 * x4 ≤ 0

    solver.addConstraint(0.5 * x1 – 12 * x2 – 0.5 * x3 + 3 * x4, EquationType_E::LessEqual_E, 0);

    //4. 添加约束: x3 ≤ 1

    solver.addConstraint(x3, EquationType_E::LessEqual_E, 1);

    //5. 设置最大化目标函数: 0.75 * x1 – 20 * x2 + 0.5 * x3 – 6 * x4

    solver.setMaxExpress(0.75 * x1 – 20 * x2 + 0.5 * x3 – 6 * x4);

    //6. 设置求解变量

    solver.setOptVar({ &x1, &x2, &x3, &x4 });

    //7. 求解

    if (solver.solve())

    {

        std::cout << "MaxExpress: " << solver.getMaxExpressValue() << "\n";

        std::cout << "x1 = " << x1.value << "\n";

        std::cout << "x2 = " << x2.value << "\n";

        std::cout << "x3 = " << x3.value << "\n";

        std::cout << "x4 = " << x4.value << "\n";

    }

    system("pause");

    return 0;

}

快速傅里叶变换

1 前言

傅里叶变换作为一种有效的数学分析工具,应用非常广泛。网上也有比较多的文章记录快速傅里叶变换的原理,但是它们大都基于二分思想导出对应公式。在这里,采用同样的思路,力争导出一种更具通用性的分治算法。也就是说,当数据长度为2、3、5、7等等的倍数时,也能进行快速傅里叶变换

 

2 离散傅里叶变换

对于序列\(x_0,x_1,…,x_N\),其离散傅里叶变换定义如下

(2-1)

$$
X_k = \sum_{j = 0}^{N-1}x_j \exp(-\frac{2\pi kj}{N}i)
$$

而傅里叶逆变换定义如下

(2-2)

$$
x_k =\frac{1}{N} \sum_{j = 0}^{N-1}X_j \exp(\frac{2\pi kj}{N}i)
$$

上面公式中\(k = 0,1,2,…, N – 1\).

 

3 变换的周期性

在正式导出快速傅里叶变换的公式前,从定义上来看一个性质

(3-1)

$$
\begin{aligned}
X_{k+N} &= \sum_{j = 0}^{N-1}x_j \exp(-\frac{2\pi (k+N)j}{N}i) 
\\
&= \sum_{j = 0}^{N-1}x_j (\exp(-2\pi j)\exp(-\frac{2\pi kj}{N}i))
\\
&= \sum_{j = 0}^{N-1}x_j \exp(-\frac{2\pi kj}{N}i)
\\
&= X_k
\end{aligned}
$$

同理,傅里叶逆变换可以得到类似的结果

(3-2)

$$
x_{k+N} =x_k
$$

 

4 快速傅里叶变换

这里假设\(N = nM\),那么,我们就可以将序列分成\(n\)份进行计算,每一份的数据长度为\(M\),即

(4-1)

$$
\begin{aligned}  
X_k &= \sum_{j = 0}^{N-1} x_j \exp(-\dfrac{2\pi kj}{N}i)
\\
&= \sum_{j = 0}^{M-1} x_{nj} \exp(-\dfrac{2\pi k(nj)}{N}i) + \sum_{j = 0}^{M-1} x_{nj+1} \exp(-\dfrac{2\pi k(nj+1)}{N}i) 
\\
&\ \ \ \ + … + \sum_{j = 0}^{M-1} x_{nj+n-1} \exp(-\dfrac{2\pi k(nj+n-1)}{N}i)
\\
&= \sum_{p = 0}^{n-1} \sum_{j = 0}^{M-1} x_{nj + p} \exp(-\dfrac{2\pi k(nj+p)}{N}i) 
\\
&= \sum_{p = 0}^{n-1} \sum_{j = 0}^{M-1} x_{nj + p} (\exp(-\dfrac{2\pi kp}{N}i) \exp(-\dfrac{2\pi knj}{N}i))
\\
&= \sum_{p = 0}^{n-1} ((\exp(-\dfrac{2\pi kp}{N}i)\sum_{j = 0}^{M-1} x_{nj + p}\exp(-\dfrac{2\pi knj}{N}i))
\\
&= \sum_{p = 0}^{n-1} ((\exp(-\dfrac{2\pi kp}{N}i)\sum_{j = 0}^{M-1} x_{nj + p}\exp(-\dfrac{2\pi kj}{M}i))
\end{aligned}
$$

这里令

(4-2)

$$
G_{p,M} = \sum_{j = 0}^{M-1} x_{nj + p}\exp(-\dfrac{2\pi kj}{M}i)
\\
H_{p,M} = \frac{1}{M} \sum_{j = 0}^{M-1} X_{nj + p}\exp(\dfrac{2\pi kj}{M}i)
$$

最终可以得到

(4-3)

$$
X_k = \sum_{p = 0}^{n-1} ((\exp(-\dfrac{2\pi kp}{N}i)G_{p,M})
$$

同理,对于傅里叶逆变换,可以得到如下结论

(4-4)

$$
x_k =\frac{1}{n} \sum_{p = 0}^{n-1} ((\exp(\dfrac{2\pi kp}{N}i)H_{p,M})
$$

最终,公式(4-3)、(4-4)就是快速傅里叶变换以及逆变换的表达式,这里表示将序列分成\(n\)份来求解。而公式(4-2)中的表达式,又是范围更小的傅里叶变换,因此可根据条件对其进行再分割,进而达到降低计算次数的目的。

特别地,当\(n = 2\)时,就上网上常见的快速傅里叶变换表达式。

 

5 验证

可点击这里下载上面实现的c++代码:快速傅里叶变换代码

在自己电脑上,伪造了长度\(N = 143325\)的序列数据,然后测试结果如下:

使用原始定义计算耗时:52128.6ms;使用快速傅里叶变换耗时:49.5361 ms。也就说,使用快速傅里叶变换,耗时变为原来的0.00095026722375倍。当然,如果\(N\)值更大,效果越明显。

 

6 二维FFT

对于二维FFT,可以直接对二维数据每一行单独做FFT,然后对新的二维数据的每一列再单独FFT。而对于FFT逆变换,则按相反操作即可。

MathSword教程1.基础操作

1. 程序主界面主要分【矩阵运算】与【算式解析】

2. 【矩阵运算】主要承担绝大部分的操作:实数矩阵变量的创建,优化问题求解,绘图等一系列操作。

3. 【算式解析】主要针对复数表达式进行计算,一些不常使用的特殊函数,也主要内置到这个模块进行计算。

4. 程序主要操作是遵循,变量赋值,然后函数调用,类似这么一种操作,所有的函数调用都有使用说明,具体可右键查看命令帮助。程序的

一些基本的操作,可看如下视频

 

MathSword教程3.非线性拟合

 

PS:推荐优先使用【启发式优化算法】工具求解,结果不理想再按如下方式求解!

1. 虽然程序提供多个函数进行非线性拟合,但是推荐使用程序界面上方的【高级优化】进行非线性拟合、非线性方程组求解、以及非线性规划问题。

2. 【高级优化】窗口里可自定义一些代码设置,可以指定优化变量取值范围,以及设置混合整数优化问题。

3. 非线性拟合整体步骤如下:

setp1. 将已知数据以矩阵变量形式存储(导入)在一个变量当中。

setp2. 打开【高级优化】窗口,在【模板】当中选择【函数拟合】

setp3. 在Code下误差函数ErrorFun里定义误差(其实就是把函数表达式写一遍)

setp4. 在Parameter下设置data关键字(这个就是step1里的变量)

setp5. 在Parameter下设置varCount关键字(即优化变量个数)

setp6. 在Parameter下设置errorCount关键字(即误差函数个数,这个根据step3当中定义的误差而定,一般就1个)

setp7. 然后点击【求解】即可

step8. 计算完毕后,回到程序命令行主窗口,输入Parameter下设置varName关键字后的名称,查看结果。

PS:当中还有很多参数,可具体看参数设置说明,不同参数,优化结果以及优化效率会有一定差异。

 

MathSword教程4.非线性方程组求解

1 简单版本

如果想快速定义非线性方程问题求解,且书写语言更加接近数学表达式形式,可参考线性/非线性规划教程

2 进阶版本

这里需要一定的编程基础,可灵活构建所需的数学问题。 优先推荐使用【启发式优化算法】工具求解,结果不理想再按如下方式求解!

非线性方程组求解步骤如下:

setp1. 打开【高级优化】窗口,在【模板】当中选择【方程求解】

setp2. 在Code下误差函数ErrorFun里定义误差(其实就是把方程组写入到里面然后构建误差)

setp3. 在Parameter下设置varCount关键字(即优化变量个数)

setp4. 在Parameter下设置errorCount关键字(即误差函数个数,这个一般就是方程个数)

setp5. 然后点击【求解】即可

step6. 计算完毕后,回到程序命令行主窗口,输入Parameter下设置varName关键字后的名称,查看结果。(结果的顺序对应在step2当中变量定义的顺序)

 

 

 

 

 

 

MathSword教程5.数学公式搜索

1. 前言

工作中,我们可能会涉及到一些数据处理的问题。在数据处理过程中,很多时候,虽然我们有一堆数据,但是我们并不清楚这些数据之间的关系。因为某些原因,我们想观察这些数据之间的关系。这时,就希望能找到某种数学表达式来把这些数据变量给关联起来。

如果对参数个数没有特别要求的话,神经网络可以作为参考的数学模型一种选择。如果对参数个数有限制要求,则就需要另一种方式来查找对应表达式了,也就是我们希望找一种表达更简洁,且变量尽量少的数学表达式来构建数据之间的关系。常规地,我们希望找到一种如下形式的公式

(1)

$$
y = F(x_1,x_2,…, x_n, a_1, a_2, …, a_m)
$$

公式(1)中的\(x_1,x_2,…,x_n\)为自变量,\(y\)为因变量,\(a_1,a_2,…,a_m\)为公式中的额外变量参数,\(F\)为我们需要查找的数学模型。在这里,我们已知的是数据集\(x_1,x_2,…,x_n, y\),而未知的是\(a_1,a_2,…,a_m,F\)

在实际情况中,如果我们通过观察\(x_1,x_2,…,x_n, y\)的规律,或者变量之间的物理意义,如果能定向构建出\(F\),再进行优化求解,那这时就会省很多事情。这里,假如,我们并不能观察出\(x_1,x_2,…,x_n, y\)的规律,这时,构建\(F\)就无外乎进行盲试,MathSword提供一些盲试的数学模型。

 

2. 应用

 

 

 

MathSword教程6.神经网络

1. 前言

有些数学模型难以查找时,我们会尝试使用神经网络来进行逼近。MathSword提供一种简单快速的神经网络搭建模式,在下面的视频例子中,可直接使用对应代码更改相应参数即达到不同模型建立的目的。使用时,注意如下参数的设置:

NetCreateConnectLayer :使用本函数快速搭建一个前馈网络

NetCreateSolver :迭代不收敛时,尝试重新用本函数设置求解算法以及学习速率

NetOptSetOptVarValue :迭代不收敛时,尝试用本函数给优化变量重新设置迭代初值

PS:迭代不收敛时,确定输入输出样本的值域是否合理(比如选择Tanh作为输出层函数, 这时对应的样本值域就应该在-1到1之间),不合理时,需要先对数据进行映射(比如映射到-1与1之间等等)!

2. 应用

2.1 代码快速构建(不推荐使用)

下载: code

下载: Mnist

 

 

 

 

 

 

2.2 图构建方式(推荐使用)