OTB  6.1.0
Orfeo Toolbox
otbStreamingInnerProductVectorImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbStreamingInnerProductVectorImageFilter_txx
23 #define otbStreamingInnerProductVectorImageFilter_txx
25 
26 #include "itkImageRegionIterator.h"
28 #include "itkNumericTraits.h"
29 #include "itkProgressReporter.h"
30 #include "otbMacro.h"
31 
32 namespace otb
33 {
34 
35 template<class TInputImage>
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  // allocate the data objects for the outputs which are
46  // just decorators around matrix types
47  typename ImageType::Pointer output1 = static_cast<ImageType*>(this->MakeOutput(0).GetPointer());
48  this->itk::ProcessObject::SetNthOutput(0, output1.GetPointer());
49  typename MatrixObjectType::Pointer output2 = static_cast<MatrixObjectType*>(this->MakeOutput(1).GetPointer());
50  this->itk::ProcessObject::SetNthOutput(1, output2.GetPointer());
51 
52  m_CenterData = true;
53 }
54 
55 template<class TInputImage>
59 {
60  switch (output)
61  {
62  case 0:
63  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
64  break;
65  case 1:
66  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
67  break;
68  default:
69  // might as well make an image
70  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
71  break;
72  }
73 }
74 
75 template<class TInputImage>
79 {
80  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
81 }
82 
83 template<class TInputImage>
87 {
88  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
89 }
90 
91 template<class TInputImage>
92 void
95 {
96  Superclass::GenerateOutputInformation();
97  if (this->GetInput())
98  {
99  this->GetOutput()->CopyInformation(this->GetInput());
100  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
101 
102  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
103  {
104  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
105  }
106  }
107 }
108 
109 template<class TInputImage>
110 void
113 {
114  // This is commented to prevent the streaming of the whole image for the first stream strip
115  // It shall not cause any problem because the output image of this filter is not intended to be used.
116  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
117  //this->GraftOutput( image );
118  // Nothing that needs to be allocated for the remaining outputs
119 }
120 
121 template<class TInputImage>
122 void
125 {
126  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
127  inputPtr->UpdateOutputInformation();
128 
129  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
130  {
131  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
132  }
133 
134  unsigned int numberOfThreads = this->GetNumberOfThreads();
135  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
136  // Set the number of training image
137  MatrixType tempMatrix;
138  tempMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages);
139  tempMatrix.fill(0);
140  m_ThreadInnerProduct = ArrayMatrixType(numberOfThreads, tempMatrix);
141 
142  MatrixType initMatrix;
143  initMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages);
144  initMatrix.fill(0);
145  this->GetInnerProductOutput()->Set(initMatrix);
146 
147 }
148 
149 template<class TInputImage>
150 void
153 {
154  // Compute Inner product Matrix
155  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
156  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
157  unsigned int numberOfThreads = this->GetNumberOfThreads();
158  MatrixType innerProduct;
159  innerProduct.set_size(numberOfTrainingImages, numberOfTrainingImages);
160  innerProduct.fill(0);
161 
162  // Concatenate threaded matrix
163  for (unsigned int thread = 0; thread < numberOfThreads; thread++)
164  {
165  innerProduct += m_ThreadInnerProduct[thread];
166  }
167 
168  //---------------------------------------------------------------------
169  // Fill the rest of the inner product matrix and make it symmetric
170  //---------------------------------------------------------------------
171  for (unsigned int band_x = 0; band_x < (numberOfTrainingImages - 1); ++band_x)
172  {
173  for (unsigned int band_y = band_x + 1; band_y < numberOfTrainingImages; ++band_y)
174  {
175  innerProduct[band_x][band_y] = innerProduct[band_y][band_x];
176  } // end band_y loop
177  } // end band_x loop
178 
179  if ((numberOfTrainingImages - 1) != 0)
180  {
181  innerProduct /= (numberOfTrainingImages - 1);
182  }
183  else
184  {
185  innerProduct.fill(0);
186  }
187 
188  // Set the output
189  this->GetInnerProductOutput()->Set(innerProduct);
190 }
191 
192 template<class TInputImage>
193 void
195 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
196 {
200  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
201  // support progress methods/callbacks
202  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
203  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
205 
206  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
207  if (m_CenterData == true)
208  {
209  it.GoToBegin();
210  // do the work
211  while (!it.IsAtEnd())
212  {
213  PixelType vectorValue = it.Get();
214  double mean(0.);
215  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
216  {
217  mean += static_cast<double>(vectorValue[i]);
218  }
219  mean /= static_cast<double>(vectorValue.GetSize());
220 
221  // Matrix iteration
222  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
223  {
224  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
225  {
226  m_ThreadInnerProduct[threadId][band_x][band_y] +=
227  (static_cast<double>(vectorValue[band_x]) - mean) * (static_cast<double>(vectorValue[band_y]) - mean);
228  } // end: band_y loop
229  } // end: band_x loop
230  ++it;
231  progress.CompletedPixel();
232  } // end: looping through the image
233  }
234  else
235  {
236  it.GoToBegin();
237  // do the work
238  while (!it.IsAtEnd())
239  {
240  PixelType vectorValue = it.Get();
241  // Matrix iteration
242  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
243  {
244  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
245  {
246  m_ThreadInnerProduct[threadId][band_x][band_y] +=
247  (static_cast<double>(vectorValue[band_x])) * (static_cast<double>(vectorValue[band_y]));
248  } // end: band_y loop
249  } // end: band_x loop
250  ++it;
251  progress.CompletedPixel();
252  } // end: looping through the image
253  }
254 }
255 
256 template <class TImage>
257 void
259 ::PrintSelf(std::ostream& os, itk::Indent indent) const
260 {
261  Superclass::PrintSelf(os, indent);
262  os << indent << "m_CenterData: " << m_CenterData << std::endl;
263  os << indent << "InnerProduct: " << this->GetInnerProductOutput()->Get() << std::endl;
264 
265 }
266 
267 } // end namespace otb
268 #endif
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE
unsigned int ThreadIdType
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) ITK_OVERRIDE
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
bool IsAtEnd(void) const
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
DataObject * GetOutput(const DataObjectIdentifierType &key)