OTB  5.0.0
Orfeo Toolbox
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>
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 
212 /********END TODO ******/
213  }
214  this->GetResultOutput()->Set(resultMatrix);
215 }
216 
217 template<class TInputImage, class TInputImage2>
218 void
220 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
221 {
225  InputImagePointer input1Ptr = const_cast<TInputImage *>(this->GetFirstInput());
226  InputImagePointer input2Ptr = const_cast<TInputImage2 *>(this->GetSecondInput());
228 
229  // support progress methods/callbacks
230  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
231  input1Ptr->SetRequestedRegion(outputRegionForThread);
232  input2Ptr->SetRequestedRegion(outputRegionForThread);
233  input1Ptr->PropagateRequestedRegion();
234  input1Ptr->UpdateOutputData();
235  input2Ptr->PropagateRequestedRegion();
236  input2Ptr->UpdateOutputData();
237 
238  itk::ImageRegionConstIterator<TInputImage> it1(input1Ptr, outputRegionForThread);
239  itk::ImageRegionConstIterator<TInputImage2> it2(input2Ptr, outputRegionForThread);
240  it1.GoToBegin();
241  it2.GoToBegin();
242 
243  // loop the second image and get one pixel a time
244  while (!it1.IsAtEnd())
245  {
246  PixelType vectorValue1 = it1.Get();
247  PixelType2 vectorValue2 = it2.Get();
248 
249  // Add a first component to vectorValue2 and vectorValue1 filled with ones.
250  if (m_UsePadFirstInput == true)
251  {
252  PixelType vectortemp1(vectorValue1.Size() + 1);
253  vectortemp1[0] = 1;
254  for (unsigned int n = 0; n < vectorValue1.Size(); ++n)
255  {
256  vectortemp1[n + 1] = vectorValue1[n];
257 
258  }
259  vectorValue1.SetSize(vectortemp1.Size());
260  vectorValue1 = vectortemp1;
261  }
262 
263  if (m_UsePadSecondInput == true)
264  {
265  PixelType2 vectortemp2(vectorValue2.Size() + 1);
266  vectortemp2[0] = 1;
267  for (unsigned int m = 0; m < vectorValue2.Size(); m++)
268  {
269  vectortemp2[m + 1] = vectorValue2[m];
270 
271  }
272  vectorValue2.SetSize(vectortemp2.Size());
273  vectorValue2 = vectortemp2;
274  }
275 
276  for (unsigned int i = 0; i < vectorValue1.Size(); ++i)
277  {
278  for (unsigned int j = 0; j < vectorValue2.Size(); ++j)
279  {
280  m_ThreadSum[threadId](i, j) += static_cast<RealType>(vectorValue1[i]) * static_cast<RealType>(vectorValue2[j]);
281  }
282 
283  }
284  ++it1;
285  ++it2;
286  progress.CompletedPixel();
287  }
288 }
289 
290 template<class TInputImage, class TInputImage2>
291 void
293 ::PrintSelf(std::ostream& os, itk::Indent indent) const
294 {
295  Superclass::PrintSelf(os, indent);
296 
297  os << indent << "Result: " << this->GetResultOutput()->Get() << std::endl;
298 }
299 
300 } // end namespace otb
301 #endif
unsigned int Rows() const
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx)
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
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
unsigned int ThreadIdType
DataObject * GetOutput(const DataObjectIdentifierType &key)