My Project
triangulate.hh
Go to the documentation of this file.
1 /* -*- mia-c++ -*-
2  *
3  * This file is part of MIA - a toolbox for medical image analysis
4  * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
5  *
6  * MIA is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #include <mia/mesh/clist.hh>
22 #include <mia/core/msgstream.hh>
23 #include <iostream>
24 
26 
37 template <class VertexVector, class Polygon>
39 {
40 public:
45  TPolyTriangulator(const VertexVector& vv);
46 
56  template <class TriangleList>
57  bool triangulate(TriangleList& output, const Polygon& poly) const;
58 
59 private:
61  typedef typename VertexVector::value_type Vector;
62 
63  Vector eval_orientation(const Polygon& poly) const;
64  bool is_convex(const typename CPoly::const_iterator& i, bool debug = false) const;
65  bool is_ear(const typename CPoly::const_iterator& p, const CPoly& cpoly, bool debug = false) const;
66 
67  bool is_inside(
68  const typename VertexVector::value_type& a,
69  const typename VertexVector::value_type& b,
70  const typename VertexVector::value_type& c,
71  const typename VertexVector::value_type& p,
72  bool debug = false ) const;
73 
74 
75  const VertexVector& m_vv;
76  mutable Vector m_orientation;
77 };
78 
79 template <class VertexVector, class Polygon>
81  m_vv(vv)
82 {
83 }
84 
85 template <class VertexVector, class Polygon>
86 template <class TriangleList>
87 bool
88 TPolyTriangulator<VertexVector, Polygon>::triangulate(TriangleList& output, const Polygon& poly) const
89 {
90  size_t poly_size = poly.size();
91 
92  if ( poly_size < 3) // no triangles at all
93  return false;
94 
96  CPoly cpoly;
97  m_orientation = eval_orientation(poly);
98  typename Polygon::const_iterator pi = poly.begin();
99  typename Polygon::const_iterator pe = poly.end();
100 
101  while (pi != pe)
102  cpoly.push_back(*pi++);
103 
104  typename CPoly::iterator p_i = cpoly.begin();
105  typename CPoly::iterator p_e = cpoly.end();
106  p_i = p_i->succ;
107 
108  while ( (p_i != p_e) && (poly_size > 3)) {
109  if (is_ear(p_i, cpoly, false)) {
110  // we have a valid triangle, store it
111  output.push_back(typename TriangleList::value_type(**p_i->succ, **p_i, **p_i->prev));
112  // set the current middle node
113  typename CPoly::iterator p_r = p_i;
114  p_i = (p_i->prev != cpoly.begin()) ? p_i->prev : p_i->succ;
115  cpoly.remove(p_r);
116  --poly_size;
117  } else
118  p_i = p_i->succ;
119  }
120 
121  if ((p_i == p_e) && (poly_size > 3)) {
122  cvdebug() << "gotcha\n";
123  p_i = cpoly.begin();
124  p_i = p_i->succ;
125 
126  while ( (p_i != p_e) && (poly_size > 3)) {
127  if (is_ear(p_i, cpoly, true)) {
128  // we have a valid triangle, store it
129  output.push_back(typename TriangleList::value_type(**p_i->succ, **p_i, **p_i->prev));
130  // set the current middle node
131  typename CPoly::iterator p_r = p_i;
132  p_i = (p_i->prev != cpoly.begin()) ? p_i->prev : p_i->succ;
133  cpoly.remove(p_r);
134  --poly_size;
135  } else
136  p_i = p_i->succ;
137  }
138  }
139 
140  output.push_back(typename TriangleList::value_type(**p_i->succ, **p_i, **p_i->prev));
141  return true;
142 }
143 
144 
145 
146 template <class VertexVector, class Polygon>
147 typename TPolyTriangulator<VertexVector, Polygon>::Vector
149 {
150  typename VertexVector::value_type result(0, 0, 0);
151  typename Polygon::const_iterator pb = poly.begin();
152  typename Polygon::const_iterator be = poly.end();
153  typename Polygon::const_iterator c1 = pb;
154  ++c1;
155  typename Polygon::const_iterator c2 = c1;
156  ++c2;
157  typename VertexVector::value_type a = m_vv[*pb];
158 
159  while (c1 != be && c2 != be) {
160  result += (m_vv[*c1++] - a) ^ (m_vv[*c2++] - a);
161  }
162 
163  return result;
164 }
165 
166 template <class VertexVector, class Polygon>
167 bool TPolyTriangulator<VertexVector, Polygon>::is_convex(const typename CPoly::const_iterator& i, bool /*debug*/) const
168 {
169  typename VertexVector::value_type a = m_vv[ **i->prev];
170  typename VertexVector::value_type b = m_vv[ **i];
171  typename VertexVector::value_type c = m_vv[ **i->succ];
172  typename VertexVector::value_type ab = a - b;
173  typename VertexVector::value_type cb = c - b;
174  typename VertexVector::value_type cross = ab ^ cb;
175  const bool result = dot(cross, m_orientation) > 0;
176  return result;
177 }
178 
179 template <class VertexVector, class Polygon>
180 bool TPolyTriangulator<VertexVector, Polygon>::is_ear(const typename CPoly::const_iterator& p, const CPoly& cpoly, bool debug) const
181 {
182  if (!is_convex(p, debug)) {
183  cvdebug() << "corner is concave\n";
184  return false;
185  }
186 
187  typename VertexVector::value_type a = m_vv[ **p->prev];
188  typename VertexVector::value_type b = m_vv[ **p];
189  typename VertexVector::value_type c = m_vv[ **p->succ];
190  cvdebug() << "check triangle" << a << b << c << " = (" << **p->prev << "," << **p << "," << **p->succ << "\n";
191  typename CPoly::const_iterator i = cpoly.begin();
192  i = i->succ;
193 
194  while (i != cpoly.end()) {
195  if (i != p && i != p->prev && i != p->succ)
196  if (!is_convex(i, debug) && is_inside(a, b, c, m_vv[ **i], debug)) {
197  cvdebug() << "point " << **i << ":" << m_vv[ **i] << " is concave and inside\n";
198  return false;
199  }
200 
201  i = i->succ;
202  }
203 
204  return true;
205 }
206 
207 template <class VertexVector, class Polygon>
209  const typename VertexVector::value_type& a,
210  const typename VertexVector::value_type& b,
211  const typename VertexVector::value_type& c,
212  const typename VertexVector::value_type& p,
213  bool /*debug*/) const
214 {
215  double abc = ((a - b) ^ (c - b)).norm() * 0.5;
216  double abp = ((a - p) ^ (b - p)).norm() * 0.5;
217  double acp = ((a - p) ^ (c - p)).norm() * 0.5;
218  double bcp = ((b - p) ^ (c - p)).norm() * 0.5;
219  const bool result = (abc >= abp + acp + bcp);
220  return result;
221 }
222 
clist::node::prev
node * prev
Definition: clist.hh:38
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
cross
T cross(const T2DVector< T > &a, const T2DVector< T > &b)
Definition: 2d/vector.hh:518
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
clist.hh
msgstream.hh
cvdebug
CDebugSink & cvdebug()
Definition: msgstream.hh:226
gsl::dot
double EXPORT_GSL dot(const gsl_vector *lhs, const gsl_vector *rhs)
clist::end
iterator end()
Definition: clist.hh:81
clist::node::succ
node * succ
Definition: clist.hh:39
clist::push_back
void push_back(T val)
Definition: clist.hh:114
TPolyTriangulator::triangulate
bool triangulate(TriangleList &output, const Polygon &poly) const
Definition: triangulate.hh:88
clist::remove
void remove(node *n)
Definition: clist.hh:97
clist::node
Definition: clist.hh:36
clist
Definition: clist.hh:30
clist::begin
iterator begin()
Definition: clist.hh:76
TPolyTriangulator::TPolyTriangulator
TPolyTriangulator(const VertexVector &vv)
Definition: triangulate.hh:80
TPolyTriangulator
class to make a triangle mesh from a closed polygon
Definition: triangulate.hh:39