Orfeo Toolbox  4.0
otbFillGapsFilter.cxx
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 
19 #include "otbFillGapsFilter.h"
20 
21 namespace otb
22 {
23 
29 {
32 
34 
35  this->ProcessObjectType::SetNthOutput(0, output);
36 
37  m_Radius = 4.0;
38  m_AngularBeam = 1.0;
39 }
40 
41 void
44 {
45  this->ProcessObjectType::SetNthInput(0,
46  const_cast<LineSpatialObjectListType *>(input));
47 }
48 
52 {
53  return static_cast<const LineSpatialObjectListType *>
54  (this->ProcessObjectType::GetInput(0));
55 }
56 
60 {
61  return static_cast<LineSpatialObjectListType *>
62  (this->ProcessObjectType::GetOutput(0));
63 }
64 
69 void
72 {
73  // Get the LineSpatialObject
74 
75  const LineSpatialObjectList * inputLine = this->GetInput();
76  LineSpatialObjectList * outputLine = this->GetOutput();
77 
78  // Get the list of points which consists of two points to represent a
79  // straight line
80 
81  LineSpatialObjectListType::const_iterator itLineListA = inputLine->begin();
82  LineSpatialObjectListType::const_iterator itLineListAEnd = inputLine->end();
83  LineSpatialObjectListType::const_iterator itLineListB;
84  LineSpatialObjectListType::const_iterator itLineListBEnd = inputLine->end();
85 
86  double x1(0.), y1(0.), x2(0.), y2(0.), x3(0.), y3(0.), x4(0.), y4(0.);
87  double xTemp(0.), yTemp(0.);
88  double R13(0.), R14(0.), R23(0.), R24(0.);
89  double CosTheta(0.);
90 
91  PointListType::const_iterator itPoints;
92 
93  CosTheta = vcl_cos(m_AngularBeam);
94  --itLineListAEnd;
95 
96  while (itLineListA != itLineListAEnd)
97  {
98  itLineListB = itLineListA;
99  ++itLineListB;
100 
101  PointListType& pointsList = (*itLineListA)->GetPoints();
102  itPoints = pointsList.begin();
103 
104  x1 = (*itPoints).GetPosition()[0];
105  y1 = (*itPoints).GetPosition()[1];
106 
107  ++itPoints;
108  x2 = (*itPoints).GetPosition()[0];
109  y2 = (*itPoints).GetPosition()[1];
110 
111  PointType point;
112  PointListType pointList;
113 
114  point.SetPosition(x1, y1);
115  pointList.push_back(point);
116  point.SetPosition(x2, y2);
117  pointList.push_back(point);
118 
119  LineSpatialObjectType::Pointer line = LineSpatialObjectType::New();
120  line->SetId(0);
121  line->SetPoints(pointList);
122  line->ComputeBoundingBox();
123  outputLine->push_back(line);
124 
125  while (itLineListB != itLineListBEnd)
126  {
127 
128  pointsList = (*itLineListB)->GetPoints();
129  itPoints = pointsList.begin();
130 
131  x3 = (*itPoints).GetPosition()[0];
132  y3 = (*itPoints).GetPosition()[1];
133 
134  ++itPoints;
135  x4 = (*itPoints).GetPosition()[0];
136  y4 = (*itPoints).GetPosition()[1];
137 
138  // Calculate the radius for each point of each line
139 
140  R13 = vcl_sqrt((x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3));
141  R14 = vcl_sqrt((x1 - x4) * (x1 - x4) + (y1 - y4) * (y1 - y4));
142  R23 = vcl_sqrt((x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3));
143  R24 = vcl_sqrt((x2 - x4) * (x2 - x4) + (y2 - y4) * (y2 - y4));
144 
145  double Rmin = m_Radius;
146 
147  // Test the lower Radius
148  if (R13 < Rmin) Rmin = R13;
149  if (R14 < Rmin) Rmin = R14;
150  if (R23 < Rmin) Rmin = R23;
151  if (R24 < Rmin) Rmin = R24;
152 
153  if (Rmin < m_Radius)
154  {
155  // Sort Points such as the radius of P2 and P3 is the smallest one.
156  if (Rmin == R24)
157  {
158  xTemp = x3;
159  yTemp = y3;
160  x3 = x4;
161  y3 = y4;
162  x4 = xTemp;
163  y4 = yTemp;
164  }
165  if (Rmin == R13)
166  {
167  xTemp = x1;
168  yTemp = y1;
169  x1 = x2;
170  y1 = y2;
171  x2 = xTemp;
172  y2 = yTemp;
173  }
174  if (Rmin == R14)
175  {
176  xTemp = x3;
177  yTemp = y3;
178  x3 = x4;
179  y3 = y4;
180  x4 = xTemp;
181  y4 = yTemp;
182  xTemp = x1;
183  yTemp = y1;
184  x1 = x2;
185  y1 = y2;
186  x2 = xTemp;
187  y2 = yTemp;
188  }
189 
190  //Estimate the norm each line
191  /* double Norm12, Norm23, Norm34;
192  Norm12 = vcl_sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
193  Norm23 = vcl_sqrt( (x2-x3)*(x2-x3) + (y2-y3)*(y2-y3) );
194  Norm34 = vcl_sqrt( (x3-x4)*(x3-x4) + (y3-y4)*(y3-y4) );
195  */
196  double Angle12_23, Angle12_34, Angle23_34;
197  //Estimate the angle between lines 12-23 and lines 12-34
198  /*Angle12_23 = (x2-x1)*(x3-x2) + (y2-y1)*(y3-y2);
199  Angle12_23 = Angle12_23 / Norm12 / Norm23;
200 
201  Angle12_34 = (x2-x1)*(x4-x3) + (y2-y1)*(y4-y3);
202  Angle12_34 = Angle12_34 / Norm12 / Norm34; */
203 
204  Angle12_23 = vcl_cos(vcl_atan2((y2 - y1), (x2 - x1)) - vcl_atan2((y3 - y2), (x3 - x2)));
205  Angle12_34 = vcl_cos(vcl_atan2((y2 - y1), (x2 - x1)) - vcl_atan2((y4 - y3), (x4 - x3)));
206  Angle23_34 = vcl_cos(vcl_atan2((y3 - y2), (x3 - x2)) - vcl_atan2((y4 - y3), (x4 - x3)));
207 
208  if ((Angle12_23 > CosTheta) && (Angle12_34 > CosTheta) && (Angle23_34 > CosTheta))
209  {
210 
211  // Store 23-segment
212  PointType point;
213  PointListType pointList;
214 
215  point.SetPosition(x2, y2);
216  pointList.push_back(point);
217  point.SetPosition(x3, y3);
218  pointList.push_back(point);
219 
220  LineSpatialObjectType::Pointer line = LineSpatialObjectType::New();
221  line->SetId(0);
222  line->SetPoints(pointList);
223  line->ComputeBoundingBox();
224  outputLine->push_back(line);
225 
226  pointList.clear();
227 
228  }
229  } // if(Rmin < m_Radius)
230 
231  ++itLineListB;
232  } // while(itLineListB != itLineListBEnd)
233  ++itLineListA;
234  }
235 
236  // Insert the last element
237 
238  PointType point;
239  PointListType pointList;
240 
241  point.SetPosition(x3, y3);
242  pointList.push_back(point);
243  point.SetPosition(x4, y4);
244  pointList.push_back(point);
245 
246  LineSpatialObjectType::Pointer line = LineSpatialObjectType::New();
247  line->SetId(0);
248  line->SetPoints(pointList);
249  line->ComputeBoundingBox();
250  outputLine->push_back(line);
251 
252  pointList.clear();
253 
254 }
255 
260 void
262 ::PrintSelf(std::ostream& os, itk::Indent indent) const
263 {
264  Superclass::PrintSelf(os, indent);
265 }
266 
267 } // end namespace otb

Generated at Sat Mar 8 2014 15:55:44 for Orfeo Toolbox with doxygen 1.8.3.1