Orfeo ToolBox  4.2
Orfeo ToolBox is not a black box
otbStreamingMatrixTransposeMatrixImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12  Some parts of this code are derived from ITK. See ITKCopyright.txt
13  for details.
14 
15 
16  This software is distributed WITHOUT ANY WARRANTY; without even
17  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
18  PURPOSE. See the above copyright notices for more information.
19 
20 =========================================================================*/
21 #ifndef __otbStreamingMatrixTransposeMatrixImageFilter_txx
22 #define __otbStreamingMatrixTransposeMatrixImageFilter_txx
23 
25 
26 #include "itkImageRegionIterator.h"
27 #include "itkNumericTraits.h"
28 #include "itkProgressReporter.h"
29 
30 namespace otb
31 {
32 
33 template<class TInputImage, class TInputImage2>
36 {
37  this->SetNumberOfRequiredInputs(2);
38 
39  // first output is a copy of the image, DataObject created by
40  // superclass
41  //
42  // allocate the data objects for the outputs which are
43  // just decorators around pixel types
44 
45  typename ImageType::Pointer output1 = static_cast<ImageType*>(this->MakeOutput(0).GetPointer());
46  this->itk::ProcessObject::SetNthOutput(0, output1.GetPointer());
47  typename MatrixObjectType::Pointer output2 = static_cast<MatrixObjectType*>(this->MakeOutput(1).GetPointer());
48  this->itk::ProcessObject::SetNthOutput(1, output2.GetPointer());
49 
50  // false means no pad added
51  m_UsePadFirstInput = false;
52  m_UsePadSecondInput = false;
53 
54  // Number of component initialization
55  m_NumberOfComponents1 = 0;
56  m_NumberOfComponents2 = 0;
57 }
58 
59 template<class TInputImage, class TInputImage2>
60 itk::DataObject::Pointer
62 ::MakeOutput(unsigned int output)
63 {
64  switch (output)
65  {
66  case 0:
67  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
68  break;
69  case 1:
70  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
71  break;
72  default:
73  // might as well make an image
74  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
75  break;
76  }
77 
78 }
79 template<class TInputImage, class TInputImage2>
83 {
84  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
85 }
86 
87 template<class TInputImage, class TInputImage2>
91 {
92  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
93 }
94 
95 template<class TInputImage, class TInputImage2>
96 void
99 {
100  Superclass::GenerateInputRequestedRegion();
101 
102  if (this->GetFirstInput() && this->GetSecondInput())
103  {
104  InputImagePointer image = const_cast<typename Superclass::InputImageType *>(this->GetFirstInput());
105  InputImagePointer image2 = const_cast<typename Superclass::InputImageType *>(this->GetSecondInput());
106  image->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
107  image2->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
108  }
109 }
110 template<class TInputImage, class TInputImage2>
111 void
114 {
115  Superclass::GenerateOutputInformation();
116  if (this->GetFirstInput())
117  {
118  this->GetOutput()->CopyInformation(this->GetFirstInput());
119  this->GetOutput()->SetLargestPossibleRegion(this->GetFirstInput()->GetLargestPossibleRegion());
120  }
121 
122  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
123  {
124  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
125  }
126 }
127 
128 template<class TInputImage, class TInputImage2>
129 void
132 {
133  // This is commented to prevent the streaming of the whole image for the first stream strip
134  // It shall not cause any problem because the output image of this filter is not intended to be used.
135  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
136  //this->GraftOutput( image );
137  // Nothing that needs to be allocated for the remaining outputs
138 }
139 
140 template<class TInputImage, class TInputImage2>
141 void
144 {
145 
146  TInputImage * inputPtr1 = const_cast<TInputImage *>(this->GetFirstInput());
147  inputPtr1->UpdateOutputInformation();
148  TInputImage2 * inputPtr2 = const_cast<TInputImage2 *>(this->GetSecondInput());
149  inputPtr2->UpdateOutputInformation();
150 
151  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
152  {
153  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
154  }
155 
156  if (inputPtr1->GetLargestPossibleRegion().GetSize() != inputPtr2->GetLargestPossibleRegion().GetSize())
157  {
158  itkExceptionMacro(<< " Can't multiply the transposed matrix of a "
159  << inputPtr1->GetLargestPossibleRegion().GetSize()
160  << " and a "
161  << inputPtr2->GetLargestPossibleRegion().GetSize()
162  << " matrix ");
163  }
164 
165  m_NumberOfComponents1 = inputPtr1->GetNumberOfComponentsPerPixel();
166  m_NumberOfComponents2 = inputPtr2->GetNumberOfComponentsPerPixel();
167  unsigned int numberOfThreads = this->GetNumberOfThreads();
168 
169  if (m_UsePadFirstInput == true)
170  {
171  m_NumberOfComponents1++;
172  }
173  if (m_UsePadSecondInput == true)
174  {
175  m_NumberOfComponents2++;
176  }
177 
178  MatrixType tempMatrix, initMatrix;
179  tempMatrix.SetSize(m_NumberOfComponents1, m_NumberOfComponents2);
181  m_ThreadSum = ArrayMatrixType(numberOfThreads, tempMatrix);
182 
183  initMatrix.SetSize(m_NumberOfComponents2, m_NumberOfComponents2);
185  this->GetResultOutput()->Set(initMatrix);
186 }
187 
188 template<class TInputImage, class TInputImage2>
189 void
192 {
193  unsigned int numberOfThreads = this->GetNumberOfThreads();
194  MatrixType resultMatrix;
195  resultMatrix.SetSize(m_NumberOfComponents1, m_NumberOfComponents2);
197 
198  for (unsigned int thread = 0; thread < numberOfThreads; thread++)
199  {
204  for (unsigned int i = 0; i < resultMatrix.Rows(); ++i)
205  {
206  for (unsigned int j = 0; j < resultMatrix.Cols(); ++j)
207  {
208  resultMatrix(i, j) += m_ThreadSum[thread](i, j);
209  }
210  }
211  /********END TODO ******/
212  }
213  this->GetResultOutput()->Set(resultMatrix);
214 }
215 
216 template<class TInputImage, class TInputImage2>
217 void
219 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
220 {
224  InputImagePointer input1Ptr = const_cast<TInputImage *>(this->GetFirstInput());
225  InputImagePointer input2Ptr = const_cast<TInputImage2 *>(this->GetSecondInput());
226 
227  // support progress methods/callbacks
228  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
229  input1Ptr->SetRequestedRegion(outputRegionForThread);
230  input2Ptr->SetRequestedRegion(outputRegionForThread);
231  input1Ptr->PropagateRequestedRegion();
232  input1Ptr->UpdateOutputData();
233  input2Ptr->PropagateRequestedRegion();
234  input2Ptr->UpdateOutputData();
235 
236  itk::ImageRegionConstIterator<TInputImage> it1(input1Ptr, outputRegionForThread);
237  itk::ImageRegionConstIterator<TInputImage2> it2(input2Ptr, outputRegionForThread);
238  it1.GoToBegin();
239  it2.GoToBegin();
240 
241  // loop the second image and get one pixel a time
242  while (!it1.IsAtEnd())
243  {
244  PixelType vectorValue1 = it1.Get();
245  PixelType2 vectorValue2 = it2.Get();
246 
247  // Add a first component to vectorValue2 and vectorValue1 filled with ones.
248  if (m_UsePadFirstInput == true)
249  {
250  PixelType vectortemp1(vectorValue1.Size() + 1);
251  vectortemp1[0] = 1;
252  for (unsigned int n = 0; n < vectorValue1.Size(); ++n)
253  {
254  vectortemp1[n + 1] = vectorValue1[n];
255 
256  }
257  vectorValue1.SetSize(vectortemp1.Size());
258  vectorValue1 = vectortemp1;
259  }
260 
261  if (m_UsePadSecondInput == true)
262  {
263  PixelType2 vectortemp2(vectorValue2.Size() + 1);
264  vectortemp2[0] = 1;
265  for (unsigned int m = 0; m < vectorValue2.Size(); m++)
266  {
267  vectortemp2[m + 1] = vectorValue2[m];
268 
269  }
270  vectorValue2.SetSize(vectortemp2.Size());
271  vectorValue2 = vectortemp2;
272  }
273 
274  for (unsigned int i = 0; i < vectorValue1.Size(); ++i)
275  {
276  for (unsigned int j = 0; j < vectorValue2.Size(); ++j)
277  {
278  m_ThreadSum[threadId](i, j) += static_cast<RealType>(vectorValue1[i]) * static_cast<RealType>(vectorValue2[j]);
279  }
280 
281  }
282  ++it1;
283  ++it2;
284  progress.CompletedPixel();
285  }
286 }
287 
288 template<class TInputImage, class TInputImage2>
289 void
291 ::PrintSelf(std::ostream& os, itk::Indent indent) const
292 {
293  Superclass::PrintSelf(os, indent);
294 
295  os << indent << "Result: " << this->GetResultOutput()->Get() << std::endl;
296 }
297 
298 } // end namespace otb
299 #endif
unsigned int Rows() const
unsigned int Cols() const
bool SetSize(unsigned int r, unsigned int c)
virtual void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId)
void Fill(const T &value)
bool IsAtEnd(void) const
PixelType Get(void) const
unsigned int ThreadIdType