ANIMA  4.0
animaSpectralClusteringFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 namespace anima
5 {
6 
7 template <class ScalarType>
8 SpectralClusteringFilter <ScalarType>
10 {
11  m_ClassesMembership.clear();
12  m_Centroids.clear();
13  m_SpectralVectors.clear();
14  m_InputData.set_size(1,1);
15 
16  m_NbClass = 0;
17  m_MaxIterations = 100;
18 
19  m_CMeansAverageType = CMeansFilterType::ApproximateSpherical;
20 
21  m_Verbose = true;
22  m_RelStopCriterion = 1.0e-4;
23  m_MValue = 2;
24 
25  m_SigmaWeighting = 1;
26 }
27 
28 template <class ScalarType>
29 void
31 ::SetInputData(MatrixType &data)
32 {
33  if (data.rows() == 0)
34  return;
35 
36  m_InputData = data;
37 }
38 
39 template <class ScalarType>
40 void
42 ::Update()
43 {
44  if (m_NbClass > m_InputData.rows())
45  throw itk::ExceptionObject(__FILE__,__LINE__,"More classes than inputs...",ITK_LOCATION);
46 
47  m_DataWeights.resize(m_InputData.rows());
48  std::fill(m_DataWeights.begin(),m_DataWeights.end(),1.0 / m_InputData.rows());
49 
50  this->ComputeSpectralVectors();
51 
52  m_MainFilter.SetNbClass(m_NbClass);
53  m_MainFilter.SetMaxIterations(m_MaxIterations);
54  m_MainFilter.SetRelStopCriterion(m_RelStopCriterion);
55  m_MainFilter.SetMValue(m_MValue);
56 
57  m_MainFilter.SetInputData(m_SpectralVectors);
58  m_MainFilter.SetDataWeights(m_DataWeights);
59  m_MainFilter.SetVerbose(m_Verbose);
60  m_MainFilter.SetFlagSpectralClustering(true);
61  m_MainFilter.SetSphericalAverageType(m_CMeansAverageType);
62 
63  m_MainFilter.Update();
64 
65  unsigned int inputSize = m_InputData.rows();
66  m_ClassesMembership.resize(inputSize);
67 
68  for (unsigned int i = 0;i < inputSize;++i)
69  m_ClassesMembership[i] = m_MainFilter.GetClassesMembership(i);
70 
71  m_Centroids.resize(m_NbClass);
72  for (unsigned int i = 0;i < m_NbClass;++i)
73  m_Centroids[i] = m_MainFilter.GetCentroid(i);
74 }
75 
76 template <class ScalarType>
77 void
79 ::InitializeSigmaFromDistances()
80 {
81  double sigmaTmp = 0;
82  unsigned int inputSize = m_InputData.rows();
83  unsigned int numPts = inputSize*(inputSize + 1)/2 - inputSize - 1;
84 
85  for (unsigned int i = 0;i < inputSize;++i)
86  for (unsigned int j = i+1;j < inputSize;++j)
87  sigmaTmp += m_InputData(i,j);
88 
89  if (sigmaTmp > 0)
90  m_SigmaWeighting = std::sqrt(sigmaTmp/numPts);
91  else
92  m_SigmaWeighting = 1;
93 }
94 
95 template <class ScalarType>
96 void
98 ::ComputeSpectralVectors()
99 {
100  unsigned int inputSize = m_InputData.rows();
101  m_WMatrix.set_size(inputSize,inputSize);
102  m_WMatrix.fill(0.0);
103  m_DValues.resize(inputSize);
104 
105  for (unsigned int i = 0;i < inputSize;++i)
106  {
107  for (unsigned int j = i+1;j < inputSize;++j)
108  {
109  m_WMatrix(i,j) = std::exp(- m_InputData(i,j) / (2.0 * m_SigmaWeighting * m_SigmaWeighting));
110  m_WMatrix(j,i) = m_WMatrix(i,j);
111  }
112  }
113 
114  for (unsigned int i = 0;i < inputSize;++i)
115  {
116  m_DValues[i] = 0;
117  for (unsigned int j = 0;j < inputSize;++j)
118  m_DValues[i] += m_WMatrix(i,j);
119 
120  m_DValues[i] = 1.0/std::sqrt(m_DValues[i]);
121  }
122 
123  for (unsigned int i = 0;i < inputSize;++i)
124  {
125  for (unsigned int j = i+1;j < inputSize;++j)
126  {
127  m_WMatrix(i,j) *= m_DValues[i] * m_DValues[j];
128  m_WMatrix(j,i) = m_WMatrix(i,j);
129  }
130  }
131 
132  // Compute eigenvalues and vectors, keep only m_NbClasses largest values and corresponding truncated vectors
133  m_EigenAnalyzer.SetDimension(inputSize);
134  m_EigenAnalyzer.SetOrder(inputSize);
135  m_EigVals.set_size(inputSize);
136  m_EigVecs.set_size(inputSize,inputSize);
137  m_EigenAnalyzer.ComputeEigenValuesAndVectors(m_WMatrix,m_EigVals,m_EigVecs);
138 
139  m_SpectralVectors.resize(inputSize);
140  m_WorkVec.resize(m_NbClass);
141  for (unsigned int i = 0;i < inputSize;++i)
142  {
143  double tmpSum = 0;
144  for (unsigned int j = 0;j < m_NbClass;++j)
145  {
146  m_WorkVec[j] = m_EigVecs.get(inputSize - j - 1,i);
147  tmpSum += m_WorkVec[j]*m_WorkVec[j];
148  }
149 
150  tmpSum = std::sqrt(tmpSum);
151  for (unsigned int j = 0;j < m_NbClass;++j)
152  m_WorkVec[j] /= tmpSum;
153 
154  m_SpectralVectors[i] = m_WorkVec;
155  }
156 }
157 
158 template <class ScalarType>
159 std::vector <unsigned int>
161 ::GetClassMembers(unsigned int i)
162 {
163  std::vector <unsigned int> resVal;
164 
165  if (m_ClassesMembership.size() != m_InputData.rows())
166  return resVal;
167 
168  unsigned int inputSize = m_InputData.rows();
169 
170  for (unsigned int j = 0;j < inputSize;++j)
171  {
172  bool isClassIMax = true;
173  for (unsigned int k = 0;k < m_NbClass;++k)
174  {
175  if (m_ClassesMembership[j][k] > m_ClassesMembership[j][i])
176  {
177  isClassIMax = false;
178  break;
179  }
180  }
181 
182  if (isClassIMax)
183  resVal.push_back(j);
184  }
185 
186  return resVal;
187 }
188 
189 template <class ScalarType>
190 double
192 ::ComputeClustersSpreading()
193 {
194  if (m_ClassesMembership.size() != m_InputData.rows())
195  return -1;
196 
197  unsigned int inputSize = m_InputData.rows();
198  std::vector <double> classDistance(m_NbClass,0);
199 
200  for (unsigned int i = 0;i < m_NbClass;++i)
201  {
202  double sumWeights = 0;
203  for (unsigned int j = 0;j < inputSize;++j)
204  {
205  sumWeights += m_ClassesMembership[j][i];
206 
207  for (unsigned int k = 0;k < m_NbClass;++k)
208  classDistance[i] += m_ClassesMembership[j][i] * std::abs(m_Centroids[i][k] - m_SpectralVectors[j][k]);
209  }
210 
211  classDistance[i] /= sumWeights;
212  }
213 
214  double resVal = 0;
215  for (unsigned int i = 0;i < m_NbClass;++i)
216  resVal += classDistance[i];
217 
218  return resVal/m_NbClass;
219 }
220 
221 } // end namespace anima
Provides an implementation of spectral clustering, as proposed in A.Y. Ng, M.I. Jordan and Y...