OTB  7.4.0
Orfeo Toolbox
otbStreamingInnerProductVectorImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2020 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_hxx
23 #define otbStreamingInnerProductVectorImageFilter_hxx
25 
26 #include "itkImageRegionIterator.h"
27 #include "itkImageRegionConstIteratorWithIndex.h"
28 #include "itkNumericTraits.h"
29 #include "itkProgressReporter.h"
30 #include "otbMacro.h"
31 
32 namespace otb
33 {
34 
35 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>
56 {
57  switch (output)
58  {
59  case 0:
60  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
61  break;
62  case 1:
63  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
64  break;
65  default:
66  // might as well make an image
67  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
68  break;
69  }
70 }
71 
72 template <class TInputImage>
74 {
75  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
76 }
77 
78 template <class TInputImage>
81 {
82  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
83 }
84 
85 template <class TInputImage>
87 {
88  Superclass::GenerateOutputInformation();
89  if (this->GetInput())
90  {
91  this->GetOutput()->CopyInformation(this->GetInput());
92  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
93 
94  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
95  {
96  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
97  }
98  }
99 }
100 
101 template <class TInputImage>
103 {
104  // This is commented to prevent the streaming of the whole image for the first stream strip
105  // It shall not cause any problem because the output image of this filter is not intended to be used.
106  // InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
107  // this->GraftOutput( image );
108  // Nothing that needs to be allocated for the remaining outputs
109 }
110 
111 template <class TInputImage>
113 {
114  TInputImage* inputPtr = const_cast<TInputImage*>(this->GetInput());
115  inputPtr->UpdateOutputInformation();
116 
117  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
118  {
119  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
120  }
121 
122  unsigned int numberOfThreads = this->GetNumberOfThreads();
123  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
124  // Set the number of training image
125  MatrixType tempMatrix;
126  tempMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages);
127  tempMatrix.fill(0);
128  m_ThreadInnerProduct = ArrayMatrixType(numberOfThreads, tempMatrix);
129 
130  MatrixType initMatrix;
131  initMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages);
132  initMatrix.fill(0);
133  this->GetInnerProductOutput()->Set(initMatrix);
134 }
135 
136 template <class TInputImage>
138 {
139  // Compute Inner product Matrix
140  TInputImage* inputPtr = const_cast<TInputImage*>(this->GetInput());
141  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
142  unsigned int numberOfThreads = this->GetNumberOfThreads();
143  MatrixType innerProduct;
144  innerProduct.set_size(numberOfTrainingImages, numberOfTrainingImages);
145  innerProduct.fill(0);
146 
147  // Concatenate threaded matrix
148  for (unsigned int thread = 0; thread < numberOfThreads; thread++)
149  {
150  innerProduct += m_ThreadInnerProduct[thread];
151  }
152 
153  //---------------------------------------------------------------------
154  // Fill the rest of the inner product matrix and make it symmetric
155  //---------------------------------------------------------------------
156  for (unsigned int band_x = 0; band_x < (numberOfTrainingImages - 1); ++band_x)
157  {
158  for (unsigned int band_y = band_x + 1; band_y < numberOfTrainingImages; ++band_y)
159  {
160  innerProduct[band_x][band_y] = innerProduct[band_y][band_x];
161  } // end band_y loop
162  } // end band_x loop
163 
164  if ((numberOfTrainingImages - 1) != 0)
165  {
166  innerProduct /= (numberOfTrainingImages - 1);
167  }
168  else
169  {
170  innerProduct.fill(0);
171  }
172 
173  // Set the output
174  this->GetInnerProductOutput()->Set(innerProduct);
175 }
176 
177 template <class TInputImage>
178 void PersistentInnerProductVectorImageFilter<TInputImage>::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
179 {
183  InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput());
184  // support progress methods/callbacks
185  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
186  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
188 
189  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
190  if (m_CenterData == true)
191  {
192  it.GoToBegin();
193  // do the work
194  while (!it.IsAtEnd())
195  {
196  PixelType vectorValue = it.Get();
197  double mean(0.);
198  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
199  {
200  mean += static_cast<double>(vectorValue[i]);
201  }
202  mean /= static_cast<double>(vectorValue.GetSize());
203 
204  // Matrix iteration
205  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
206  {
207  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
208  {
209  m_ThreadInnerProduct[threadId][band_x][band_y] +=
210  (static_cast<double>(vectorValue[band_x]) - mean) * (static_cast<double>(vectorValue[band_y]) - mean);
211  } // end: band_y loop
212  } // end: band_x loop
213  ++it;
214  progress.CompletedPixel();
215  } // end: looping through the image
216  }
217  else
218  {
219  it.GoToBegin();
220  // do the work
221  while (!it.IsAtEnd())
222  {
223  PixelType vectorValue = it.Get();
224  // Matrix iteration
225  for (unsigned int band_x = 0; band_x < numberOfTrainingImages; ++band_x)
226  {
227  for (unsigned int band_y = 0; band_y <= band_x; ++band_y)
228  {
229  m_ThreadInnerProduct[threadId][band_x][band_y] += (static_cast<double>(vectorValue[band_x])) * (static_cast<double>(vectorValue[band_y]));
230  } // end: band_y loop
231  } // end: band_x loop
232  ++it;
233  progress.CompletedPixel();
234  } // end: looping through the image
235  }
236 }
237 
238 template <class TImage>
239 void PersistentInnerProductVectorImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
240 {
241  Superclass::PrintSelf(os, indent);
242  os << indent << "m_CenterData: " << m_CenterData << std::endl;
243  os << indent << "InnerProduct: " << this->GetInnerProductOutput()->Get() << std::endl;
244 }
245 
246 } // end namespace otb
247 #endif
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
itk::SimpleDataObjectDecorator< MatrixType > MatrixObjectType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.