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

Generated at Sat Aug 30 2014 16:27:09 for Orfeo Toolbox with doxygen 1.8.3.1