OTB  5.7.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 
232  itk::ImageRegionConstIterator<TInputImage> it1(input1Ptr, outputRegionForThread);
233  itk::ImageRegionConstIterator<TInputImage2> it2(input2Ptr, outputRegionForThread);
234  it1.GoToBegin();
235  it2.GoToBegin();
236 
237  // loop the second image and get one pixel a time
238  while (!it1.IsAtEnd())
239  {
240  PixelType vectorValue1 = it1.Get();
241  PixelType2 vectorValue2 = it2.Get();
242 
243  // Add a first component to vectorValue2 and vectorValue1 filled with ones.
244  if (m_UsePadFirstInput == true)
245  {
246  PixelType vectortemp1(vectorValue1.Size() + 1);
247  vectortemp1[0] = 1;
248  for (unsigned int n = 0; n < vectorValue1.Size(); ++n)
249  {
250  vectortemp1[n + 1] = vectorValue1[n];
251 
252  }
253  vectorValue1.SetSize(vectortemp1.Size());
254  vectorValue1 = vectortemp1;
255  }
256 
257  if (m_UsePadSecondInput == true)
258  {
259  PixelType2 vectortemp2(vectorValue2.Size() + 1);
260  vectortemp2[0] = 1;
261  for (unsigned int m = 0; m < vectorValue2.Size(); m++)
262  {
263  vectortemp2[m + 1] = vectorValue2[m];
264 
265  }
266  vectorValue2.SetSize(vectortemp2.Size());
267  vectorValue2 = vectortemp2;
268  }
269 
270  for (unsigned int i = 0; i < vectorValue1.Size(); ++i)
271  {
272  for (unsigned int j = 0; j < vectorValue2.Size(); ++j)
273  {
274  m_ThreadSum[threadId](i, j) += static_cast<RealType>(vectorValue1[i]) * static_cast<RealType>(vectorValue2[j]);
275  }
276 
277  }
278  ++it1;
279  ++it2;
280  progress.CompletedPixel();
281  }
282 }
283 
284 template<class TInputImage, class TInputImage2>
285 void
287 ::PrintSelf(std::ostream& os, itk::Indent indent) const
288 {
289  Superclass::PrintSelf(os, indent);
290 
291  os << indent << "Result: " << this->GetResultOutput()->Get() << std::endl;
292 }
293 
294 } // end namespace otb
295 #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)