OTB  5.0.0
Orfeo Toolbox
otbLineDetectorImageFilterBase.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 __otbLineDetectorImageFilterBase_txx
19 #define __otbLineDetectorImageFilterBase_txx
20 
22 
23 #include "itkDataObject.h"
24 
28 #include "itkImageRegionIterator.h"
30 #include "itkProgressReporter.h"
31 #include "otbMacro.h"
32 #include "otbMath.h"
33 
34 namespace otb
35 {
36 
40 template <class TInputImage, class TOutputImage, class TOutputImageDirection, class InterpolatorType>
43 {
44  this->SetNumberOfRequiredOutputs(2);
45  this->SetNumberOfRequiredOutputs(1);
46  m_Radius.Fill(1);
47  m_LengthLine = 1;
48  m_WidthLine = 0;
49  m_Threshold = 0;
50  m_NumberOfDirections = 4;
51  m_FaceList.Fill(0);
52 // THOMAS : donc contructeur de base
53 // this->SetNthOutput(0, OutputImageType::New());
54 // this->SetNthOutput(1, OutputImageType::New());
55 }
57 
58 template <class TInputImage, class TOutputImage, class TOutputImageDirection, class InterpolatorType>
59 void
61 ::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
62  {
63  // call the superclass' implementation of this method
64  Superclass::GenerateInputRequestedRegion();
65 
66  // get pointers to the input and output
67  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
68  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
69 
70  if (!inputPtr || !outputPtr)
71  {
72  return;
73  }
74 
75  // get a copy of the input requested region (should equal the output
76  // requested region)
77  typename TInputImage::RegionType inputRequestedRegion;
78  inputRequestedRegion = inputPtr->GetRequestedRegion();
79 
80  // Define the size of the region by the radius
81  m_Radius[1] = static_cast<unsigned int>(3 * (2 * m_WidthLine + 1) + 2);
82  m_Radius[0] = 2 * m_LengthLine + 1;
83 
84  // Define the size of the facelist by taking into account the rotation of the region
85  m_FaceList[0] =
86  static_cast<unsigned int>(vcl_sqrt(static_cast<double>(
87  (m_Radius[0] * m_Radius[0]) + (m_Radius[1] * m_Radius[1]))) + 1);
88  m_FaceList[1] = m_FaceList[0];
89 
90  otbMsgDevMacro(<< "Radius : " << m_Radius[0] << " " << m_Radius[1]);
91  otbMsgDevMacro(<< "-> FaceList : " << m_FaceList[0] << " " << m_FaceList[1]);
92 
93  // pad the input requested region by the operator radius
94  inputRequestedRegion.PadByRadius(m_FaceList);
95 
96  // crop the input requested region at the input's largest possible region
97  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
98  {
99  inputPtr->SetRequestedRegion(inputRequestedRegion);
100  return;
101  }
102  else
103  {
104  // Couldn't crop the region (requested region is outside the largest
105  // possible region). Throw an exception.
106 
107  // store what we tried to request (prior to trying to crop)
108  inputPtr->SetRequestedRegion(inputRequestedRegion);
109 
110  // build an exception
111  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
112  std::ostringstream msg;
113  msg << static_cast<const char *>(this->GetNameOfClass())
114  << "::GenerateInputRequestedRegion()";
115  e.SetLocation(msg.str().c_str());
116  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
117  e.SetDataObject(inputPtr);
118  throw e;
119  }
120  }
121 
122 /*
123  * Set up state of filter before multi-threading.
124  * InterpolatorType::SetInputImage is not thread-safe and hence
125  * has to be set up before ThreadedGenerateData
126  */
127 template <class TInputImage, class TOutputImage, class TOutputImageDirection, class InterpolatorType>
128 void
131 {
132  typename OutputImageType::Pointer output = this->GetOutput();
133  output->FillBuffer(0);
134  typename OutputImageType::Pointer outputDirection = this->GetOutputDirection();
135  outputDirection->FillBuffer(0);
136 }
137 
138 template <class TInputImage, class TOutputImage, class TOutputImageDirection, class InterpolatorType>
139 void
142  const OutputImageRegionType& outputRegionForThread,
143  itk::ThreadIdType threadId
144  )
145 {
146 
147  typename InputImageType::ConstPointer input = this->GetInput();
148 
149  InterpolatorPointer interpolator2 = InterpolatorType::New();
150  interpolator2->SetInputImage(input);
151 
157 
158  // Allocate output
159  typename OutputImageType::Pointer output = this->GetOutput();
160  typename OutputImageType::Pointer outputDir = this->GetOutputDirection();
161 
162  // Find the data-set boundary "faces"
165 
167  faceList = bC(input, outputRegionForThread, m_FaceList);
168 
169  // support progress methods/callbacks
170  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
171 
172  typename TInputImage::IndexType bitIndex;
173  typename InterpolatorType::ContinuousIndexType Index;
174 
175  // --------------------------------------------------------------------------
176 
177  // Number of direction
178  const unsigned int NB_DIR = this->GetNumberOfDirections();
179  // Number of zone
180  const int NB_ZONE = 3;
181  // Definition of the 4 directions
182  //double Theta[NB_DIR];
183  //ROMAIN
184  double* Theta = new double[NB_DIR];
185 
186  // La rotation nulle correspond a un contour horizontal -> 0 !!
187  for (unsigned int i = 0; i < NB_DIR; ++i)
188  {
189  Theta[i] = (CONST_PI * (i / double(NB_DIR)));
190  /* if(Theta[i]>CONST_PI)
191  Theta[i] = Theta[i]-CONST_PI;
192  if((i/double(NB_DIR))==0.5)
193  Theta[i]=0.; */
194  }
195 
196  // Number of the zone
197  unsigned int zone;
198 
199  // Intensity of the linear feature
200  double R;
201 
202  // Direction of detection
203  double Direction = 0.;
204 
205  // Pixel location in the input image
206  int X, Y;
207 
208  // Pixel location after rotation in the system axis of the region
209  double xout, yout;
210 
211  // location of the central pixel in the input image
212  int Xc, Yc;
213 
214  // location of the central pixel between zone 1 and 2 and between zone 1 and 3
215  int Yc12, Yc13;
216 
217  //---------------------------------------------------------------------------
218  otbMsgDevMacro(<< "Theta : " << Theta[0] << " " << Theta[1] << " " << Theta[2] << " " << Theta[3]);
219 
220  // Process each of the boundary faces. These are N-d regions which border
221  // the edge of the buffer.
222 
223  //bool interiorFace = true;
224 
225  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
226  {
227  bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
228  cit = itk::ConstNeighborhoodIterator<InputImageType>(m_FaceList, input, *fit);
229 
231  itdir = itk::ImageRegionIterator<OutputImageType>(outputDir, *fit);
232 
233  bit.OverrideBoundaryCondition(&nbc);
234  bit.GoToBegin();
235  cit.OverrideBoundaryCondition(&nbc);
236  cit.GoToBegin();
237 
238  otbMsgDevMacro(<< " ------------------- FaceList --------------------------");
239 
240  while ((!bit.IsAtEnd()) && (!cit.IsAtEnd()))
241  {
242  InterpolatorPointer interpolator = InterpolatorType::New();
243  // Location of the central pixel of the region
244  off.Fill(0);
245  bitIndex = bit.GetIndex(off);
246  Xc = bitIndex[0];
247  Yc = bitIndex[1];
248 
249  // JULIEN : If the processed region is the center face
250  // the input image can be used for the interpolation
251  //if (interiorFace)
252  // {
253  // interpolator->SetInputImage(input);
254  // }
255  // else we must feed the interpolator with a partial image corresponding
256  // to the boundary conditions
257  //else
258  // {
259  typename InputImageType::RegionType tempRegion;
260  typename InputImageType::SizeType tempSize;
261  tempSize[0] = 2 * m_FaceList[0] + 1;
262  tempSize[1] = 2 * m_FaceList[1] + 1;
263  tempRegion.SetSize(tempSize);
265  tempIndex[0] = off[0] - m_FaceList[0];
266  tempIndex[1] = off[1] - m_FaceList[1];
267  tempRegion.SetIndex(cit.GetIndex(tempIndex));
268  typename InputImageType::Pointer tempImage = InputImageType::New();
269  tempImage->SetRegions(tempRegion);
270  tempImage->Allocate();
271 
272  for (unsigned int p = 0; p <= 2 * m_FaceList[0]; ++p)
273  {
274  for (unsigned int q = 0; q <= 2 * m_FaceList[1]; q++)
275  {
277  index[0] = p - m_FaceList[0];
278  index[1] = q - m_FaceList[1];
279  tempImage->SetPixel(cit.GetIndex(index), cit.GetPixel(index));
280  }
281  }
282  interpolator->SetInputImage(tempImage);
283  // }
284 
285  // Location of the central pixel between zone 1 and zone 2
286  Yc12 = Yc - m_WidthLine - 1;
287 
288  // Location of the central pixel between zone 1 and zone 3
289  Yc13 = Yc + m_WidthLine + 1;
290 
291  // Contains for the 4 directions the the pixels belonging to each zone
292  //std::vector<double> PixelValues[NB_DIR][NB_ZONE];
293  // ROMAIN
294  std::vector<double>** PixelValues = NULL;
295  PixelValues = new std::vector<double>*[NB_DIR];
296  for (unsigned int i = 0; i < NB_DIR; ++i)
297  {
298  PixelValues[i] = NULL;
299  PixelValues[i] = new std::vector<double>[NB_ZONE];
300  }
301  //otbMsgDevMacro( << "\tCentre Xc/Yc="<<Xc<<" "<<Yc<<" Yc12/Yc13="<<Yc12<<" "<<Yc13);
302  // Loop on the region
303  for (unsigned int i = 0; i < m_Radius[0]; ++i)
304  for (unsigned int j = 0; j < m_Radius[1]; ++j)
305  {
306 
307  off[0] = i - m_Radius[0] / 2;
308  off[1] = j - m_Radius[1] / 2;
309 
310  bitIndex = bit.GetIndex(off);
311  X = bitIndex[0];
312  Y = bitIndex[1];
313 
314  // We determine in the horizontal direction with which zone the pixel belongs.
315  if (Y < Yc12) zone = 1;
316  else if ((Yc12 < Y) && (Y < Yc13)) zone = 0;
317  else if (Y > Yc13) zone = 2;
318  else continue;
319  //otbMsgDevMacro( << "\t\tPoint traite (i, j)=("<<i<<","<<j<<") -> X, Y="<<X<<","<<Y<<" zone="<<zone);
320  // Loop on the directions
321  for (unsigned int dir = 0; dir < NB_DIR; ++dir)
322  {
323  //ROTATION( (X-Xc), (Y-Yc), Theta[dir], xout, yout);
324 
325  xout = (X - Xc) * vcl_cos(Theta[dir]) - (Y - Yc) * vcl_sin(Theta[dir]);
326  yout = (X - Xc) * vcl_sin(Theta[dir]) + (Y - Yc) * vcl_cos(Theta[dir]);
327 
328  Index[0] = static_cast<CoordRepType>(xout + Xc);
329  Index[1] = static_cast<CoordRepType>(yout + Yc);
330 
331  PixelValues[dir][zone].push_back(static_cast<double>(interpolator->EvaluateAtContinuousIndex(Index)));
332  }
333  } // end of the loop on the pixels of the region
334 
335  R = 0.;
336  Direction = 0.;
337 
338  // Loop on the 4 directions
339 
340  for (unsigned int dir = 0; dir < NB_DIR; ++dir)
341  {
342 
343  double Rtemp = this->ComputeMeasure(&PixelValues[dir][0], &PixelValues[dir][1], &PixelValues[dir][2]);
344 
345  if (Rtemp > R)
346  {
347  R = Rtemp;
348  Direction = Theta[dir];
349  }
350 
351  } // end of the loop on the directions
352 
353  //otbMsgDevMacro( << "\t\tR, Direction : "<<R<<","<<Direction);
354  if (R >= this->GetThreshold())
355  {
356 
357  // Assignment of this value to the output pixel
358  it.Set(static_cast<OutputPixelType>(R));
359 
360  // Assignment of this value to the "outputdir" pixel
361  itdir.Set(static_cast<OutputPixelType>(Direction));
362  }
363  else
364  {
365 
367 
368  itdir.Set(static_cast<OutputPixelType>(0));
369  }
370  ++bit;
371  ++cit;
372  ++it;
373  ++itdir;
374  //interiorFace = false;
375  progress.CompletedPixel();
376 
377  // ROMAIN
378  for (unsigned int i = 0; i < NB_DIR; ++i)
379  {
380  delete[] PixelValues[i];
381  PixelValues[i] = NULL;
382  }
383  delete[] PixelValues;
384  PixelValues = NULL;
385  }
386 
387  }
388  delete[] Theta;
389 }
390 
391 template <class TInputImage, class TOutputImage, class TOutputImageDirection, class InterpolatorType>
392 double
394 ::ComputeMeasure(std::vector<double>* itkNotUsed(m1), std::vector<double>* itkNotUsed(m2), std::vector<double>* itkNotUsed(m3))
395 {
396  return 0;
397 }
398 
402 template <class TInputImage, class TOutputImage, class TOutputImageDirection, class InterpolatorType>
403 void
405 ::PrintSelf(std::ostream& os, itk::Indent indent) const
406 {
407  Superclass::PrintSelf(os, indent);
408  os << indent << "Length: " << m_LengthLine << std::endl;
409  os << indent << "Width: " << m_WidthLine << std::endl;
411 
412 }
413 
414 } // end namespace otb
415 
416 #endif
virtual IndexType GetIndex(void) const
const double CONST_PI
Definition: otbMath.h:45
virtual bool IsAtEnd() const
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
void Set(const PixelType &value) const
virtual void OverrideBoundaryCondition(const ImageBoundaryConditionPointerType i)
InterpolatorType::CoordRepType CoordRepType
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
void Fill(OffsetValueType value)
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId)
void PrintSelf(std::ostream &os, itk::Indent indent) const
virtual double ComputeMeasure(std::vector< double > *m1, std::vector< double > *m2, std::vector< double > *m3)
OutputImageType::RegionType OutputImageRegionType
#define NULL
unsigned int ThreadIdType
#define otbMsgDevMacro(x)
Definition: otbMacro.h:95