diff --git a/Code/Hyperspectral/otbMDMDNMFImageFilter.h b/Code/Hyperspectral/otbMDMDNMFImageFilter.h index 3005a6de4b990e3ba94ab47353c2cf15f935244b..401a2afdbbd08d68023d58e54dc3d915c7a1c042 100644 --- a/Code/Hyperspectral/otbMDMDNMFImageFilter.h +++ b/Code/Hyperspectral/otbMDMDNMFImageFilter.h @@ -34,46 +34,46 @@ namespace otb * works: * K. S. F.J. Theis and T. Tanaka, First results on uniqueness of sparse non-negative matrix factorisation. - * M. G. A. Huck and J. Blanc-Talon, IEEE TGRS, vol. 48, no. 6,pp. 2590-2602, 2010. + * M. G. A. Huck and J. Blanc-Talon, IEEE TGRS, vol. 48, no. 6, pp. 2590-2602, 2010. * A. Huck and M. Guillaume, in WHISPERS, 2010, Grenoble. - * + * * Let $\matR$ be the matrix of the hyperspectral data, whose $I$ columns are the * spectral pixels and the $L$ rows are the vectorial spectral band * images. The linear mixing model can be written as follow : - * \begin{equation} + * \begin{equation} * \matR=\matA \matS + \matN= \matX + \matN - * \end{equation} + * \end{equation} * The $I$ columns of $\matR$ contain the spectral pixels * and the $I$ columns of $\matS$ hold their respective sets of abundance * fractions. The $J$ rows of $\matS$ are the abundance maps * corresponding to the respective end-members. The $J$ columns of * $\matA$ are the end members spectra, and $\matX$ is the signal * matrix. Both $\matA$ and $\matS$ are unknown. - * + * * The basic NMF formulation is to find two matrices $\hat{\matA}$ and - * $\hat{ \matS}$ such as: - * \begin{equation} - * \matX\simeq \hat{\matA} \hat{\matS} - * \end{equation} + * $\hat{ \matS}$ such as: + * \begin{equation} + * \matX\simeq \hat{\matA} \hat{\matS} + * \end{equation} * NMF based algorithms consider the - * properties of the dual spaces $span^+(\matA')$ and $span^+(\matS)$,in + * properties of the dual spaces $span^+(\matA')$ and $span^+(\matS)$, in * which $span^+(\mathbf m^1 ...\mathbf m^d)=\{\mathbf v=\sum_i \mathbf * m^i\mathbf a_i|\mathbf a\in \matR _+^d\}$. The * positiveness is then a fundamental assumption and is exploited to * restrict the admissible solutions set. - * + * * A common used solution is to minimize the reconstruction quadratic - * error $RQE({\matA},{\matS})=\|\matR-{\matA} {\matS}\|^2_F$. In order to + * error $RQE({\matA}, {\matS})=\|\matR-{\matA} {\matS}\|^2_F$. In order to * satisfy the sum-to-one constraint for hyperspectral data, a * regularization term $STU(\matS)$ can be added to the objective * function. - * + * * A generic expression for the optimized function is $$ * f(\matA,\matS)=\|\matA \matS-\matR\|_{norm}+\sum_i \lambda_i * D_i(\matA) + \sum_j \lambda_j D_j(\matS)$$ in which $\|\matA * \matS-\matR\|_{norm}$ is the connection-to-the-data term, and - * $\lambda_{\{i,j\}}$ are regularization parameters for end members and - * abundances constraints $D_{\{i,j\}}$. + * $\lambda_{\{i, j\}}$ are regularization parameters for end members and + * abundances constraints $D_{\{i, j\}}$. * In \cite{Huck2010a}, they * propose an other regularization term, * $D_A(\matA)=Tr(\matA^T\matA)-\frac{1}{L}Tr\left ( \matA^T \1_{LL}\matA @@ -89,7 +89,7 @@ namespace otb * minimizes the following function $ f(\matA,\matS) =RQE(\matA, * \matS)+\delta.STU(\matS)+\lambda_A D_A(\matA)-\lambda_S D_S(\matS)$, * $STU$ the sum-to-one constraint. - * + * * In the literature, NMF based optimization algorithms are mainly based * on multiplicative rules, or else alternate gradient descent * iterations, or else on alternate least square methods. In MDMD-NMF, the update rules @@ -110,7 +110,7 @@ namespace otb * \delta\cdot\1_{1I}\end{array}\right],\enspace \bar\matA=\left[ * \begin{array}{c} \matA \\ * \delta\cdot\1_{1J}\end{array}\right]\enspace$. - * + * * \ingroup ImageFilters */ template <class TInputImage, class TOutputImage> @@ -198,7 +198,7 @@ private: MDMDNMFImageFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented - static void AddOneRowOfOnes(const MatrixType & m,MatrixType & M); + static void AddOneRowOfOnes(const MatrixType & m, MatrixType & M); static double Criterion(const MatrixType & X, const MatrixType & A, diff --git a/Code/Hyperspectral/otbMDMDNMFImageFilter.txx b/Code/Hyperspectral/otbMDMDNMFImageFilter.txx index 1445246df1187536a8b134c567c93695610d4aaf..3bbd4f0b4282b38156b980e3669bb0c15396854e 100644 --- a/Code/Hyperspectral/otbMDMDNMFImageFilter.txx +++ b/Code/Hyperspectral/otbMDMDNMFImageFilter.txx @@ -56,7 +56,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { this->GetOutput()->SetNumberOfComponentsPerPixel(m_Endmembers.columns()); } - else + else { throw itk::ExceptionObject(__FILE__, __LINE__, "Endmembers matrix columns size required to know the output size", @@ -119,7 +119,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> E1 = Xsu - Asu*S; // Bloc 2 - E2 = S - 1./ ((double) nbEndmembers);// * ones - S; + E2 = S - 1./ ((double) nbEndmembers); // * ones - S; // Computing trace(transpose(A)*A) double trAtA = 0; @@ -147,7 +147,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { for (int j2=0; j2<nbEndmembers; ++j2) { - E3(j1,j2) = sumColsOfA(j1)*sumColsOfA(j2); + E3(j1, j2) = sumColsOfA(j1)*sumColsOfA(j2); } } */ @@ -168,7 +168,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> const double &m_LambdS, MatrixType & gradS) { - // Calculus of: gradS = 2 * Asu' * (Asu*S-Xsu) - lambd * 2 * (S - 1/J*ones(J,I)); + // Calculus of: gradS = 2 * Asu' * (Asu*S-Xsu) - lambd * 2 * (S - 1/J*ones(J, I)); MatrixType Xsu, Asu, ones; Xsu.set_size(X.rows()+1, X.cols()); @@ -192,8 +192,8 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> MatrixType &gradA) { // Compute gradA - // = (A*S-X) * (transpose(S)) + m_LambdD*(A-1/nbBands*ones(L,L)*A) - // = (A*S-X) * (transpose(S)) + m_LambdD*A- m_LambdD*/nbBands*ones(L,L)*A) + // = (A*S-X) * (transpose(S)) + m_LambdD*(A-1/nbBands*ones(L, L)*A) + // = (A*S-X) * (transpose(S)) + m_LambdD*A- m_LambdD*/nbBands*ones(L, L)*A) MatrixType onesA; VectorType sumColulmnsOfA; @@ -310,8 +310,8 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { for (unsigned int j=0; j<M.cols(); ++j) { - if (M(i,j)<0) - M(i,j) = 0; + if (M(i, j)<0) + M(i, j) = 0; } } } @@ -328,7 +328,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { for (unsigned int j=0; j<M.cols(); ++j) { - M(i,j) = M1(i,j) * M2(i,j); + M(i, j) = M1(i, j) * M2(i, j); } } return M; @@ -435,7 +435,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> MatrixType A(this->m_Endmembers); MatrixType S; - S.set_size(nbEndmembers,nbPixels); + S.set_size(nbEndmembers, nbPixels); S.fill(1.); //std::cout << "S " << S.cols() << std::endl; //----------- Declaration of useful variables -----------// @@ -477,19 +477,19 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> std::cout << "Iteration = " << counter << std::endl; std::cout << "Criterion = " << Criterion(X, A, S, m_Delt, m_LambdS, m_LambdD) << std::endl; std::cout << "statGradS = " << gradS.fro_norm() << std::endl; - std::cout << "gradS(0,0) = " << gradS(0,0) << std::endl; + std::cout << "gradS(0, 0) = " << gradS(0, 0) << std::endl; std::cout << "alphS = " << alphS << std::endl; std::cout << "normS = " << S.fro_norm() << std::endl; - std::cout << "S(0,0) = " << S(0,0) << std::endl; + std::cout << "S(0, 0) = " << S(0, 0) << std::endl; } - ProjGradOneStep(X, A, gradS, sig, bet, m_Delt,m_LambdS, m_LambdD, S, alphS, false); + ProjGradOneStep(X, A, gradS, sig, bet, m_Delt, m_LambdS, m_LambdD, S, alphS, false); if (counter%divisorParam == 0) { std::cout << "alphS = " << alphS << std::endl; std::cout << "normS = " << S.fro_norm() << std::endl; - std::cout << "S(0,0) = " << S(0,0) << std::endl; + std::cout << "S(0, 0) = " << S(0, 0) << std::endl; } //---------------- Update A -----------------// @@ -499,7 +499,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> if (counter%divisorParam == 0) { - std::cout << "gradA(0,0) = " << gradA(0,0) << std::endl; + std::cout << "gradA(0, 0) = " << gradA(0, 0) << std::endl; } @@ -507,7 +507,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { std::cout << "alphA = " << alphA << std::endl; std::cout << "normA = " << A.fro_norm() << std::endl; - std::cout << "A(0,0) = " << A(0,0) << std::endl; + std::cout << "A(0, 0) = " << A(0, 0) << std::endl; } ProjGradOneStep(X, S, gradA, sig, bet, m_Delt, m_LambdS, m_LambdD, A, alphA, true); @@ -515,7 +515,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { std::cout << "alphA = " << alphA << std::endl; std::cout << "normA = " << A.fro_norm() << std::endl; - std::cout << "A(0,0) = " << A(0,0) << std::endl; + std::cout << "A(0, 0) = " << A(0, 0) << std::endl; } //------------ crit evaluation --------------// @@ -523,7 +523,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> critA = Adiff.absolute_value_max(); Sdiff = Sold - S; critS = Sdiff.absolute_value_max(); - crit = std::max(critA,critS); + crit = std::max(critA, critS); if (counter%divisorParam == 0) { @@ -538,7 +538,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> } //--- Putting the rows of in the bands of the output vector image ---// - //TODO + //TODO // Could be impoved choosing an imageList for the abundance maps // and a vector list for the endmember spectra (columns of A). @@ -553,7 +553,7 @@ MDMDNMFImageFilter<TInputImage, TOutputImage> { for (unsigned int j=0; j<nbEndmembers; ++j) { - vectorPixel.SetElement(j, S(j,i)); + vectorPixel.SetElement(j, S(j, i)); } outputIt.Set(vectorPixel); ++i; diff --git a/Code/Simulation/otbReduceSpectralResponse.h b/Code/Simulation/otbReduceSpectralResponse.h index b50bd1583b0f027945af820e137e5251956b20c5..1922668fae208b760fddad9ebb97887da33f228e 100644 --- a/Code/Simulation/otbReduceSpectralResponse.h +++ b/Code/Simulation/otbReduceSpectralResponse.h @@ -75,20 +75,20 @@ public: typedef typename std::vector<ValuePrecisionType> ReduceSpectralResponseVectorType; /** Standard macros */ itkNewMacro(Self) - ; itkTypeMacro(ReduceSpectralResponse, DataObject) - ; +; itkTypeMacro(ReduceSpectralResponse, DataObject) +; itkGetConstObjectMacro(InputSatRSR, InputRSRType) - ; itkSetObjectMacro(InputSatRSR, InputRSRType) - ; +; itkSetObjectMacro(InputSatRSR, InputRSRType) +; itkGetConstObjectMacro(InputSpectralResponse, InputSpectralResponseType) - ; itkSetObjectMacro(InputSpectralResponse, InputSpectralResponseType) - ; +; itkSetObjectMacro(InputSpectralResponse, InputSpectralResponseType) +; /** The GetReduceResponse method gives the output. The first value in the pair is the central wavelength of the band (see SpectralResponse). */ itkGetObjectMacro(ReduceResponse, InputSpectralResponseType) - ; +; /** Clear the vector data */ virtual bool Clear(); @@ -119,7 +119,7 @@ protected: virtual ~ReduceSpectralResponse() { } - ; +; /** PrintSelf method */ //void PrintSelf(std::ostream& os, itk::Indent indent) const; diff --git a/Code/Simulation/otbReduceSpectralResponse.txx b/Code/Simulation/otbReduceSpectralResponse.txx index 9c7676f8c89d12f516de91d1428671cd73bbe290..7058060d0b877694b7a2c704623ba8ba134688ea 100644 --- a/Code/Simulation/otbReduceSpectralResponse.txx +++ b/Code/Simulation/otbReduceSpectralResponse.txx @@ -64,7 +64,6 @@ ReduceSpectralResponse<TSpectralResponse , TRSR> PrecisionType lambda2; - typename VectorPairType::const_iterator it; VectorPairType pairs = (m_InputSatRSR->GetRSR())[numBand]->GetResponse(); it = pairs.begin(); diff --git a/Code/Simulation/otbSatelliteRSR.h b/Code/Simulation/otbSatelliteRSR.h index d3dfbb17ce3ce412e17c3e6f82220d5657d61d1c..a0c01707b4392a30e71c5390ec8806226b3a59c6 100644 --- a/Code/Simulation/otbSatelliteRSR.h +++ b/Code/Simulation/otbSatelliteRSR.h @@ -54,17 +54,17 @@ public: /** Standard macros */ itkNewMacro(Self) - ; itkTypeMacro(SatelliteRSR, DataObject) - ; +; itkTypeMacro(SatelliteRSR, DataObject) +; /** Set the number of band of the satellite from an ASCII file * Need to parse first all the file to determine the number of columns */ itkGetConstMacro(NbBands, unsigned int) - ; itkSetMacro(NbBands, unsigned int) - ; +; itkSetMacro(NbBands, unsigned int) +; itkSetMacro(SortBands, bool) - ; +; /** Template parameters typedef */ typedef TPrecision PrecisionType; @@ -141,7 +141,7 @@ protected: virtual ~SatelliteRSR() { } - ; +; bool m_SortBands; diff --git a/Code/Simulation/otbSpectralResponse.h b/Code/Simulation/otbSpectralResponse.h index 500a4f89e2c271311a88eef8ca7ffd6689f62c26..64856b52be112364e5be86cae9829e5ef5b878fd 100644 --- a/Code/Simulation/otbSpectralResponse.h +++ b/Code/Simulation/otbSpectralResponse.h @@ -119,7 +119,6 @@ public: inline ValuePrecisionType operator()(const PrecisionType & lambda); - /** Operator for comparing Pair Lambda/Response * Pairs are ordered by wavelength */ @@ -162,7 +161,7 @@ protected: virtual ~SpectralResponse() { } - ; + ; /** PrintSelf method */ //void PrintSelf(std::ostream& os, itk::Indent indent) const; diff --git a/Code/Simulation/otbSpectralResponse.txx b/Code/Simulation/otbSpectralResponse.txx index 1b8c1d01dcbb78c86d76bf69b0ba1e41025ee29b..9307dde59739234338b064f9fc696da749378000 100644 --- a/Code/Simulation/otbSpectralResponse.txx +++ b/Code/Simulation/otbSpectralResponse.txx @@ -88,7 +88,6 @@ unsigned int SpectralResponse<TPrecision, TValuePrecision>::Size() const } - template<class TPrecision, class TValuePrecision> void SpectralResponse<TPrecision, TValuePrecision>::SetPosGuessMin(const PrecisionType & lambda) { @@ -114,7 +113,6 @@ void SpectralResponse<TPrecision, TValuePrecision>::SetPosGuessMin(const Precisi } - template<class TPrecision, class TValuePrecision> inline typename SpectralResponse<TPrecision, TValuePrecision>::ValuePrecisionType SpectralResponse<TPrecision, TValuePrecision>::operator()(const PrecisionType & lambda) diff --git a/Testing/Code/Hyperspectral/otbMDMDNMFImageFilter.cxx b/Testing/Code/Hyperspectral/otbMDMDNMFImageFilter.cxx index c33ff78c4e9074c70510b289aad66455d1c3b55d..f5a7efeedd70cb91b7f6808296a174836f1575bf 100644 --- a/Testing/Code/Hyperspectral/otbMDMDNMFImageFilter.cxx +++ b/Testing/Code/Hyperspectral/otbMDMDNMFImageFilter.cxx @@ -56,10 +56,10 @@ int otbMDMDNMFImageFilterTest(int argc, char * argv[]) MDMDNMFImageFilterType::MatrixType A; A.set_size(readerImage->GetOutput()->GetNumberOfComponentsPerPixel(), 5); A.fill(100.); - A.set_column(1,200.); - A.set_column(2,300.); - A.set_column(3,400.); - A.set_column(4,500.); + A.set_column(1, 200.); + A.set_column(2, 300.); + A.set_column(3, 400.); + A.set_column(4, 500.); unmixer->SetEndmembersMatrix(A); unmixer->SetMaxIter(maxIter);