OTB  6.1.0
Orfeo Toolbox
otbISRAUnmixingImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbISRAUnmixingImageFilter_txx
22 #define otbISRAUnmixingImageFilter_txx
23 
25 #include <algorithm>
26 
27 namespace otb
28 {
29 
30 namespace Functor
31 {
32 
33 template <class TInput, class TOutput, class TPrecision>
36  : m_OutputSize(0),
37  m_MaxIteration(100)
38 {
39 }
40 
41 template <class TInput, class TOutput, class TPrecision>
44 {
45 }
46 
47 template <class TInput, class TOutput, class TPrecision>
48 unsigned int
51 {
52  return m_OutputSize;
53 }
54 
55 template <class TInput, class TOutput, class TPrecision>
56 bool
58 ::operator != (const Self& itkNotUsed(other)) const
59 {
60  return true;
61 }
62 
63 template <class TInput, class TOutput, class TPrecision>
64 bool
66 ::operator == (const Self& itkNotUsed(other)) const
67 {
68  return false;
69 }
70 
71 template <class TInput, class TOutput, class TPrecision>
72 void
75 {
76  m_U = U;
77  m_OutputSize = m_U.cols();
78  m_Svd.reset( new SVDType(m_U) );
79 }
80 
81 
82 template <class TInput, class TOutput, class TPrecision>
86 {
87  return m_U;
88 }
89 
90 template <class TInput, class TOutput, class TPrecision>
93 ::operator ()(const InputType& in) const
94 {
95  // TODO : support different types between input and output ?
96  VectorType inVector(in.Size());
97  for (unsigned int i = 0; i < in.GetSize(); ++i )
98  {
99  inVector[i] = in[i];
100  }
101 
102  // Initialize with Unconstrained Least Square solution
103  VectorType outVector = m_Svd->solve(inVector);
104 
105  unsigned int nbEndmembers = m_OutputSize;
106  unsigned int nbBands = in.Size();
107 
108  // Apply ISRA iterations
109  for (unsigned int i = 0; i < m_MaxIteration; ++i)
110  {
111 
112  // Use a temporary storage since it is used
113  // inside the iterations
114  VectorType outVectorNew = outVector;
115  for (unsigned int e = 0; e < nbEndmembers; ++e)
116  {
117  PrecisionType numerator = 0;
118  PrecisionType denominator = 0;
119 
120  for (unsigned int b = 0; b < nbBands; ++b)
121  {
122  numerator += in[b] * m_U(b, e);
123 
124  PrecisionType dot = 0;
125  for (unsigned int s = 0; s < nbEndmembers; ++s)
126  {
127  // Use outVector from previous iteration here
128  dot += m_U(b, s) * outVector[s];
129  }
130  denominator += dot * m_U(b, e);
131  }
132 
133  outVectorNew[e] *= (numerator/denominator);
134  }
135 
136  // Prepare for next iteration
137  outVector = outVectorNew;
138  }
139 
140  OutputType out(outVector.size());
141  for (unsigned int i = 0; i < out.GetSize(); ++i )
142  {
143  out[i] = outVector[i];
144  }
145  return out;
146 }
147 
148 }
149 
150 template <class TInputImage, class TOutputImage, class TPrecision>
153 {
154 }
155 
156 template <class TInputImage, class TOutputImage, class TPrecision>
159 {
160 }
161 
162 template <class TInputImage, class TOutputImage, class TPrecision>
163 void
166 {
167  this->GetFunctor().SetEndmembersMatrix(m);
168  this->Modified();
169 }
170 
171 template <class TInputImage, class TOutputImage, class TPrecision>
175 {
176  return this->GetFunctor().GetEndmembersMatrix();
177 }
178 
179 template <class TInputImage, class TOutputImage, class TPrecision>
180 void
182 ::PrintSelf(std::ostream& os, itk::Indent indent) const
183 {
184  Superclass::PrintSelf(os, indent);
185 }
186 
187 } // end namespace
188 
189 #endif
void SetEndmembersMatrix(const MatrixType &m)
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
bool operator==(const ISRAUnmixingFunctor &other) const
const MatrixType & GetEndmembersMatrix() const
OutputType operator()(const InputType &in) const
bool operator!=(const ISRAUnmixingFunctor &other) const
const MatrixType & GetEndmembersMatrix(void) const