Levenberg-Marquardt算法与透视变换矩阵优化

2024-01-05 13:40
文章标签 算法 优化 矩阵 透视 变换 levenberg marquardt

本文主要是介绍Levenberg-Marquardt算法与透视变换矩阵优化，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

一般来说我们利用牛顿法使用来求f(x)=0的解。求解方法如下： 
先对f(x)一阶泰勒展开得 
                                \LARGE f(x+\triangle ) = f(x) +f'(x)\triangle

所以我们有
                    \LARGE \triangle = x -x_{0} = -\frac{f(x_{0})}{f'(x_{0}}，即\LARGE \LARGE x = x_{0}-\frac{f(x_{0})}{f'(x_{0})}

因此也就得到了我们的牛顿迭代公式： 
                               \LARGE x_{n} = x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})}
求解最优化问题
                                \LARGE min f(x)

牛顿法首先则是将问题转化为求
                                \LARGE f'(x) = 0
这个方程的根。 
        一阶展开：
                                \LARGE f'(x) \approx f'(x_{_{0}}) +(x - x_{_{0}})f''(x_{_{0}})                    

令
                               \LARGE f'(x_{_{0}}) +(x - x_{_{0}})f''(x_{_{0}}) = 0

         求解得到\LARGE x，相比于\LARGE x_{0}，\LARGE f'(x) < f'(x_{0})
 

高斯牛顿法

在讲牛顿法的时候，我们举的例子\LARGE x是一维的，若如果我们遇到多维的x该如何办呢？这时我们就可以利用雅克比，海赛矩阵之类的来表示高维求导函数了。 
比如
                             \LARGE f(X) = 0有,其中\LARGE X = [x_{0},x_{1},...,x_{n}]
 

所以我们有雅克比矩阵： 
                                      \LARGE J_{f} = \begin{bmatrix} \frac{\partial f_{1} }{\partial x_{0}} &... &\frac{\partial f_{1} }{\partial x_{n}}\\ &. .. & \\ \frac{\partial f_{n} }{\partial x_{0}}& ... & \frac{\partial f_{n} }{\partial x_{n}} \end{bmatrix}

有海赛矩阵： 
                                     \LARGE H_{f} = \begin{bmatrix} \frac{\partial^2 f_{1} }{\partial x_{0}^2}& \frac{\partial^2 f_{1} }{\partial x_{0}\partial x_{1}}&...&\frac{\partial^2 f_{1} }{\partial x_{0}\partial x_{n}}\\ ...& ... & ...&...\\ \frac{\partial^2 f_{n} }{\partial x_{0}^2}& \frac{\partial^2 f_{n} }{\partial x_{0}\partial x_{1}}&...&\frac{\partial^2 f_{n} }{\partial x_{0}\partial x_{n}} \end{bmatrix}
所以高维牛顿法解最优化问题又可写成： 
                                   \LARGE X_{n} = X_{n-1} - H_{f}\triangle
 

梯度 代替了低维情况中的一阶导 
Hessian矩阵代替了二阶导 
求逆代替了除法 
例：不妨设目标函数为： 
                               \LARGE s(x) = \sum_{i=0}^{n} f^{2}(x_{i})

所以梯度向量在方向上的分量： 
                              \LARGE g_{j} = 2\sum_{i=0}^{n} f_{i}\frac{\partial f_{i}}{\partial x_{j}}

Hessian 矩阵的元素则直接在梯度向量的基础上求导： 
                             \LARGE H_{jk} = 2\sum_{i=0}^{n} (\frac{\partial f_{i}}{\partial x_{j}}\frac{\partial f_{i}}{\partial x_{k}} + \frac{\partial^2 f_{i}}{\partial x_{j}\partial x_{k}})
 

高斯牛顿法的一个小技巧是，将二次偏导省略，于是： 
                             \LARGE H_{jk} \approx 2\sum_{i=0}^{n} \frac{\partial f_{i}}{\partial x_{j}}\frac{\partial f_{i}}{\partial x_{k}} = 2\sum_{i=0}^{n} J_{ij}J_{ik}

其中\LARGE J_{ij}J_{ik}为雅克比矩阵中的第i行j列元素 
改写成 矩阵相乘形式： 
                                 \LARGE H\approx 2J_{f}^{T}J_{f}  

代入牛顿法高维迭代方程的基本形式，得到高斯牛顿法迭代方程： 
                           \LARGE X_{n+1} = X_{n} - \triangle,其中\LARGE \triangle = -(J_{f}^{T}J_{f})^{-1}J_{f}^{T}f

Levenberg-Marquardt算法

引用维基百科的一句话就是：
莱文贝格－马夸特方法（Levenberg–Marquardt algorithm）能提供数非线性最小化（局部最小）的数值解。此算法能借由执行时修改参数达到结合高斯-牛顿算法以及梯度下降法的优点，并对两者之不足作改善（比如高斯-牛顿算法之反矩阵不存在或是初始值离局部极小值太远）

在我看来，就是在高斯牛顿基础上修改了一点。 
在高斯牛顿迭代法中，我们已经知道 
                            \LARGE \triangle = -(J_{f}^{T}J_{f})^{-1}J_{f}^{}

在莱文贝格－马夸特方法算法中则是 
                \LARGE \triangle = -(J_{f}^{T}J_{f} + \lambda I)^{-1}J            

在我看来好像就这点区别。至少我看的维基百科是这样的。 
然后Levenberg-Marquardt方法的好处就是在于可以调节: 
如果下降太快，使用较小的λ，使之更接近高斯牛顿法 
如果下降太慢，使用较大的λ，使之更接近梯度下降法

 

H矩阵本身已经最优，手动改变了一点测试LM。

close all;
clear all;
clc;% 返回值 H 是一个3*3的矩阵
% pts1 和 pts2是2*n的坐标矩阵对应特征点的(x,y)坐标
CordReal = load('TestReal.txt');
CordImg = load('TestImg.txt');
pts1 = CordReal(:,1:2)';
pts2 = CordImg(:,1:2)';
n = size(pts1,2);
A = zeros(2*n,9);
A(1:2:2*n,1:2) = pts1';
A(1:2:2*n,3) = 1;
A(2:2:2*n,4:5) = pts1';
A(2:2:2*n,6) = 1;
x1 = pts1(1,:)';
y1 = pts1(2,:)';
x2 = pts2(1,:)';
y2 = pts2(2,:)';
A(1:2:2*n,7) = -x2.*x1;
A(2:2:2*n,7) = -y2.*x1;
A(1:2:2*n,8) = -x2.*y1;
A(2:2:2*n,8) = -y2.*y1;
A(1:2:2*n,9) = -x2;
A(2:2:2*n,9) = -y2;[evec,D] = eig(A'*A);
H = reshape(evec(:,1),[3,3])';
H = H/H(end); % make H(3,3) = 1TestReal = load('TestReal.txt');
TestImg = load('TestImg.txt');
N = size(TestReal,1);
TestRealQici = [TestReal(:,1:2),ones(N,1)];
TestImgQici = [TestImg(:,1:2),ones(N,1)];
Result = zeros(N,3);
ImgInvH = zeros(N,3);
for i=1:NResult(i,1)=((H(1,1)+0.01)*x1(i)+H(1,2)*y1(i)+H(1,3))./(H(3,1)*x1(i)+H(3,2)*y1(i)+H(3,3));Result(i,2)=(H(2,1)*x1(i)+H(2,2)*y1(i)+H(2,3))./(H(3,1)*x1(i)+H(3,2)*y1(i)+H(3,3));%     Result(i,:)=H*TestRealQici(i,:)';
%     ImgInvH(i,:)=H\TestImgQici(i,:)';
endD = TestImg(:,1:2)-Result(:,1:2);
xlswrite('Result.xls',[TestImgQici Result])%%
width = round(max(TestReal(:,1)))+50;
height = round(max(TestReal(:,2)))+50;
RealGen = 255*ones(height,width);
imshow(RealGen);
hold on;
plot(TestReal(:,1),TestReal(:,2),'k+');%黑十字
hold on;
plot(ImgInvH(:,1),ImgInvH(:,2),'r.');%红点%___________
figure();
width = round(max(TestImg(:,1)))+30;
height = round(max(TestImg(:,2)))+30;
ImgGen = 255*ones(height,width);
imshow(ImgGen);
hold on;
plot(TestImg(:,1),TestImg(:,2),'k+');%黑十字
%___________
hold on;
plot(Result(:,1),Result(:,2),'R.');%红点%%
%H的LM优化
% % 计算函数f的雅克比矩阵，是解析式
% syms h1 h2 h3 h4 h5 h6 h7 h8 h9 xi yi xo yo real;
% % f=(xout-(h1*xin+h2*yin+h3)./(h7*xin+h8*yin+h9)).^2+(yout-(h4*xin+h5*yin+h6)./(h7*xin+h8*yin+h9)).^2; %自己定义误差函数为各点的误差和
% f1=(h1*xi+h2*yi+h3)./(h7*yi+h8*yi+h9);
% f2=(h4*xi+h5*yi+h6)./(h7*yi+h8*yi+h9);
% Jsym = jacobian([f1;f2],[h1 h2 h3 h4 h5 h6 h7 h8 h9])% 拟合用数据为x1,y1,x2,y2。% 2. LM算法
% 以近似解H作为迭代初始值
h1_0 = H(1,1);
h2_0 = H(1,2);
h3_0 = H(1,3);
h4_0 = H(2,1);
h5_0 = H(2,2);
h6_0 = H(2,3);
h7_0 = H(3,1);
h8_0 = H(3,2);
h9_0 = H(3,3);% 数据个数
Ndata=length(x1);
% 参数维数
Nparams=9;
% 迭代最大次数
n_iters=50;
% LM算法的阻尼系数初值
lamda=0.01;
lamda1=0.0001;% step1: 变量赋值
updateJ=1;
h1_est = h1_0+0.01;
h2_est = h2_0;
h3_est = h3_0;
h4_est = h4_0;
h5_est = h5_0;
h6_est = h6_0;
h7_est = h7_0;
h8_est = h8_0;
h9_est = h9_0;
x_Cord_est = x2;
y_Cord_est = y2;% step2: 迭代
% step2: 迭代
it = 1;     % 迭代计数
id = 1;     % 有效结果计数
repeat = 1; % 若误差未减小的情况连续出现次数
while repeat < 10   % 连续出现10次if updateJ==1% 根据当前估计值，计算雅克比矩阵J=zeros(2,Nparams);Jit = zeros(2*length(x1),Nparams);DeltaF = zeros(Nparams,2 );Fh0 =0;for i=1:length(x1)Jit(i,:)=[x1(i)./(h7_est*x1(i)+h8_est*y1(i)+h9_est),y1(i)./(h7_est*x1(i)+h8_est*y1(i)+h9_est),1./(h7_est*x1(i)+h8_est*y1(i)+h9_est),0,0,0,(-x1(i)*(h1_est*x1(i)+h2_est*y1(i)+h3_est))./((h7_est*x1(i)+h8_est*y1(i)+h9_est)^2),(-y1(i)*(h1_est*x1(i)+h2_est*y1(i)+h3_est))./((h7_est*x1(i)+h8_est*y1(i)+h9_est)^2),(-(h1_est*x1(i)+h2_est*y1(i)+h3_est))./((h7_est*x1(i)+h8_est*y1(i)+h9_est)^2)];Jit(i+1,:)=[0,0,0,x1(i)./(h7_est*x1(i)+h8_est*y1(i)+h9_est),y1(i)./(h7_est*x1(i)+h8_est*y1(i)+h9_est),1./(h7_est*x1(i)+h8_est*y1(i)+h9_est),-(x1(i)*(h4_est*x1(i)+h5_est*y1(i)+h6_est))./((h7_est*x1(i)+h8_est*y1(i)+h9_est)^2),(-y1(i)*(h4_est*x1(i)+h5_est*y1(i)+h6_est))./((h7_est*x1(i)+h8_est*y1(i)+h9_est)^2),(-(h4_est*x1(i)+h5_est*y1(i)+h6_est))./((h7_est*x1(i)+h8_est*y1(i)+h9_est)^2)];end f1_est = (h1_est*x1+h2_est*y1+h3_est)./(h7_est*x1+h8_est*y1+h9_est);f2_est = (h4_est*x1+h5_est*y1+h6_est)./(h7_est*x1+h8_est*y1+h9_est);dx = x_Cord_est- f1_est;dy = y_Cord_est- f2_est;d = cat(1,dx,dy);Hess=Jit'*Jit;if repeat==1e=dot(d,d);end            end H_lm=Hess+(lamda*eye(Nparams,Nparams));    % 计算步长dp，并根据步长计算新的可能的参数估计值if(rank(H_lm) == 9)break;enddp=(H_lm)\(Jit'*d(:));   h1_lm=h1_est+dp(1);h2_lm=h2_est+dp(2);h3_lm=h3_est+dp(3);h4_lm=h4_est+dp(4);h5_lm=h5_est+dp(5);h6_lm=h6_est+dp(6);h7_lm=h7_est+dp(7);h8_lm=h8_est+dp(8);h9_lm=h9_est+dp(9);f1_est_lm = (h1_lm*x1+h2_lm*y1+h3_lm)./(h7_lm*x1+h8_lm*y1+h9_lm);f2_est_lm = (h4_lm*x1+h5_lm*y1+h6_lm)./(h7_lm*x1+h8_lm*y1+h9_lm);dx_lm = x_Cord_est- f1_est_lm ;dy_lm = y_Cord_est- f2_est_lm;d_lm = cat(1,dx_lm ,dy_lm );e_lm=dot(d_lm,d_lm);if e_lm<e
%         if e_lm<10
%             break;
%         elselamda=lamda1/10;h1_est=h1_lm;h2_est=h2_lm;h3_est=h3_lm;h4_est=h4_lm;h5_est=h5_lm;h6_est=h6_lm;h7_est=h7_lm;h8_est=h8_lm;h9_est=h9_lm;e=e_lm;disp(e);updateJ=1;repeat=repeat+1;
%          endelseupdateJ=0;lamda=lamda*5;end
end
%显示优化的结果
H(1,1)=h1_est;
H(1,2)=h2_est;
H(1,3)=h3_est;
H(2,1)=h4_est;
H(2,2)=h5_est;
H(2,3)=h6_est;
H(3,1)=h7_est;
H(3,2)=h8_est;
H(3,3)=h9_est;TestReal = load('TestReal.txt');
TestImg = load('TestImg.txt');
N = size(TestReal,1);
TestRealQici = [TestReal(:,1:2),ones(N,1)];
TestImgQici = [TestImg(:,1:2),ones(N,1)];
Result = zeros(N,3);
ImgInvH = zeros(N,3);
for i=1:NResult(i,1)=(H(1,1)*x1(i)+H(1,2)*y1(i)+H(1,3))./(H(3,1)*x1(i)+H(3,2)*y1(i)+H(3,3));Result(i,2)=(H(2,1)*x1(i)+H(2,2)*y1(i)+H(2,3))./(H(3,1)*x1(i)+H(3,2)*y1(i)+H(3,3));%Result(i,:)=H\TestRealQici(i,:)';%ImgInvH(i,:)=H*TestImgQici(i,:)';
endxlswrite('Result.xls',[TestImgQici Result])%%
width = round(max(TestReal(:,1)))+50;
height = round(max(TestReal(:,2)))+50;
RealGen = 255*ones(height,width);
imshow(RealGen);
hold on;
plot(TestReal(:,1),TestReal(:,2),'k+');%黑十字
hold on;
plot(ImgInvH(:,1),ImgInvH(:,2),'r.');%红点%___________
figure();
width = round(max(TestImg(:,1)))+30;
height = round(max(TestImg(:,2)))+30;
ImgGen = 255*ones(height,width);
imshow(ImgGen);
hold on;
plot(TestImg(:,1),TestImg(:,2),'k+');%黑十字
%___________
hold on;
plot(Result(:,1),Result(:,2),'R.');%红点


 

优化后结果

 

这篇关于Levenberg-Marquardt算法与透视变换矩阵优化的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/573043。

相关文章

正则表达式高级应用与性能优化记录

正则表达式高级应用与性能优化记录

《正则表达式高级应用与性能优化记录》本文介绍了正则表达式的高级应用和性能优化技巧,包括文本拆分、合并、XML/HTML解析、数据分析、以及性能优化方法,通过这些技巧,可以更高效地利用正则表达式进行复杂... 目录第6章：正则表达式的高级应用6.1 模式匹配与文本处理6.1.1 文本拆分6.1.2 文本合并6
阅读更多...
Vue3 的 shallowRef 和 shallowReactive：优化性能

Vue3 的 shallowRef 和 shallowReactive：优化性能

大家对 Vue3 的 ref 和 reactive 都很熟悉，那么对 shallowRef 和 shallowReactive 是否了解呢？ 在编程和数据结构中，“shallow”（浅层）通常指对数据结构的最外层进行操作，而不递归地处理其内部或嵌套的数据。这种处理方式关注的是数据结构的第一层属性或元素，而忽略更深层次的嵌套内容。 1. 浅层与深层的对比 1.1 浅层（Shallow） 定义
阅读更多...
不懂推荐算法也能设计推荐系统

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例，告诉我们不懂推荐算法的产品，也能从产品侧出发， 设计出一款不错的推荐系统。 相信很多新手产品，看到算法二字，多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法（注：传统是指科班出身的产品都会接触过）。但对于推荐算法，多数产品对着网上搜到的资源，都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后，更是加大了理解的难度。 但，不了解推荐算法，就无法做推荐系
阅读更多...
HDFS—存储优化（纠删码）

HDFS—存储优化（纠删码）

纠删码原理 HDFS 默认情况下，一个文件有3个副本，这样提高了数据的可靠性，但也带来了2倍的冗余开销。 Hadoop3.x 引入了纠删码，采用计算的方式，可以节省约50％左右的存储空间。 此种方式节约了空间，但是会增加 cpu 的计算。 纠删码策略是给具体一个路径设置。所有往此路径下存储的文件，都会执行此策略。 默认只开启对 RS-6-3-1024k
阅读更多...
康拓展开（hash算法中会用到）

康拓展开（hash算法中会用到）

康拓展开是一个全排列到一个自然数的双射（也就是某个全排列与某个自然数一一对应） 公式： X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。（a[i]在不同应用中的含义不同）； 典型应用： 计算当前排列在所有由小到大全排列中的顺序，也就是说求当前排列是第
阅读更多...
使用opencv优化图片（画面变清晰）

使用opencv优化图片（画面变清晰）

文章目录 需求影响照片清晰度的因素 实现降噪测试代码 锐化空间锐化Unsharp Masking频率域锐化对比测试 对比度增强常用算法对比测试 需求 对图像进行优化，使其看起来更清晰，同时保持尺寸不变，通常涉及到图像处理技术如锐化、降噪、对比度增强等 影响照片清晰度的因素 影响照片清晰度的因素有很多，主要可以从以下几个方面来分析 1. 拍摄设备 相机传感器：相机传
阅读更多...
csu 1446 Problem J Modified LCS （扩展欧几里得算法的简单应用）

csu 1446 Problem J Modified LCS （扩展欧几里得算法的简单应用）

这是一道扩展欧几里得算法的简单应用题，这题是在湖南多校训练赛中队友ac的一道题，在比赛之后请教了队友，然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意：把题意转变下就变成了：求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点：求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个
阅读更多...
综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动，这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升：减少抖动，提高图像的清晰度和细节表现力，使得监控画面更加真实可信。 - 细节增强：在低光条件下，抖
阅读更多...
【数据结构】——原来排序算法搞懂这些就行，轻松拿捏

【数据结构】——原来排序算法搞懂这些就行，轻松拿捏

前言：快速排序的实现最重要的是找基准值，下面让我们来了解如何实现找基准值 基准值的注释：在快排的过程中，每一次我们要取一个元素作为枢纽值，以这个数字来将序列划分为两部分。 在此我们采用三数取中法，也就是取左端、中间、右端三个数，然后进行排序，将中间数作为枢纽值。 快速排序实现主框架： //快速排序 void QuickSort(int* arr, int left, int rig
阅读更多...
poj 3974 and hdu 3068 最长回文串的O（n）解法(Manacher算法)

poj 3974 and hdu 3068 最长回文串的O（n）解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大，用原来的O（n^2）会爆。 小白上的O(n^2)解法代码：TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2