OTB  6.7.0
Orfeo Toolbox
otbStreamingStatisticsMapFromLabelImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2019 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 otbStreamingStatisticsMapFromLabelImageFilter_hxx
23 #define otbStreamingStatisticsMapFromLabelImageFilter_hxx
25 
27 #include "itkImageRegionIterator.h"
28 #include "itkProgressReporter.h"
29 #include "otbMacro.h"
30 #include <cmath>
31 
32 namespace otb
33 {
34 
35 template<class TInputVectorImage, class TLabelImage>
38  : m_UseNoDataValue()
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  typename PixelValueMapObjectType::Pointer output
46  = static_cast<PixelValueMapObjectType*>(this->MakeOutput(1).GetPointer());
48 
49  this->Reset();
50 
51 }
52 
53 template<class TInputVectorImage, class TLabelImage>
57 {
58  return static_cast<itk::DataObject*>(PixelValueMapObjectType::New().GetPointer());
59 }
60 
61 template<class TInputVectorImage, class TLabelImage>
62 void
65 {
66  // Process object is not const-correct so the const_cast is required here
68  const_cast< LabelImageType * >( input ) );
69 
70 }
71 
72 template<class TInputVectorImage, class TLabelImage>
76 {
77  return static_cast< const TLabelImage * >
79 }
80 
81 template<class TInputVectorImage, class TLabelImage>
85 {
86  return m_MeanRadiometricValue;
87 }
88 
89 template<class TInputVectorImage, class TLabelImage>
93 {
94  return m_StDevRadiometricValue;
95 }
96 
97 template<class TInputVectorImage, class TLabelImage>
101 {
102  return m_MinRadiometricValue;
103 }
104 
105 template<class TInputVectorImage, class TLabelImage>
109 {
110  return m_MaxRadiometricValue;
111 }
112 
113 template<class TInputVectorImage, class TLabelImage>
117 {
118  return m_LabelPopulation;
119 }
120 
121 template<class TInputVectorImage, class TLabelImage>
122 void
125 {
126  Superclass::GenerateOutputInformation();
127 
128  if (this->GetInput())
129  {
130  this->GetOutput()->CopyInformation(this->GetInput());
131  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
132 
133  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
134  {
135  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
136  }
137  }
138 }
139 
140 template<class TInputVectorImage, class TLabelImage>
141 void
144 {
145  // This is commented to prevent the streaming of the whole image for the first stream strip
146  // It shall not cause any problem because the output image of this filter is not intended to be used.
147  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
148  //this->GraftOutput( image );
149  // Nothing that needs to be allocated for the remaining outputs
150 }
151 
152 template<class TInputVectorImage, class TLabelImage>
153 void
156  {
157  // Update temporary accumulator
158  AccumulatorMapType outputAcc;
159  auto endAcc = outputAcc.end();
160 
161  for (auto const& threadAccMap: m_AccumulatorMaps)
162  {
163  for(auto const& it: threadAccMap)
164  {
165  auto label = it.first;
166  auto itAcc = outputAcc.find(label);
167  if (itAcc == endAcc)
168  {
169  outputAcc.emplace(label, it.second);
170  }
171  else
172  {
173  itAcc->second.Update(it.second);
174  }
175  }
176  }
177 
178  // Publish output maps
179  for(auto& it: outputAcc)
180  {
181  const LabelPixelType label = it.first;
182  const auto &bandCount = it.second.GetBandCount();
183  const auto &sum = it.second.GetSum();
184  const auto &sqSum = it.second.GetSqSum();
185 
186  // Count
187  m_LabelPopulation[label] = it.second.GetCount();
188 
189  // Mean & stdev
191  RealVectorPixelType std (sqSum);
192  for (unsigned int band = 0 ; band < mean.GetSize() ; band++)
193  {
194  // Number of valid pixels in band
195  auto count = bandCount[band];
196  // Mean
197  mean[band] /= count;
198 
199  // Unbiased standard deviation (not sure unbiased is usefull here)
200  const double variance = (sqSum[band] - (sum[band] * mean[band])) / (count - 1);
201  std[band] = std::sqrt(variance);
202  }
203  m_MeanRadiometricValue[label] = mean;
204  m_StDevRadiometricValue[label] = std;
205 
206  // Min & max
207  m_MinRadiometricValue[label] = it.second.GetMin();
208  m_MaxRadiometricValue[label] = it.second.GetMax();
209  }
210 
211  }
212 
213 template<class TInputVectorImage, class TLabelImage>
214 void
217 {
218  m_AccumulatorMaps.clear();
219 
220  m_MeanRadiometricValue.clear();
221  m_StDevRadiometricValue.clear();
222  m_MinRadiometricValue.clear();
223  m_MaxRadiometricValue.clear();
224  m_LabelPopulation.clear();
225  m_AccumulatorMaps.resize(this->GetNumberOfThreads());
226 
227 }
228 
229 template<class TInputVectorImage, class TLabelImage>
230 void
233 {
234  // The Requested Regions of all the inputs are set to their Largest Possible Regions
236 
237  // Iteration over all the inputs of the current filter (this)
238  for( itk::InputDataObjectIterator it( this ); !it.IsAtEnd(); it++ )
239  {
240  // Check whether the input is an image of the appropriate dimension
241  // dynamic_cast of all the input images as itk::ImageBase objects
242  // in order to pass the if ( input ) test whatever the inputImageType (vectorImage or labelImage)
243  ImageBaseType * input = dynamic_cast< ImageBaseType *>( it.GetInput() );
244 
245  if ( input )
246  {
247  // Use the function object RegionCopier to copy the output region
248  // to the input. The default region copier has default implementations
249  // to handle the cases where the input and output are the same
250  // dimension, the input a higher dimension than the output, and the
251  // input a lower dimension than the output.
252  InputImageRegionType inputRegion;
253  this->CallCopyOutputRegionToInputRegion( inputRegion, this->GetOutput()->GetRequestedRegion() );
254  input->SetRequestedRegion(inputRegion);
255  }
256  }
257 }
258 
259 template<class TInputVectorImage, class TLabelImage>
260 void
262 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId )
263 {
267  InputVectorImagePointer inputPtr = const_cast<TInputVectorImage *>(this->GetInput());
268  LabelImagePointer labelInputPtr = const_cast<TLabelImage *>(this->GetInputLabelImage());
270 
271  itk::ImageRegionConstIterator<TInputVectorImage> inIt(inputPtr, outputRegionForThread);
272  itk::ImageRegionConstIterator<TLabelImage> labelIt(labelInputPtr, outputRegionForThread);
273  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
274 
275  auto &acc = m_AccumulatorMaps[threadId];
276  auto endAcc = acc.end();
277 
278  // do the work
279  for (inIt.GoToBegin(), labelIt.GoToBegin();
280  !inIt.IsAtEnd() && !labelIt.IsAtEnd();
281  ++inIt, ++labelIt)
282  {
283  const auto &value = inIt.Get();
284  auto label = labelIt.Get();
285 
286  // Update the accumulator
287  auto itAcc = acc.find(label);
288  if (itAcc == endAcc)
289  {
290  acc.emplace(label, AccumulatorType(this->GetNoDataValue(), this->GetUseNoDataValue(), value));
291  }
292  else
293  {
294  itAcc->second.Update(value);
295  }
296 
297  progress.CompletedPixel();
298  }
299 }
300 
301 template<class TInputVectorImage, class TLabelImage>
302 void
304 ::PrintSelf(std::ostream& os, itk::Indent indent) const
305 {
306  Superclass::PrintSelf(os, indent);
307 }
308 
309 } // end namespace otb
310 #endif
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
ObjectType * GetPointer() const
virtual void GenerateInputRequestedRegion()
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
unsigned int ThreadIdType
DataObject * GetInput(const DataObjectIdentifierType &key)
std::unordered_map< LabelPixelType, RealVectorPixelType > PixelValueMapType
virtual void SetNthInput(DataObjectPointerArraySizeType num, DataObject *input)
bool IsAtEnd(void) const
unsigned int GetSize(void) const noexcept
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
virtual void SetRequestedRegion(const RegionType &region)