7 template <
class ScalarType>
8 SpectralClusteringFilter <ScalarType>
11 m_ClassesMembership.clear();
13 m_SpectralVectors.clear();
14 m_InputData.set_size(1,1);
17 m_MaxIterations = 100;
19 m_CMeansAverageType = CMeansFilterType::ApproximateSpherical;
22 m_RelStopCriterion = 1.0e-4;
28 template <
class ScalarType>
39 template <
class ScalarType>
44 if (m_NbClass > m_InputData.rows())
45 throw itk::ExceptionObject(__FILE__,__LINE__,
"More classes than inputs...",ITK_LOCATION);
47 m_DataWeights.resize(m_InputData.rows());
48 std::fill(m_DataWeights.begin(),m_DataWeights.end(),1.0 / m_InputData.rows());
50 this->ComputeSpectralVectors();
52 m_MainFilter.SetNbClass(m_NbClass);
53 m_MainFilter.SetMaxIterations(m_MaxIterations);
54 m_MainFilter.SetRelStopCriterion(m_RelStopCriterion);
55 m_MainFilter.SetMValue(m_MValue);
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);
63 m_MainFilter.Update();
65 unsigned int inputSize = m_InputData.rows();
66 m_ClassesMembership.resize(inputSize);
68 for (
unsigned int i = 0;i < inputSize;++i)
69 m_ClassesMembership[i] = m_MainFilter.GetClassesMembership(i);
71 m_Centroids.resize(m_NbClass);
72 for (
unsigned int i = 0;i < m_NbClass;++i)
73 m_Centroids[i] = m_MainFilter.GetCentroid(i);
76 template <
class ScalarType>
79 ::InitializeSigmaFromDistances()
82 unsigned int inputSize = m_InputData.rows();
83 unsigned int numPts = inputSize*(inputSize + 1)/2 - inputSize - 1;
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);
90 m_SigmaWeighting = std::sqrt(sigmaTmp/numPts);
95 template <
class ScalarType>
98 ::ComputeSpectralVectors()
100 unsigned int inputSize = m_InputData.rows();
101 m_WMatrix.set_size(inputSize,inputSize);
103 m_DValues.resize(inputSize);
105 for (
unsigned int i = 0;i < inputSize;++i)
107 for (
unsigned int j = i+1;j < inputSize;++j)
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);
114 for (
unsigned int i = 0;i < inputSize;++i)
117 for (
unsigned int j = 0;j < inputSize;++j)
118 m_DValues[i] += m_WMatrix(i,j);
120 m_DValues[i] = 1.0/std::sqrt(m_DValues[i]);
123 for (
unsigned int i = 0;i < inputSize;++i)
125 for (
unsigned int j = i+1;j < inputSize;++j)
127 m_WMatrix(i,j) *= m_DValues[i] * m_DValues[j];
128 m_WMatrix(j,i) = m_WMatrix(i,j);
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);
139 m_SpectralVectors.resize(inputSize);
140 m_WorkVec.resize(m_NbClass);
141 for (
unsigned int i = 0;i < inputSize;++i)
144 for (
unsigned int j = 0;j < m_NbClass;++j)
146 m_WorkVec[j] = m_EigVecs.get(inputSize - j - 1,i);
147 tmpSum += m_WorkVec[j]*m_WorkVec[j];
150 tmpSum = std::sqrt(tmpSum);
151 for (
unsigned int j = 0;j < m_NbClass;++j)
152 m_WorkVec[j] /= tmpSum;
154 m_SpectralVectors[i] = m_WorkVec;
158 template <
class ScalarType>
159 std::vector <unsigned int>
161 ::GetClassMembers(
unsigned int i)
163 std::vector <unsigned int> resVal;
165 if (m_ClassesMembership.size() != m_InputData.rows())
168 unsigned int inputSize = m_InputData.rows();
170 for (
unsigned int j = 0;j < inputSize;++j)
172 bool isClassIMax =
true;
173 for (
unsigned int k = 0;k < m_NbClass;++k)
175 if (m_ClassesMembership[j][k] > m_ClassesMembership[j][i])
189 template <
class ScalarType>
192 ::ComputeClustersSpreading()
194 if (m_ClassesMembership.size() != m_InputData.rows())
197 unsigned int inputSize = m_InputData.rows();
198 std::vector <double> classDistance(m_NbClass,0);
200 for (
unsigned int i = 0;i < m_NbClass;++i)
202 double sumWeights = 0;
203 for (
unsigned int j = 0;j < inputSize;++j)
205 sumWeights += m_ClassesMembership[j][i];
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]);
211 classDistance[i] /= sumWeights;
215 for (
unsigned int i = 0;i < m_NbClass;++i)
216 resVal += classDistance[i];
218 return resVal/m_NbClass;
Provides an implementation of spectral clustering, as proposed in A.Y. Ng, M.I. Jordan and Y...
vnl_matrix< double > MatrixType