基于光流法相位提取算法---MATLAB实现

基于光流法相位提取算法原理的参考文献为:点击打开链接

%% *******************************************
%% **********************************************
clear;close all;clc;
N = 512;
xmax = 1;
ymax = 1;
delta = [0,pi/3];
x = linspace(-xmax,xmax,N);
y = linspace(-ymax,ymax,N);
[X,Y] = meshgrid(x,y);
%参数设置,背景光强,调制幅度,物体相位
A = 0.2*exp(-1.8*(X.^2+Y.^2));
B = 0.2*exp(-0.2*(X.^2+Y.^2));
phi = pi*5*(X.^2 + Y.^2);
figure;mesh(phi);
phi_w = mod(phi,2*pi);
figure;imshow(phi_w,[]);
M = 2;
I = zeros(N,N,M);


for k = 1:M
    I(:,:,k) = A + B.*cos(phi+delta(k)) + 0.05*rand(N,N);
        [NR, NC]=size(I(:,:,1));
        [u,v]=meshgrid(1:NC, 1:NR);
        %Temporal variables
        u0=floor(NC/2)+1; v0=floor(NR/2)+1;           
        u=u-u0; v=v-v0; 
        %High pass Fourier filtering by Gaussian filter with sigma=freq         
        freq = 3;
        H=1-exp(-(u.^2+v.^2)/(2*(freq)^2));
        C=fft2(I(:,:,k));    
        CH=C.*ifftshift(H);    
        I(:,:,k)=real(ifft2(CH));  
end
phi_cor = rofDemod(I(:,:,1),I(:,:,2));
figure;imshow(phi_cor,[]);
程序中rofDemod函数见: rofDemo

  • 1
    点赞
  • 14
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
Lucas-Kanade光流算法matlab实现: ```matlab function [u,v] = lucaskanade(frame1,frame2,window_size) % convert frames to grayscale if size(frame1, 3) == 3 frame1 = rgb2gray(frame1); end if size(frame2, 3) == 3 frame2 = rgb2gray(frame2); end % calculate gradients fx = conv2(frame1,[-1 1; -1 1],'valid'); fy = conv2(frame1,[-1 -1; 1 1],'valid'); ft = conv2(frame2, ones(2), 'valid') + conv2(frame1, -ones(2), 'valid'); u = zeros(size(frame1)); v = zeros(size(frame2)); % iterate through each pixel in the frame for i = window_size+1:size(fx,1)-window_size-1 for j = window_size+1:size(fx,2)-window_size-1 % extract the window around the pixel Ix = fx(i-window_size:i+window_size, j-window_size:j+window_size); Iy = fy(i-window_size:i+window_size, j-window_size:j+window_size); It = ft(i-window_size:i+window_size, j-window_size:j+window_size); % flatten the windows into vectors Ix = Ix(:); Iy = Iy(:); b = -It(:); A = [Ix Iy]; % solve the equation Ax = b if rank(A'*A) >= 2 nu = pinv(A)*b; else nu = [0;0]; end u(i,j)=nu(1); v(i,j)=nu(2); end end end ``` 基于Kalman滤波的光流算法matlab实现: ```matlab function [u,v] = kalmanflow(frame1,frame2,window_size) % convert frames to grayscale if size(frame1, 3) == 3 frame1 = rgb2gray(frame1); end if size(frame2, 3) == 3 frame2 = rgb2gray(frame2); end % calculate gradients fx = conv2(frame1,[-1 1; -1 1],'valid'); fy = conv2(frame1,[-1 -1; 1 1],'valid'); ft = conv2(frame2, ones(2), 'valid') + conv2(frame1, -ones(2), 'valid'); u = zeros(size(frame1)); v = zeros(size(frame2)); % iterate through each pixel in the frame for i = window_size+1:size(fx,1)-window_size-1 for j = window_size+1:size(fx,2)-window_size-1 % extract the window around the pixel Ix = fx(i-window_size:i+window_size, j-window_size:j+window_size); Iy = fy(i-window_size:i+window_size, j-window_size:j+window_size); It = ft(i-window_size:i+window_size, j-window_size:j+window_size); % flatten the windows into vectors Ix = Ix(:); Iy = Iy(:); b = -It(:); A = [Ix Iy]; % initial guess for velocity nu = pinv(A)*b; % initialize Kalman filter x = nu; P = eye(2); Q = 0.1*eye(2); R = 0.1*eye(2); % iterate through each frame for k = 1:10 % predict next state x = A*x; P = A*P*A' + Q; % update state based on measurement z = [fx(i+k,j+k); fy(i+k,j+k)]; K = P*A'*inv(A*P*A' + R); x = x + K*(z - A*x); P = (eye(2) - K*A)*P; end u(i,j)=x(1); v(i,j)=x(2); end end end ``` 基于Horn-Schunck光流算法matlab实现: ```matlab function [u, v] = hornschunck(frame1, frame2, alpha, iterations) % convert frames to grayscale if size(frame1, 3) == 3 frame1 = rgb2gray(frame1); end if size(frame2, 3) == 3 frame2 = rgb2gray(frame2); end % calculate gradients fx = conv2(frame1, [-1 1; -1 1], 'valid'); fy = conv2(frame1, [-1 -1; 1 1], 'valid'); ft = conv2(frame2, ones(2), 'valid') + conv2(frame1, -ones(2), 'valid'); % initialize velocity vectors u = zeros(size(frame1)); v = zeros(size(frame2)); % iterate through each pixel in the frame for i = 1:iterations % calculate average velocity for each pixel u_avg = conv2(u, 0.25*[1 1; 1 1], 'same'); v_avg = conv2(v, 0.25*[1 1; 1 1], 'same'); % calculate intermediate variables rho = fx.*u_avg + fy.*v_avg + ft; delta = alpha^2 + fx.^2 + fy.^2; % update velocity vectors u = u_avg - fx.*rho./delta; v = v_avg - fy.*rho./delta; end end ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值