Using a complex steerable filter on an image

I'm trying to use a complex steerable filter on an image to give me different phase-images for a project. The whole process is documented in Pintea, et al "Hand tremor frequency estimation in videos". To do this I used the first derivative of the Gaussian and took the hilbert transform of it to give me a complex filter. However, after filtering it and trying to reconstruct it, I only get a blank gray image. I've tried just filtering the image with the first derivative of the Gaussian and it works. However, taking the hilbert transform of the Gaussian first derivative doesn't. I tried to get the phase of the complex filter and reconstrucintg using ifft2(). I've also tried to normalize the values from 0 to 255 after reconstructing using mat2gray(), which didn't help either. My code is shown below. Any help would be much appreciated!

% Define orientations and scales
orientations = [0, pi/4, pi/2, 3*pi/4]; % angles in radians
sigma = [1 0.5 0.25]; % standard deviation of Gaussian/scales
filter_size = 2*ceil(2*sigma)+1; % filter size

% Load input image 
input_image = imread('surfing.jpg');
targetSize = [300 600];
r = centerCropWindow2d(size(input_image),targetSize);
input_image = imcrop(input_image,r);
input_image = rgb2gray(input_image);

% Perform cropping or select the region for the hand image
region = input_image;

% Initialize complex steerable pyramid
complex_pyramid = cell(length(orientations) * length(sigma), 1);


index = 1;
for i = 1:length(orientations)
    for j = 1:length(sigma)
        % Create complex steerable filters G and H for each orientation and scale

        % Create a Gaussian filter
        G = fspecial('gaussian', [filter_size(j), filter_size(j)], sigma(j));
        % Compute the gradient of the Gaussian filter
        [Gx, Gy] = gradient(G);

        G1 = cos(orientations(i)) * Gx + sin(orientations(i)) * Gy; % Combine gradients based on orientation
        rotated_filter = imrotate(G1, -orientations(i)*180/pi, 'loose'); % Negative angle for correct orientation

        % Hilbert transform of Gaussian
        h = hilbert(rotated_filter);
        complex_filt{i+j-1} = h;

    end
end

% testing complex gaussian steerability
filtered_images = cell(1, numel(orientations));
for i = 1:numel(orientations)
    for j = 1:length(sigma)
        conv_image = conv2(complex_filt{i},input_image);
        
        % get phase only to normalize image
        conv_image = double(conv_image);
        phase = angle(conv_image);
        complex_phase = exp(1j*phase);
        
        % reconstruct image
        reconstruct = ifft2(complex_phase);
        filtered_images{i+j-1} = reconstruct
        filtered_images{i+j-1} = conv_image;
    end
end
figure;
imagesc(input_image); colormap gray;

% visualize
for i = 1:numel(orientations)
    for j = 1:length(sigma)
        figure;
        imagesc(mat2gray(real(filtered_images{i}))*255); colormap gray;
        title(['Orientation: ', num2str(orientations(i)), ', Scale: ', num2str(sigma(j))]);

        % Display the imaginary and real parts of complex filters for debugging
        figure;
        subplot(1, 2, 1);
        imagesc(real(complex_filt{i+j-1})); colormap gray;
        title(['Real Part - Orientation: ', num2str(orientations(i)), ', Scale: ', num2str(sigma(j))]);
        
        subplot(1, 2, 2);
        imagesc(imag(complex_filt{i+j-1})); colormap gray;
        title(['Imaginary Part - Orientation: ', num2str(orientations(i)), ', Scale: ', num2str(sigma(j))]);
    end
end

Answer 1

[Shubh]

Hi Raymond,

It looks like there are a few issues with your code that could be causing the problems you're experiencing. Let's go through each part and make the necessary corrections.

1. Complex Filter Creation: The process of creating the complex filter using the Hilbert transform seems correct. However, you need to ensure that the size of the Gaussian filter (filter_size) is properly computed for each scale (sigma). The current calculation might not be accurate.
2. Filtering Image with Complex Filter: When you convolve the image with the complex filter, it's important to ensure that the convolution is done correctly. MATLAB's "conv2" function, by default, does a full convolution, which might not be what you want. You should consider using the 'same' option to get an output of the same size as the input image.
3. Phase Extraction and Reconstruction: The phase extraction seems fine, but the reconstruction using "ifft2" might not be appropriate here. The inverse Fourier transform is used to convert frequency-domain data back to the spatial domain. However, in your case, you're working directly in the spatial domain, so using "ifft2" on the phase image might not be correct. Instead, you should be visualizing the phase image directly.
4. Normalization using "mat2gray": The way you're using "mat2gray" is not standard. "mat2gray" automatically scales the data to fit in the [0, 1] range. Multiplying by 255 afterwards is not necessary and might be causing issues.

Here's a revised version of your code with these points considered:

% Define orientations and scales
orientations = [0, pi/4, pi/2, 3*pi/4]; % angles in radians
sigma = [1, 0.5, 0.25]; % standard deviation of Gaussian/scales

% Load input image 
input_image = imread('surfing.jpg');
targetSize = [300, 600];
r = centerCropWindow2d(size(input_image),targetSize);
input_image = imcrop(input_image,r);
input_image = rgb2gray(input_image);

% Perform cropping or select the region for the hand image
region = input_image;

% Initialize complex steerable pyramid
complex_filt = cell(length(orientations) * length(sigma), 1);

for i = 1:length(orientations)
    for j = 1:length(sigma)
        % Adjust filter size for each scale
        filter_size = 2 * ceil(2 * sigma(j)) + 1;

        % Create a Gaussian filter
        G = fspecial('gaussian', [filter_size, filter_size], sigma(j));
        
        % Compute the gradient of the Gaussian filter
        [Gx, Gy] = gradient(G);
        G1 = cos(orientations(i)) * Gx + sin(orientations(i)) * Gy;
        
        % Hilbert transform of the oriented Gaussian
        h = hilbert(G1);
        complex_filt{i+j-1} = h;
    end
end

% Filtering the image with complex filters
filtered_images = cell(numel(orientations), length(sigma));

for i = 1:numel(orientations)
    for j = 1:length(sigma)
        conv_image = conv2(double(region), double(real(complex_filt{i+j-1})), 'same');
        phase = angle(conv_image);
        
        % Storing the phase images
        filtered_images{i, j} = phase;
    end
end

% Visualizing the phase images
for i = 1:numel(orientations)
    for j = 1:length(sigma)
        figure;
        imagesc(mat2gray(filtered_images{i, j})); colormap gray;
        title(['Orientation: ', num2str(orientations(i)), ', Scale: ', num2str(sigma(j))]);
    end
end

This revised code should help you visualize the phase images correctly. Remember to adjust the path to your image file as needed. Hope this helps!