Square wave grating does not display clear white and black bars in output

Pre-warning: Uncomfortable visual stimuli included in this question

I'm trying to create square wave grating patterns at spatial frequencies of 3 and 14 cycles per degree inside an ellipse shape. First, I have a working example when I create the 3 cycles per degree grating:

% Clear background
close all
clc
clear
% Parameters
visual_frequency_cpd = 6;        % Visual frequency in cycles per degree (visual_frequency_cpd/2 for correct value when rescaled in the lab)
viewing_distance_cm = 80;       % Viewing distance in centimeters
grating_size = 499.2;             % Size of the grating in pixels (adjust to change overall size)
num_cycles = visual_frequency_cpd * (grating_size / viewing_distance_cm);  % Number of cycles in the grating
pixels_per_cycle = 50;          % Number of pixels per cycle (adjust to change spatial frequency)
% Image resolution (pixels per millimeter)
resolution = 10;
% Convert ellipse dimensions from millimeters to pixels
ellipse_width_mm = 198;
ellipse_height_mm = 154;
ellipse_width_pixels = ellipse_width_mm * resolution;
ellipse_height_pixels = ellipse_height_mm * resolution;
% Calculate the total number of bars
num_bars = num_cycles * 2;
% Calculate the height of each bar
bar_height = ellipse_height_pixels / num_bars;
% Create the grating matrix initialised in white
grating = ones(ellipse_height_pixels, ellipse_width_pixels);
% Calculate the center coordinates of the ellipse
center_x = (ellipse_width_pixels + 1) / 2;
center_y = (ellipse_height_pixels + 1) / 2;
% Calculate the radii of the ellipse
ellipse_radius_x = (ellipse_width_pixels - 1) / 2;
ellipse_radius_y = (ellipse_height_pixels - 1) / 2;
% Calculate the phase shift needed to start and finish with the same color
phase_shift = 0.5 * bar_height;
for x = 1:ellipse_width_pixels
    for y = 1:ellipse_height_pixels
        % Calculate the distance from the center of the ellipse
        distance_x = x - center_x;
        distance_y = y - center_y;
        
        % Check if the pixel is inside the ellipse
        inside_ellipse = ((distance_x / ellipse_radius_x) ^ 2 + (distance_y / ellipse_radius_y) ^ 2) <= 1;
        
        % Set the intensity value for the current pixel
        if inside_ellipse
            % Calculate the phase for the current position
            phase = mod((y - phase_shift) / bar_height, 2);
            if phase >= 1
                grating(y, x) = 0.25;  % Set black bars inside the ellipse
            end
        else
            % Set the pixel to grey if it's outside the ellipse
            grating(y, x) = 0.5;
        end
    end
end
% Scale the grating to the desired size
grating = imresize(grating, grating_size/ellipse_height_pixels);
% Display the grating with preserved aspect ratio
imshow(grating, 'Border', 'tight');
colormap(gray);

Output:

I have a similar code for drawing the 14 cycles per degree grating here:

% Clear background
close all
clc
clear

% Parameters
visual_frequency_cpd = 28;        % Visual frequency in cycles per degree (visual_frequency_cpd/2 for correct value when rescaled in the lab)
viewing_distance_cm = 80;       % Viewing distance in centimeters
grating_size = 499.4;             % Size of the grating in pixels (adjust to change overall size)

num_cycles = visual_frequency_cpd * (grating_size / viewing_distance_cm);  % Number of cycles in the grating

% Image resolution (pixels per millimeter)
resolution = 10;

% Convert ellipse dimensions from millimeters to pixels
ellipse_width_mm = 198;
ellipse_height_mm = 154;
ellipse_width_pixels = ellipse_width_mm * resolution;
ellipse_height_pixels = ellipse_height_mm * resolution;

% Calculate the total number of bars
num_bars = num_cycles * 2;

% Calculate the height of each bar
bar_height = ellipse_height_pixels / num_bars;

% Create the grating matrix initialised in white
grating = ones(ellipse_height_pixels, ellipse_width_pixels);

% Calculate the center coordinates of the ellipse
center_x = (ellipse_width_pixels + 1) / 2;
center_y = (ellipse_height_pixels + 1) / 2;

% Calculate the radii of the ellipse
ellipse_radius_x = (ellipse_width_pixels - 1) / 2;
ellipse_radius_y = (ellipse_height_pixels - 1) / 2;

% Calculate the phase shift needed to start and finish with the same color
phase_shift = 0.5 * bar_height;

for x = 1:ellipse_width_pixels
    for y = 1:ellipse_height_pixels
        % Calculate the distance from the center of the ellipse
        distance_x = x - center_x;
        distance_y = y - center_y;
        
        % Check if the pixel is inside the ellipse
        inside_ellipse = ((distance_x / ellipse_radius_x) ^ 2 + (distance_y / ellipse_radius_y) ^ 2) <= 1;
        
        % Set the intensity value for the current pixel
        if inside_ellipse
            % Calculate the phase for the current position
            phase = mod((y - phase_shift) / bar_height, 2);
            if phase >= 1
                grating(y, x) = 0.25;  % Set black bars inside the ellipse
            end
        else
            % Set the pixel to grey if it's outside the ellipse
            grating(y, x) = 0.5;
        end
    end
end

% Scale the grating to the desired size
grating = imresize(grating, grating_size/ellipse_height_pixels);

% Display the grating with preserved aspect ratio
imshow(grating, 'Border', 'tight');
colormap(gray);

Output:

As you can see, I don't get the clear divisions between black and white lines for the 14 cycles per degree output. Does anybody know how to draw a grating at 14 cycles per degree with clear divisions between the black and white lines?

[KALYAN ACHARJYA]

The lines are clearly visible. Please zoom in for a closer look

[Haydn Farrelly]

Hi @KALYAN ACHARJYA, as I say it does not have clearly defined black and white lines, as would be expected from a square wave grating. These are various shades of gray.

Answer 1

[Voss]

imresize introduces grey stripes by default (default method is 'bicubic'), but you can change the method imresize uses and see if you get something you like better.

The following is the result from your code with visual_frequency_cpd = 28;, zoomed in to the bottom 5% of the image. You can see the grey stripes in with the black and white stripes. In fact they all seem to be various shades of grey, with few if any actually black or actually white.

% get_grating() is your code, made into a function that takes
% visual_frequency_cpd and a resize method for imresize() as its inputs,
% and returns grating (which is the same as your variable grating) and 
% grating_original, which is the image data before resizing via imresize().
[grating,grating_o] = get_grating(28,'bicubic');

figure
imshow(grating, 'Border', 'tight');
colormap(gray);
set(gca(),'DataAspectRatioMode','auto');
ylim(gca().YLim(2)*[0.95 1])

What follows here is the same result but before imresize is applied, i.e., it's your variable grating but the version before imresize not the final version, again zoomed in to the bottom 5%. You can see the original contains only black and white stripes. This is just to illustrate that it is in fact imresize introducing the grey.

figure
imshow(grating_o, 'Border', 'tight');
colormap(gray);
set(gca(),'DataAspectRatioMode','auto');
ylim(gca().YLim(2)*[0.95 1])

I assume that you need the final result to be a certain size, and that's why you're using imresize. If you also want the final result to contain only black and white (no grey) stripes, then one thing you can do is to change the method used by imresize to 'nearest'. This is the result of doing that:

grating = get_grating(28,'nearest');

figure
imshow(grating, 'Border', 'tight');
colormap(gray);
set(gca(),'DataAspectRatioMode','auto');
ylim(gca().YLim(2)*[0.95 1])

Or, without zooming in:

figure
% Display the grating with preserved aspect ratio
imshow(grating, 'Border', 'tight');
colormap(gray);

% get_grating() is your code, made into a function that takes
% visual_frequency_cpd and a resize method for imresize() as its inputs,
% and returns grating (which is the same as your variable grating) and 
% grating_original, which is the image data before resizing via imresize().

function [grating,grating_original] = get_grating(visual_frequency_cpd,resize_method)

% Parameters
% visual_frequency_cpd = 2;        % Visual frequency in cycles per degree (visual_frequency_cpd/2 for correct value when rescaled in the lab)
viewing_distance_cm = 80;       % Viewing distance in centimeters
grating_size = 499.4;             % Size of the grating in pixels (adjust to change overall size)

num_cycles = visual_frequency_cpd * (grating_size / viewing_distance_cm);  % Number of cycles in the grating

% Image resolution (pixels per millimeter)
resolution = 10;

% Convert ellipse dimensions from millimeters to pixels
ellipse_width_mm = 198;
ellipse_height_mm = 154;
ellipse_width_pixels = ellipse_width_mm * resolution;
ellipse_height_pixels = ellipse_height_mm * resolution;

% Calculate the total number of bars
num_bars = num_cycles * 2;

% Calculate the height of each bar
bar_height = ellipse_height_pixels / num_bars;

% Create the grating matrix initialised in white
grating = ones(ellipse_height_pixels, ellipse_width_pixels);

% Calculate the center coordinates of the ellipse
center_x = (ellipse_width_pixels + 1) / 2;
center_y = (ellipse_height_pixels + 1) / 2;

% Calculate the radii of the ellipse
ellipse_radius_x = (ellipse_width_pixels - 1) / 2;
ellipse_radius_y = (ellipse_height_pixels - 1) / 2;

% Calculate the phase shift needed to start and finish with the same color
phase_shift = 0.5 * bar_height;

for x = 1:ellipse_width_pixels
    for y = 1:ellipse_height_pixels
        % Calculate the distance from the center of the ellipse
        distance_x = x - center_x;
        distance_y = y - center_y;
        
        % Check if the pixel is inside the ellipse
        inside_ellipse = ((distance_x / ellipse_radius_x) ^ 2 + (distance_y / ellipse_radius_y) ^ 2) <= 1;
        
        % Set the intensity value for the current pixel
        if inside_ellipse
            % Calculate the phase for the current position
            phase = mod((y - phase_shift) / bar_height, 2);
            if phase >= 1
                grating(y, x) = 0.25;  % Set black bars inside the ellipse
            end
        else
            % Set the pixel to grey if it's outside the ellipse
            grating(y, x) = 0.5;
        end
    end
end

grating_original = grating;

% Scale the grating to the desired size
grating = imresize(grating, grating_size/ellipse_height_pixels,resize_method);

end

[Haydn Farrelly]

Thank you, that makes a lot of sense! The only problem with that output is the irregulairty of the width of the bars which would affect the visual frequency, do you know how I could restore the regular bar height?

[Voss]

You're welcome!

Rather than have the code construct a matrix of some size and then imresize it to the desired final size, I would try to have the code construct the matrix to be the final size from the start.