| Re: Plane fit [message #39094] |
Thu, 22 April 2004 22:43  |
MKatz843
Messages: 98
Registered: March 2002
|
Member |
|
|
> I was wondering if there is any function/routine to compute a 2
> dimensional fit (i.e. a surface) that fits a series of points spread
> in 3-dim space.
Here's what I use. My own home-cooked simple plane fit based on least
squares. It couldn't be simpler. Well I suppose without all the
keywords it could, but what fun is that?
; M. Katz 1/26/04
; IDL function to perform a least-squares fit a plane, based on
; Ax + By + C = z
;
; ABC = plane_fit(x, y, z, error=error)
;
function plane_fit, x, y, z, error=error, noerror=noerror,
noshow=noshow
tx2 = total(x^2)
ty2 = total(y^2)
txy = total(x*y)
tx = total(x)
ty = total(y)
N = n_elements(x)
A = [[tx2, txy, tx], $
[txy, ty2, ty], $
[tx, ty, N ]]
b = [total(z*x), total(z*y), total(z)]
out = invert(a) # b
if not keyword_set(noshow) then begin
print, 'Plane Fit: Ax + By + C = z'
print, 'A = ', out(0)
print, 'B = ', out(1)
print, 'C = ', out(2)
endif
if not keyword_set(noerror) then begin
error = stdev(out(0)*x + out(1)*y + out(2) - z)
if not keyword_set(noshow) $
then print, 's = ', error
endif
return, out ;--- [A,B,C]
end
Unless you set the noerror keyword, the program also calculates the
RMS error and returns it as a scalar in the error keyword. The noshow
keyword suppresses printing of the A,B,C values to the screen.
|
|
|
|
| Re: Plane fit [message #39149 is a reply to message #39094] |
Tue, 20 April 2004 08:55  |
Dick Jackson
Messages: 347
Registered: August 1998
|
Senior Member |
|
|
Hi Steve,
"Steve" <Steve.Morris@libero.it> wrote in message
news:606669c3.0404200355.71209718@posting.google.com...
> Hi!
> I was wondering if there is any function/routine to compute a 2
> dimensional fit (i.e. a surface) that fits a series of points spread
> in 3-dim space.
> Cheers,
> S.
This is what I use, based on Craig Markwardt's excellent MPFIT routines
(you'll need at least mpfit.pro and mpfit2dfun.pro from
http://cow.physics.wisc.edu/~craigm/idl/fitting.html):
;----------------------------------------
FUNCTION Plane2DFunction, x, y, parms, dParms
; parms define a plane: [c00, c01, c10] (1st-order polynomial in 2
vbles)
Return, parms[0] + parms[1]*x + parms[2]*y
END
;----------------------------------------
;; Then, if data is a (3, n) array of xyz points...
fitParms = MPFIT2DFUN('Plane2DFunction', $
data[0, *], data[1, *], data[2, *], $
0D, $ ; Error measure, ignored with Weights...
[Mean(data[2, *]), 0D, 0D], $ ; Estimate
Weights=Replicate(1D, nPts), /Quiet)
;; Then, for any x and y values...
fitZ = Plane2DFunction(x, y, fitParms)
;-----------------------------------------
Hope this helps!
Cheers,
--
-Dick
Dick Jackson / dick@d-jackson.com
D-Jackson Software Consulting / http://www.d-jackson.com
Calgary, Alberta, Canada / +1-403-242-7398 / Fax: 241-7392
|
|
|
|
| Re: Plane fit [message #39189 is a reply to message #39149] |
Tue, 27 April 2004 14:18 |
Steve.Morris
Messages: 17
Registered: November 2002
|
Junior Member |
|
|
Hi!
Thanks both, the function plane_fit is really straightforward and work
excellent!
I have tried to get more confident with the powerfull mpfit package,
but I'm still stuck at the first step :(
Here is what I do (where x,y,z have the following format array[n] and
I have no errors for my data ... gedankexperiment ;) !)
FitParms = MPFIT2DFUN('Plane2DFunction',$
x ,y, z,$
start_params=[0.D,0.D,0.D],$
/WEIGHTS, /Quiet)
where the fitting function is defined as follows
FUNCTION Plane2DFunction, x, y, P
Return, P[0] + P[1]*x + P[2]*y
END
If I use the keyword Quiet, and query(help,FitParms) about FitParms,
idl reply me
FITPARMS DOUBLE = NaN
If I drop the keyword Quiet, I have the following message
MPFIT2DFUN: ERROR: must pass parameters in P or PARINFO
and FitParms is still NaN ...
Where do I wrong? Probably I'm doing something really stupid
Thanks!Cheers,
S.
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
