OTB  5.0.0
Orfeo Toolbox
otbStreamingInnerProductVectorImageFilter.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 __otbStreamingInnerProductVectorImageFilter_txx
22 #define __otbStreamingInnerProductVectorImageFilter_txx
24 
25 #include "itkImageRegionIterator.h"
27 #include "itkNumericTraits.h"
28 #include "itkProgressReporter.h"
29 #include "otbMacro.h"
30 
31 namespace otb
32 {
33 
34 template<class TInputImage>
37 {
38  // first output is a copy of the image, DataObject created by
39  // superclass
40  //
41  // allocate the data objects for the outputs which are
42  // just decorators around pixel types
43 
44  // allocate the data objects for the outputs which are
45  // just decorators around matrix types
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  m_CenterData = true;
52 }
53 
54 template<class TInputImage>
58 {
59  switch (output)
60  {
61  case 0:
62  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
63  break;
64  case 1:
65  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
66  break;
67  default:
68  // might as well make an image
69  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
70  break;
71  }
72 }
73 
74 template<class TInputImage>
78 {
79  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
80 }
81 
82 template<class TInputImage>
86 {
87  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
88 }
89 
90 template<class TInputImage>
91 void
94 {
95  Superclass::GenerateOutputInformation();
96  if (this->GetInput())
97  {
98  this->GetOutput()->CopyInformation(this->GetInput());
99  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
100 
101  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
102  {
103  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
104  }
105  }
106 }
107 
108 template<class TInputImage>
109 void
112 {
113  // This is commented to prevent the streaming of the whole image for the first stream strip
114  // It shall not cause any problem because the output image of this filter is not intended to be used.
115  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
116  //this->GraftOutput( image );
117  // Nothing that needs to be allocated for the remaining outputs
118 }
119 
120 template<class TInputImage>
121 void
124 {
125  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
126  inputPtr->UpdateOutputInformation();
127 
128  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
129  {
130  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
131  }
132 
133  unsigned int numberOfThreads = this->GetNumberOfThreads();
134  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
135  // Set the number of training image
136  MatrixType tempMatrix;
137  tempMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages);
138  tempMatrix.fill(0);
139  m_ThreadInnerProduct = ArrayMatrixType(numberOfThreads, tempMatrix);
140 
141  MatrixType initMatrix;
142  initMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages);
143  initMatrix.fill(0);
144  this->GetInnerProductOutput()->Set(initMatrix);
145 
146 }
147 
148 template<class TInputImage>
149 void
152 {
153  // Compute Inner product Matrix
154  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
155  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
156  unsigned int numberOfThreads = this->GetNumberOfThreads();
157  MatrixType innerProduct;
158  innerProduct.set_size(numberOfTrainingImages, numberOfTrainingImages);
159  innerProduct.fill(0);
160 
161  // Concatenate threaded matrix
162  for (unsigned int thread = 0; thread < numberOfThreads; thread++)
163  {
164  innerProduct += m_ThreadInnerProduct[thread];
165  }
166 
167  //---------------------------------------------------------------------
168  // Fill the rest of the inner product matrix and make it symmetric
169  //---------------------------------------------------------------------
170  for (unsigned int band_x = 0; band_x < (numberOfTrainingImages - 1); ++band_x)
171  {
172  for (unsigned int band_y = band_x + 1; band_y < numberOfTrainingImages; ++band_y)
173  {
174  innerProduct[band_x][band_y] = innerProduct[band_y][band_x];
175  } // end band_y loop
176  } // end band_x loop
177 
178  if ((numberOfTrainingImages - 1) != 0)
179  {
180  innerProduct /= (numberOfTrainingImages - 1);
181  }
182  else
183  {
184  innerProduct.fill(0);
185  }
186 
187  // Set the output
188  this->GetInnerProductOutput()->Set(innerProduct);
189 }
190 
191 template<class TInputImage>
192 void
194 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
195 {
199  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
200  // support progress methods/callbacks
201  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
202  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
204 
205  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
206  if (m_CenterData == true)
207  {
208  it.GoToBegin();
209  // do the work
210  while (!it.IsAtEnd())
211  {
212  PixelType vectorValue = it.Get();
213  double mean(0.);
214  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
215  {
216  mean += static_cast<double>(vectorValue[i]);
217  }
218  mean /= static_cast<double>(vectorValue.GetSize());
219 
220  // Matrix iteration
221  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
222  {
223  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
224  {
225  m_ThreadInnerProduct[threadId][band_x][band_y] +=
226  (static_cast<double>(vectorValue[band_x]) - mean) * (static_cast<double>(vectorValue[band_y]) - mean);
227  } // end: band_y loop
228  } // end: band_x loop
229  ++it;
230  progress.CompletedPixel();
231  } // end: looping through the image
232  }
233  else
234  {
235  it.GoToBegin();
236  // do the work
237  while (!it.IsAtEnd())
238  {
239  PixelType vectorValue = it.Get();
240  // Matrix iteration
241  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
242  {
243  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
244  {
245  m_ThreadInnerProduct[threadId][band_x][band_y] +=
246  (static_cast<double>(vectorValue[band_x])) * (static_cast<double>(vectorValue[band_y]));
247  } // end: band_y loop
248  } // end: band_x loop
249  ++it;
250  progress.CompletedPixel();
251  } // end: looping through the image
252  }
253 }
254 
255 template <class TImage>
256 void
258 ::PrintSelf(std::ostream& os, itk::Indent indent) const
259 {
260  Superclass::PrintSelf(os, indent);
261  os << indent << "m_CenterData: " << m_CenterData << std::endl;
262  os << indent << "InnerProduct: " << this->GetInnerProductOutput()->Get() << std::endl;
263 
264 }
265 
266 } // end namespace otb
267 #endif
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId)
virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx)
virtual void PrintSelf(std::ostream &os, itk::Indent indent) const
bool IsAtEnd(void) const
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
unsigned int ThreadIdType
DataObject * GetOutput(const DataObjectIdentifierType &key)