| 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.
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
