OTB  6.7.0
Orfeo Toolbox
otbLabelImageSmallRegionMergingFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 
22 #ifndef otbLabelImageSmallRegionMergingFilter_hxx
23 #define otbLabelImageSmallRegionMergingFilter_hxx
24 
27 #include "itkProgressReporter.h"
28 
29 namespace otb
30 {
31 template <class TInputLabelImage >
34 {
35 }
36 
37 template <class TInputLabelImage >
38 void
40 ::SetLabelPopulation( LabelPopulationType const & labelPopulation )
41 {
42  m_LabelPopulation = labelPopulation;
43 
44  // Initialize m_CorrespondingMap to the identity (i.e. m[label] = label)
45  for (auto label : m_LabelPopulation)
46  {
47  m_LUT[label.first] = label.first;
48  }
49 }
50 
51 template <class TInputLabelImage >
53 ::LabelPopulationType const &
56 {
57  return m_LabelPopulation;
58 }
59 
60 
61 template <class TInputLabelImage >
62 void
64 ::SetLabelStatistic( LabelStatisticType const & labelStatistic )
65 {
66  m_LabelStatistic = labelStatistic;
67 }
68 
69 template <class TInputLabelImage >
71 ::LabelStatisticType const &
74 {
75  return m_LabelStatistic;
76 }
77 
78 template <class TInputLabelImage >
80 ::LUTType const &
82 ::GetLUT() const
83 {
84  return m_LUT;
85 }
86 
87 template <class TInputLabelImage >
88 void
91 {
92  m_NeighboursMapsTmp.clear();
93  m_NeighboursMapsTmp.resize( this->GetNumberOfThreads() );
94 }
95 
96 template <class TInputLabelImage >
97 void
100 {
101  NeigboursMapType neighboursMap;
102 
103  // Merge the neighbours maps from all threads
104  for( unsigned int threadId = 0; threadId < this->GetNumberOfThreads();
105  threadId++)
106  {
107  for (auto const & neighbours : m_NeighboursMapsTmp[threadId])
108  {
109  neighboursMap[ neighbours.first ].insert
110  ( neighbours.second.begin(), neighbours.second.end() );
111  }
112  }
113 
114  // For each label of the label map, find the "closest" connected label,
115  // according to the euclidian distance between the corresponding
116  // m_labelStatistic elements.
117  for (auto const & neighbours : neighboursMap)
118  {
119  double proximity = itk::NumericTraits<double>::max();
120  InputLabelType label = neighbours.first;
121  InputLabelType closestNeighbour = label;
122  auto const& statsLabel = m_LabelStatistic[ label ];
123 
124  for (auto const & neighbour : neighbours.second)
125  {
126  auto const & statsNeighbour = m_LabelStatistic[ neighbour ];
127  assert( statsLabel.Size() == statsNeighbour.Size() );
128 
129  double distance = (statsLabel - statsNeighbour).GetSquaredNorm();
130  if (distance < proximity)
131  {
132  proximity = distance;
133  closestNeighbour = neighbour;
134  }
135  }
136 
137  auto curLabelLUT = FindCorrespondingLabel(label);
138  auto adjLabelLUT = FindCorrespondingLabel(closestNeighbour);
139 
140  // Keep the smallest label (this prevents infinite loop in the LUT
141  // (like LUT[i]=j and LUT[j]=i)
142  if(curLabelLUT < adjLabelLUT)
143  {
144  m_LUT[adjLabelLUT] = curLabelLUT;
145  }
146  else
147  {
148  m_LUT[curLabelLUT] = adjLabelLUT;
149  }
150  }
151 
152  // Update the LUT
153  for (auto & label : m_LUT)
154  {
155  label.second = FindCorrespondingLabel( label.first );
156  }
157 
158  // Update statistics : for each newly merged segments, sum the population, and
159  // recompute the mean.
160  for (auto const & label : m_LUT)
161  {
162  if ((label.second != label.first) && (m_LabelPopulation[label.first]!=0))
163  {
164  // Cache values to reduce number of lookups
165  auto const & populationFirst = m_LabelPopulation[label.first];
166  auto const & populationSecond = m_LabelPopulation[label.second];
167  auto const & statisticFirst = m_LabelStatistic[label.first];
168  auto const & statisticSecond = m_LabelStatistic[label.second];
169 
170  m_LabelStatistic[ label.second ] =
171  ( (statisticFirst * populationFirst)
172  + (statisticSecond * populationSecond) )
173  / (populationFirst + populationSecond);
174 
175  m_LabelPopulation[ label.second ] += populationFirst;
176 
177  // Do not use this label anymore
178  m_LabelPopulation[ label.first ] = 0;
179  }
180  }
181 }
182 
183 template <class TInputLabelImage >
185 ::InputLabelType
188  < TInputLabelImage >::InputLabelType label)
189 {
190  auto correspondingLabel = m_LUT[label];
191  while (label != correspondingLabel)
192  {
193  label = correspondingLabel;
194  correspondingLabel = m_LUT[correspondingLabel];
195  }
196  return correspondingLabel;
197 }
198 
199 template <class TInputLabelImage >
200 void
203 {
204  // call the superclass' implementation of this method
205  Superclass::GenerateInputRequestedRegion();
206 
207  // get pointers to the input
208  auto inputPtr = const_cast<TInputLabelImage *>(this->GetInput());
209 
210 
211  // get a copy of the input requested region
212  auto inputRequestedRegion = inputPtr->GetRequestedRegion();
213 
214  // pad the input requested region by the operator radius
215  inputRequestedRegion.PadByRadius(1);
216  // crop the input requested region at the input's largest possible region
217  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
218  {
219  inputPtr->SetRequestedRegion(inputRequestedRegion);
220  return;
221  }
222  else
223  {
224  // Couldn't crop the region (requested region is outside the largest
225  // possible region). Throw an exception.
226 
227  // store what we tried to request (prior to trying to crop)
228  inputPtr->SetRequestedRegion(inputRequestedRegion);
229 
230  // build an exception
231  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
232  e.SetLocation(ITK_LOCATION);
233  e.SetDescription("Requested region is (at least partially) outside the "
234  "largest possible region.");
235  e.SetDataObject(inputPtr);
236  throw e;
237  }
238 
239 }
240 
241 template <class TInputLabelImage >
242 void
244 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
245  itk::ThreadIdType threadId )
246 {
248  using NeighborhoodIteratorType =
250 
251  typename NeighborhoodIteratorType::RadiusType radius;
252  radius.Fill(1);
253 
254  auto labelImage = this->GetInput();
255 
256  IteratorType it(labelImage, outputRegionForThread);
257  NeighborhoodIteratorType itN(radius, labelImage, outputRegionForThread);
258 
259  // 4 connected Neighborhood (top, bottom, left and right)
260  typename IteratorType::OffsetType top = {{0,-1}};
261  itN.ActivateOffset(top);
262  typename IteratorType::OffsetType bottom = {{0,1}};
263  itN.ActivateOffset(bottom);
264  typename IteratorType::OffsetType right = {{1,0}};
265  itN.ActivateOffset(right);
266  typename IteratorType::OffsetType left = {{-1,0}};
267  itN.ActivateOffset(left);
268 
269  for (it.GoToBegin(); ! it.IsAtEnd(); ++it, ++itN)
270  {
271  assert( !itN.IsAtEnd() );
272  int currentLabel = m_LUT[ it.Get() ];
273 
274  if ( m_LabelPopulation[currentLabel] == m_Size )
275  {
276  for (auto ci = itN.Begin() ; !ci.IsAtEnd(); ci++)
277  {
278  int neighbourLabel = m_LUT[ ci.Get() ];
279  if (neighbourLabel != currentLabel)
280  m_NeighboursMapsTmp[threadId][ currentLabel ].insert(neighbourLabel);
281  }
282  }
283  }
284 }
285 
286 template <class TInputLabelImage >
287 void
289 ::PrintSelf(std::ostream& os, itk::Indent indent) const
290 {
291  Superclass::PrintSelf(os, indent);
292 }
293 
294 template <class TInputLabelImage >
297 {
298  m_SmallRegionMergingFilter = LabelImageSmallRegionMergingFilterType::New();
299 }
300 
301 template <class TInputLabelImage >
302 void
304 ::SetInputLabelImage( const TInputLabelImage * labelImage )
305 {
306  m_SmallRegionMergingFilter->GetFilter()->SetInput( labelImage );
307 }
308 
309 template <class TInputLabelImage >
310 void
312 ::SetLabelPopulation( LabelPopulationType const & labelPopulation )
313 {
314  m_SmallRegionMergingFilter->GetFilter()
315  ->SetLabelPopulation(labelPopulation);
316 }
317 
318 template <class TInputLabelImage >
320 ::LabelPopulationType const &
323 {
324  return m_SmallRegionMergingFilter->GetFilter()->GetLabelPopulation();
325 }
326 
327 template <class TInputLabelImage >
328 void
330 ::SetLabelStatistic( LabelStatisticType const & labelStatistic )
331 {
332  m_SmallRegionMergingFilter->GetFilter()->SetLabelStatistic(labelStatistic);
333 }
334 
335 template <class TInputLabelImage >
337 ::LabelStatisticType const &
340 {
341  return m_SmallRegionMergingFilter->GetFilter()->GetLabelStatistic();
342 }
343 
344 template <class TInputLabelImage >
346 ::LUTType const &
348 ::GetLUT() const
349 {
350  return m_SmallRegionMergingFilter->GetFilter()->GetLUT();
351 }
352 
353 template <class TInputLabelImage>
354 void
357 {
358  this->GenerateData();
359 }
360 
361 template <class TInputLabelImage >
362 void
365 {
366  this->SetProgress(0.0);
367 
368  // Update the filter for all sizes.
369  for (unsigned int size = 1; size < m_MinSize; size++)
370  {
371  m_SmallRegionMergingFilter->GetFilter()->SetSize( size) ;
372  m_SmallRegionMergingFilter->Update();
373  this->UpdateProgress(static_cast<double>(size+1)/m_MinSize);
374  }
375 }
376 
377 
378 } // end namespace otb
379 
380 #endif
void SetLabelStatistic(LabelStatisticType const &labelStatistic)
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
void PrintSelf(std::ostream &os, itk::Indent indent) const override
static ITK_CONSTEXPR_FUNC T max(const T &)
PersistentLabelImageSmallRegionMergingFilterType::LabelPopulationType LabelPopulationType
void SetLabelPopulation(LabelPopulationType const &labelPopulation)
void SetLabelStatistic(LabelStatisticType const &labelStatistic)
unsigned int ThreadIdType
std::unordered_map< InputLabelType, RealVectorPixelType > LabelStatisticType
std::unordered_map< InputLabelType, double > LabelPopulationType
void SetInputLabelImage(const TInputLabelImage *labelImage)
PersistentLabelImageSmallRegionMergingFilterType::LabelStatisticType LabelStatisticType
void SetLabelPopulation(LabelPopulationType const &labelPopulation)
std::unordered_map< InputLabelType, std::set< InputLabelType > > NeigboursMapType