Need help with implementing a 2D elliptical Gaussian function

I'm trying to implement a 2D gaussian function, which has an elliptical shape rather than circular. For example if I have a standard Gaussian fuction such as: f=a1*exp(-(((x-b1)/c1).^2+((y-b1)/c1).^2))

where [x,y]=meshgrid(xmin:spacing:xmax,ymin:spacing,yman); I can just f=a1*exp(-((r-b1)/c).^2) where r=sqrt(x^2+y^2);

however if my c1 is not the same, I will have an ellipse. I'm not sure how to implement this. Any help will be highly appreciated.

Answer 1

[David Young]

You can make an elliptical filter aligned with the axes by replacing the second occurrence of c1 by a different variable, say c2, and giving c1 and c2 different values. Here's some code that demonstrates what I mean; I've taken the opportunity to change b1 to b2 as well, so the centre of the filter can be moved to an arbitrary point:

% Set up mesh
xmin = -100;
xmax = 100;
ymin = -100;
ymax = 100;
spacing = 1;
xvals = xmin:spacing:xmax+spacing/2;
yvals = ymin:spacing:ymax+spacing/2;
[x,y] = meshgrid(xvals, yvals);

% parameters for the gaussian
a1 = 1;
b1 = 20;
b2 = 40;
c1 = 10;
c2 = 40;

% Compute the filter and display it
f=a1.*exp(-(((x-b1)./c1).^2+((y-b2)./c2).^2));
contour(x, y, f);

If you want the ellipse to be oriented in an arbitrary direction, you need to rotate the axes before the computation. This involves multiplying x and y by a rotation matrix. Please say if you also need help with this.

[MJ HL]

Hello David, ... As you guessed I'm from that kind of people that need more help in rotating this filter :) . How should I do this? I need to rotate this filter in an arbitrary direction... thanks

Answer 2

[Ronni]

Thanks David. I tried that and it seems to be working. I also had something similar what the problem that I had was that my function was sum of multiple functions. That I still can't get it working. Hopefully I can explain here what I mean.

I have this 2D data, which looks like a combination of gaussians. So since it was centered around zero, to fit this 2D data, I just took 1D profile across the center and fitted it with just using x variable. I assumed I can use the same parameters for y since for my initial test it was just a circular distribution. Later, I will be tweeking it so the FWHM of the added of function of one side is longer than the other. Hence, it will turn into an elliptical multi-gaussian function rather than just a circular mult-gaussian function.

This is what I have written, but the contour looks weird:

       xgrid=-2:0.05:2;
       ygrid=-2:0.05:2;
       [x,y]=meshgrid(xgrid,ygrid);

       f=zeros(size(x));

       a1 =       57.27;
       b1 =      0.5494;
       c1 =      0.3432;
       a2 =       58.35;
       b2 =     -0.5444;
       c2 =      0.3461;
       a3 =       33.19;
       b3 =      0.8118;
       c3 =      0.2174 ;
       a4 =       33.38 ;
       b4 =     -0.8114;
       c4 =      0.2188;
       a5 =       90.86;
       b5 =    0.005504;
       c5 =      0.5079;

       b_r =      -10.81;
       a =        9.4e+005;
       b_l =      10.81;

       f_l_35=    a*exp(b_l*x+b_l*y);        
       f_m_35=    a1*exp(-(((x-b1)/c1).^2+((y-b1)/c1).^2)) ...
                  + a2*exp(-(((x-b2)/c2).^2+((y-b2)/c2).^2)) ...
                  + a3*exp(-(((x-b3)/c3).^2+((y-b3)/c3).^2)) ...
                  + a4*exp(-(((x-b4)/c4).^2+((y-b4)/c4).^2))...
                  + a5*exp(-(((x-b5)/c5).^2+((y-b5)/c5).^2));

       f_r_35=    a*exp(b_r*x+b_r*y);

%I'm fitting the tails to the exponential function.

          f(find(x<=-1))= f_l_35(find(x<=-1));
          f(find(y<=-1))= f_l_35(find(y<=-1));

%Middle part to mult-gaussian

          f(find(x<1 & x>-1))= f_m_35(find(x<1 & x>-1));
          f(find(y<1 & y>-1))= f_m_35(find(y<1 & y>-1));

%again my right tail to exponential

          f(find(x>=1))= f_r_35(find(x>=1));
          f(find(y>=1))= f_r_35(find(y>=1));

Again, I got these coefficients by fitting central 1D profile. I kind of know why the distribution looks weird. I think when I did 1D profile, I did it at y=0 and so I also assumed b1_y (written as b1) to be zero. When I plot in 2D, it's shifting the centeral profile I fitted by b1/b2/b3 etc.

The distribution should really be looking like

    f_m_35_revised=  a1*exp(-(((r-b1)/c1).^2)) ...
                  + a2*exp(-(((r-b2)/c2).^2))...
                  + a3*exp(-(((r-b3)/c3).^2))...
                  + a4*exp(-(((r-b4)/c4).^2))...
                  + a5*exp(-(((r-b5)/c5).^2));

where r = sqrt(x.^2+y.^2);

Since I know my data is center about zero. May be I will change my fitting model and change the b parameters to 0. so more like f_m_35(x)= a1*exp(-(((x)/c1).^2)) ... + a2*exp(-(((x)/c2).^2))... + a3*exp(-(((x)/c3).^2))... + a4*exp(-(((x)/c4).^2))... + a5*exp(-(((x)/c5).^2)); I'll see how that goes...