How to do a Gauss2D*3 fitting
Dahair
I'm analyzing a images which contains 3 bright disks very close to each other in y direction and I need to analyze the intensity. If these 3 disks are isolated, I can easily do 3 Gauss2D fitting on them. But they are too close to each other that Gauss2D fitting will results in the results yWidth twice as large as xWidth.
So I wondered whether I can do a Gauss2D*3.
I saw in Igor build-in functions there is one Gauss2D*2 which can solve similar condition in which 2 disks are close to each other and the function has a form of:
f(x,y) = z0+(A*exp((-1/(2*(1-(cor1)^2)))*((((x-x1)/xw1)^2)+(((y-y1)/yw1)^2)-((2*cor1*(x-x1)*(y-y1))/(xw1*yw1)))))+ (B*exp((-1/(2*(1-(cor2)^2)))*((((x-x2)/xw2)^2)+(((y-y2)/yw2)^2)-((2*cor2*(x-x2)*(y-y2))/(xw2*yw2))))). There is no term that involves x1 and x2 (or y1 and y2) together. So I wondered if I want to do a Gauss2D*3, can I just add another term ((C*exp((-1/(2*(1-(cor3)^2)))*((((x-x3)/xw3)^2)+(((y-y3)/yw3)^2)-((2*cor3*(x-x3)*(y-y3))/(xw3*yw3)))))) to this function?
But here's another way to do it:
Wave pw
Variable xx, yy
Make/D/FREE/N=7 onepeakcoefs
onepeakcoefs = pw[p]
Variable result = Gauss2D(onepeakcoefs, xx, yy)
onepeakcoefs = pw[p+7]
onepeakcoefs[0] = 0
result += Gauss2D(onepeakcoefs, xx, yy)
onepeakcoefs = pw[p+14]
onepeakcoefs[0] = 0
result += Gauss2D(onepeakcoefs, xx, yy)
return result
end
This one can be used for any number of peaks:
Wave pw
Variable xx, yy
Make/D/N=7/FREE onepeakcoefs
Variable i
Variable npeaks = (numpnts(pw)-1)/6
Variable result = pw[0]
onepeakcoefs[0] = 0
for (i = 0; i < npeaks; i += 1)
onepeakcoefs[1,] = pw[p + 6*i]
result += Gauss2D(onepeakcoefs, xx, yy)
endfor
return result
end
In both cases, the first element of the coefficient wave is a vertical offset for the entire group of peaks.
Then coefficients 1-6 are x1, xw1, y1, yw1, and cor1. Coefficients 7-12 are x2, xw2, y2, yw2 and cor2. Etc...
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
August 27, 2014 at 02:08 pm - Permalink
Thanks John! I'm new in Igor so I'm not pretty sure each line's meaning. In Function Gauss2d3peaks(pw, xx, yy) : FitFunc;
1. Is pw[p] used to store coeffs for all peaks? What's the value of "P", there is no initializing.
2. What's the meaning of this line: Variable result = Gauss2D(onepeakcoefs, xx, yy)? I suppose result doesn't have any value (like 1 or 1.5) here. Does it store the equation of "z0+(A*exp((-1/(2*(1-(cor1)^2)))*((((x-x1)/xw1)^2)+(((y-y1)/yw1)^2)-((2*cor1*(x-x1)*(y-y1))/(xw1*yw1)))))"?
3. So "onepeakcoefs = pw[p+7]
onepeakcoefs[0] = 0
result += Gauss2D(onepeakcoefs, xx, yy)"
and "onepeakcoefs = pw[p+14]
onepeakcoefs[0] = 0
result += Gauss2D(onepeakcoefs, xx, yy)"
is adding the second and third peak part expression to the equation, right? And finally, we can get the whole equation. And about the coeffs, are coeffs stored in pw[] or onepeakcoefs? I think it should be pw[], right? But how many elements in pw[]? Is it 21 with 2 element = 0, or these two 0 are not written into pw[] and pw[] has the size of 19?
Thanks!
August 28, 2014 at 11:37 am - Permalink
Hey John,
Another question. By saying "Then coefficients 1-6 are x1, xw1, y1, yw1, and cor1. Coefficients 7-12 are x2, xw2, y2, yw2 and cor2. Etc...", what the positions of A and B?
Thanks
August 28, 2014 at 12:39 pm - Permalink
The function p represents the current point number in the implied loop in a wave assignment. It's fundamental to understanding how Igor works- perhaps you should do at least the first half of the Guided Tour:
DisplayHelpTopic "Getting Started"
And learn more about wave assignments:
DisplayHelpTopic "Waveform Arithmetic and Assignment"
The function Gauss2D returns the value of a 2D Gaussian peak (that monstrous equation that you wrote out) at given values of X and Y. The code line you are asking about creates a local variable called Result and it puts the value of the first peak into it. The variable is later used with the += operator to add the second and third peaks. That operator does two things at once- it reads the contents of the left hand side, adds the right hand side to it, and stores it back into the left hand side. So this line:
A += B
is the same as
A = A + B
And to understand that you have to remember that in programming the equals sign is the assignment operator, which is quite different from the mathematical equals sign that asserts equality of two expressions.
Yes.
The input wave pw contains coefficients for the entire fit, all three peaks. You can fit only one Y offset coefficient, so pw[0] contains Y0 for the entire fit. The wave onepeakcoefs contains the coefficients for a single peak; the Gauss2D function requires a wave containing coefficients for a single peak. The assignment of zero to onepeakcoefs[0] sets the y0 value for peaks 2 and 3 to zero.
So if you followed the previous bit, you should now realize that there are 19 coefficients for the three-peak fit. There are 6 coefficients for each peak, plus a Y0 coefficient for all three peaks.
Woops. That should read "Then coefficients 1-6 are A, x1, xw1, y1, yw1, and cor1. Coefficients 7-12 are B, x2, xw2, y2, yw2 and cor2. Etc..."
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
August 29, 2014 at 09:45 am - Permalink