Umasoft
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
FelkelComponents.h
1 //
2 // FelkelComponents.h
3 //
4 // Copyright (c) 2003-2011 Virtual Terrain Project
5 // Free for all uses, see license.txt for details.
6 //
7 // Straight skeleton algorithm and original implementation
8 // courtesy of Petr Felkel and Stepan Obdrzalek (petr.felkel@tiani.com)
9 // Re-implemented for the Virtual Terrain Project (vterrain.org)
10 // by Roger James (www.beardandsandals.co.uk)
11 //
12 
13 #ifndef FELKELCOMPONENTSH
14 #define FELKELCOMPONENTSH
15 
16 #include <math.h>
17 #include <assert.h>
18 #include <vector>
19 #include <list>
20 
21 #include "vtdata/vtLog.h"
22 #include "vtdata/MathTypes.h"
23 
24 // forward declarations for CEdge constructor
25 class vtString;
26 class vtMaterial;
27 
28 #ifdef _MSC_VER
29 #pragma warning(disable: 4786) // prevent common warning about templates
30 #include <limits>
31 #endif
32 
33 using namespace std;
34 
35 #define CN_PI ((CNumber ) 3.14159265358979323846)
36 #ifdef _MSC_VER
37 //#define CN_INFINITY ((CNumber ) std::numeric_limits<double>::infinity())
38 #define CN_INFINITY ((CNumber ) DBL_MAX)
39 #else
40 #define CN_INFINITY ((CNumber ) 1.797693E+308)
41 #endif
42 #define CN_SLOPE_MAX ((CNumber ) CN_PI * 89 /180)
43 #define CN_SLOPE_MIN ((CNumber ) CN_PI / 180)
44 
45 #define MIN_DIFF 0.00005 // For vtp this gives about a 5 cm resolution
46 #define MIN_ANGLE_DIFF 0.00005
47 
48 #define SIMILAR(a,b) ((a)-(b) < MIN_DIFF && (b)-(a) < MIN_DIFF)
49 #define ANGLE_SIMILAR(a,b) (a.NormalizedAngle() - b.NormalizedAngle() < MIN_ANGLE_DIFF && b.NormalizedAngle() - a.NormalizedAngle() < MIN_ANGLE_DIFF)
50 
51 
52 class CNumber
53 {
54 public:
55  // Constructor
56  CNumber (double x = 0.0);
57  // Operator overloads
58  // operator const double& (void) const { return n; }
59  CNumber& operator = (const CNumber &x) { m_n = x.m_n; return *this; }
60  CNumber& operator = (const double &x) { m_n = x; return *this; }
61  operator double& (void) const { return (double &) m_n; }
62  bool operator == (const CNumber &x) const { return SIMILAR (m_n, x.m_n); }
63  bool operator != (const CNumber &x) const { return !SIMILAR (m_n, x.m_n); }
64  bool operator <= (const CNumber &x) const { return m_n < x.m_n || *this == x; }
65  bool operator >= (const CNumber &x) const { return m_n > x.m_n || *this == x; }
66  bool operator < (const CNumber &x) const { return m_n < x.m_n && *this != x; }
67  bool operator > (const CNumber &x) const { return m_n > x.m_n && *this != x; }
68 
69  bool operator == (const double x) const { return *this == CNumber (x); }
70  bool operator != (const double x) const { return *this != CNumber (x); }
71  bool operator <= (const double x) const { return *this <= CNumber (x); }
72  bool operator >= (const double x) const { return *this >= CNumber (x); }
73  bool operator < (const double x) const { return *this < CNumber (x); }
74  bool operator > (const double x) const { return *this > CNumber (x); }
75  // Functions
76  CNumber &NormalizeAngle();
77  CNumber NormalizedAngle();
78 private:
79  double m_n;
80 };
81 
82 // This is not really a 3D class beware !!!!!!
83 // Most things only work on x and z co-ords
84 
85 // This is not really a 3D class beware !!!!!!
86 // Most things only work on x and z co-ords
87 
88 class C3DPoint
89 {
90 public:
91  C3DPoint (void) : m_x(0), m_y(0), m_z(0) { }
92  C3DPoint (CNumber X, CNumber Y, CNumber Z) : m_x (X), m_y (Y), m_z(Z) { }
93  bool operator == (const C3DPoint &p) const
94  {
95  if (m_x == p.m_x && m_z == p.m_z && m_y == p.m_y)
96  return true;
97  else
98  {
99 #ifdef FELKELDEBUG
100  if (m_x == p.m_x && m_z == p.m_z)
101  VTLOG("C3DPoint - 2D equality not maintained y1 %e y2 %e\n", m_y, p.m_y);
102 #endif
103  return false;
104  }
105  }
106  bool operator != (const C3DPoint &p) const
107  {
108  if (m_x != p.m_x || m_z != p.m_z || m_y != p.m_y)
109  {
110 #ifdef FELKELDEBUG
111  if (m_x == p.m_x && m_z == p.m_z)
112  VTLOG("C3DPoint - 2D inequality not maintained y1 %e y2 %e\n", m_y, p.m_y);
113 #endif
114  return true;
115  }
116  else
117  return false;
118  }
119  C3DPoint operator - (const C3DPoint &p) const {return C3DPoint(m_x - p.m_x, m_y, m_z - p.m_z);}
120  C3DPoint operator + (const C3DPoint &p) const {return C3DPoint(m_x + p.m_x, m_y, m_z + p.m_z);}
121  C3DPoint operator * (const CNumber &n) const { return C3DPoint (n*m_x, m_y, n*m_z); }
122  bool IsInfiniteXZ (void) { return *this == C3DPoint (CN_INFINITY, CN_INFINITY, CN_INFINITY) ? true : false; }
123  CNumber DistXZ(const C3DPoint &p) const {return sqrt ((m_x-p.m_x)*(m_x-p.m_x) + (m_z - p.m_z)*(m_z - p.m_z));}
124  CNumber LengthXZ() {return sqrt (m_x * m_x + m_z * m_z);}
125  CNumber DotXZ(const C3DPoint &p) const {return m_x * p.m_x + m_z * p.m_z; }
126  CNumber CrossXZ(const C3DPoint &p) {return (m_x * p.m_z - m_z * p.m_x);} // 2D pseudo cross
127 
128  CNumber m_x, m_y, m_z;
129 };
130 
131 class CEdge
132 {
133 public:
134  CEdge (CNumber X, CNumber Y, CNumber Z, CNumber Slope,
135  const vtString *pMaterial, RGBi Color) :
136  m_Point(X, Y, Z), m_Slope(Slope), m_pMaterial(pMaterial), m_Color(Color) {};
137  bool operator == (const CEdge &p) const { return m_Point == p.m_Point; }
138  C3DPoint m_Point;
139  CNumber m_Slope;
140  const vtString *m_pMaterial;
141  RGBi m_Color;
142 };
143 
144 typedef vector <CEdge> Contour;
145 typedef vector <Contour> ContourVector;
146 
148 {
149 public:
150  CRidgeLine(const C3DPoint &p = C3DPoint (0, 0, 0), const C3DPoint &q = C3DPoint(0, 0, 0), const CNumber &Slope = -1, const bool IsRidgeLine = false);
151  CRidgeLine(const C3DPoint &p, const CNumber &a, const CNumber &Slope, const bool IsRidgeLine = true) : m_Origin (p), m_Angle (a), m_Slope(Slope), m_IsRidgeLine(IsRidgeLine)
152  {
153  m_Angle.NormalizeAngle();
154  };
155  CRidgeLine Opaque(void) const { return CRidgeLine (m_Origin, m_Angle + CN_PI, m_Slope, m_IsRidgeLine); }
156  static CRidgeLine AngleAxis (const C3DPoint &b, const C3DPoint &a, const CNumber &sa, const C3DPoint &c, const CNumber &sc);
157  C3DPoint Intersection(const CRidgeLine &a);
158  C3DPoint IntersectionAnywhere (const CRidgeLine& a) const;
159  bool Colinear(const CRidgeLine &a) const;
160  inline bool PointOnRidgeLine (const C3DPoint &p) const
161  {
162  return (p == m_Origin || CRidgeLine(m_Origin, p, -1).m_Angle == m_Angle) ? true : false;
163  }
164  inline bool FacingTowards (const CRidgeLine &a) const
165  {
166 // return (a.PointOnRidgeLine(m_Origin) && PointOnRidgeLine(a.m_Origin) && !(m_Origin == a.m_Origin)) ? true : false;
167  if (m_IsRidgeLine != a.m_IsRidgeLine)
168  VTLOG("%s %d m_IsRidgeLine %d\n", __FILE__, __LINE__, m_IsRidgeLine);
169  if (m_IsRidgeLine)
170  {
171  if (a.PointOnRidgeLine(m_Origin) && PointOnRidgeLine(a.m_Origin) && !(m_Origin == a.m_Origin) && SIMILAR(m_Slope, 0.0) && SIMILAR(m_Slope, 0.0))
172  {
173  VTLOG("Forcing return false\n");
174  }
175  return false;
176  }
177  else
178  return (a.PointOnRidgeLine(m_Origin) && PointOnRidgeLine(a.m_Origin) && !(m_Origin == a.m_Origin)) ? true : false;
179  }
180  CNumber Dist(const C3DPoint &p) const;
181 
182 
183  C3DPoint m_Origin;
184  CNumber m_Angle;
185  CNumber m_Slope;
186  bool m_IsRidgeLine;
187 };
188 
189 class CSkeletonLine;
190 class CVertexList;
191 class CIntersection;
192 
193 class CVertex
194 {
195 public:
196  CVertex (void) : m_ID (-1) { };
197  CVertex (const C3DPoint &p, const C3DPoint &prev = C3DPoint (), const CNumber &prevslope = -1, const C3DPoint &next = C3DPoint (), const CNumber &nextslope = -1)
198  : m_point (p), m_axis (CRidgeLine::AngleAxis (p, prev, prevslope, next, nextslope)), m_leftLine (p, prev, prevslope), m_rightLine (p, next, nextslope), m_higher (NULL),
199  m_leftVertex (NULL), m_rightVertex (NULL), m_nextVertex (NULL), m_prevVertex (NULL), m_done (false), m_ID (-1),
200  m_leftSkeletonLine (NULL), m_rightSkeletonLine (NULL), m_advancingSkeletonLine (NULL) { }
201  CVertex (const C3DPoint &p, CVertex &left, CVertex &right);
202  CVertex *Highest (void) { return m_higher ? m_higher -> Highest () : this; }
203  bool AtContour (void) const { return m_leftVertex == this && m_rightVertex == this; }
204  bool operator == (const CVertex &v) const { return m_point == v.m_point; }
205  bool operator < (const CVertex &) const { VTLOG("%s %d Assert failed\n", __FILE__, __LINE__); return false; }
206  C3DPoint CoordinatesOfAnyIntersectionOfTypeB(const CVertex &left, const CVertex &right);
207  C3DPoint IntersectionOfTypeB(const CVertex &left, const CVertex &right);
208  CNumber NearestIntersection (CVertexList &vl, CVertex **left, CVertex **right, C3DPoint &p);
209  bool InvalidIntersection(CVertexList &vl, const CIntersection &is);
210  bool VertexInCurrentContour(CVertex& Vertex);
211  // data
212  C3DPoint m_point;
213  CRidgeLine m_axis; // the axis (ridgeline) for this vertex
214  CRidgeLine m_leftLine, m_rightLine; // vectors for the original left and right edge contour lines
215  CVertex *m_leftVertex, *m_rightVertex; // Current contour chain (List of active vertices)
216  CVertex *m_nextVertex, *m_prevVertex; // Overall vertex list
217  CVertex *m_higher; // chain to next higher point in the skeleton (set when intersection is applied)
218  bool m_done;
219  int m_ID;
220 
221  // The following two fields are used to fill in the skeleton line structure when an intersection has been computed
222  CSkeletonLine *m_leftSkeletonLine;
223  CSkeletonLine *m_rightSkeletonLine;
224  CSkeletonLine *m_advancingSkeletonLine; // This field is filled in when an intersection is computed but is not used
225 };
226 
227 class CVertexList : public list <CVertex>
228 {
229 public:
230  CVertexList (void) { }
231  iterator prev (const iterator &i) { iterator tmp (i); if (tmp == begin ()) tmp = end (); tmp --; return tmp; }
232  iterator next (const iterator &i) { iterator tmp (i); tmp ++; if (tmp == end ()) tmp = begin (); return tmp; }
233  void push_back (const CVertex& x)
234  {
235 #if VTDEBUG
236  if (!(x.m_prevVertex == NULL || x.m_leftLine.FacingTowards (x.m_prevVertex -> m_rightLine)))
237  VTLOG("%s %d Assert failed\n", __FILE__, __LINE__);
238  if (!(x.m_nextVertex == NULL || x.m_rightLine.FacingTowards (x.m_nextVertex -> m_leftLine)))
239  VTLOG("%s %d Assert failed\n", __FILE__, __LINE__);
240 #endif
241  ((CVertex &)x).m_ID = size (); // automatic ID generation
242  list <CVertex> :: push_back (x);
243 #ifdef FELKELDEBUG
244  VTLOG("Vertex %d x %f y %f z %f ridge angle %f ridge slope %f\n left %d right %d prev %d next %d added to list\n",
245  ((CVertex &)x).m_ID,
246  ((CVertex &)x).m_point.m_x, ((CVertex &)x).m_point.m_y, ((CVertex &)x).m_point.m_z,
247  ((CVertex &)x).m_axis.m_Angle, ((CVertex &)x).m_axis.m_Slope,
248  (((CVertex &)x).m_leftVertex == NULL) ? -999 : ((CVertex &)x).m_leftVertex->m_ID,
249  (((CVertex &)x).m_rightVertex == NULL) ? -999 : ((CVertex &)x).m_rightVertex->m_ID,
250  (((CVertex &)x).m_prevVertex == NULL) ? -999 : ((CVertex &)x).m_prevVertex->m_ID,
251  (((CVertex &)x).m_nextVertex == NULL) ? -999 : ((CVertex &)x).m_nextVertex->m_ID);
252 #endif
253  }
254  void Dump ()
255  {
256  VTLOG("Dumping Vertex list\n");
257  for (iterator i = begin(); i != end(); i++)
258  {
259  VTLOG("Vertex %d x %e y %e z %e ridge angle %e ridge slope %e left %d right %d prev %d next\n",
260  (*i).m_ID,
261  (double) (*i).m_point.m_x, (double) (*i).m_point.m_y, (double) (*i).m_point.m_z,
262  (double) (*i).m_axis.m_Angle, (double) (*i).m_axis.m_Slope,
263  ((*i).m_leftVertex == NULL) ? -999 : (*i).m_leftVertex->m_ID,
264  ((*i).m_rightVertex == NULL) ? -999 : (*i).m_rightVertex->m_ID,
265  ((*i).m_prevVertex == NULL) ? -999 : (*i).m_prevVertex->m_ID,
266  ((*i).m_nextVertex == NULL) ? -999 : (*i).m_nextVertex->m_ID);
267  }
268  VTLOG("Vertex list dump complete\n");
269  }
270 };
271 
272 class CSegment
273 {
274 public:
275  CSegment (const C3DPoint &p = C3DPoint(), const C3DPoint &q = C3DPoint()) : m_a(p), m_b(q) { };
276  CSegment (CNumber x1, CNumber z1, CNumber x2, CNumber z2) : m_a (x1, 0, z1), m_b (x2, 0, z2) { }
277  CNumber Dist(const C3DPoint &p);
278 
279 private:
280  C3DPoint m_a;
281  C3DPoint m_b;
282 };
283 
285 {
286 public:
287  CSkeletonLine (void) : m_ID (-1) { }
288  CSkeletonLine (const CVertex &l, const CVertex &h) : m_lower (l), m_higher (h), m_ID (-1) { };
289  operator CSegment (void) { return CSegment (m_lower.m_vertex -> m_point, m_higher.m_vertex -> m_point); }
290 
292  {
293  SkeletonPoint (const CVertex &v = CVertex (), CSkeletonLine *l = NULL, CSkeletonLine *r = NULL) : m_vertex (&v), m_left (l), m_right (r) { }
294  const CVertex *m_vertex;
295  CSkeletonLine *m_left, *m_right;
296  int LeftID (void) const { if (!m_left) return -1; return m_left -> m_ID; }
297  int RightID (void) const { if (!m_right) return -1; return m_right -> m_ID; }
298  int VertexID (void) const { if (!m_vertex) return -1; return m_vertex -> m_ID; }
299  } m_lower, m_higher;
300  bool operator == (const CSkeletonLine &s) const
301  {
302  return m_higher.m_vertex -> m_ID == s.m_higher.m_vertex -> m_ID && m_lower.m_vertex -> m_ID == s.m_lower.m_vertex -> m_ID ;
303  }
304  bool operator < (const CSkeletonLine &) const { VTLOG("%s %d Assert failed\n", __FILE__, __LINE__); return false; }
305  int m_ID;
306 };
307 
308 
309 class CSkeleton : public list <CSkeletonLine>
310 {
311 public:
312  void push_back (const CSkeletonLine &x)
313  {
314  ((CSkeletonLine &)x).m_ID = size (); // automatically assign the ID number
315 #ifdef FELKELDEBUG
316  {
317  VTLOG("New skeleton line %d lower %d higher %d\n",
318  ((CSkeletonLine &)x).m_ID,
319  ((CSkeletonLine &)x).m_lower.m_vertex->m_ID,
320  ((CSkeletonLine &)x).m_higher.m_vertex->m_ID);
321  }
322 #endif
323  list <CSkeletonLine> :: push_back (x);
324  }
325 };
326 
327 #endif // FELKELCOMPONENTSH