1 package org.lcsim.spacegeom;
2
3
4
5
6
7
8
9
10
11
12 public class PrincipalAxesLineFitter
13 {
14
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
24
25
26 public PrincipalAxesLineFitter()
27 {
28 }
29
30
31
32
33
34 public void fit(double[][] points)
35 {
36 int npoints = points[0].length;
37
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
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
83 Eigensystem es = _tensor.eigensystem();
84 es.eigensort();
85
86
87 Matrix eigenvalues = es.eigenvalues();
88
89
90
91
92 Matrix EVEC = es.eigenvectors();
93
94
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
103
104
105
106
107
108
109
110
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
122
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
133
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
144
145
146
147 public Matrix dirCovariance()
148 {
149 return _dircov;
150 }
151
152
153
154
155
156 public double phi()
157 {
158 return _phi;
159 }
160
161
162
163
164
165 public double theta()
166 {
167 return _theta;
168 }
169 }
170
171