Orfeo Toolbox  4.0
otbTouziEdgeDetectorImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbTouziEdgeDetectorImageFilter_txx
19 #define __otbTouziEdgeDetectorImageFilter_txx
20 
22 
23 #include "itkDataObject.h"
26 #include "itkImageRegionIterator.h"
29 #include "itkProgressReporter.h"
30 #include "otbMacro.h"
31 #include "otbMath.h"
32 
33 namespace otb
34 {
35 
39 template <class TInputImage, class TOutputImage, class TOutputImageDirection>
41 {
42  m_Radius.Fill(1);
43 }
44 
45 template <class TInputImage, class TOutputImage, class TOutputImageDirection>
46 void TouziEdgeDetectorImageFilter<TInputImage, TOutputImage,
47  TOutputImageDirection>::GenerateInputRequestedRegion() throw (
48  itk::InvalidRequestedRegionError)
49  {
50  // call the superclass' implementation of this method
51  Superclass::GenerateInputRequestedRegion();
52 
53  // get pointers to the input and output
54  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
55  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
56 
57  if (!inputPtr || !outputPtr)
58  {
59  return;
60  }
61 
62  // get a copy of the input requested region (should equal the output
63  // requested region)
64  typename TInputImage::RegionType inputRequestedRegion;
65  inputRequestedRegion = inputPtr->GetRequestedRegion();
66 
67  // pad the input requested region by the operator radius
68  inputRequestedRegion.PadByRadius(m_Radius);
69 
70  // crop the input requested region at the input's largest possible region
71  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
72  {
73  inputPtr->SetRequestedRegion(inputRequestedRegion);
74  return;
75  }
76  else
77  {
78  // Couldn't crop the region (requested region is outside the largest
79  // possible region). Throw an exception.
80 
81  // store what we tried to request (prior to trying to crop)
82  inputPtr->SetRequestedRegion(inputRequestedRegion);
83 
84  // build an exception
85  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
86  std::ostringstream msg;
87  msg << static_cast<const char *>(this->GetNameOfClass())
88  << "::GenerateInputRequestedRegion()";
89  e.SetLocation(msg.str().c_str());
90  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
91  e.SetDataObject(inputPtr);
92  throw e;
93  }
94  }
95 
101 template <class TInputImage, class TOutputImage, class TOutputImageDirection>
102 void
105 {
106 
107  typename OutputImageDirectionType::RegionType region;
108  typename OutputImageType::Pointer output = this->GetOutput();
109 
110  OutputImageDirectionType * direction = this->GetOutputDirection();
111 
112  region.SetSize(output->GetRequestedRegion().GetSize());
113  region.SetIndex(output->GetRequestedRegion().GetIndex());
114  direction->SetRegions(region);
115  direction->SetOrigin(output->GetOrigin());
116  direction->SetSpacing(output->GetSpacing());
117  direction->Allocate();
118 }
119 
120 template <class TInputImage, class TOutputImage, class TOutputImageDirection>
121 void
124  const OutputImageRegionType& outputRegionForThread,
125  itk::ThreadIdType threadId
126  )
127 {
128  unsigned int i;
133 
134  // Allocate output
135  typename OutputImageType::Pointer output = this->GetOutput();
136  typename InputImageType::ConstPointer input = this->GetInput();
137  typename OutputImageDirectionType::Pointer outputDir = this->GetOutputDirection();
138 
139  // Find the data-set boundary "faces"
142 
144  faceList = bC(input, outputRegionForThread, m_Radius);
145 
146  // support progress methods/callbacks
147  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
148 
149  typename TInputImage::IndexType bitIndex;
150 
151  // Initializations
152  // ---------------
153  // Number of direction
154  const int NB_DIR = 4;
155  // Number of region of the filter
156  const int NB_REGION = 2;
157  // Definition of the 4 directions
158  double Theta[NB_DIR];
159 
160  Theta[0] = 0.;
161  Theta[1] = CONST_PI_4;
162  Theta[2] = CONST_PI_2;
163  Theta[3] = 3 * CONST_PI / 4.;
164 
165  // contains for the 4 directions the sum of the pixels belonging to each region
166  double Sum[NB_DIR][NB_REGION];
167  // Mean of region 1
168  double M1;
169  // Mean of region 2
170  double M2;
171  // Result of the filter for each direction
172  double R_theta[NB_DIR];
173  double Sum_R_theta = 0.;
174  // Intensity of the contour
175  double R_contour;
176  // Direction of the contour
177  double Dir_contour = 0.;
178  // sign of the contour
179  int sign;
180 
181  // Pixel location in the input image
182  int x;
183  int y;
184  // Location of the central pixel in the input image
185  int xc;
186  int yc;
187 
188  int cpt = 0;
189 
190  // Process each of the boundary faces. These are N-d regions which border
191  // the edge of the buffer.
192  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
193  {
194 
195  cpt += 1;
196 
197  bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
198  unsigned int neighborhoodSize = bit.Size();
199 
201  it_dir = itk::ImageRegionIterator<OutputImageDirectionType>(outputDir, *fit);
202 
203  bit.OverrideBoundaryCondition(&nbc);
204  bit.GoToBegin();
205 
206  while (!bit.IsAtEnd())
207  {
208 
209  // Location of the pixel central
210  bitIndex = bit.GetIndex();
211 
212  xc = bitIndex[0];
213  yc = bitIndex[1];
214 
215  // Initializations
216  for (int dir = 0; dir < NB_DIR; ++dir)
217  {
218  for (int m = 0; m < NB_REGION; m++)
219  Sum[dir][m] = 0.;
220  }
221 
222  R_contour = -1;
223  Dir_contour = 0.;
224  Sum_R_theta = 0.;
225 
226  // Loop on pixels of the filter
227  for (i = 0; i < neighborhoodSize; ++i)
228  {
229 
230  bitIndex = bit.GetIndex(i);
231  x = bitIndex[0];
232  y = bitIndex[1];
233 
234  // We determine for each direction with which region the pixel belongs.
235 
236  // Horizontal direction
237  if (y < yc) Sum[0][0] += static_cast<double>(bit.GetPixel(i));
238  else if (y > yc) Sum[0][1] += static_cast<double>(bit.GetPixel(i));
239 
240  // Diagonal direction 1
241  if ((y - yc) < (x - xc)) Sum[1][0] += static_cast<double>(bit.GetPixel(i));
242  else if ((y - yc) > (x - xc)) Sum[1][1] += static_cast<double>(bit.GetPixel(i));
243 
244  // Vertical direction
245  if (x > xc) Sum[2][0] += static_cast<double>(bit.GetPixel(i));
246  else if (x < xc) Sum[2][1] += static_cast<double>(bit.GetPixel(i));
247 
248  // Diagonal direction 2
249  if ((y - yc) > -(x - xc)) Sum[3][0] += static_cast<double>(bit.GetPixel(i));
250  else if ((y - yc) < -(x - xc)) Sum[3][1] += static_cast<double>(bit.GetPixel(i));
251 
252  } // end of the loop on pixels of the filter
253 
254  // Loop on the 4 directions
255  for (int dir = 0; dir < NB_DIR; ++dir)
256  {
257  // Calculation of the mean of the 2 regions
258  M1 = Sum[dir][0] / static_cast<double>(m_Radius[0] * (2 * m_Radius[0] + 1));
259  M2 = Sum[dir][1] / static_cast<double>(m_Radius[0] * (2 * m_Radius[0] + 1));
260 
261  // Calculation of the intensity of the contour
262  if ((M1 != 0) && (M2 != 0)) R_theta[dir] = static_cast<double>(1 - std::min((M1 / M2), (M2 / M1)));
263  else R_theta[dir] = 0.;
264 
265  // Determination of the maximum intensity of the contour
266  R_contour = static_cast<double>(std::max(R_contour, R_theta[dir]));
267 
268  // Determination of the sign of contour
269  if (M2 > M1) sign = +1;
270  else sign = -1;
271 
272  Dir_contour += sign * Theta[dir] * R_theta[dir];
273  Sum_R_theta += R_theta[dir];
274 
275  } // end of the loop on the directions
276 
277  // Assignment of this value to the output pixel
278  it.Set(static_cast<OutputPixelType>(R_contour));
279 
280  // Determination of the direction of the contour
281  if (Sum_R_theta != 0.) Dir_contour = Dir_contour / Sum_R_theta;
282 
283  // Assignment of this value to the "outputdir" pixel
284  it_dir.Set(static_cast<OutputPixelDirectionType>(Dir_contour));
285 
286  ++bit;
287  ++it;
288  ++it_dir;
289  progress.CompletedPixel();
290 
291  }
292 
293  }
294 }
295 
299 template <class TInputImage, class TOutputImage, class TOutputImageDirection>
300 void
302 ::PrintSelf(std::ostream& os, itk::Indent indent) const
303 {
304  Superclass::PrintSelf(os, indent);
305  os << indent << "Radius: " << m_Radius << std::endl;
306 }
307 
308 } // end namespace otb
309 
310 #endif

Generated at Sat Mar 8 2014 16:22:47 for Orfeo Toolbox with doxygen 1.8.3.1