View Javadoc

1   /*
2    * CentralMomentsCalculator.java
3    *
4    * Created on April 5, 2006, 11:25 AM
5    *
6    * $Id: CentralMomentsCalculator.java,v 1.1.1.1 2010/11/30 21:32:00 jeremy Exp $
7    */
8   
9   package org.lcsim.math.moments;
10  
11  import Jama.EigenvalueDecomposition;
12  import Jama.Matrix;
13  
14  /**
15   * calculates rotational and translational invariant moments of spatial distributions
16   * @author Norman Graf
17   */
18  public class CentralMomentsCalculator
19  {
20      // 0
21      private double m000;
22      
23      // 1
24      private double m100, m010, m001;
25      
26      // 2
27      private double m110, m101, m011;
28      private double m200, m020, m002;
29      
30      // centroids
31      private double xc, yc, zc;
32      
33      // invariants
34      private double J1, J2, J3;
35      
36      
37      
38      private Matrix _tensor = new Matrix(3,3);
39      
40      // eigenvalues
41      private double[] _eigenvalues = new double[3];
42      
43      // eigenvectors
44      private Matrix _eigenvectors;
45      
46      // direction cosines for the cluster direction
47      private double[] _dircos = new double[3];
48  
49      
50      /**
51       * Creates a new instance of CentralMomentsCalculator
52       */
53      public CentralMomentsCalculator()
54      {
55      }
56      
57      public void calculateMoments(double[] x, double[] y, double[] z, double[] w)
58      {
59          reset();
60          //TODO make this more efficient.
61          int n = x.length;
62          // TODO check that all arrays are the same size.
63          
64          // calculate centroids
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          // on to the higher moments wrt centroid
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          // normalize
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         // now for eigenvalues, eigenvectors
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 //        System.out.println("eigenvalues: "+_eigenvalues[0]+" "+_eigenvalues[1]+" "+_eigenvalues[2]);
121         _eigenvectors = eig.getV();
122 //        System.out.println("eigenvectors:");
123 //        _eigenvectors.print(4,4);
124         
125         // direction cosines are the eigenvector elements corresponding to the lowest eigenvalue
126         // note that eigenvalues are sorted in ascending order by EigenvalueDecomposition
127         _dircos[0] = -_eigenvectors.get(0,0);
128         _dircos[1] = -_eigenvectors.get(1,0);
129         _dircos[2] = -_eigenvectors.get(2,0);
130 //        System.out.println("dircos: "+_dircos[0]+" "+_dircos[1]+" "+_dircos[2]);
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 }