基于Lucas-Kanade算法的三维光流提取matlab仿真

目录

1.算法描述

2.仿真效果预览

3.MATLAB核心程序

4.完整MATLAB

---

1.算法描述

光流的概念:(Optical flow or optic flow)
       它是一种运动模式，这种运动模式指的是一个物体、表面、边缘在一个视角下由一个观察者（比如眼睛、摄像头等）和背景之间形成的明显移动。光流技术，如运动检测和图像分割，时间碰撞，运动补偿编码，三维立体视差，都是利用了这种边缘或表面运动的技术。

二维图像的移动相对于观察者而言是三维物体移动的在图像平面的投影。
有序的图像可以估计出二维图像的瞬时图像速率或离散图像转移。

光流算法：
      它评估了两幅图像的之间的变形，它的基本假设是体素和图像像素守恒。它假设一个物体的颜色在前后两帧没有巨大而明显的变化。基于这个思路，我们可以得到图像约束方程。不同的光流算法解决了假定了不同附加条件的光流问题。

Lucas–Kanade算法：
      这个算法是最常见，最流行的。它计算两帧在时间t 到t + δt之间每个每个像素点位置的移动。 由于它是基于图像信号的泰勒级数，这种方法称为差分，这就是对于空间和时间坐标使用偏导数。
图像约束方程可以写为I (x ,y ,z ,t ) = I (x + δx ,y + δy ,z + δz ,t + δt )
I(x, y,z, t) 为在（x,y,z）位置的体素。

我们假设移动足够的小，那么对图像约束方程使用泰勒公式，我们可以得到：

H.O.T. 指更高阶，在移动足够小的情况下可以忽略。从这个方程中我们可以得到：

2.仿真效果预览

matlab2022a仿真如下：

3.MATLAB核心程序

function [ux,uy,uz]=LK3D( image1, image2, r )
%This function estimates deformations between two subsequent 3-D images
 
%   June 2012
 
%  Default parameter
if nargin==2
    r=2;
end
 
[height,width,depth]=size(image1); 
image1=im2double(image1);
image2=im2double(image2);
 
% Initializing flow vectors
ux=zeros(size(image1)); uy=ux; uz=ux;
 
% Computing image derivatives
[Ix,Iy,Iz,It]=imageDerivatives3D(image1,image2);
 
for i=(r+1):(height-r)
    for j=(r+1):(width-r)
        for k=(r+1):(depth-r)
        
        blockofIx=Ix(i-r:i+r,j-r:j+r,k-r:k+r);
        blockofIy=Iy(i-r:i+r,j-r:j+r,k-r:k+r);
        blockofIz=Iz(i-r:i+r,j-r:j+r,k-r:k+r);
        blockofIt=It(i-r:i+r,j-r:j+r,k-r:k+r);
 
               
        A=zeros(3,3);
        B=zeros(3,1);
        
        A(1,1)=sum(sum(sum(blockofIx.^2)));
        A(1,2)=sum(sum(sum(blockofIx.*blockofIy)));
        A(1,3)=sum(sum(sum(blockofIx.*blockofIz)));
        
        A(2,1)=sum(sum(sum(blockofIy.*blockofIx)));
        A(2,2)=sum(sum(sum(blockofIy.^2)));
        A(2,3)=sum(sum(sum(blockofIy.*blockofIz)));
 
        A(3,1)=sum(sum(sum(blockofIz.*blockofIx)));
        A(3,2)=sum(sum(sum(blockofIz.*blockofIy)));
        A(3,3)=sum(sum(sum(blockofIz.^2)));
       
        B(1,1)=sum(sum(sum(blockofIx.*blockofIt)));
        B(2,1)=sum(sum(sum(blockofIy.*blockofIt)));
        B(3,1)=sum(sum(sum(blockofIz.*blockofIt)));
        
        invofA=pinv(A);
        
        V=invofA*(-B);
        ux(i,j,k)=V(1,1);
        uy(i,j,k)=V(2,1);
        uz(i,j,k)=V(3,1);
        end
    end
end
end
a92

4.完整MATLAB

V