NeighborhoodIterators6.cxx 7.93 KB
Newer Older
1
/*
Julien Michel's avatar
Julien Michel committed
2
 * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
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's avatar
Jordi Inglada committed
20

21

Jordi Inglada's avatar
Jordi Inglada committed
22 23 24
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
25
#include "itkUnaryFunctorImageFilter.h"
Jordi Inglada's avatar
Jordi Inglada committed
26 27
#include "itkRescaleIntensityImageFilter.h"
#include "itkNeighborhoodIterator.h"
28
#include "itkImageRegionIterator.h"
Jordi Inglada's avatar
Jordi Inglada committed
29 30 31
#include "itkFastMarchingImageFilter.h"
#include "itkAddImageFilter.h"

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's avatar
Jordi Inglada committed
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
50
// such as these can be efficiently written using the random access
Jordi Inglada's avatar
Jordi Inglada committed
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.

67
int main(int argc, char* argv[])
Jordi Inglada's avatar
Jordi Inglada committed
68
{
OTB Bot's avatar
OTB Bot committed
69
  if (argc < 4)
70
  {
Jordi Inglada's avatar
Jordi Inglada committed
71 72
    std::cerr << "Missing parameters. " << std::endl;
    std::cerr << "Usage: " << std::endl;
73
    std::cerr << argv[0] << " outputImageFile startX startY" << std::endl;
Jordi Inglada's avatar
Jordi Inglada committed
74
    return -1;
75
  }
Jordi Inglada's avatar
Jordi Inglada committed
76

77 78 79
  using PixelType                = float;
  using ImageType                = otb::Image<PixelType, 2>;
  using NeighborhoodIteratorType = itk::NeighborhoodIterator<ImageType>;
Jordi Inglada's avatar
Jordi Inglada committed
80

81
  using FastMarchingFilterType = itk::FastMarchingImageFilter<ImageType, ImageType>;
Jordi Inglada's avatar
Jordi Inglada committed
82 83 84

  FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();

85 86
  using NodeContainer = FastMarchingFilterType::NodeContainer;
  using NodeType      = FastMarchingFilterType::NodeType;
Jordi Inglada's avatar
Jordi Inglada committed
87 88

  NodeContainer::Pointer seeds = NodeContainer::New();
89

OTB Bot's avatar
OTB Bot committed
90
  ImageType::IndexType seedPosition;
91

92 93
  seedPosition[0]              = 128;
  seedPosition[1]              = 128;
Jordi Inglada's avatar
Jordi Inglada committed
94 95 96 97
  const double initialDistance = 1.0;

  NodeType node;

OTB Bot's avatar
OTB Bot committed
98
  const double seedValue = -initialDistance;
Jordi Inglada's avatar
Jordi Inglada committed
99 100

  ImageType::SizeType size = {{256, 256}};
101

OTB Bot's avatar
OTB Bot committed
102 103
  node.SetValue(seedValue);
  node.SetIndex(seedPosition);
Jordi Inglada's avatar
Jordi Inglada committed
104
  seeds->Initialize();
OTB Bot's avatar
OTB Bot committed
105
  seeds->InsertElement(0, node);
Jordi Inglada's avatar
Jordi Inglada committed
106

OTB Bot's avatar
OTB Bot committed
107 108
  fastMarching->SetTrialPoints(seeds);
  fastMarching->SetSpeedConstant(1.0);
Jordi Inglada's avatar
Jordi Inglada committed
109

110
  itk::AddImageFilter<ImageType, ImageType, ImageType>::Pointer adder = itk::AddImageFilter<ImageType, ImageType, ImageType>::New();
111

112
  // Allocate the noise image
113
  ImageType::Pointer    noise = ImageType::New();
114 115 116 117
  ImageType::RegionType noiseRegion;
  noiseRegion.SetSize(size);
  noise->SetRegions(noiseRegion);
  noise->Allocate();
118

119 120 121 122 123 124
  // Fill the noise image
  itk::ImageRegionIterator<ImageType> itNoise(noise, noiseRegion);
  itNoise.GoToBegin();

  // Random number seed
  unsigned int sample_seed = 12345;
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<double>(sample_seed) / 2147483711UL;
    rnd         = (1.0 - u) * dMin + u * dMax;

    itNoise.Set((PixelType)rnd);
137
    ++itNoise;
138
  }
139 140

  adder->SetInput1(noise);
Jordi Inglada's avatar
Jordi Inglada committed
141
  adder->SetInput2(fastMarching->GetOutput());
142

Jordi Inglada's avatar
Jordi Inglada committed
143
  try
144
  {
OTB Bot's avatar
OTB Bot committed
145
    fastMarching->SetOutputSize(size);
Jordi Inglada's avatar
Jordi Inglada committed
146 147 148
    fastMarching->Update();

    adder->Update();
149
  }
OTB Bot's avatar
OTB Bot committed
150
  catch (itk::ExceptionObject& excep)
151
  {
Jordi Inglada's avatar
Jordi Inglada committed
152 153
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
154
  }
Jordi Inglada's avatar
Jordi Inglada committed
155 156 157

  ImageType::Pointer input = adder->GetOutput();

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's avatar
Jordi Inglada committed
161 162 163 164 165

  ImageType::IndexType index;
  index[0] = ::atoi(argv[2]);
  index[1] = ::atoi(argv[3]);

166
  // Next we create the neighborhood iterator and position it at the seed point.
Jordi Inglada's avatar
Jordi Inglada committed
167 168 169 170 171 172 173

  NeighborhoodIteratorType::RadiusType radius;
  radius.Fill(1);
  NeighborhoodIteratorType it(radius, input, input->GetRequestedRegion());

  it.SetLocation(index);

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's avatar
Jordi Inglada committed
181 182

  bool flag = true;
OTB Bot's avatar
OTB Bot committed
183
  while (flag == true)
184
  {
Jordi Inglada's avatar
Jordi Inglada committed
185 186 187 188
    NeighborhoodIteratorType::OffsetType nextMove;
    nextMove.Fill(0);

    flag = false;
189

Jordi Inglada's avatar
Jordi Inglada committed
190
    PixelType min = it.GetCenterPixel();
191
    for (unsigned i = 0; i < it.Size(); ++i)
192
    {
OTB Bot's avatar
OTB Bot committed
193
      if (it.GetPixel(i) < min)
194 195
      {
        min      = it.GetPixel(i);
Jordi Inglada's avatar
Jordi Inglada committed
196
        nextMove = it.GetOffset(i);
197
        flag     = true;
Jordi Inglada's avatar
Jordi Inglada committed
198
      }
199
    }
OTB Bot's avatar
OTB Bot committed
200
    it.SetCenterPixel(255.0);
Jordi Inglada's avatar
Jordi Inglada committed
201
    it += nextMove;
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's avatar
Jordi Inglada committed
219

220 221 222
  using WritePixelType = unsigned char;
  using WriteImageType = otb::Image<WritePixelType, 2>;
  using WriterType     = otb::ImageFileWriter<WriteImageType>;
223

224
  using RescaleFilterType = itk::RescaleIntensityImageFilter<ImageType, WriteImageType>;
225

Jordi Inglada's avatar
Jordi Inglada committed
226
  RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
227

OTB Bot's avatar
OTB Bot committed
228 229 230
  rescaler->SetOutputMinimum(0);
  rescaler->SetOutputMaximum(255);
  rescaler->SetInput(input);
231

Jordi Inglada's avatar
Jordi Inglada committed
232
  WriterType::Pointer writer = WriterType::New();
OTB Bot's avatar
OTB Bot committed
233 234
  writer->SetFileName(argv[1]);
  writer->SetInput(rescaler->GetOutput());
Jordi Inglada's avatar
Jordi Inglada committed
235
  try
236
  {
Jordi Inglada's avatar
Jordi Inglada committed
237
    writer->Update();
238
  }
OTB Bot's avatar
OTB Bot committed
239
  catch (itk::ExceptionObject& err)
240
  {
Jordi Inglada's avatar
Jordi Inglada committed
241 242 243
    std::cout << "ExceptionObject caught !" << std::endl;
    std::cout << err << std::endl;
    return -1;
244
  }
245
  return EXIT_SUCCESS;
Jordi Inglada's avatar
Jordi Inglada committed
246
}