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