Orfeo Toolbox  4.0
otbPixelSuppressionByDirectionImageFilter.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 __otbPixelSuppressionByDirectionImageFilter_txx
19 #define __otbPixelSuppressionByDirectionImageFilter_txx
20 
22 
23 #include "itkDataObject.h"
24 #include "itkMacro.h"
27 #include "itkImageRegionIterator.h"
30 #include "itkOffset.h"
31 #include "itkProgressReporter.h"
32 #include "otbMath.h"
33 
34 namespace otb
35 {
36 
40 template <class TInputImage, class TOutputImage>
42 {
43 
44  m_Radius.Fill(1);
45  m_AngularBeam = static_cast<double>(0.);
46 
47 }
48 
49 template <class TInputImage, class TOutputImage>
50 void
53 {
54  this->SetInput(0, input);
55 }
56 
57 template <class TInputImage, class TOutputImage>
58 void
61 {
62  this->SetInput(1, input);
63 }
64 
65 template <class TInputImage, class TOutputImage>
66 const
70 {
71  if (this->GetNumberOfInputs() < 1)
72  {
73  return 0;
74  }
75 
76  return static_cast<const TInputImage *>
77  (this->GetInput(0));
78 }
79 
80 template <class TInputImage, class TOutputImage>
81 const
85 {
86  if (this->GetNumberOfInputs() < 1)
87  {
88  return 0;
89  }
90 
91  return static_cast<const TInputImage *>
92  (this->GetInput(1));
93 }
94 
95 template <class TInputImage, class TOutputImage>
97  itk::InvalidRequestedRegionError)
98  {
99  // call the superclass' implementation of this method
100  Superclass::GenerateInputRequestedRegion();
101 
102  // get pointers to the input and output
103  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInputImageDirection());
104  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
105 
106  if (!inputPtr || !outputPtr)
107  {
108  return;
109  }
110 
111  // get a copy of the input requested region (should equal the output
112  // requested region)
113  typename TInputImage::RegionType inputRequestedRegion;
114  inputRequestedRegion = inputPtr->GetRequestedRegion();
115 
116  // pad the input requested region by the operator radius
117  inputRequestedRegion.PadByRadius(m_Radius);
118 
119  // crop the input requested region at the input's largest possible region
120  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
121  {
122  inputPtr->SetRequestedRegion(inputRequestedRegion);
123  return;
124  }
125  else
126  {
127  // Couldn't crop the region (requested region is outside the largest
128  // possible region). Throw an exception.
129 
130  // store what we tried to request (prior to trying to crop)
131  inputPtr->SetRequestedRegion(inputRequestedRegion);
132 
133  // build an exception
134  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
135  std::ostringstream msg;
136  msg << static_cast<const char *>(this->GetNameOfClass())
137  << "::GenerateInputRequestedRegion()";
138  e.SetLocation(msg.str().c_str());
139  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
140  e.SetDataObject(inputPtr);
141  throw e;
142  }
143  }
144 
145 template<class TInputImage, class TOutputImage>
147  const OutputImageRegionType& outputRegionForThread,
148  itk::ThreadIdType threadId
149  )
150 {
151 
153  const InputPixelType cvalue = 255;
154  cbc.SetConstant(cvalue);
155 
159 
160  // Allocate output
161  typename OutputImageType::Pointer output = this->GetOutput();
162  typename InputImageType::ConstPointer input = this->GetInputImage();
163  typename InputImageType::ConstPointer inputDirection = this->GetInputImageDirection();
164 
165  // Find the data-set boundary "faces"
168 
170  faceList = bC(inputDirection, outputRegionForThread, m_Radius);
171 
172  // support progress methods/callbacks
173  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
174 
175  //typename TInputImage::IndexType bitIndex;
176 
177  //---------------------------------------------------------------------------
178 
179  InputPixelType PixelValue;
180 
181  // Location of the central pixel in the input image
182 // int Xc, Yc;
183 
184  // Pixel location in the system axis of the region
185  int x, y;
186 
187  // Pixel location in the system axis of the region after rotation of theta
188  // where theta is the direction of the cantral pixel
189 
190  // Pixel Direction
191  double ThetaXcYc, Thetaxtyt;
192 
193  //---------------------------------------------------------------------------
194 
195  // Process each of the boundary faces. These are N-d regions which border
196  // the edge of the buffer.
197  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
198  {
199  bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, inputDirection, *fit);
200 
202  itout = itk::ImageRegionIterator<OutputImageType>(output, *fit);
203 
204  bit.OverrideBoundaryCondition(&cbc);
205  bit.GoToBegin();
206 
207  while (!bit.IsAtEnd())
208  {
209 
210  /*// Location of the central pixel of the region in the input image
211  bitIndex = bit.GetIndex();
212 
213  Xc = bitIndex[0];
214  Yc = bitIndex[1];
215  */
216  // Get Pixel Direction from the image of directions
217  ThetaXcYc = static_cast<double>(bit.GetCenterPixel());
218 
219  // Pixel intensity in the input image
220  PixelValue = itin.Get();
221 
222  bool IsLine = false;
223 
225  // Loop on the region
226  for (unsigned int i = 0; i < 2 * m_Radius[0] + 1; ++i)
227  for (unsigned int j = 0; j < 2 * m_Radius[1] + 1; ++j)
228  {
229 
230  off[0] = i - m_Radius[0];
231  off[1] = j - m_Radius[1];
232 
233  x = off[0];
234  y = off[1];
235 
236  // No calculation on the central pixel
237  if ((x == 0) && (y == 0)) continue;
238 
239  Thetaxtyt = vcl_atan2(static_cast<double>(y), static_cast<double>(x)); //result is [-PI, PI]
240  while (Thetaxtyt < 0)
241  Thetaxtyt = CONST_PI + Thetaxtyt; // Theta is now [0, PI] as is
242  // the result of detectors
243  while (Thetaxtyt > CONST_PI_2)
244  Thetaxtyt = Thetaxtyt - CONST_PI; // Theta is now [-PI/2, PI/2]
245 
246  if ((vcl_abs(vcl_cos(Thetaxtyt - ThetaXcYc)) >= vcl_cos(m_AngularBeam)) // this
247  // pixel
248  // is
249  // in
250  // the
251  // angular beam
252  && (vcl_abs(vcl_cos(bit.GetPixel(off) - ThetaXcYc)) >= vcl_cos(m_AngularBeam))) //and
253  //its
254  //direction
255  //is
256  //also
257  //in
258  //the beam
259  {
260  IsLine = true;
261  continue;
262  }
263 
264  }
265 
266  // end of the loop on the pixels of the region
267 
268  // Assignment of this value to the output pixel
269  if (IsLine == true)
270  {
271  itout.Set(static_cast<OutputPixelType>(PixelValue));
272  }
273  else
274  {
275  itout.Set(static_cast<OutputPixelType>(0.));
276  }
277 
278  ++bit;
279  ++itin;
280  ++itout;
281  progress.CompletedPixel();
282 
283  }
284 
285  }
286 }
287 
291 template <class TInputImage, class TOutput>
292 void
294 {
295  Superclass::PrintSelf(os, indent);
296  os << indent << "Radius: " << m_Radius << std::endl;
297 }
298 
299 } // end namespace otb
300 
301 #endif

Generated at Sat Mar 8 2014 16:13:10 for Orfeo Toolbox with doxygen 1.8.3.1