當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

HS光流法及其代码示例

發布時間：2023/12/15 编程问答 32 豆豆

生活随笔收集整理的這篇文章主要介紹了 HS光流法及其代码示例小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

前言

光流是指空間運動物體在觀測成像面上的像素運動的瞬時速度，它利用圖像序列像素強度的時域變化和相關性來確定各自像素位置的“運動”，即反映圖像灰度在時間上的變化與景物中物體結構與其運動的關系。將二維圖像平面特定坐標上的灰度瞬時變化率定義為光流矢量。

光流的物理意義

當人的眼睛觀察運動物體時，物體的景象在人眼的視網膜上形成一系列連續變化的圖像，這一系列連續變化的信息不斷“流過”視網膜（即圖像平面），好像一種光的“流”，故稱之為光流。光流表達了圖像的變化，由于它包含了目標運動的信息，因此可被觀察者用來確定目標的運動情況。

光流場

在空間中，運動可以用運動場描述，而在一個圖像平面上，物體的運動往往是通過圖像序列中不同圖像灰度分布的不同體現的，從而，空間中的運動場轉移到圖像上就表示為光流場（optical flow field）。
　　光流場是一個二維矢量場，它反映了圖像上每一點灰度的變化趨勢，可看成是帶有灰度的像素點在圖像平面上運動而產生的瞬時速度場。它包含的信息即是各像點的瞬時運動速度矢量信息。
　　研究光流場的目的就是為了從序列圖像中近似計算不能直接得到的運動場。光流場在理想情況下，光流場對應于運動場。

三維空間的矢量場及其在二維平面內的投影

光流場計算的基本原理

設在t 時刻，像素點 (x,y) 處的灰度值為I(x,y,t)；在(t+ $Δ\Delta$ t)，該點運動到新的位置，它在圖像上的位置變為(x+ $Δ\Delta$ x, y+ $Δ\Delta$ y)，灰度值記為I(x+ $Δ\Delta$ x, y+ $Δ\Delta$ y)。根據圖像一致性假設，即圖像沿著運動軌跡的亮度保持不變，滿足 $dI(x,y,t)dt\frac{dI(x,y,t)}{dt}$ =0，則： $(1)I(x,y,t)=I(x+Δx,y+Δy+t+Δt)I(x,y,t)=I(x+\Delta x, y+\Delta y + t+\Delta t)\tag {1}$ 　　設u和 $υ\upsilon$ 分別為該點的光流矢量沿x和y方向的兩個分量，且 u= $dxdt\frac{dx}{dt}$ ， $υ\upsilon$ = $dydt\frac{dy}{dt}$ ，將式(1)等號右邊用泰勒公式展開看，得到： $(2)I(x+Δx,y+Δy,t+Δt)=I(x,y,t)+?I?xΔx+?I?yΔy+?I?tΔt+εI(x+\Delta x, y+\Delta y, t+ \Delta t) = I(x,y,t)+\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t+\varepsilon \tag {2}$ 　　忽略二階以上的高次項，則有： $(3)?I?xΔx+?I?yΔy+?I?tΔt=0\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t=0 \tag{3}$ 　　由于 $Δt→0\Delta t\to 0$ ，于是有： $(4)?I?xdxdt+?I?ydydt+?I?t=0\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}+\frac{\partial I}{\partial t} = 0\tag{4}$ 　　也即： $(5)Ixu+Iyυ+It=0I_xu+I_y\upsilon+I_t= 0\tag{5}$ 　　式（5）是光流基本等式。設 $I_x、I_y和I_t$ 分別為參考點像素的灰度值沿x,y,t這三個方向的偏導數。
　　Horn-Schunck算法引用的附加約束條件的基本思想是：在求解光流時，要求光流本身盡可能地平滑，即引入對光流的整體平滑性約束求解光流方程病態問題。所謂平滑，就是在給定的領域內 $(?2u+?2υ)(\nabla^2u+\nabla^2\upsilon)$ 應盡量地小，這就是求條件極值時的約束條件。對 $u、υu、\upsilon$ 的附加條件如下： $(6)min{[?u?x]2+[?u?y]2]+[?υ?x]2+[?υ?y]2}min\{[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2]+[\frac{\partial \upsilon}{\partial x}]^2+[\frac{\partial \upsilon}{\partial y}]^2\}\tag{6}$ 　　式中， $?2u=[?u?x]2+[?u?y]2\nabla ^2u=[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2$ 是 $u$ 的拉普拉斯算子， $?2υ=[?u?x]2+[?u?y]2\nabla ^2\upsilon=[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2$ 是 $υ\upsilon$ 的拉普拉斯算子。綜合式(5)和式(6)，Horn-Schunck算法將光流 $u、υu、\upsilon$ 的計算歸結為如下問題： $(7)min{?(Ixu+Iyυ+It)2+α2[[?u?x]2+[?u?y]2+[?υ?x]2+[?υ?y]2]}min\{\iint (I_xu+I_y\upsilon + I_t)^2 + \alpha ^ 2 [[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2+[\frac{\partial \upsilon}{\partial x}]^2+[\frac{\partial \upsilon}{\partial y}]^2]\}\tag{7}$ 　　因而，可以得到其相應的歐拉-拉格朗日方程，并利用高斯-塞德爾方法進行求解，得到圖像每個次置第一次至第(n+1)次迭代估計 $(un+1,υn+1)(u^{n+1}, \upsilon^{n+1})$ 為：

$(8)un+1=un￣?Ix2￣Ix￣un￣+Iy￣υn￣+Itα2+Ix￣2+Iy￣2u^{n+1} = \overline{u^n}-\overline{I_x^2}\frac{\overline{I_x}\overline{u^n}+\overline{I_y}\overline{\upsilon^n}+I_t}{\alpha^2 +\overline{I_x}^2+ \overline{I_y}^2}\tag{8}$ $(9)υn+1=un￣?Iy2￣Ix￣un￣+Iy￣υn￣+Itα2+Ix￣2+Iy￣2\upsilon^{n+1} = \overline{u^n}-\overline{I_y^2}\frac{\overline{I_x}\overline{u^n}+\overline{I_y}\overline{\upsilon^n}+I_t}{\alpha^2 +\overline{I_x}^2+ \overline{I_y}^2}\tag{9}$

求解過程要得到穩定的解通常需要上百次的迭代。整個迭代過程既與圖像尺寸有關，又與每次的傳遞量（速度的改變量）有關。由迭代公式可以發現，在一些缺乏特征且較為平坦的區域（梯度為0或較小），其速度由迭代公式的第一項決定，該點的速度信息需要從特征較為豐富的區域傳遞過來。為了加快算法的收斂速度，一方面可以用金字塔的層次結構來減小圖像的尺寸擴散，另一方面可以采用增加擴散量的方法來加速算法。

Horn-Schunck算法MATLAB實現程序如下：

1.求導函數 computeDerivatices

function [fx, fy, ft] = computeDerivatives(im1, im2) % 功能：求輸入圖像參考像素點的像素值沿三軸方向的偏導數 % 輸入： % im1-輸入圖像1 % im2-輸入圖像2 % 輸出： % fx-參考像素點的灰度值沿x方向的偏導數 % fy-參考像素點的灰度值沿y方向的偏導數 % fz-參考像素點的灰度值沿z方向的偏導數 if size(im2,1)==0im2=zeros(size(im1)); end % 利用標準模板求得式（5）中的偏導數Ix, Iy, It fx = conv2(im1,0.25* [-1 1; -1 1],'same') + conv2(im2, 0.25*[-1 1; -1 1],'same'); fy = conv2(im1, 0.25*[-1 -1; 1 1], 'same') + conv2(im2, 0.25*[-1 -1; 1 1], 'same'); ft = conv2(im1, 0.25*ones(2),'same') + conv2(im2, -0.25*ones(2),'same');

2.高斯濾波函數 gaussFilter

function G=gaussFilter(segma,kSize) % 功能：實現高斯濾波 % 輸入： % sigma-高斯分布的概率密度函數的方差 % kSize-高斯向量的模板尺寸大小 % 輸出： % G-方差為segma，大小為kSize的一維高斯向量模板 if nargin<1segma=1; end if nargin<2kSize=2*(segma*3); end % 利用均值為0，方差為segma高斯分布概率密度函數求解一維高斯向量模板 G=(1/(sqrt(2*pi)*segma)) * exp (-(x.^2)/(2*segma^2));

3.平滑性約束條件函數 smoothImg

function smoothedImg=smoothImg(img,segma) % 功能：實現平滑性約束條件 % 輸入： % img-數字圖像 % sigma-高斯分布的方差 % 輸出： % smoothedImg-經高斯濾波的圖像矩陣if nargin<2segma=1; end %調用高斯濾波函數 G=gaussFilter(segma); %根據（6）對圖像進行平滑約束 smoothedImg=conv2(img,G,'same'); % 二次平滑 smoothedImg=conv2(smoothedImg,G','same');

4.求解光流場函數 HS

function [u, v] = HS(im1, im2, alpha, ite, uInitial, vInitial, displayFlow, displayImg) % 功能：求解光流場 % 輸入： % im1-輸入圖像1 % im2-輸入圖像2 % alpha-反映HS光流算法的平滑性約束條件的參數 % ita-（8）和（9）式中的迭代次數 % uInitial-光流橫向分量初始值 % vInitial-光流縱向分量初始值 % displayFlow-光流場顯示參數，其值為1時顯示，為0時不顯示 % displayImg-顯示光流場的指定圖像，如果為空矩陣，則無指定圖像輸出 % 輸出： % u-橫向光流矢量 % v-縱向光流矢量% 初始化參數 if nargin<1 || nargin<2im1=imread('HS1.tif');im2=imread('HS2.tif'); end if nargin<3alpha=1; end if nargin<4ite=100; end if nargin<5 || nargin<6uInitial = zeros(size(im1(:,:,1)));vInitial = zeros(size(im2(:,:,1))); elseif size(uInitial,1) ==0 || size(vInitial,1)==0uInitial = zeros(size(im1(:,:,1)));vInitial = zeros(size(im2(:,:,1))); end if nargin<7displayFlow=1; end if nargin<8displayImg=im1; end % 將RGB圖像轉化為灰度圖像 if size(size(im1),2)==3im1=rgb2gray(im1); end if size(size(im2),2)==3im2=rgb2gray(im2); end im1=double(im1); im2=double(im2); % 調用平滑性約束函數對圖像進行平滑 im1=smoothImg(im1,1); im2=smoothImg(im2,1); tic; % 為光流矢量設置初始值 u = uInitial; v = vInitial; [fx, fy, ft] = computeDerivatives(im1, im2); % 調用求導函數對時間分量和空間分量進行求導 kernel_1=[1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12]; % 均值模板 % 根據式（3.15.9）迭代求解，迭代次數為100 for i=1:ite% 計算光流矢量的局部均值uAvg=conv2(u,kernel_1,'same');vAvg=conv2(v,kernel_1,'same');% 根據式（3.15.9）用迭代法求解光流矢量u= uAvg - ( fx .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha^2 + fx.^2 + fy.^2); v= vAvg - ( fy .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha^2 + fx.^2 + fy.^2); end u(isnan(u))=0; v(isnan(v))=0; % 畫圖 if displayFlow==1plotFlow(u, v, displayImg, 5, 5); % 調用畫圖函數 end

5.畫圖函數 plotFlow

function plotFlow(u, v, imgOriginal, rSize, scale) % 功能：繪制光流場圖 % 輸入：u-橫向光流矢量 v-縱向光流矢量 imgOriginal-光流場顯示的圖像 % rSize-可見光流矢量區域的尺寸 scale-光流場規模 figure(); if nargin>2if sum(sum(imgOriginal))~=0imshow(imgOriginal,[0 255]);hold on;end end if nargin<4rSize=5; end if nargin<5scale=3; end for i=1:size(u,1)for j=1:size(u,2)if floor(i/rSize)~=i/rSize || floor(j/rSize)~=j/rSizeu(i,j)=0;v(i,j)=0;endend end quiver(u, v, scale, 'color', 'b', 'linewidth', 2); set(gca,'YDir','reverse');

運行結果如下所示：

參考
https://blog.csdn.net/qq_41368247/article/details/82562165
《現代數字圖像-處理技術提高及應用案例詳解》趙小川

總結

以上是生活随笔為你收集整理的HS光流法及其代码示例的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。