OTB  9.0.0
Orfeo Toolbox
otbRegionImageToRectangularPathListFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 // Adapted from otbImageToPathListAlignFilter.hxx
23 
24 #ifndef otbRegionImageToRectangularPathListFilter_hxx
25 #define otbRegionImageToRectangularPathListFilter_hxx
26 
27 #include <iostream>
28 #include <exception>
29 #include "otbRegionImageToRectangularPathListFilter.h"
30 #include "itkImageRegionIterator.h"
31 #include "itkImageRegionIteratorWithIndex.h"
32 #include "itkConstNeighborhoodIterator.h"
33 #include "itkPathIterator.h"
34 #include "itkImageLinearIteratorWithIndex.h"
35 #include "otbMath.h"
36 #include "itkNeighborhoodBinaryThresholdImageFunction.h"
37 #include "itkFloodFilledImageFunctionConditionalIterator.h"
38 #include "itkConstNeighborhoodIterator.h"
39 #include "itkConstantBoundaryCondition.h"
40 
41 namespace otb
42 {
43 
45 template <class TInputImage, class TOutputPath>
46 RegionImageToRectangularPathListFilter<TInputImage, TOutputPath>::RegionImageToRectangularPathListFilter()
47 {
48  this->SetNumberOfRequiredInputs(1);
49  m_MinimumFit = 0.0;
50  m_MinimumSize = 0.0;
51 }
53 
55 template <class TInputImage, class TOutputPath>
56 RegionImageToRectangularPathListFilter<TInputImage, TOutputPath>::~RegionImageToRectangularPathListFilter()
57 {
58 }
59 
60 /* Algorithm */
61 
62 template <class TInputImage, class TOutputPath>
63 void RegionImageToRectangularPathListFilter<TInputImage, TOutputPath>::GenerateData(void)
64 {
65  typename InputImageType::SizeType Taille;
66  // this->DebugOn();
67  // itkDebugMacro(<< "RegionImageToRectangularPathListFilter::GenerateData() called");
68 
69  // Get the input and output pointers
70  const InputImageType* InputImage = this->GetInput();
71  OutputPathListType* OutputPath = this->GetOutput();
72  // Generate the image
73 
74  /* Filter algorithm */
75 
76  Taille = InputImage->GetLargestPossibleRegion().GetSize();
77  itkDebugMacro(<< "Input image size : " << Taille);
78 
79  typename InputImageType::PointType origin;
80  typename InputImageType::SpacingType spacing;
81  origin = InputImage->GetOrigin();
82  spacing = InputImage->GetSignedSpacing();
83  std::cout << "Image origin : " << origin << std::endl;
84  std::cout << "Image spacing : " << spacing << std::endl;
85 
86  typedef itk::ImageRegionConstIterator<InputImageType> IteratorType;
87  IteratorType it(InputImage, InputImage->GetLargestPossibleRegion());
88  int pixelCount = 0;
89  int regionCount = 0;
90  int selectedRegionCount = 0;
91  int pixelDebugNumber = 0; // 20;
92  int regionDebugNumber = 0; // 5000; //20;
93 
94  typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
95  typename NeighborhoodIteratorType::RadiusType radius;
96  radius.Fill(1); // square one-pixel neighborhood
97  // 2 std::cout << "radius : " << radius << std::endl;
98  NeighborhoodIteratorType nit(radius, InputImage, InputImage->GetLargestPossibleRegion());
99  NeighborhoodIteratorType nit2(radius, InputImage, InputImage->GetLargestPossibleRegion());
100  // Set up the boundary condition to be zero padded (used on input image)
101  itk::ConstantBoundaryCondition<TInputImage> BC;
102  // BC.SetConstant(itk::NumericTraits<PixelType>::Zero);
103  BC.SetConstant(itk::NumericTraits<PixelType>::max());
104  nit.OverrideBoundaryCondition(&BC); // assign the boundary condition
105  nit2.OverrideBoundaryCondition(&BC); // assign the boundary condition
106 
107  // Build a temporary image of chars for use in the flood algorithm
108  markerImage = MarkerImageType::New();
109  typename MarkerImageType::RegionType markerRegion = InputImage->GetLargestPossibleRegion();
110 
111  markerImage->SetLargestPossibleRegion(markerRegion);
112  markerImage->SetBufferedRegion(markerRegion);
113  markerImage->SetRequestedRegion(markerRegion);
114  markerImage->Allocate();
115  unsigned char maxValue = itk::NumericTraits<typename MarkerImageType::PixelType>::max();
116  markerImage->FillBuffer(maxValue);
117  unsigned char zeroValue = itk::NumericTraits<typename MarkerImageType::PixelType>::Zero;
118  unsigned char regionValue;
119  int neighbor;
120  // markerImage->FillBuffer(itk::NumericTraits<typename MarkerImageType::PixelType>::max());
121 
122  // typedef itk::ImageRegionIterator<MarkerImageType> MarkerIteratorType;
123  typedef itk::ImageRegionIteratorWithIndex<MarkerImageType> MarkerIteratorType;
124  MarkerIteratorType mit(markerImage, markerRegion);
125 
126  typedef typename TInputImage::IndexType IndexType;
127  std::vector<IndexType> regionContainer; // Pb for growing from within loop
128  typename std::vector<IndexType>::iterator regionIterator;
129 
130  regionContainer.reserve(Taille[0] * Taille[1]); // to avoid growth problems
131  IndexType explorerIndex; // position whose neighbors are to be checked for inclusion in current region
132 
133  try
134  {
135  OutputPath->Clear();
136 
137  for (nit.GoToBegin(), pixelCount = 0; !nit.IsAtEnd(); ++nit)
138  {
139  pixelCount++;
140  if (pixelCount <= pixelDebugNumber)
141  {
142  std::cout << "Pixel #" << pixelCount << " : " << nit.GetCenterPixel() << std::endl;
143  for (neighbor = 1; neighbor <= 7; neighbor += 2) // 4 neighbors : above, left, right, below
144  std::cout << "Neighbor #" << neighbor << " : " << nit.GetPixel(neighbor) << std::endl;
145  }
146  mit.SetIndex(nit.GetIndex());
147  if (mit.Get() == maxValue) // pixel not yet processed
148  {
149  mit.Set(zeroValue);
150  regionContainer.clear();
151  regionContainer.push_back(nit.GetIndex());
152  regionValue = nit.GetCenterPixel();
153  regionCount++;
154  if (pixelCount <= pixelDebugNumber)
155  {
156  std::cout << "Starting new region at " << nit.GetIndex() << " with value " << (unsigned short)regionValue << std::endl;
157  }
158  // Collect all pixels in same region (same value), reachable in 4-connectivity
159  for (regionIterator = regionContainer.begin(); regionIterator != /*<*/ regionContainer.end(); ++regionIterator /* regionContainer grows within loop */)
160  {
161  // add neighbors not yet processed, until whole region has been collected (no "new" neighbor)
162  explorerIndex = *regionIterator;
163  nit2.SetLocation(explorerIndex);
164  if (pixelCount <= pixelDebugNumber)
165  std::cout << "Exploring neighbors of " << explorerIndex << std::endl;
166 
167  for (neighbor = 1; neighbor <= 7; neighbor += 2)
168  if (nit2.GetPixel(neighbor) == regionValue) // ZZ and not yet processed...
169  {
170  mit.SetIndex(nit2.GetIndex(neighbor));
171  if (mit.Get() == maxValue) // pixel not yet processed
172  {
173  mit.Set(zeroValue);
174  regionContainer.push_back(nit2.GetIndex(neighbor));
175  if (pixelCount <= pixelDebugNumber)
176  {
177  std::cout << "Adding " << nit2.GetIndex(neighbor) << std::endl;
178  std::cout << "Added " << regionContainer.back() << std::endl;
179  }
180  }
181  }
182  }
183  if (pixelCount <= pixelDebugNumber)
184  {
185  std::cout << "Region queue (" << regionContainer.size() << " elements) : " << std::endl;
186  for (regionIterator = regionContainer.begin(); regionIterator != regionContainer.end(); ++regionIterator)
187  std::cout << *regionIterator;
188  std::cout << std::endl;
189  }
190 
191  // Compute variance-covariance matrix of x-y coordinates of region pixels
192  double sumX = 0, sumY = 0, sumX2 = 0, sumY2 = 0, sumXY = 0;
193  double avgX, avgY, varX, varY, covarXY;
194  int n = regionContainer.size();
195  for (regionIterator = regionContainer.begin(); regionIterator != regionContainer.end(); ++regionIterator)
196  {
197  explorerIndex = *regionIterator;
198  sumX += explorerIndex[0];
199  sumY += explorerIndex[1];
200  sumX2 += explorerIndex[0] * explorerIndex[0];
201  sumY2 += explorerIndex[1] * explorerIndex[1];
202  sumXY += explorerIndex[0] * explorerIndex[1];
203  }
204  avgX = (double)sumX / n;
205  avgY = (double)sumY / n;
206  varX = (double)sumX2 / n - avgX * avgX;
207  varY = (double)sumY2 / n - avgY * avgY;
208  covarXY = (double)sumXY / n - avgX * avgY;
209 
210  // Compute average deviation (or mean absolute deviation) matrix instead, because rectangles resulting from above statistics are too elongated
211  double sumAX = 0, sumAY = 0, crossTermAXY = 0;
212  double adevX, adevY, adevXY;
213  double ax, ay;
214  for (regionIterator = regionContainer.begin(); regionIterator != regionContainer.end(); ++regionIterator)
215  {
216  explorerIndex = *regionIterator;
217  ax = std::abs(explorerIndex[0] - avgX);
218  ay = std::abs(explorerIndex[1] - avgY);
219  sumAX += ax;
220  sumAY += ay;
221  crossTermAXY += ax * ay;
222  }
223  adevX = sumAX / n;
224  adevY = sumAY / n;
225  adevXY = std::sqrt(crossTermAXY) / n;
226 
227  // Compute eigenvalues and eigenvectors of variance-covariance matrix (for DIRECTION)
228  double delta, l1, l2, /* eigenvalues */
229  y1, y2, /* eigenvectors y coordinate, for x = 1*/
230  x1 = 1, /* first eigenvector x coordinate */
231  x2 = 1, /* second eigenvector x coordinate, 1 except in special case when covarXY == 0 */
232  alpha /* main direction */;
233  delta = (varX - varY) * (varX - varY) + 4 * covarXY * covarXY;
234  l1 = (varX + varY + std::sqrt(delta)) / 2;
235  l2 = (varX + varY - std::sqrt(delta)) / 2;
236  if (covarXY != 0.0) // ZZ or larger than a small epsilon ? (eg 10^(-15)...)
237  {
238  y1 = (l1 - varX) / covarXY; // for x1 = 1
239  y2 = (l2 - varX) / covarXY; // for x2 = 1
240  }
241  else // matrix was already diagonal
242  {
243  y1 = 0;
244  x2 = 0;
245  y2 = 1;
246  }
247 
248  // Compute eigenvalues and eigenvectors of absolute mean deviation matrix (for PROPORTIONS)
249  double adelta, al1, al2, /* eigenvalues */
250  ay1, ay2, /* eigenvectors y coordinate, for x = 1*/
251  ax1 = 1, /* first eigenvector x coordinate */
252  ax2 = 1; /* second eigenvector x coordinate, 1 except in special case when covarXY == 0 */
253  adelta = (adevX - adevY) * (adevX - adevY) + 4 * adevXY * adevXY;
254  al1 = (adevX + adevY + std::sqrt(adelta)) / 2;
255  al2 = (adevX + adevY - std::sqrt(adelta)) / 2;
256  if (adevXY != 0.0) // ZZ or larger than a small epsilon ? (eg 10^(-15)...)
257  {
258  ay1 = (al1 - adevX) / adevXY; // for x1 = 1
259  ay2 = (al2 - adevX) / adevXY; // for x2 = 1
260  }
261  else // matrix was already diagonal
262  {
263  ay1 = 0;
264  ax2 = 0;
265  ay2 = 1;
266  }
267 
268  if (y1 != 0)
269  alpha = std::atan(1 / y1) * 180 / vnl_math::pi;
270  else
271  alpha = 90;
272  if (alpha < 0)
273  alpha += 180; // Conventionnaly given as a value between 0 and 180 degrees
274 
275  // Compute equivalent length and width (based on equal area criterion)
276  double length, width;
277  if (al2 != 0)
278  {
279  length = std::sqrt(std::abs(al1 / al2) * n);
280  // length = std::sqrt(l1 / l2 * n);
281  if (al1 != 0)
282  width = std::abs(al2 / al1) * length;
283  else // l1 == 0 and l2 == 0
284  {
285  length = width = std::sqrt(static_cast<double>(n)); // should happen only when n == 1 anyway
286  }
287  }
288  else
289  {
290  length = n; // Arbitrary representation for degenerate case
291  width = 1;
292  }
293 
294  // Normalize eigenvectors (for following tests)
295  double norm;
296  norm = std::sqrt(x1 * x1 + y1 * y1);
297  assert(norm != 0); //(by construction of eigenvectors)
298  if (norm != 0)
299  {
300  x1 /= norm;
301  y1 /= norm;
302  }
303  norm = std::sqrt(x2 * x2 + y2 * y2);
304  assert(norm != 0); //(by construction of eigenvectors)
305  if (norm != 0)
306  {
307  x2 /= norm;
308  y2 /= norm;
309  }
310 
311  // Normalize eigenvectors (for following tests) (No : used only for printing values for debugging)
312  norm = std::sqrt(ax1 * ax1 + ay1 * ay1);
313  assert(norm != 0); //(by construction of eigenvectors)
314  if (norm != 0)
315  {
316  ax1 /= norm;
317  ay1 /= norm;
318  }
319  norm = std::sqrt(ax2 * ax2 + ay2 * ay2);
320  assert(norm != 0); //(by construction of eigenvectors)
321  if (norm != 0)
322  {
323  ax2 /= norm;
324  ay2 /= norm;
325  }
326 
327  // Count the proportion of region pixels contained within rectangle model, to evaluate rectangular fit, or "rectangularity"
328  // Rectangle model uses [x1, y1] and [x2, y2] for direction (angle), and size derived from adev matrix
329  double vx, vy; // x-y coordinates relative to rectangle center
330  double halfLength = length / 2, halfWidth = width / 2; // to write formulas more easily
331  int countWithin = 0; // number of pixels contained within rectangle
332  for (regionIterator = regionContainer.begin(); regionIterator != regionContainer.end(); ++regionIterator)
333  {
334  explorerIndex = *regionIterator;
335  vx = explorerIndex[0] - avgX;
336  vy = explorerIndex[1] - avgY;
337  if (std::abs(vx * x1 + vy * y1) <= halfLength && std::abs(vx * x2 + vy * y2) <= halfWidth)
338  countWithin++;
339  }
340 
341  if (regionCount <= regionDebugNumber)
342  {
343  std::cout << std::endl << "Region " << regionCount << " (area = " << n << " pixels)" << std::endl;
344  std::cout << "sumX = " << sumX << "; sumY = " << sumY << "; sumX2 = " << sumX2 << "; sumY2 = " << sumY2 << "; sumXY = " << sumXY << std::endl;
345  std::cout << "avgX = " << avgX << "; avgY = " << avgY << std::endl;
346  std::cout << "varX = " << varX << "; varY = " << varY << "; covarXY = " << covarXY << std::endl;
347  std::cout << "adevX = " << adevX << "; adevY = " << adevY << "; adevXY = " << adevXY << std::endl;
348  std::cout << "crossTermAXY = " << crossTermAXY << std::endl;
349  std::cout << "eigenvalue 1 = " << l1 << "; eigenvalue 2 = " << l2 << std::endl;
350  std::cout << "eigenvector 1 = [" << x1 << ", " << y1 << "]; eigenvector 2 = [" << x2 << ", " << y2 << "]" << std::endl;
351  std::cout << "A-eigenvalue 1 = " << al1 << "; A-eigenvalue 2 = " << al2 << std::endl;
352  std::cout << "A-eigenvector 1 = [" << ax1 << ", " << ay1 << "]; A-eigenvector 2 = [" << ax2 << ", " << ay2 << "]" << std::endl;
353  std::cout << "length = " << length << "; width = " << width << std::endl;
354  std::cout << "main direction = " << alpha << " degree" << std::endl;
355  std::cout << "rectangular fit = " << (float)countWithin / n << std::endl;
356  }
357 
358  // Build rectangle list, converting image coordinates into physical coordinates
359  typedef typename OutputPathType::ContinuousIndexType ContinuousIndexType;
360 
361  ContinuousIndexType point;
362 
363  OutputPathPointerType path = OutputPathType::New();
364 
365  path->Initialize();
366  point[0] = (avgX + x1 * halfLength + x2 * halfWidth) * spacing[0] + origin[0];
367  point[1] = (avgY + y1 * halfLength + y2 * halfWidth) * spacing[1] + origin[1];
368  path->AddVertex(point);
369  if (regionCount <= regionDebugNumber)
370  std::cout << "corner 1 : [" << point[0] << ", " << point[1] << "]" << std::endl;
371  point[0] = (avgX - x1 * halfLength + x2 * halfWidth) * spacing[0] + origin[0];
372  point[1] = (avgY - y1 * halfLength + y2 * halfWidth) * spacing[1] + origin[1];
373  path->AddVertex(point);
374  if (regionCount <= regionDebugNumber)
375  std::cout << "corner 2 : [" << point[0] << ", " << point[1] << "]" << std::endl;
376  point[0] = (avgX - x1 * halfLength - x2 * halfWidth) * spacing[0] + origin[0];
377  point[1] = (avgY - y1 * halfLength - y2 * halfWidth) * spacing[1] + origin[1];
378  path->AddVertex(point);
379  if (regionCount <= regionDebugNumber)
380  std::cout << "corner 3 : [" << point[0] << ", " << point[1] << "]" << std::endl;
381  point[0] = (avgX + x1 * halfLength - x2 * halfWidth) * spacing[0] + origin[0];
382  point[1] = (avgY + y1 * halfLength - y2 * halfWidth) * spacing[1] + origin[1];
383  path->AddVertex(point);
384  if (regionCount <= regionDebugNumber)
385  std::cout << "corner 4 : [" << point[0] << ", " << point[1] << "]" << std::endl;
386  point[0] = (avgX + x1 * halfLength + x2 * halfWidth) * spacing[0] + origin[0];
387  point[1] = (avgY + y1 * halfLength + y2 * halfWidth) * spacing[1] + origin[1];
388  path->AddVertex(point);
389 
390  path->SetValue((double)countWithin / n);
391 
392  if ((float)countWithin / n >= m_MinimumFit // keep only rectangles with fit larger than minimumFit
393  && n >= m_MinimumSize) // and size larger than minimumSize
394  {
395  OutputPath->PushBack(path);
396  selectedRegionCount++;
397  }
398  }
399  }
400  }
401  catch (std::exception& e)
402  {
403  std::cerr << "Caught exception " << e.what() << std::endl;
404  }
405  std::cout << pixelCount << " pixels seen." << std::endl;
406  std::cout << regionCount << " regions processed, " << selectedRegionCount << " regions selected." << std::endl;
407 
408  itkDebugMacro(<< "ImageToPathListAlignFilter::GenerateData() finished");
409 
410 } // end update function
411 
412 template <class TInputImage, class TOutputPath>
413 void RegionImageToRectangularPathListFilter<TInputImage, TOutputPath>::PrintSelf(std::ostream& os, itk::Indent indent) const
414 {
415  Superclass::PrintSelf(os, indent);
416 }
417 
418 } // end namespace otb
419 
420 #endif
otb::ObjectList::Clear
void Clear(void)
Definition: otbObjectList.hxx:201
otb::ObjectList::PushBack
void PushBack(ObjectType *element)
Definition: otbObjectList.hxx:82
otb::RegionImageToRectangularPathListFilter::RegionImageToRectangularPathListFilter
RegionImageToRectangularPathListFilter()
Definition: otbRegionImageToRectangularPathListFilter.hxx:46
otbMath.h
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::RegionImageToRectangularPathListFilter::~RegionImageToRectangularPathListFilter
~RegionImageToRectangularPathListFilter() override
Definition: otbRegionImageToRectangularPathListFilter.hxx:56
otb::ImageToPathListFilter::InputImageType
TInputImage InputImageType
Definition: otbImageToPathListFilter.h:54
otb::RegionImageToRectangularPathListFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbRegionImageToRectangularPathListFilter.hxx:413
otb::ImageToPathListFilter::OutputPathPointerType
Superclass::OutputPathPointerType OutputPathPointerType
Definition: otbImageToPathListFilter.h:66
otbRegionImageToRectangularPathListFilter.h
otb::RegionImageToRectangularPathListFilter::GenerateData
void GenerateData() override
Definition: otbRegionImageToRectangularPathListFilter.hxx:63
otb::ObjectList
This class is a generic all-purpose wrapping around an std::vector<itk::SmartPointer<ObjectType> >.
Definition: otbObjectList.h:40
otb::Image::RegionType
Superclass::RegionType RegionType
Definition: otbImage.h:154
Generated at Wed Feb 28 2024 15:23:35 for OTB with   doxygen 1.8.17