Path: news.mathworks.com!newsfeed-00.mathworks.com!newscon06.news.prodigy.com!newscon02.news.prodigy.com!prodigy.net!news.glorb.com!news.ks.uiuc.edu!not-for-mail
From: Allen <ajhalldeleteme@gmaildot.com>
Newsgroups: comp.soft-sys.matlab
Subject: Help with 2d or 3d data... (plane fit/lsplane.m)
Date: Tue, 21 Jun 2005 01:46:42 -0500
Organization: Usenet @ UIUC - http://www.killfile.org/usenet/
Lines: 102
Message-ID: <d98d4i$7ao$1@news.ks.uiuc.edu>
NNTP-Posting-Host: cis-mac1.mrl.uiuc.edu
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
X-Trace: news.ks.uiuc.edu 1119336402 7512 130.126.103.251 (21 Jun 2005 06:46:42 GMT)
X-Complaints-To: news+abuse@ks.uiuc.edu
NNTP-Posting-Date: Tue, 21 Jun 2005 06:46:42 +0000 (UTC)
User-Agent: Mozilla Thunderbird 1.0 (Macintosh/20041206)
X-Accept-Language: en-us, en
Xref: news.mathworks.com comp.soft-sys.matlab:285817


Hi everyone,

I'm extremely new to matlab, so I'm very sorry if I'm missing some 
simple things... I've been trying to learn as much as possible... I need 
some help thinking about this current problem.

I have height data taken from an AFM (atomic force microscope), each z 
height (nanometers) is taken at some x,y position (also in nanometers). 
  Currently I have imported the z array (256x256 double), and have also 
created an x and y array(s) corresponding to the positions along x and y 
of each point.  The data can be specified in 3d, in that it has a 
position in x, y, and z.  I would like to plane-fit the data with say, 
the lsplane.m file of I.M. Smith (included below for those who don't 
have it from the repository).  [or other plane-fitting routine you may 
suggest]

I have some confusion about my data and the array required for 
lsplane.m.  Here's the input array information (quoting):  "Array [x y 
z] where x=vector of x-coords, y=vector of y-coords, z=vector of z 
coords, Dimension mx3."

Now, how are my X, Y, Z matrices related to the array Lsplane.m is 
requesting?  I can't seem to merge the data together into a 3d array. 
[in a spread sheet, the x,y would simply be the index for the z data... 
x's size of course is 1x256 (so is Y, but I can transpose it), and Z's 
size is 256x256.

Any help appreciated in dealing with my silly confusion... as I said, 
I'm painfully new to all of this.  :)  Always fun.  I guess I'll be 
learning a lot.  ;)

Thanks very much for any help and suggestions.  I greatly appreciate it!!
-Allen

ps- please feel free to e-mail me regarding the question, but I will 
assume I need to look here so that the answer will be helpful to all who 
come after me.  ;)  I did do a lot of looking in the google archives and 
help files etc., first before asking, but perhaps not enough... [for 
e-mail demunge the address]

---------function: lsplane.m included from repository-------------

   function [x0, a, d, normd] = lsplane(X)
% ---------------------------------------------------------------------
% LSPLANE.M   Least-squares plane (orthogonal distance
%             regression).
%
% Version 1.0
% Last amended   I M Smith 27 May 2002.
% Created        I M Smith 08 Mar 2002
% ---------------------------------------------------------------------
% Input
% X        Array [x y z] where x = vector of x-coordinates,
%          y = vector of y-coordinates and z = vector of
%          z-coordinates.
%          Dimension: m x 3.
%
% Output
% x0       Centroid of the data = point on the best-fit plane.
%          Dimension: 3 x 1.
%
% a        Direction cosines of the normal to the best-fit
%          plane.
%          Dimension: 3 x 1.
%
% <Optional...
% d        Residuals.
%          Dimension: m x 1.
%
% normd    Norm of residual errors.
%          Dimension: 1 x 1.
% ...>
%
% [x0, a <, d, normd >] = lsplane(X)
% ---------------------------------------------------------------------
% check number of data points
   m = size(X, 1);
   if m < 3
     error('At least 3 data points required: ' )
   end
%
% calculate centroid
   x0 = mean(X)';
%
% form matrix A of translated points
   A = [(X(:, 1) - x0(1)) (X(:, 2) - x0(2)) (X(:, 3) - x0(3))];
%
% calculate the SVD of A
   [U, S, V] = svd(A, 0);
%
% find the smallest singular value in S and extract from V the
% corresponding right singular vector
   [s, i] = min(diag(S));
   a = V(:, i);
%
% calculate residual distances, if required
   if nargout > 2
     d = U(:, i)*s;
     normd = norm(d);
   end
% ---------------------------------------------------------------------
% End of LSPLANE.M.