Orfeo Toolbox  4.0
otbVectorDataToRightAngleVectorDataFilter.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 __otbVectorDataToRightAngleVectorDataFilter_txx
19 #define __otbVectorDataToRightAngleVectorDataFilter_txx
20 
22 #include "otbVectorData.h"
23 
24 namespace otb
25 {
26 template <class TVectorData>
29 {
30  this->SetNumberOfRequiredInputs(1);
31  this->SetNumberOfRequiredOutputs(1);
32 
33  m_DistanceThreshold = 20.;
34  m_AngleThreshold = CONST_PI / 30.; //we want a threshold at 6 degrees
35 }
36 
37 template <class TVectorData>
38 void
41 {
42  // Get the input segments
43  typename VectorDataType::Pointer vData = const_cast<VectorDataType *>(this->GetInput());
44 
45  // Output
46  this->GetOutput(0)->SetMetaDataDictionary(vData->GetMetaDataDictionary());
47 
48  // Retrieving root node
49  typename DataNodeType::Pointer root = this->GetOutput(0)->GetDataTree()->GetRoot()->Get();
50  // Create the document node
51  typename DataNodeType::Pointer document = DataNodeType::New();
52  document->SetNodeType(otb::DOCUMENT);
53  // Adding the layer to the data tree
54  this->GetOutput(0)->GetDataTree()->Add(document, root);
55  // Create the folder node
56  typename DataNodeType::Pointer folder = DataNodeType::New();
57  folder->SetNodeType(otb::FOLDER);
58  this->GetOutput(0)->GetDataTree()->Add(folder, document);
59  this->GetOutput(0)->SetProjectionRef(vData->GetProjectionRef());
60 
61  // Itterate on the vector data
62  TreeIteratorType itVectorRef(vData->GetDataTree()); // Reference
63  itVectorRef.GoToBegin();
64 
65  while (!itVectorRef.IsAtEnd())
66  {
67  if (!itVectorRef.Get()->IsLineFeature())
68  {
69  ++itVectorRef;
70  continue; // do not process if it's not a line
71  }
72  TreeIteratorType itVectorCur = itVectorRef; // Current
73 
74  while (!itVectorCur.IsAtEnd())
75  {
76  if (!itVectorCur.Get()->IsLineFeature())
77  {
78  ++itVectorCur;
79  continue; // do not process if it's not a line
80  }
81  // Compute the angle formed by the two segments
82  double Angle = this->ComputeAngleFormedBySegments(itVectorRef.Get()->GetLine(),
83  itVectorCur.Get()->GetLine());
84 
85  // Check if the angle is a right one
86  if (vcl_abs(Angle - CONST_PI_2) <= m_AngleThreshold)
87  {
88  // Right angle coordinate
89  PointType RightAngleCoordinate;
90  RightAngleCoordinate = this->ComputeRightAngleCoordinate(itVectorRef.Get()->GetLine(),
91  itVectorCur.Get()->GetLine());
92 
93  // Compute the distance between the two segments and the right angle formed by this segments
94  double dist1_2 = this->ComputeDistanceFromPointToSegment(RightAngleCoordinate,
95  itVectorRef.Get()->GetLine());
96  double dist2_2 = this->ComputeDistanceFromPointToSegment(RightAngleCoordinate,
97  itVectorCur.Get()->GetLine());
98 
99  double threshold_2 = m_DistanceThreshold * m_DistanceThreshold;
100  if (dist1_2 < threshold_2 && dist2_2 < threshold_2)
101  {
102  // If Right Angle & not so far from segments: Add to the output
103  typename DataNodeType::Pointer CurrentGeometry = DataNodeType::New();
104  CurrentGeometry->SetNodeId("FEATURE_POINT");
105  CurrentGeometry->SetNodeType(otb::FEATURE_POINT);
106  CurrentGeometry->SetPoint(RightAngleCoordinate);
107  this->GetOutput(0)->GetDataTree()->Add(CurrentGeometry, folder);
108  }
109  }
110 
111  ++itVectorCur;
112  }
113  ++itVectorRef;
114  }
115 }
116 
117 template <class TVectorData>
118 double
121 {
122  VertexListType * vertexList = const_cast<VertexListType *>(line->GetVertexList());
123 
124  VertexType v0 = vertexList->GetElement(0);
125  VertexType v1 = vertexList->GetElement(1);
126 
127  // Case the extremities of the segment are the same
128  double l2 = v0.SquaredEuclideanDistanceTo(v1);
129 
130  // Is the projection of rAngle on the segment inside (0<u<1) or
131  // inside the segment bounds
132  double u = ((rAngle[0] - v0[0]) *(v1[0] - v0[0] ) +
133  (rAngle[1] - v0[1]) *(v1[1] - v0[1])) / l2;
134 
135  if( u < 1e-10 ) u = 0.;
136  if( u -1. > 1e-10 ) u = 1.;
137 
138  double x = v0[0] + u *(v1[0] - v0[0] );
139  double y = v0[1] + u *(v1[1] - v0[1] );
140 
141  double dx = x - rAngle[0];
142  double dy = y - rAngle[1];
143 
144  return dx*dx + dy*dy;
145 }
146 
147 template <class TVectorData>
148 double
151 {
152  double oriDst = this->ComputeOrientation(lineDst);
153  double oriSrc = this->ComputeOrientation(lineSrc);
154 
155  return vcl_abs(oriDst - oriSrc);
156 }
157 
158 template <class TVectorData>
159 double
162 {
163  VertexListType * vertexList = const_cast<VertexListType *>(line->GetVertexList());
164 
165  double Xp1 = vertexList->GetElement(0)[0];
166  double Yp1 = vertexList->GetElement(0)[1];
167 
168  double Xp2 = vertexList->GetElement(1)[0];
169  double Yp2 = vertexList->GetElement(1)[1];
170 
171  //Compute the orientation
172  double dx = Xp1 - Xp2;
173  double dy = Yp1 - Yp2;
174  double orientation = vcl_atan2(dy, dx);
175  if (orientation < 0) orientation += CONST_PI;
176 
177  return orientation;
178 }
179 
180 template <class TVectorData>
182 ::PointType
185 {
186  PointType P;
187 
188  VertexListType * vertexListSrc = const_cast<VertexListType *>(lineSrc->GetVertexList());
189  double Xp1 = vertexListSrc->GetElement(0)[0];
190  double Yp1 = vertexListSrc->GetElement(0)[1];
191  double Xp2 = vertexListSrc->GetElement(1)[0];
192  double Yp2 = vertexListSrc->GetElement(1)[1];
193 
194  VertexListType * vertexListDst = const_cast<VertexListType *>(lineDst->GetVertexList());
195  double Xq1 = vertexListDst->GetElement(0)[0];
196  double Yq1 = vertexListDst->GetElement(0)[1];
197  double Xq2 = vertexListDst->GetElement(1)[0];
198  double Yq2 = vertexListDst->GetElement(1)[1];
199 
200  double numU = (Xq2 - Xq1)*(Yp1-Yq1) - (Yq2-Yq1)*(Xp1-Xq1);
201  double denumU = (Yq2-Yq1)*(Xp2-Xp1)-(Xq2-Xq1)*(Yp2-Yp1);
202 
203  double u = numU/ denumU;
204 
205  P[0] = Xp1 + u*(Xp2-Xp1);
206  P[1] = Yp1 + u * (Yp2-Yp1);
207 
208  return P;
209 }
210 
211 template <class TVectorData>
212 void
214 ::PrintSelf(std::ostream& os, itk::Indent indent) const
215 {
216  Superclass::PrintSelf(os, indent);
217 }
218 
219 } // end namespace otb
220 
221 #endif

Generated at Sat Mar 8 2014 16:24:44 for Orfeo Toolbox with doxygen 1.8.3.1