2d gaussian function

I need to plot a 2d gaussian function, where x and y corresponds to the image pixels, my code uses a nested for loop which makes my program run extremely slow, is there a way to write this in a more faster way?

(right now it takes about 8-10 sec to run on 1920*1080 size matrix and i need to produce 173,340 2d gaussian functions images which is too much time...)

thanks in advance for any help.

this is my current code:

function mat = gauss2d(mat, sigma, center)
gsize = size(mat);
for r=1:gsize(1)
    for c=1:gsize(2)
        mat(r,c) = gaussC(r,c, sigma, center);
    end
end

function val = gaussC(x, y, sigma, center)
xc = center(1);
yc = center(2);
exponent = ((x-xc).^2 + (y-yc).^2)./(2*sigma);
val       = (exp(-exponent));

[Janak]

I found a small error in gaussC function. Instead of 2*sigma in exponent, it should be 2*sigma^2.

[Alan]

The val should be changed to "val = 1./2*pi*sigma^2*(exp(-exponent))".

[Image Analyst]

Why? Neither your nor ital's equation is a Gaussian (normal) distribution function. A Gaussian would be

exponent = ((x-xc).^2 + (y-yc).^2)./(2*sigma^2);
amplitude = 1 / (sigma * sqrt(2*pi));  
% The above is very much different than Alan's "1./2*pi*sigma^2"
% which is the same as pi*Sigma^2 / 2.
val       = amplitude  * exp(-exponent);

[Walter Roberson]

I know that I do not have any code for coherence length gaussian beam on multibeam interference.

[Rik]

It sounds like you are getting a lesson in planning then. Why is your bad planning an excuse to expect other to drop everything and help you? This sounds like a relatively big project. Did you make a plan how to break it up in solvable pieces?

[Walter Roberson]

If you need a program that large written in that short of a time, you need to hire a team of about 10 people experienced in rush programming. I would expect that the project cost would be at least a quarter million dollars.

[Walter Roberson]

It is already Saturday in Atlantic Time and everywhere east of that. There is zero chance that anyone is going to write the code for you in the next few hours before you meet with your professor.

However, if you ask specific questions about MATLAB then there is a possibility that someone might answer in the next few hours.

Answer 1

[Walter Roberson]

function mat = gauss2d(mat, sigma, center)
gsize = size(mat);
[R,C] = ndgrid(1:gsize(1), 1:gsize(2));
mat = gaussC(R,C, sigma, center);

[lital]

thank you so much!!!

[BBB]

Hi, can you please help me create a function for the multivariate gaussian 2-D distribution for the following given equation:

.

The function has to give a final plot of the gaussian bump using the imagesc in 2D.

[Walter Roberson]

1/2 .* pi .* sigma_x .* sigma(y) .* exp( - (x.^2)./2 .* sigma(x).^2 - (y.^2)./2 .* sigma(y).^2)

However you need to check whether the is intended to be or is intended to be and likewise if is intended to be or

[BBB]

Hi Walter,

It's perfect and working! Thanks :)

Answer 2

[Image Analyst]

Is that 8-10 seconds for all 173,340 Gaussians? That's not bad. I have a demo that randomly places Gaussians in a larger image using fspecial() and indexing, not one pixel at a time like you did. It takes about 7.4 seconds on an old computer for 173,340 randomly placed Gaussians in a 1920x1080 image. Here's the code, in case you're interested:

% Demo to randomly place Gaussians in an image.
% By ImageAnalyst
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
tic;
% Set up some parameters.
fontSize = 20;
backgroundGrayLevel = 128;
windowSize = 50; % Could be random if you want.
sigma = 10; % Could be random if you want.
numberOfGaussians = 173340;
rows = 1080;
columns = 1920;
% Create one Gaussian.
g = fspecial('gaussian', windowSize, sigma);
grayImage = backgroundGrayLevel * ones(rows, columns);
% Create random signs so that the Gaussians are
% randomly brighter or darker than the background.
s = 2*randi(2, [1 numberOfGaussians])-3;
% Note: g and grayImage are floating point images, not uint8,
% though you could modify the program to have them be uint8 if you wanted.
% Get a list of random locations.
randomRow = randi(rows-windowSize+1, [1 numberOfGaussians]);
randomCol = randi(columns-windowSize+1, [1 numberOfGaussians]);
% Place the Gaussians on the image at those random locations.
for k = 1 : numberOfGaussians
  grayImage(randomRow(k):randomRow(k)+windowSize-1, randomCol(k):randomCol(k)+windowSize-1) = ...
    grayImage(randomRow(k):randomRow(k)+windowSize-1, randomCol(k):randomCol(k)+windowSize-1) + ...
    s(k) * g;
end
toc;
% Display the final image.
imshow(grayImage, []);
caption = sprintf('%d Gaussians, Randomly Placed', numberOfGaussians);
title(caption, 'FontSize', fontSize);
axis on;

% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1])
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off')

[Muhammed Sameed]

This is a really good piece of code. But how would you do the reverse i.e. from a noisy image with N randomly placed Gaussians, extract the centroid location and sigma of each Gaussian? I am interested in N=2 for the moment, but a generalized function would be helpful.

[Walter Roberson]

Hmmm, would that be like finding the values of x(n), y(n), width(n), height(n) such that the sum of the gaussians generated by those parameters is everywhere less than (image + 1/2*EPS(image)), where EPS(image) is 1 for integer-valued images and eps(image) for floating point? Which is to say, that the sum of the gaussians can exceed the value of the image at any given location, but must not exceed it by enough that the total would become the next representable number ? For example on an integer image, if an image location was 42, and the sum of gaussians predicted 42.42 there, then that would be okay because uint8(42.42) would be 42, but predicting 42.52 would not be okay because uint8(42.52) would be 43.

Answer 3

[Mahdi]

I faced the same problem, just so others know you can use fspecial('gaussian', hsize, sigma) intrinsic function.

[Image Analyst]

Thanks for calling it out specially - it was kind of buried in my demo code and hard to notice.

Answer 4

[Logan]

where do you get the gaussC function?

[Walter Roberson]

The code for gaussC is shown in the original Question by lital.

Answer 5

[Karbala'a Unvi. Science]

Dear Sir, I am interested about the code that you wrote about the 2D Gaussian. I have a problem that I want to an image data to be distributed in another image ( image A is the Original, image B is the data one) so that when you see image A you find that there is a noise in it ( where that noise is image B)... I hope that is a good information to help me in building the code... Help will be appreciated.. Thank you in advance

[Image Analyst]

This is not an answer. Who are you talking to? You should have put it as a comment under their answer. Why don't you just try your best and then post your code as a new question? And explain it better there. I don't even understand what you want. I don't know if B is a "data" image or a "noise" image and if you can just add B to A.

Answer 6

[rusgu fcf]

what does the variable 'center' signify in the code of gaussian2d by lital ?

[Image Analyst]

It's the location of the center of the Gaussian - where the peak is.

[Rym Benchaabane]

is center a real number ie 0? or a two element array ie [0,0]?

Thank you

[Image Analyst]

It's a 2-element vector, as you can see by looking inside his gaussC() function.

I recommend you use my code or Walter's code since lital's had problems.