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