OTB  5.11.0
Orfeo Toolbox
otbStreamingMatrixTransposeMatrixImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbStreamingMatrixTransposeMatrixImageFilter_txx
23 #define otbStreamingMatrixTransposeMatrixImageFilter_txx
24 
26 
27 #include "itkImageRegionIterator.h"
28 #include "itkNumericTraits.h"
29 #include "itkProgressReporter.h"
30 
31 namespace otb
32 {
33 
34 template<class TInputImage, class TInputImage2>
37 {
38  this->SetNumberOfRequiredInputs(2);
39 
40  // first output is a copy of the image, DataObject created by
41  // superclass
42  //
43  // allocate the data objects for the outputs which are
44  // just decorators around pixel types
45 
46  typename ImageType::Pointer output1 = static_cast<ImageType*>(this->MakeOutput(0).GetPointer());
47  this->itk::ProcessObject::SetNthOutput(0, output1.GetPointer());
48  typename MatrixObjectType::Pointer output2 = static_cast<MatrixObjectType*>(this->MakeOutput(1).GetPointer());
49  this->itk::ProcessObject::SetNthOutput(1, output2.GetPointer());
50 
51  // false means no pad added
52  m_UsePadFirstInput = false;
53  m_UsePadSecondInput = false;
54 
55  // Number of component initialization
56  m_NumberOfComponents1 = 0;
57  m_NumberOfComponents2 = 0;
58 }
59 
60 template<class TInputImage, class TInputImage2>
64 {
65  switch (output)
66  {
67  case 0:
68  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
69  break;
70  case 1:
71  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
72  break;
73  default:
74  // might as well make an image
75  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
76  break;
77  }
78 
79 }
80 template<class TInputImage, class TInputImage2>
84 {
85  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
86 }
87 
88 template<class TInputImage, class TInputImage2>
92 {
93  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
94 }
95 
96 template<class TInputImage, class TInputImage2>
97 void
100 {
101  Superclass::GenerateInputRequestedRegion();
102 
103  if (this->GetFirstInput() && this->GetSecondInput())
104  {
105  InputImagePointer image = const_cast<typename Superclass::InputImageType *>(this->GetFirstInput());
106  InputImagePointer image2 = const_cast<typename Superclass::InputImageType *>(this->GetSecondInput());
107  image->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
108  image2->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
109  }
110 }
111 template<class TInputImage, class TInputImage2>
112 void
115 {
116  Superclass::GenerateOutputInformation();
117  if (this->GetFirstInput())
118  {
119  this->GetOutput()->CopyInformation(this->GetFirstInput());
120  this->GetOutput()->SetLargestPossibleRegion(this->GetFirstInput()->GetLargestPossibleRegion());
121  }
122 
123  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
124  {
125  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
126  }
127 }
128 
129 template<class TInputImage, class TInputImage2>
130 void
133 {
134  // This is commented to prevent the streaming of the whole image for the first stream strip
135  // It shall not cause any problem because the output image of this filter is not intended to be used.
136  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
137  //this->GraftOutput( image );
138  // Nothing that needs to be allocated for the remaining outputs
139 }
140 
141 template<class TInputImage, class TInputImage2>
142 void
145 {
146 
147  TInputImage * inputPtr1 = const_cast<TInputImage *>(this->GetFirstInput());
148  inputPtr1->UpdateOutputInformation();
149  TInputImage2 * inputPtr2 = const_cast<TInputImage2 *>(this->GetSecondInput());
150  inputPtr2->UpdateOutputInformation();
151 
152  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
153  {
154  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
155  }
156 
157  if (inputPtr1->GetLargestPossibleRegion().GetSize() != inputPtr2->GetLargestPossibleRegion().GetSize())
158  {
159  itkExceptionMacro(<< " Can't multiply the transposed matrix of a "
160  << inputPtr1->GetLargestPossibleRegion().GetSize()
161  << " and a "
162  << inputPtr2->GetLargestPossibleRegion().GetSize()
163  << " matrix ");
164  }
165 
166  m_NumberOfComponents1 = inputPtr1->GetNumberOfComponentsPerPixel();
167  m_NumberOfComponents2 = inputPtr2->GetNumberOfComponentsPerPixel();
168  unsigned int numberOfThreads = this->GetNumberOfThreads();
169 
170  if (m_UsePadFirstInput == true)
171  {
172  m_NumberOfComponents1++;
173  }
174  if (m_UsePadSecondInput == true)
175  {
176  m_NumberOfComponents2++;
177  }
178 
179  MatrixType tempMatrix, initMatrix;
180  tempMatrix.SetSize(m_NumberOfComponents1, m_NumberOfComponents2);
182  m_ThreadSum = ArrayMatrixType(numberOfThreads, tempMatrix);
183 
184  initMatrix.SetSize(m_NumberOfComponents2, m_NumberOfComponents2);
186  this->GetResultOutput()->Set(initMatrix);
187 }
188 
189 template<class TInputImage, class TInputImage2>
190 void
193 {
194  unsigned int numberOfThreads = this->GetNumberOfThreads();
195  MatrixType resultMatrix;
196  resultMatrix.SetSize(m_NumberOfComponents1, m_NumberOfComponents2);
198 
199  for (unsigned int thread = 0; thread < numberOfThreads; thread++)
200  {
205  for (unsigned int i = 0; i < resultMatrix.Rows(); ++i)
206  {
207  for (unsigned int j = 0; j < resultMatrix.Cols(); ++j)
208  {
209  resultMatrix(i, j) += m_ThreadSum[thread](i, j);
210  }
211  }
212 
213 /********END TODO ******/
214  }
215  this->GetResultOutput()->Set(resultMatrix);
216 }
217 
218 template<class TInputImage, class TInputImage2>
219 void
221 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
222 {
226  InputImagePointer input1Ptr = const_cast<TInputImage *>(this->GetFirstInput());
227  InputImagePointer input2Ptr = const_cast<TInputImage2 *>(this->GetSecondInput());
229 
230  // support progress methods/callbacks
231  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
232 
233  itk::ImageRegionConstIterator<TInputImage> it1(input1Ptr, outputRegionForThread);
234  itk::ImageRegionConstIterator<TInputImage2> it2(input2Ptr, outputRegionForThread);
235  it1.GoToBegin();
236  it2.GoToBegin();
237 
238  // loop the second image and get one pixel a time
239  while (!it1.IsAtEnd())
240  {
241  PixelType vectorValue1 = it1.Get();
242  PixelType2 vectorValue2 = it2.Get();
243 
244  // Add a first component to vectorValue2 and vectorValue1 filled with ones.
245  if (m_UsePadFirstInput == true)
246  {
247  PixelType vectortemp1(vectorValue1.Size() + 1);
248  vectortemp1[0] = 1;
249  for (unsigned int n = 0; n < vectorValue1.Size(); ++n)
250  {
251  vectortemp1[n + 1] = vectorValue1[n];
252 
253  }
254  vectorValue1.SetSize(vectortemp1.Size());
255  vectorValue1 = vectortemp1;
256  }
257 
258  if (m_UsePadSecondInput == true)
259  {
260  PixelType2 vectortemp2(vectorValue2.Size() + 1);
261  vectortemp2[0] = 1;
262  for (unsigned int m = 0; m < vectorValue2.Size(); m++)
263  {
264  vectortemp2[m + 1] = vectorValue2[m];
265 
266  }
267  vectorValue2.SetSize(vectortemp2.Size());
268  vectorValue2 = vectortemp2;
269  }
270 
271  for (unsigned int i = 0; i < vectorValue1.Size(); ++i)
272  {
273  for (unsigned int j = 0; j < vectorValue2.Size(); ++j)
274  {
275  m_ThreadSum[threadId](i, j) += static_cast<RealType>(vectorValue1[i]) * static_cast<RealType>(vectorValue2[j]);
276  }
277 
278  }
279  ++it1;
280  ++it2;
281  progress.CompletedPixel();
282  }
283 }
284 
285 template<class TInputImage, class TInputImage2>
286 void
288 ::PrintSelf(std::ostream& os, itk::Indent indent) const
289 {
290  Superclass::PrintSelf(os, indent);
291 
292  os << indent << "Result: " << this->GetResultOutput()->Get() << std::endl;
293 }
294 
295 } // end namespace otb
296 #endif
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE
unsigned int Rows() const
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) ITK_OVERRIDE
unsigned int Cols() const
bool SetSize(unsigned int r, unsigned int c)
unsigned int ThreadIdType
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
void Fill(const T &value)
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
DataObject * GetOutput(const DataObjectIdentifierType &key)