船载视频稳定和校正的地平线跟踪方法研究（Matlab代码实现）

  💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

 ⛳️赠与读者

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、数据、文章下载

---

 ⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

     或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

参考文献：

摘要：提出了一种用于稳定和校正浮动平台数字视频的算法。该方法依赖于一种水平跟踪技术，该技术在北太平洋研究航行期间，在各种照明和海况条件下，对12天内的48小时视频数据进行了测试。在这个数据集中，地平线在92%的帧中得到了正确的标记。理想化的摄像机模型假设纯粹的俯仰和横滚运动，平坦的海面和畅通无阻的地平线。俯仰和横滚是沿着摄像机观察方向定义的，而不是在传统的船舶坐标系中定义的，因此该方法可用于相对于船舶的任何航向。俯仰和横滚的不确定性是根据地平线发现方法的不确定性估计的。这些误差在横摇方向上为0.68，在纵摇方向上为0.38。校正误差主要由摄像机高度的不确定性引起，这可能会随着浮动平台的垂荡运动而变化。这些误差的传播被证明为破碎波分布L(c)。
关键词：吸积 吸积盘 恒星：个体：SS 433 V1343 天鹰座超巨星 技术：光谱学 X射线：双星

1)使用霍夫变换检测地平线(为此需要图像处理或计算机视觉工具箱)。
2) 根据地平线方向计算精确的相机俯仰和横滚。
3)稳定从不同相机角度拍摄的多幅图像。
4)使用已知的相机高度将图像校正为平坦的海面。
5)如果已知相机方向的不确定性，则估计稳定或校正中的误差。

📚2 运行结果

文献：

复现结果：

主函数代码：

clc
clearvars
close all

% Set camera pose and parameters
H = 10.7;
inc = 75*pi/180;
roll = 0*pi/180;
azi = 0*pi/180;
load('intrinsicMatrix.mat')

% Establish domain
tMin = 0;  % sec
tMax = 60*60*1;  % sec, use long time for smooth Lambda Results
T = tMax-tMin;
xMin = -10;  % m
xMax = 10;  % m
yMin = 30;  % m
yMax = 50;  % m
A = (yMax-yMin)*(xMax-xMin);  % m^2
 
% Set breaker characteristics
cAvg = 3;   % m/s
cStd = 1;   % m/s
LAvg = 3;   % breaker crest length, m (constant) 
TAvg = 2;   % breaking duration, sec (constant) 
breakerDirX = true;  % direction of breaking: x = true, y = false

% Set Lambda probability density function
dc = 0.25;  % m/s
cVec = -2:dc:8;
lambdaPdf = normpdf(cVec,cAvg,cStd);

% Determine number of breakers needed for chosen breaking rate (or just set
% numBreakers directly)
brkRate = 120/3600; % Hz
numBreakers = floor(brkRate*A*T*dc/(TAvg*LAvg*sum(lambdaPdf.*cVec*dc)));
lTotal = TAvg*LAvg*numBreakers*sum(lambdaPdf*dc)/(A*T);
LambdaInput = lTotal*lambdaPdf;

% Randomly sample position, time, and speed of each breaker
xCenter = xMin + (xMax-xMin)*rand(numBreakers,1);
yCenter = yMin + (yMax-yMin)*rand(numBreakers,1);
tBegin = tMin + T*rand(numBreakers,1);
cRand = cAvg + cStd*randn(numBreakers,1);

% Simulate breaker moving forward in time
dt = 1/15;  % "Frame Rate" of observations
distTraveled = cRand*(0:dt:TAvg);   % vector of distance traveled by breaker
numTimeSteps = TAvg/dt + 1;

% Calculate positions of breaker end points in time
if breakerDirX
    x1 = repmat(xCenter,[1,numTimeSteps])+distTraveled;
    x2 = repmat(xCenter,[1,numTimeSteps])+distTraveled;
    y1 = repmat(yCenter-LAvg/2,[1,numTimeSteps]);
    y2 = repmat(yCenter+LAvg/2,[1,numTimeSteps]);
else
    x1 = repmat(xCenter-LAvg/2,[1,numTimeSteps]);
    x2 = repmat(xCenter+LAvg/2,[1,numTimeSteps]);
    y1 = repmat(yCenter,[1,numTimeSteps])+distTraveled;
    y2 = repmat(yCenter,[1,numTimeSteps])+distTraveled;
end

% Set camera motion frequency and amplitude, then find height for each
% time during breaker duration
timeElapsed = repmat(0:dt:TAvg,[numBreakers,1]);
tMatrix = repmat(tBegin,[1,numTimeSteps])+timeElapsed;
amp = 1;    % m (remember wave height is twice wave amplitude)
freq = 1/8;  % Hz
radFreq = 2*pi*freq;
HMatrix = amp*cos(tMatrix*radFreq);

% Determine falsely observed position of breaker endpoints due to camera
% motion
x1Obs = zeros(numBreakers,numTimeSteps);
y1Obs = zeros(numBreakers,numTimeSteps);
x2Obs = zeros(numBreakers,numTimeSteps);
y2Obs = zeros(numBreakers,numTimeSteps);
for i = 1:numBreakers
    for j= 1:numTimeSteps
        [dx1,dy1,~,~] = ProjectionErrorsXY(x1(i,j),y1(i,j),H,inc,roll,azi,K,HMatrix(i,j),0,0,0);
        x1Obs(i,j) = x1(i,j)+dx1;
        y1Obs(i,j) = y1(i,j)+dy1;
        [dx2,dy2,~,~] = ProjectionErrorsXY(x2(i,j),y2(i,j),H,inc,roll,azi,K,HMatrix(i,j),0,0,0);
        x2Obs(i,j) = x2(i,j)+dx2;
        y2Obs(i,j) = y2(i,j)+dy2;
    end
end

% calculate true and observed breaker speeds and lengths, based on x and y
% coordinates of breaker endpoints
if breakerDirX
    cMatrixTrue = (diff(x1,1,2)/dt+diff(x2,1,2)/dt)/2;
    LMatrixTrue = y2-y1;
    cMatrixObs = (diff(x1Obs,1,2)/dt+diff(x2Obs,1,2)/dt)/2;
    LMatrixObs = y2Obs-y1Obs;
else
    cMatrixTrue = (diff(y1,1,2)/dt+diff(y2,1,2)/dt)/2;
    LMatrixTrue = x2-x1;
    cMatrixObs = (diff(y1Obs,1,2)/dt+diff(y2Obs,1,2)/dt)/2;
    LMatrixObs = x2Obs-x1Obs;
end
LMatrixTrue = LMatrixTrue(:,1:(numTimeSteps-1));
LMatrixObs = LMatrixObs(:,1:(numTimeSteps-1));

% calculate true and observed lambda
LambdaTrue = nan(length(cVec),1);
LambdaObs = nan(length(cVec),1);
for i=1:length(cVec)
    ind = cMatrixTrue>(cVec(i)-dc/2) & cMatrixTrue<(cVec(i)+dc/2);
    LambdaTrue(i) = sum(LMatrixTrue(ind))/(A*dc*T/dt);
    indObs = cMatrixObs>(cVec(i)-dc/2) & cMatrixObs<(cVec(i)+dc/2);
    LambdaObs(i) = sum(LMatrixObs(indObs))/(A*dc*T/dt);
end

% show 10 example crests
figure(1)
hold on
set(gca,'xlim',[xMin xMax],'ylim',[yMin yMax])
for i=1:10
    for j=1:numTimeSteps
        plot([x1(i,j),x2(i,j)],[y1(i,j),y2(i,j)],'-r')
        plot([x1Obs(i,j),x2Obs(i,j)],[y1Obs(i,j),y2Obs(i,j)],'-b')
        pause(.01)
    end
end

% Plot input, true, and observed Lambdas
figure(2)
plot(cVec,LambdaTrue,'-r')
hold on
plot(cVec,LambdaInput,'-k')
plot(cVec,LambdaObs,'-b')
hold off

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

[1] Schwendeman M , Thomson J .A Horizon-Tracking Method for Shipboard Video Stabilization and Rectification*[J].Journal of Atmospheric & Oceanic Technology, 2015, 32(1):164-176.DOI:10.1175/JTECH-D-14-00047.1.

🌈4 Matlab代码、数据、文章下载

资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取