1
2
3
4
5
6
7
8
9 package org.lcsim.math.moments;
10
11 import Jama.EigenvalueDecomposition;
12 import Jama.Matrix;
13
14
15
16
17
18 public class CentralMomentsCalculator
19 {
20
21 private double m000;
22
23
24 private double m100, m010, m001;
25
26
27 private double m110, m101, m011;
28 private double m200, m020, m002;
29
30
31 private double xc, yc, zc;
32
33
34 private double J1, J2, J3;
35
36
37
38 private Matrix _tensor = new Matrix(3,3);
39
40
41 private double[] _eigenvalues = new double[3];
42
43
44 private Matrix _eigenvectors;
45
46
47 private double[] _dircos = new double[3];
48
49
50
51
52
53 public CentralMomentsCalculator()
54 {
55 }
56
57 public void calculateMoments(double[] x, double[] y, double[] z, double[] w)
58 {
59 reset();
60
61 int n = x.length;
62
63
64
65 for(int i=0; i<n; ++i)
66 {
67 m000 += w[i];
68 m100 += x[i]*w[i];
69 m010 += y[i]*w[i];
70 m001 += z[i]*w[i];
71 }
72
73 xc = m100/m000;
74 yc = m010/m000;
75 zc = m001/m000;
76
77
78 double xa, ya, za;
79 for(int i=0; i<n; ++i)
80 {
81 xa = x[i]-xc;
82 ya = y[i]-yc;
83 za = z[i]-zc;
84
85 m110 += xa*ya*w[i];
86 m101 += xa*za*w[i];
87 m011 += ya*za*w[i];
88
89 m200 += xa*xa*w[i];
90 m020 += ya*ya*w[i];
91 m002 += za*za*w[i];
92 }
93
94
95 m110 /= m000;
96 m101 /= m000;
97 m011 /= m000;
98
99 m200 /= m000;
100 m020 /= m000;
101 m002 /= m000;
102
103 J1 = m200 + m020 + m002;
104 J2 = m200*m020 + m200*m002 + m020*m002 - m110*m110 - m101*m101 - m011*m011;
105 J3 = m200*m020*m002 + 2.*m110*m101*m011 - m002*m110*m110 - m020*m101*m101 - m200*m011*m011;
106
107
108 _tensor.set(0, 0, m020 + m002);
109 _tensor.set(1,0, - m110);
110 _tensor.set(2,0, - m101);
111 _tensor.set(0,1, - m110);
112 _tensor.set(1,1, + m200 + m002);
113 _tensor.set(2,1, - m011);
114 _tensor.set(0,2, - m101);
115 _tensor.set(1,2, - m011);
116 _tensor.set(2,2, + m200 + m020);
117
118 EigenvalueDecomposition eig = _tensor.eig();
119 _eigenvalues = eig.getRealEigenvalues();
120
121 _eigenvectors = eig.getV();
122
123
124
125
126
127 _dircos[0] = -_eigenvectors.get(0,0);
128 _dircos[1] = -_eigenvectors.get(1,0);
129 _dircos[2] = -_eigenvectors.get(2,0);
130
131
132 }
133
134 public double[] eigenvalues()
135 {
136 return _eigenvalues;
137 }
138
139 public double[] directionCosines()
140 {
141 return _dircos;
142 }
143
144 public double[] centroid()
145 {
146 return new double[] {xc, yc, zc};
147 }
148
149 public double[] centroidVariance()
150 {
151 return new double[] {m200, m020, m002, m110, m101, m011};
152 }
153
154 public double[] invariants()
155 {
156
157 return new double[] {J1, J2, J3};
158 }
159
160 void reset()
161 {
162 m000 = 0.;
163 m100 = 0.;
164 m010 = 0.;
165 m001 = 0.;
166 m110 = 0.;
167 m101 = 0.;
168 m011 = 0.;
169 m200 = 0.;
170 m020 = 0.;
171 m002 = 0.;
172 }
173
174 }