1 package org.lcsim.spacegeom;
2 /**
3 * Class for determining straight line parameters for 3D spacepoints
4 * by solving for the principal axes of the moment of inertia tensor.
5 * This version uses unweighted cells.
6 * Array centroid is a point on the line.
7 * Array dircos contains the direction cosines for the line.
8 *
9 *@version $Id: PrincipalAxesLineFitter.java,v 1.1.1.1 2010/12/01 00:15:57 jeremy Exp $
10 *
11 */
12 public class PrincipalAxesLineFitter
13 {
14 // private boolean _fit = false;
15 private Matrix _tensor = new Matrix(3,3);
16 private double[] _centroid = new double[3];
17 private double[] _dircos = new double[3];
18 private Matrix _dircov = new Matrix(3,3);
19 private double _phi=0.;
20 private double _theta=0.;
21
22 /**
23 * Default Constructor
24 *
25 */
26 public PrincipalAxesLineFitter()
27 {
28 }
29
30 /**
31 * derive the line parameters given the moment of inertia tensor.
32 * @param points the moment of inertia tensor
33 */
34 public void fit(double[][] points)
35 {
36 int npoints = points[0].length;
37 // Calculate the centroid of the N points
38
39 for(int i=0; i<3; ++i)
40 {
41 for(int j=0; j<npoints; ++j)
42 {
43 _centroid[i] += points[i][j];
44 }
45 _centroid[i]/=(double)npoints;
46 }
47
48 // Accumulate the tensor
49
50 double dx=0.;
51 double dy=0.;
52 double dz=0.;
53 double sumXX=0.;
54 double sumXY=0.;
55 double sumXZ=0.;
56 double sumYY=0.;
57 double sumYZ=0.;
58 double sumZZ=0.;
59 for(int i = 0; i<npoints; ++i)
60 {
61 dx = points[0][i] - _centroid[0];
62 dy = points[1][i] - _centroid[1];
63 dz = points[2][i] - _centroid[2];
64
65 sumXX += dx*dx;
66 sumXY += dx*dy;
67 sumXZ += dx*dz;
68 sumYY += dy*dy;
69 sumYZ += dy*dz;
70 sumZZ += dz*dz;
71 }
72 _tensor.set(0,0, + sumYY+sumZZ);
73 _tensor.set(1,0, - sumXY);
74 _tensor.set(2,0, - sumXZ);
75 _tensor.set(0,1, - sumXY);
76 _tensor.set(1,1, + sumXX+sumZZ);
77 _tensor.set(2,1, - sumYZ);
78 _tensor.set(0,2, - sumXZ);
79 _tensor.set(1,2, - sumYZ);
80 _tensor.set(2,2, + sumXX+sumYY);
81
82 // Calculate eigenvalues and eigenvectors
83 Eigensystem es = _tensor.eigensystem();
84 es.eigensort();
85
86 //eigenvalues...
87 Matrix eigenvalues = es.eigenvalues();
88 // System.out.println("Sorted (descending) eigenvalues of K:");
89 // eigenvalues.print(2);
90
91 // eigenvectors...
92 Matrix EVEC = es.eigenvectors();
93 // System.out.println("Eigenvectors of K sorted by descending eigenvalues:");
94 // EVEC.print(6);
95
96 _phi = Math.atan2(-EVEC.at(1,2),-EVEC.at(0,2));
97 _theta = Math.acos(-EVEC.at(2,2));
98
99 _dircos[0] = -EVEC.at(0,2);
100 _dircos[1] = -EVEC.at(1,2);
101 _dircos[2] = -EVEC.at(2,2);
102 // calculate the errors
103 // Strandlie, Fruehwirth, NIM A 480 (2002) 734-740
104 // Mardia, et al. Multivariate Analysis, Academic press, 6th printing, 1997
105 //
106 // Covariance matrix for normalized eigenvector is:
107 //
108 // V_i = (lambda_i/N) Sum(j!=i)[(lambda_j/(lambda_j-lambda_i)**2)gamma_j*gamma_j^T]
109 //
110 // Matrix cov = new Matrix(3,3);
111 for(int i = 0; i<2; ++i)
112 {
113 Matrix tmp = EVEC.column(i).times(EVEC.column(i).transposed());
114 tmp.times(eigenvalues.at(i,0)/(eigenvalues.at(i,0)-eigenvalues.at(2,0)));
115 _dircov.plusEqual(tmp);
116 }
117 _dircov.times(eigenvalues.at(2,0)/npoints);
118 }
119
120 /**
121 * Returns the (x, y, z) centroid of the set of points.
122 * @return The centroid of the points
123 */
124 public double[] centroid()
125 {
126 double[] res = new double[3];
127 System.arraycopy(_centroid, 0, res, 0, 3);
128 return res;
129 }
130
131 /**
132 * Returns the direction cosines of the set of points.
133 * @return the direction cosines
134 */
135 public double[] dircos()
136 {
137 double[] res = new double[3];
138 System.arraycopy(_dircos, 0, res, 0, 3);
139 return res;
140 }
141
142 /**
143 * Returns the uncertainty on the direction cosines.
144 * This has not been tested
145 * @return the uncertainty on the direction cosines
146 */
147 public Matrix dirCovariance()
148 {
149 return _dircov;
150 }
151
152 /**
153 * The phi angle of the centroid
154 * @return the phi angle of the centroid
155 */
156 public double phi()
157 {
158 return _phi;
159 }
160
161 /**
162 * The theta angle of the centroid
163 * @return the theta angle of the centroid
164 */
165 public double theta()
166 {
167 return _theta;
168 }
169 }
170
171