matlab - 使用riesz变换分析后的相移图像内容

标签 matlab image-processing motion multiscaleimage

我正在从事一个涉及视频运动放大算法的项目。目前我正在尝试使用 riesz 金字塔来理解基于相位的运动放大。我的主要信息来源是这个文档:

Riesz Pyramids for Fast Phase-Based Video Magnification \

我已执行以下步骤来尝试重现论文中的一些结果:

  1. 使用为 riesz 金字塔提供的 matlab 代码将图像分解为多个尺度

  2. 使用论文中介绍的近似 riesz 变换,通过将金字塔的一个子带与 [-0.5, 0, 0.5] 和 [-0.5, 0, 0.5]' 进行卷积,生成图像 Riesz1 和 Riesz2。

  3. 通过计算atan(R2/R1) 确定子带每个像素中的主导局部方向。该计算源自论文中的方程 3。

  4. 将变换引导到主要局部方向并计算所得的正交对

  5. 使用正交对生成复数 (I + iQ),然后使用其相位 I 来确定特定像素中的局部相位。

这是我创建的 Matlab 代码:

%Generate a circle image
img = zeros(512, 512);
img(:) = 128;
rad = 180;
for i = size(img, 1)/2 - rad : size(img,1)/2 + rad
    for j = size(img, 2)/2 - rad : size(img,2)/2 + rad
        deltaX = abs(size(img, 1)/2 - i);
        deltaY = abs(size(img, 2)/2 - j);
        if (sqrt(deltaX^2+deltaY^2) <= rad)
           img(i, j) = 255;
        end
    end
end

%build riesz pyramid
[pyr, pind] = buildNewPyr(img);

%extract band2 from pyramid (no orientation information yet)
I = pyrBand(pyr,pind,3);

%convolve band2 with approximate riesz filter for first quadrature pair
%element
R1 = conv2(I, [0.5, 0, -0.5], 'same');

%convolve band2 with approximate riesz filter (rotated by 90°) for second
%quadrature pair element
R2 = conv2(I, [0.5, 0, -0.5]', 'same');

% show the resulting image containing orientation information!
% imshow(band2_r2, []);

%To extract the phase, we have to steer the pyramid to its dominant local
%orientation. Orientation is calculated as atan(R2/R1)
theta = atan(R2./R1);
theta(isnan(theta) | isinf(theta)) = 0;
%imshow(theta, []);

% create quadrature pair
Q = zeros(size(theta, 1), size(theta, 2));

for i = 1:size(theta, 1)
    for j = 1:size(theta, 1)
        if theta(i, j) ~= 0
            %create rotation matrix
            rot_mat = ...
                [cos(theta(i, j)), sin(theta(i, j));...
                -sin(theta(i, j)) cos(theta(i, j))];

            %steer to dominant local orientation(theta) and set Q
            resultPair = rot_mat*[R1(i, j), R2(i,j)]';
            Q(i,j) = resultPair(1);
        end 
    end
end

% create amplitude and phase image
A = abs(complex(I, Q));
Phi = angle(complex(I, Q));

生成的图像如下所示:

Generated Images

现在我的问题:

  1. 使用 atan(R2/R1) 计算 theta 时,结果中出现很多伪影(参见图像“主导方向”)。我在这里错过了什么明显的事情/做错了吗?

  2. 假设到目前为止我的结果是正确的。为了放大运动,我不仅需要能够确定局部相位,还需要能够改变它。我似乎错过了一些明显的事情,但我该怎么做呢?我是否需要以某种方式改变金字塔子带像素的相位,然后折叠金字塔?如果是,怎么办?

我(显然)对这个主题相当陌生,并且对图像处理只有初步的了解。我将非常感谢任何答案,无论是解决我的问题还是只是引用其他有用的信息来源。

真诚的

最佳答案

我已经得到了该算法的有效实现。以下是我使用此方法成功运动放大视频所采取的步骤。

这些步骤应该应用于您的视频序列的每个 channel (我已经在 RGB 视频中尝试过,您可能可以在 YUV 视频中仅针对亮度执行此操作)。

  1. 创建每帧的图像金字塔。原始论文有一个推荐的金字塔结构,以允许更大的放大值,尽管它与拉普拉斯金字塔配合得相当好。

  2. 对于每个视频 channel 的每个金字塔级别,计算 Riesz 变换(有关变换的概述,请参阅 The Riesz transform and simultaneous representations of phase, energy and orientation in spatial vision,并参阅原始问题中的论文以了解有效的近似实现)。

  3. 使用 Riesz 变换,计算每个视频帧的每个金字塔级别的每个像素的局部幅度、方向和相位。以下 Matlab 代码将使用近似 Riesz 变换计算( double 格式)图像的局部方向、相位和幅度:

    function [orientation, phase, amplitude] = riesz(image)
    
    [imHeight, imWidth] = size(image);
    
    %approx riesz, from Riesz Pyramids for Fast Phase-Based Video Magnification
    
    dxkernel = zeros(size(image));
    dxkernel(1, 2)=-0.5;
    dxkernel(1,imWidth) = 0.5;
    
    
    dykernel = zeros(size(image));
    dykernel(2, 1) = -0.5;
    dykernel(imHeight, 1) = 0.5;
    
    R1 = ifft2(fft2(image) .* fft2(dxkernel));
    R2 = ifft2(fft2(image) .* fft2(dykernel));
    
    
    orientation = zeros(imHeight, imWidth);
    phase = zeros(imHeight, imWidth);
    
    orientation = (atan2(-R2, R1));
    
    phase = ((unwrap(atan2(sqrt(R1.^2 + R2.^2) , image))));
    
    amplitude = sqrt(image.^2 + R1.^2 + R2.^2);
    
    end   
    
  4. 对于每个金字塔级别,使用设置为适合您希望放大的运动的频率的带通滤波器临时过滤每个像素的相位值。请注意,这会消除相位值中的直流分量。

  5. 计算放大相位值

    amplifiedPhase = phase + (requiredGain * filteredPhase);
    
  6. 使用放大相位计算每个金字塔级别的新像素值

    amplifiedSequence = amplitude .* cos(amplifiedPhase);
    
  7. 折叠金字塔以生成新的、放大的视频 channel 。

  8. 将放大的 channel 重新组合到新的视频帧中。

原始论文中还有一些其他步骤可以提高噪声性能,但上面的序列可以很好地生成运动放大视频。

关于matlab - 使用riesz变换分析后的相移图像内容,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/29445339/

相关文章:

linux - MATLAB 中的最大波特率是多少?

arrays - Matlab:使用逻辑索引删除 IF-ELSE 条件语句

algorithm - 从散点评估/拟合椭圆

matlab - 从视频重建 3D 轨迹(由单个相机拍摄)

javascript - 使用 Webkit Transitions 在弧形中移动

javascript - 如何在 Javascript 中使图像无限重复滚动?

c - 如何同时调试C和Matlab代码?

matlab - "Out of Memory"Matlab

html - 与 CSS 相比,使用 ImageMagick 调整图像大小有什么好处?

iphone - 如何将 void(摇动手势)放入 IBAction(按钮)中?