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
37template <class VertexVector, class Polygon>
39{
40public:
45 TPolyTriangulator(const VertexVector& vv);
46
56 template <class TriangleList>
57 bool triangulate(TriangleList& output, const Polygon& poly) const;
58
59private:
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
79template <class VertexVector, class Polygon>
81 m_vv(vv)
82{
83}
84
85template <class VertexVector, class Polygon>
86template <class TriangleList>
87bool
88TPolyTriangulator<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
146template <class VertexVector, class Polygon>
147typename 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
166template <class VertexVector, class Polygon>
167bool 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
179template <class VertexVector, class Polygon>
180bool 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
207template <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
T cross(const T2DVector< T > &a, const T2DVector< T > &b)
Definition 2d/vector.hh:518
T dot(const T2DVector< T > &a, const T2DVector< T > &b)
Definition 2d/vector.hh:437
class to make a triangle mesh from a closed polygon
bool triangulate(TriangleList &output, const Polygon &poly) const
TPolyTriangulator(const VertexVector &vv)
Definition clist.hh:30
void push_back(T val)
Definition clist.hh:114
iterator end()
Definition clist.hh:81
iterator begin()
Definition clist.hh:76
void remove(node *n)
Definition clist.hh:97
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition defines.hh:33
#define NS_MIA_END
conveniance define to end the mia namespace
Definition defines.hh:36
CDebugSink & cvdebug()
Definition msgstream.hh:226
node * succ
Definition clist.hh:39
node * prev
Definition clist.hh:38