Orfeo Toolbox  4.2
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>
57 ::MakeOutput(unsigned int output)
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();
203 
204  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
205  if (m_CenterData == true)
206  {
207  it.GoToBegin();
208  // do the work
209  while (!it.IsAtEnd())
210  {
211  PixelType vectorValue = it.Get();
212  double mean(0.);
213  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
214  {
215  mean += static_cast<double>(vectorValue[i]);
216  }
217  mean /= static_cast<double>(vectorValue.GetSize());
218 
219  // Matrix iteration
220  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
221  {
222  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
223  {
224  m_ThreadInnerProduct[threadId][band_x][band_y] +=
225  (static_cast<double>(vectorValue[band_x]) - mean) * (static_cast<double>(vectorValue[band_y]) - mean);
226  } // end: band_y loop
227  } // end: band_x loop
228  ++it;
229  progress.CompletedPixel();
230  } // end: looping through the image
231  }
232  else
233  {
234  it.GoToBegin();
235  // do the work
236  while (!it.IsAtEnd())
237  {
238  PixelType vectorValue = it.Get();
239  // Matrix iteration
240  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
241  {
242  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
243  {
244  m_ThreadInnerProduct[threadId][band_x][band_y] +=
245  (static_cast<double>(vectorValue[band_x])) * (static_cast<double>(vectorValue[band_y]));
246  } // end: band_y loop
247  } // end: band_x loop
248  ++it;
249  progress.CompletedPixel();
250  } // end: looping through the image
251  }
252 }
253 
254 template <class TImage>
255 void
257 ::PrintSelf(std::ostream& os, itk::Indent indent) const
258 {
259  Superclass::PrintSelf(os, indent);
260  os << indent << "m_CenterData: " << m_CenterData << std::endl;
261  os << indent << "InnerProduct: " << this->GetInnerProductOutput()->Get() << std::endl;
262 
263 }
264 
265 } // end namespace otb
266 #endif

Generated at Sat Oct 11 2014 16:26:18 for Orfeo Toolbox with doxygen 1.8.3.1