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).
