matlab image-processing motion multiscaleimage

Phase shifting image content after analyzing with riesz transform

I am working on a project involving video motion magnification algorithms. Currently I am trying to understand phase based motion magnification using a riesz pyramid. My main source of information is this document:

Riesz Pyramids for Fast Phase-Based Video Magnification \

I have performed the following steps to attempt to reproduce some of the results in the paper:

Decompose an image into multiple scales using the provided matlab code for the riesz pyramid
Generate the images Riesz1 and Riesz2 by convolving one subband of the pyramid with [-0.5, 0, 0.5] and [-0.5, 0, 0.5]' using the approximate riesz transform introduced in the paper.
Determine the dominant local orientation in every pixel of the subband by calculating atan(R2/R1). This calculation is derived from equation 3 in the paper.
Steer the transform to the dominant local orientation and calculate the resulting quadrature pair
Use the quadrature pair to generate a complex number (I + iQ) whose phase I then used to determine the local phase in the specific pixel.

This is the Matlab code I created:

%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));

The generated images look like this:

Generated Images

Now my questions:

When calculating theta using atan(R2/R1) I get a lot of artifacts in the result (see image "dominant orientation"). Is there something obvious I miss here/do wrong?
Assuming what my results are correct thus far. To magnify motion I need to not only be able to determine the local phase, but also to change it. I seem to miss something obvious, but how would I go about that? Do I need to somehow change the phase of the pyramid subband pixels and then collapse the pyramid? If yes, how?

I am (obviously) quite new to this topic and only have a rudimentary understanding of image processing. I would be very thankful for any answer, be it a solution to my problems or just a referral to an other useful source of information.

Sincerely

Solution

I have got a functioning implementation of this algorithm. Here are the steps I took to successfully motion-magnify a video using this method.

These steps should be applied to each channel of a video sequence that you (I have tried it for RGB video, you could probably get away with doing it for just luminance, in a YUV video).

Create an image pyramid of each frame. The original paper has a recommended pyramid structure to allow greater magnification values, although it works fairly well with a Laplacian pyramid.
For each pyramid level of each video channel, calculate the Riesz transform (see The Riesz transform and simultaneous representations of phase, energy and orientation in spatial vision for an overview of the transform, and see the paper in the original question for an efficient approximate implementation).

Using the Riesz transform, calculate the local amplitude, orientation and phase for each pixel of each pyramid level of each video frame. The following Matlab code will calculate the local orientation, phase and amplitude of a (double format) image using the approximate Riesz transform:

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

For each pyramid level, temporally filter the phase values of each pixel using a bandpass filter set to a frequency appropriate for the motion that you wish to magnify. Note that this removes the DC component from the phase value.

Calculate the amplified phase value by

amplifiedPhase = phase + (requiredGain * filteredPhase);

Use the amplified phase to calculate new pixel values for each pyramid level by
```
amplifiedSequence = amplitude .* cos(amplifiedPhase);
```
Collapse the pyramids to generate the new, amplified video channel.
Recombine your amplified channels into a new video frame.

There are some other steps in the original paper to improve noise performance, but the sequence above produces motion amplified video quite nicely.