_________________________

		LEAST SQUARES PLANE (LSQP)
		_________________________

Contents	* DEFINITION OF PLANE
		* DEFINITION OF LEAST SQUARES PLANE
		* REFERENCES
		* PROGRAM (FORTRAN)



DEFINITION OF PLANE
___________________

A Plane is often defined by

	AX + BY + CZ + D = 0			(1)
or sometimes as
	AX + BY + CZ = D			(2)

The distance, Ri, of a point i from the plane (2) is

	Ri = AXi + BYi + CZi - D

assuming that this is a normalised plane i.e. A^2 + B^2 + C^2 = 1,
otherwise Ri = (AXi + BYi + CZi - D) / sqrt(A^2 + B^2 + C^2).

Its worth noting that Ri can be +ve or -ve. To reverse the sign of
Ri we can always multiply equation (2) by -1. The usual convention
is such that

	Distance of a point to the plane should be -ve if the plane
	lies between the point and the origin.

One may have to multiply (A,B,C,D) by -1 to adhere to above convention.


DEFINITION OF LEAST SQUARES PLANE
_________________________________

Given a set of N points (a point i has coordinates xi,yi,zi)
we wish to minimize, Q,
      
       Q = S [ (Ri)^2 ]

Where S (symbol for sigma sign) is the sum over all N points
(i.e. i = 1 to N). Using (2)
       
       Q = S [ (AXi + BYi + CZi - D)^2 ]	(3)

From (3)

	dQ/dA = S [ 2Xi(AXi + BYi + CZi - D) ] = 0  (4)
	dQ/dB = S [ 2Yi(AXi + BYi + CZi - D) ] = 0  (5)
	dQ/dC = S [ 2Zi(AXi + BYi + CZi - D) ] = 0  (6)
	dQ/dD = S [  -2(AXi + BYi + CZi - D) ] = 0  (7)

Equation (7), noting that S[D] = ND, can be rewritten as

	    D = AXo + BYo + CZo   		(8)
Where
	   Xo = S[Xi]/N , Yo = S[Yi]/N , Zo = S[Zi]/N
 
Equation (8) shows that the best plane passes through the
centre of mass. Substracting the centre of mass from each
point and substituting into (4-6) we get a set of simultanious
equations which can be written as

	   WP = 0				(9)

Where (we have dropped the index i for brevity)

               | S[(X-Xo)^2]      S[(X-Xo)(Y-Yo)]  S[(X-Xo)(Z-Zo)] |
           W = | S[(X-Xo)(Y-Yo)]  S[(Y-Yo)^2]      S[(Y-Yo)(Z-Zo)] |
               | S[(X-Xo)(Z-Zo)]  S[(Y-Yo)(Z-Zo)]  S[(Z-Zo)^2]     |

	       |A|
	   P = |B|
	       |C|

We need to solve equation (9) to obtain the LSQP. Equation (9)
can always have the trivial solution A=B=C=0. To avoid getting
the trivial solution we need to impose conditions on the coefficients
of the plane. The most common condition is

	   A^2 + B^2 + C^2 = 1		         (10)

With the condition (10), it can be shown that the solution to (9)
becomes an Eigenvalue problems i.e.

	   WE = VE				 (11)

Where V are the eigenvalues & E the eigenvectors. Equation (11) will
return three eigenvalues and a 3x3 matrix (eigenvectors). The eigenvalues
will be the sum of squares of the distances and the eigenvectors will
provide three sets (a set is a column of the matrix E) of A,B,C values.
One then simply needs to choose the set of A,B,C associated with the
smallest eigenvalue. It is of interest to note:

	(a) the three sets of A,B,C contained in E are orthogonal to
            each other. The three sets represent the best, intermediate
	    and worst planes respectively

	(b) If the eigenvalues are similar/same (within an error
	    margin) then this suggests that the data coordinates
	    are degenerate/symmetric

	(c) One set of A,B,C also represents the best Least Squares
	    Line


In summary, to calculate the LSQP we need to obtain the eigenvalues
and eigenvectors of W.


REFERENCES
__________

1. Pearson, Karl
   On lines & planes of closest fit to systems of points in space
   Philosophical Magazine pp 559-572 S. 6. Vol. 2 No. 11 Nov 1901 


2. Schomaker,V, Waser, J, Marsh, R.E. & Bergman, G.
   Acta Cryst (1959) 12, 600

3. Blow, DM
   To fit a plane to a set of points by least squares
   Acta Cryst (1960) 13, 168

Refs 2 & 3 do not refer to ref 1, and I can only assume
that they were unaware of it, and so end up re-inventing
the wheel.


PROGRAM
_______


Below is a fortran routine for computing the best plane. It should
be simple to convert to C, if you have your own eigenvalue routines.


c--------------------------------------------------------
	subroutine lsqplane(n,x,eigenvalue,eigenvector)
c--------------------------------------------------------
c	Dr Suhail A Islam
c	Structural Bioinformatics Group
c	Biochemistry Building
c	Department of Biological Sciences
c	Imperial College, London SW7 2AY, U.K.
c	email: suhail.islam@ic.ac.uk
c	www:   http://www.sbg.bio.ic.ac.uk
c	Tel:   +44-(0)20-7594-5713
c	Fax:   +44-(0)20-7594-5789

c	First version based on method of Blow - Kings College 1982
c	Adapt to use eigenvalues - ICRF 1989

c  	INPUT
c	=====
c	n = number of atoms/points
c	x = x,y,z [x(1,i),x(2,i),x(3,i) respectively for point i]
c	w = weight of atom i (NOT USED HERE, but easy to ammend)

c 	OUTPUT
c	======
c	eigenvalue = contains the sum of distances squared from plane 
c	eigenvector = columns contain direction cosines

c  	NOTES
c    	=====
c	The best plane is defined by the smallest eigenvalue, say i,
c	where the RMS distance is given by
c		RMS = sqrt(eigenvalue(i) / n)
c	The best plane (for the smallest eigenvalue, i) is given by
c		A.x + B.y + C.z = D
c		A = eigenvector(1,i)
c		B = eigenvector(2,i)
c		C = eigenvector(3,i)
c		D = centroid(1)*eigenvector(1,i) +
c	    	    centroid(2)*eigenvector(2,i) +
c		    centroid(3)*eigenvector(3,i)
c	Distance of a point (x,y,z) from this plane is
c		Distance = A.x + B.y + C.z - D
c	CONVENTION: Distance should be -VE if plane is between
c	point & origin. One may have to multiply (A,B,C,D) by -1
c	to adhere to above convention.
c	The remaining 2 eigenvalues will represent an intermediate &
c	worst least squares planes through the centroid. These are
c	useful as they are orthogonal can be used to display e.g.
c	inertial ellipsoids. Alos worth noting that the eigenvectors
c	also contain the least squares lines.
c	YOU MUST CHECK ALL 3 EIGENVALUES IN CASE EIGENVALUES ARE
c	SIMILAR/SAME i.e. COORDINATES ARE DEGENERATE/SYMMETRIC
c	Also worth noting:
c ---	TRANSFORM THIS TO THE REAL n-body SPACE AS FOLLOWS:
c	x' = xo * eigenvector(1,1) + yo * eigenvector(2,1) +
c	     zo * eigenvector(3,1)
c	y' = xo * eigenvector(1,2) + yo * eigenvector(2,2) +
c	     zo * eigenvector(3,2)
c	z' = xo * eigenvector(1,3) + yo * eigenvector(2,3) +
c	     zo * eigenvector(3,3)

       	integer n
	real x(3,n)
	real eigenvalue(3),eigenvector(3,3)
	integer iopt,nrot,i,j
	real centroid(3),tensor(3,3)
	real constant,a,b,c
	real almn(3),blmn(3),xyzo(3)

c ---	INITALISE DATA

	do i = 1 , 3
	centroid(i) = 0.0
	eigenvalue(i) = 0.0
	end do
	do i = 1 , 3
	do j = 1 , 3
	eigenvector(i,j) = 0.0
	tensor(i,j) = 0.0
	end do
	end do

c ---	CALCULATE CENTROID OF THE N POINTS

	do i = 1 , 3
	do j = 1 , n
	centroid(i) = centroid(i) + x(i,j)
	end do
	centroid(i) = centroid(i) / float(n)
	end do

c ---	FORM TENSOR

	iopt = 0
	do i = 1 , n
	dx = x(1,i) - centroid(1)
	dy = x(2,i) - centroid(2)
	dz = x(3,i) - centroid(3)
c ---	TENSOR FOR LSQ PLANE
	if(iopt.le.0) then
	tensor(1,1) = tensor(1,1) + dx*dx
	tensor(2,1) = tensor(2,1) + dx*dy
	tensor(3,1) = tensor(3,1) + dx*dz
	tensor(1,2) = tensor(1,2) + dx*dy
	tensor(2,2) = tensor(2,2) + dy*dy
	tensor(3,2) = tensor(3,2) + dy*dz
	tensor(1,3) = tensor(1,3) + dx*dz
	tensor(2,3) = tensor(2,3) + dy*dz
	tensor(3,3) = tensor(3,3) + dz*dz
	else
c ---	INERTIAL MASS TENSOR FOR INERTIAL ELLIPSOIDS
	tensor(1,1) = tensor(1,1) + dy*dy+dz*dz
	tensor(2,1) = tensor(2,1) - dx*dy
	tensor(3,1) = tensor(3,1) - dx*dz
	tensor(1,2) = tensor(1,2) - dx*dy
	tensor(2,2) = tensor(2,2) + dx*dx+dz*dz
	tensor(3,2) = tensor(3,2) - dy*dz
	tensor(1,3) = tensor(1,3) - dx*dz
	tensor(2,3) = tensor(2,3) - dy*dz
	tensor(3,3) = tensor(3,3) + dx*dx+dy*dy
	end if
	end do                 

c ---	GET EIGENVALUES & EIGENVECTORS

	call jacobi(tensor,3,3,eigenvalue,eigenvector,nrot,ierr)
	if(ierr.ne.0) return

c ---	SORT EIGENVALUES FROM LARGEST TO SMALLEST
c ---	THIS IS NOT NECESSARY, JUST CONVENIENT

	call eigsrt(eigenvalue,eigenvector,3,3)

	return
	end
c ---------------------------------------------------------------
c ---------------------------------------------------------------
c ---------------------------------------------------------------
c	ROUTINES BELOW ARE FROM NUMERICAL RECIPIES (PRESS ET AL)
c	NOT SURE ABOUT COPYRIGHT OF PASSING THIS ON !
	subroutine jacobi(a,n,np,d,v,nrot,ierr)
	parameter (nmax=100)
	dimension a(np,np),d(np),v(np,np),b(nmax),z(nmax)
	ierr = 0
	do ip=1,n
        do iq=1,n
	v(ip,iq)=0.
	end do
        v(ip,ip)=1.
	end do
	do ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.
	end do
	nrot=0
	do 24 i=1,50
	sm=0.
	do ip=1,n-1
	do iq=ip+1,n
	sm=sm+abs(a(ip,iq))
	end do
	end do
        if(sm.eq.0.)return
        if(i.lt.4)then
	tresh=0.2*sm/n**2
        else
	tresh=0.
        endif
        do 22 ip=1,n-1
	do 21 iq=ip+1,n
	g=100.*abs(a(ip,iq))
	if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip)))
     *	.and.(abs(d(iq))+g.eq.abs(d(iq))))then
	a(ip,iq)=0.
	else if(abs(a(ip,iq)).gt.tresh)then
	h=d(iq)-d(ip)
	if(abs(h)+g.eq.abs(h))then
	t=a(ip,iq)/h
	else
	theta=0.5*h/a(ip,iq)
	t=1./(abs(theta)+sqrt(1.+theta**2))
	if(theta.lt.0.)t=-t
	endif
	c=1./sqrt(1.+t**2)
	s=t*c
	tau=s/(1.+c)
	h=t*a(ip,iq)
	z(ip)=z(ip)-h
	z(iq)=z(iq)+h
	d(ip)=d(ip)-h
	d(iq)=d(iq)+h
	a(ip,iq)=0.
	do 16 j=1,ip-1
	g=a(j,ip)
	h=a(j,iq)
	a(j,ip)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
16	continue
	do 17 j=ip+1,iq-1
        g=a(ip,j)
	h=a(j,iq)
	a(ip,j)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
17	continue
	do j=iq+1,n
	g=a(ip,j)
	h=a(iq,j)
	a(ip,j)=g-s*(h+g*tau)
	a(iq,j)=h+s*(g-h*tau)
	end do
	do j=1,n
	g=v(j,ip)
	h=v(j,iq)
	v(j,ip)=g-s*(h+g*tau)
	v(j,iq)=h+s*(g-h*tau)
	end do
	nrot=nrot+1
	endif
21      continue
22      continue
        do 23 ip=1,n
	b(ip)=b(ip)+z(ip)
	d(ip)=b(ip)
	z(ip)=0.
23      continue
24	continue
c	ierr = 1
	return
	end
c ---------------------------------------------------------------
	subroutine eigsrt(d,v,n,np)
	dimension d(np),v(np,np)
	do 13 i=1,n-1
	k=i
	p=d(i)
	do 11 j=i+1,n
	if(d(j).ge.p)then
	k=j
	p=d(j)
	endif
11      continue
        if(k.ne.i)then
	d(k)=d(i)
	d(i)=p
	do 12 j=1,n
	p=v(j,i)
	v(j,i)=v(j,k)
	v(j,k)=p
12	continue
	endif
13	continue
	return
	end
c ======================== END OF PROGRAM ==========================