NeighborhoodIterators6.cxx 7.93 KB
 Sébastien Dinot committed Mar 08, 2017 1 /*  Julien Michel committed Jan 14, 2019 2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)  Sébastien Dinot committed Mar 08, 2017 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19  * * This file is part of Orfeo Toolbox * * https://www.orfeo-toolbox.org/ * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */  Jordi Inglada committed Jul 04, 2006 20   Emmanuel Christophe committed May 07, 2011 21   Jordi Inglada committed Jul 04, 2006 22 23 24 #include "otbImage.h" #include "otbImageFileReader.h" #include "otbImageFileWriter.h"  Guillaume Pasero committed May 04, 2015 25 #include "itkUnaryFunctorImageFilter.h"  Jordi Inglada committed Jul 04, 2006 26 27 #include "itkRescaleIntensityImageFilter.h" #include "itkNeighborhoodIterator.h"  Julien Malik committed Nov 15, 2013 28 #include "itkImageRegionIterator.h"  Jordi Inglada committed Jul 04, 2006 29 30 31 #include "itkFastMarchingImageFilter.h" #include "itkAddImageFilter.h"  Victor Poughon committed Feb 20, 2019 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 /* Example usage: ./NeighborhoodIterators6 Output/NeighborhoodIterators6a.png 100 100 */ /* Example usage: ./NeighborhoodIterators6 Output/NeighborhoodIterators6b.png 50 150 */ /* Example usage: ./NeighborhoodIterators6 Output/NeighborhoodIterators6c.png 150 50 */  Jordi Inglada committed Jul 04, 2006 47 48 49 // Some image processing routines do not need to visit every pixel in an // image. Flood-fill and connected-component algorithms, for example, only // visit pixels that are locally connected to one another. Algorithms  Emmanuel Christophe committed Dec 06, 2008 50 // such as these can be efficiently written using the random access  Jordi Inglada committed Jul 04, 2006 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 // capabilities of the neighborhood iterator. // // The following example finds local minima. Given a seed point, we can search // the neighborhood of that point and pick the smallest value $m$. While $m$ // is not at the center of our current neighborhood, we move in the direction // of $m$ and repeat the analysis. Eventually we discover a local minimum and // stop. This algorithm is made trivially simple in ND using an ITK // neighborhood iterator. // // To illustrate the process, we create an image that descends everywhere to a // single minimum: a positive distance transform to a point. The details of // creating the distance transform are not relevant to the discussion of // neighborhood iterators, but can be found in the source code of this // example. Some noise has been added to the distance transform image for // additional interest.  Victor Poughon committed Feb 20, 2019 67 int main(int argc, char* argv[])  Jordi Inglada committed Jul 04, 2006 68 {  OTB Bot committed Apr 01, 2010 69  if (argc < 4)  Victor Poughon committed Feb 20, 2019 70  {  Jordi Inglada committed Jul 04, 2006 71 72  std::cerr << "Missing parameters. " << std::endl; std::cerr << "Usage: " << std::endl;  Victor Poughon committed Feb 20, 2019 73  std::cerr << argv[0] << " outputImageFile startX startY" << std::endl;  Jordi Inglada committed Jul 04, 2006 74  return -1;  Victor Poughon committed Feb 20, 2019 75  }  Jordi Inglada committed Jul 04, 2006 76   Victor Poughon committed May 16, 2019 77 78 79  using PixelType = float; using ImageType = otb::Image; using NeighborhoodIteratorType = itk::NeighborhoodIterator;  Jordi Inglada committed Jul 04, 2006 80   Victor Poughon committed May 16, 2019 81  using FastMarchingFilterType = itk::FastMarchingImageFilter;  Jordi Inglada committed Jul 04, 2006 82 83 84  FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();  Victor Poughon committed May 16, 2019 85 86  using NodeContainer = FastMarchingFilterType::NodeContainer; using NodeType = FastMarchingFilterType::NodeType;  Jordi Inglada committed Jul 04, 2006 87 88  NodeContainer::Pointer seeds = NodeContainer::New();  Emmanuel Christophe committed Dec 06, 2008 89   OTB Bot committed Apr 01, 2010 90  ImageType::IndexType seedPosition;  Emmanuel Christophe committed Dec 06, 2008 91   Victor Poughon committed Feb 20, 2019 92 93  seedPosition[0] = 128; seedPosition[1] = 128;  Jordi Inglada committed Jul 04, 2006 94 95 96 97  const double initialDistance = 1.0; NodeType node;  OTB Bot committed Apr 01, 2010 98  const double seedValue = -initialDistance;  Jordi Inglada committed Jul 04, 2006 99 100  ImageType::SizeType size = {{256, 256}};  Emmanuel Christophe committed Dec 06, 2008 101   OTB Bot committed Apr 01, 2010 102 103  node.SetValue(seedValue); node.SetIndex(seedPosition);  Jordi Inglada committed Jul 04, 2006 104  seeds->Initialize();  OTB Bot committed Apr 01, 2010 105  seeds->InsertElement(0, node);  Jordi Inglada committed Jul 04, 2006 106   OTB Bot committed Apr 01, 2010 107 108  fastMarching->SetTrialPoints(seeds); fastMarching->SetSpeedConstant(1.0);  Jordi Inglada committed Jul 04, 2006 109   Victor Poughon committed Feb 20, 2019 110  itk::AddImageFilter::Pointer adder = itk::AddImageFilter::New();  Julien Malik committed Nov 15, 2013 111   Julien Malik committed Nov 15, 2013 112  // Allocate the noise image  Victor Poughon committed Feb 20, 2019 113  ImageType::Pointer noise = ImageType::New();  Julien Malik committed Nov 15, 2013 114 115 116 117  ImageType::RegionType noiseRegion; noiseRegion.SetSize(size); noise->SetRegions(noiseRegion); noise->Allocate();  Julien Malik committed Nov 15, 2013 118   Julien Malik committed Nov 15, 2013 119 120 121 122 123 124  // Fill the noise image itk::ImageRegionIterator itNoise(noise, noiseRegion); itNoise.GoToBegin(); // Random number seed unsigned int sample_seed = 12345;  Victor Poughon committed Feb 20, 2019 125 126 127 128 129 130 131 132 133 134 135 136  double u = 0.; double rnd = 0.; double dMin = -.7; double dMax = .8; while (!itNoise.IsAtEnd()) { sample_seed = (sample_seed * 16807) % 2147483647L; u = static_cast(sample_seed) / 2147483711UL; rnd = (1.0 - u) * dMin + u * dMax; itNoise.Set((PixelType)rnd);  Julien Malik committed Nov 15, 2013 137  ++itNoise;  Victor Poughon committed Feb 20, 2019 138  }  Julien Malik committed Nov 15, 2013 139 140  adder->SetInput1(noise);  Jordi Inglada committed Jul 04, 2006 141  adder->SetInput2(fastMarching->GetOutput());  Emmanuel Christophe committed Dec 06, 2008 142   Jordi Inglada committed Jul 04, 2006 143  try  Victor Poughon committed Feb 20, 2019 144  {  OTB Bot committed Apr 01, 2010 145  fastMarching->SetOutputSize(size);  Jordi Inglada committed Jul 04, 2006 146 147 148  fastMarching->Update(); adder->Update();  Victor Poughon committed Feb 20, 2019 149  }  OTB Bot committed Apr 01, 2010 150  catch (itk::ExceptionObject& excep)  Victor Poughon committed Feb 20, 2019 151  {  Jordi Inglada committed Jul 04, 2006 152 153  std::cerr << "Exception caught !" << std::endl; std::cerr << excep << std::endl;  Victor Poughon committed Feb 20, 2019 154  }  Jordi Inglada committed Jul 04, 2006 155 156 157  ImageType::Pointer input = adder->GetOutput();  Victor Poughon committed Feb 20, 2019 158 159 160  // The variable \code{input} is the pointer to the distance transform image. // The local minimum algorithm is initialized with a seed point read from the // command line.  Jordi Inglada committed Jul 04, 2006 161 162 163 164 165  ImageType::IndexType index; index[0] = ::atoi(argv[2]); index[1] = ::atoi(argv[3]);  Victor Poughon committed Feb 20, 2019 166  // Next we create the neighborhood iterator and position it at the seed point.  Jordi Inglada committed Jul 04, 2006 167 168 169 170 171 172 173  NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); NeighborhoodIteratorType it(radius, input, input->GetRequestedRegion()); it.SetLocation(index);  Victor Poughon committed Feb 20, 2019 174 175 176 177 178 179 180  // Searching for the local minimum involves finding the minimum in the current // neighborhood, then shifting the neighborhood in the direction of that // minimum. The \code{for} loop below records the \doxygen{itk}{Offset} of the // minimum neighborhood pixel. The neighborhood iterator is then moved using // that offset. When a local minimum is detected, \code{flag} will remain // false and the \code{while} loop will exit. Note that this code is // valid for an image of any dimensionality.  Jordi Inglada committed Jul 04, 2006 181 182  bool flag = true;  OTB Bot committed Apr 01, 2010 183  while (flag == true)  Victor Poughon committed Feb 20, 2019 184  {  Jordi Inglada committed Jul 04, 2006 185 186 187 188  NeighborhoodIteratorType::OffsetType nextMove; nextMove.Fill(0); flag = false;  Emmanuel Christophe committed Dec 06, 2008 189   Jordi Inglada committed Jul 04, 2006 190  PixelType min = it.GetCenterPixel();  Emmanuel Christophe committed Jul 22, 2011 191  for (unsigned i = 0; i < it.Size(); ++i)  Victor Poughon committed Feb 20, 2019 192  {  OTB Bot committed Apr 01, 2010 193  if (it.GetPixel(i) < min)  Victor Poughon committed Feb 20, 2019 194 195  { min = it.GetPixel(i);  Jordi Inglada committed Jul 04, 2006 196  nextMove = it.GetOffset(i);  Victor Poughon committed Feb 20, 2019 197  flag = true;  Jordi Inglada committed Jul 04, 2006 198  }  Victor Poughon committed Feb 20, 2019 199  }  OTB Bot committed Apr 01, 2010 200  it.SetCenterPixel(255.0);  Jordi Inglada committed Jul 04, 2006 201  it += nextMove;  Victor Poughon committed Feb 20, 2019 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218  } // Figure~\ref{fig:NeighborhoodExample6} shows the results of the algorithm // for several seed points. The white line is the path of the iterator from // the seed point to the minimum in the center of the image. The effect of the // additive noise is visible as the small perturbations in the paths. // // \begin{figure} \centering // \includegraphics[width=0.3\textwidth]{NeighborhoodIterators6a.eps} // \includegraphics[width=0.3\textwidth]{NeighborhoodIterators6b.eps} // \includegraphics[width=0.3\textwidth]{NeighborhoodIterators6c.eps} // \itkcaption[Finding local minima]{Paths traversed by the neighborhood // iterator from different seed points to the local minimum. // The true minimum is at the center // of the image. The path of the iterator is shown in white. The effect of // noise in the image is seen as small perturbations in each path. } // \protect\label{fig:NeighborhoodExample6} \end{figure}  Jordi Inglada committed Jul 04, 2006 219   Victor Poughon committed May 16, 2019 220 221 222  using WritePixelType = unsigned char; using WriteImageType = otb::Image; using WriterType = otb::ImageFileWriter;  Emmanuel Christophe committed Dec 06, 2008 223   Victor Poughon committed May 16, 2019 224  using RescaleFilterType = itk::RescaleIntensityImageFilter;  Emmanuel Christophe committed Dec 06, 2008 225   Jordi Inglada committed Jul 04, 2006 226  RescaleFilterType::Pointer rescaler = RescaleFilterType::New();  Emmanuel Christophe committed Dec 06, 2008 227   OTB Bot committed Apr 01, 2010 228 229 230  rescaler->SetOutputMinimum(0); rescaler->SetOutputMaximum(255); rescaler->SetInput(input);  Emmanuel Christophe committed Dec 06, 2008 231   Jordi Inglada committed Jul 04, 2006 232  WriterType::Pointer writer = WriterType::New();  OTB Bot committed Apr 01, 2010 233 234  writer->SetFileName(argv[1]); writer->SetInput(rescaler->GetOutput());  Jordi Inglada committed Jul 04, 2006 235  try  Victor Poughon committed Feb 20, 2019 236  {  Jordi Inglada committed Jul 04, 2006 237  writer->Update();  Victor Poughon committed Feb 20, 2019 238  }  OTB Bot committed Apr 01, 2010 239  catch (itk::ExceptionObject& err)  Victor Poughon committed Feb 20, 2019 240  {  Jordi Inglada committed Jul 04, 2006 241 242 243  std::cout << "ExceptionObject caught !" << std::endl; std::cout << err << std::endl; return -1;  Victor Poughon committed Feb 20, 2019 244  }  Emmanuel Christophe committed Nov 17, 2008 245  return EXIT_SUCCESS;  Jordi Inglada committed Jul 04, 2006 246 }