apply function with two variables

I would like to apply the 2D Gaussian function for several values of x and y, for example x = 1:1280 and y = 1:1024. The other input arguments are scalar.

function [ intensity_value ] = twoDimensionalGaussian( x, y, x_centre, y_centre, sigma_x, sigma_y, sigma_xy, amplitude, background)

x_normalised = (x - x_centre)/sigma_x;
y_normalised = (y - y_centre)/sigma_y;

intensity_value = background + amplitude*exp(-(x_normalised^2 + y_normalised^2 + 2*sigma_xy/(sigma_x*sigma_y)*x_normalised*y_normalised));

end

I could write the following for loop:

nb_pixel_horiz = 1280;
nb_pixel_verti = 1024;
I = zeros(nb_pixel_horiz,nb_pixel_verti);
for x = 1:nb_pixel_horiz
  for y = 1:nb_pixel_verti
   intensity_value = twoDimensionalGaussian( x, y, x_centre, y_centre, sigma_x, sigma_y, sigma_xy, amplitude, background);
  I(x,y) = intensity_value;
  end
end

Would this solution be the most appropriate or is there another one that is faster or easier?

Thank you!

Answer 1

[Chris Turnes]

You should be able to make your code significantly faster by vectorizing it (see this article on vectorization) instead of using two for loops.

In your case, since you have two vector inputs x and y and you want to evaluate a 2D Gaussian function at each pair of points, you can use the meshgrid function to create a grid for you (see the documentation for more information) .

Using meshgrid, a vectorized form of your function might look like

function [ intensity_value ] = twoDimensionalGaussian( x, y, x_centre, y_centre, sigma_x, sigma_y, sigma_xy, amplitude, background)

% create a meshgrid for the normalized x and y values
[x_normalised, y_normalised] = meshgrid((x - x_centre)/sigma_x, (y - y_centre)/sigma_y);

% compute the Gaussian values
intensity_value = background + amplitude*exp(-(x_normalised.^2 + y_normalised.^2 + 2*sigma_xy/(sigma_x*sigma_y)*x_normalised.*y_normalised));

end

Notice in the code above that the * and ^ operators were replaced with the .* and .^ operators, which work element-wise (see the differences between the documentations of the .* operator and the * operator if the distinction is not clear).

Then instead of looping over the possible x and y values, you could call your function as

I = twoDimensionalGaussian(1:1280, 1:1024, x_centre, y_centre, sigma_x, sigma_y, amplitude, background);

If you have the Statistics Toolbox, you might also consider using the built-in mvnpdf function instead, which can be used to compute values of a 2D Gaussian function (see the documentation page for more information).

[Guillaume]

Thank you, this is what I needed.