triangulate.hh
Go to the documentation of this file.
/* -*- mia-c++ -*-
*
* This file is part of MIA - a toolbox for medical image analysis
* Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
*
* MIA is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MIA; if not, see <http://www.gnu.org/licenses/>.
*
*/
#include <mia/mesh/clist.hh>
#include <mia/core/msgstream.hh>
#include <iostream>
NS_MIA_BEGIN
template <class VertexVector, class Polygon>
class TPolyTriangulator
{
public:
TPolyTriangulator(const VertexVector& vv);
template <class TriangleList>
bool triangulate(TriangleList& output, const Polygon& poly) const;
private:
typedef clist<typename Polygon::value_type> CPoly;
typedef typename VertexVector::value_type Vector;
Vector eval_orientation(const Polygon& poly) const;
bool is_convex(const typename CPoly::const_iterator& i, bool debug = false) const;
bool is_ear(const typename CPoly::const_iterator& p, const CPoly& cpoly, bool debug = false) const;
bool is_inside(
const typename VertexVector::value_type& a,
const typename VertexVector::value_type& b,
const typename VertexVector::value_type& c,
const typename VertexVector::value_type& p,
bool debug = false ) const;
const VertexVector& m_vv;
mutable Vector m_orientation;
};
template <class VertexVector, class Polygon>
TPolyTriangulator<VertexVector, Polygon>::TPolyTriangulator(const VertexVector& vv):
m_vv(vv)
{
}
template <class VertexVector, class Polygon>
template <class TriangleList>
bool
TPolyTriangulator<VertexVector, Polygon>::triangulate(TriangleList& output, const Polygon& poly) const
{
size_t poly_size = poly.size();
if ( poly_size < 3) // no triangles at all
return false;
typedef clist<typename Polygon::value_type> CPoly;
CPoly cpoly;
m_orientation = eval_orientation(poly);
typename Polygon::const_iterator pi = poly.begin();
typename Polygon::const_iterator pe = poly.end();
while (pi != pe)
cpoly.push_back(*pi++);
typename CPoly::iterator p_i = cpoly.begin();
typename CPoly::iterator p_e = cpoly.end();
p_i = p_i->succ;
while ( (p_i != p_e) && (poly_size > 3)) {
if (is_ear(p_i, cpoly, false)) {
// we have a valid triangle, store it
output.push_back(typename TriangleList::value_type(**p_i->succ, **p_i, **p_i->prev));
// set the current middle node
typename CPoly::iterator p_r = p_i;
p_i = (p_i->prev != cpoly.begin()) ? p_i->prev : p_i->succ;
cpoly.remove(p_r);
--poly_size;
} else
p_i = p_i->succ;
}
if ((p_i == p_e) && (poly_size > 3)) {
cvdebug() << "gotcha\n";
p_i = cpoly.begin();
p_i = p_i->succ;
while ( (p_i != p_e) && (poly_size > 3)) {
if (is_ear(p_i, cpoly, true)) {
// we have a valid triangle, store it
output.push_back(typename TriangleList::value_type(**p_i->succ, **p_i, **p_i->prev));
// set the current middle node
typename CPoly::iterator p_r = p_i;
p_i = (p_i->prev != cpoly.begin()) ? p_i->prev : p_i->succ;
cpoly.remove(p_r);
--poly_size;
} else
p_i = p_i->succ;
}
}
output.push_back(typename TriangleList::value_type(**p_i->succ, **p_i, **p_i->prev));
return true;
}
template <class VertexVector, class Polygon>
typename TPolyTriangulator<VertexVector, Polygon>::Vector
TPolyTriangulator<VertexVector, Polygon>::eval_orientation(const Polygon& poly) const
{
typename VertexVector::value_type result(0, 0, 0);
typename Polygon::const_iterator pb = poly.begin();
typename Polygon::const_iterator be = poly.end();
typename Polygon::const_iterator c1 = pb;
++c1;
typename Polygon::const_iterator c2 = c1;
++c2;
typename VertexVector::value_type a = m_vv[*pb];
while (c1 != be && c2 != be) {
result += (m_vv[*c1++] - a) ^ (m_vv[*c2++] - a);
}
return result;
}
template <class VertexVector, class Polygon>
bool TPolyTriangulator<VertexVector, Polygon>::is_convex(const typename CPoly::const_iterator& i, bool /*debug*/) const
{
typename VertexVector::value_type a = m_vv[ **i->prev];
typename VertexVector::value_type b = m_vv[ **i];
typename VertexVector::value_type c = m_vv[ **i->succ];
typename VertexVector::value_type ab = a - b;
typename VertexVector::value_type cb = c - b;
typename VertexVector::value_type cross = ab ^ cb;
const bool result = dot(cross, m_orientation) > 0;
return result;
}
template <class VertexVector, class Polygon>
bool TPolyTriangulator<VertexVector, Polygon>::is_ear(const typename CPoly::const_iterator& p, const CPoly& cpoly, bool debug) const
{
if (!is_convex(p, debug)) {
cvdebug() << "corner is concave\n";
return false;
}
typename VertexVector::value_type a = m_vv[ **p->prev];
typename VertexVector::value_type b = m_vv[ **p];
typename VertexVector::value_type c = m_vv[ **p->succ];
cvdebug() << "check triangle" << a << b << c << " = (" << **p->prev << "," << **p << "," << **p->succ << "\n";
typename CPoly::const_iterator i = cpoly.begin();
i = i->succ;
while (i != cpoly.end()) {
if (i != p && i != p->prev && i != p->succ)
if (!is_convex(i, debug) && is_inside(a, b, c, m_vv[ **i], debug)) {
cvdebug() << "point " << **i << ":" << m_vv[ **i] << " is concave and inside\n";
return false;
}
i = i->succ;
}
return true;
}
template <class VertexVector, class Polygon>
bool TPolyTriangulator<VertexVector, Polygon>::is_inside(
const typename VertexVector::value_type& a,
const typename VertexVector::value_type& b,
const typename VertexVector::value_type& c,
const typename VertexVector::value_type& p,
bool /*debug*/) const
{
double abc = ((a - b) ^ (c - b)).norm() * 0.5;
double abp = ((a - p) ^ (b - p)).norm() * 0.5;
double acp = ((a - p) ^ (c - p)).norm() * 0.5;
double bcp = ((b - p) ^ (c - p)).norm() * 0.5;
const bool result = (abc >= abp + acp + bcp);
return result;
}
NS_MIA_END
cross
T cross(const T2DVector< T > &a, const T2DVector< T > &b)
Definition 2d/vector.hh:518
dot
T dot(const T2DVector< T > &a, const T2DVector< T > &b)
Definition 2d/vector.hh:437
TPolyTriangulator::triangulate
bool triangulate(TriangleList &output, const Polygon &poly) const
Definition triangulate.hh:88
TPolyTriangulator::TPolyTriangulator
TPolyTriangulator(const VertexVector &vv)
Definition triangulate.hh:80
clist
Definition clist.hh:30
clist::push_back
void push_back(T val)
Definition clist.hh:114
clist< typename Polygon::value_type >::iterator
node * iterator
Definition clist.hh:57
clist::end
iterator end()
Definition clist.hh:81
clist::begin
iterator begin()
Definition clist.hh:76
clist< typename Polygon::value_type >::const_iterator
const node * const_iterator
Definition clist.hh:58
clist::remove
void remove(node *n)
Definition clist.hh:97
clist.hh
CCmdOptionFlags::output
@ output
Definition cmdoptionflags.hh:32
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition defines.hh:33
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition defines.hh:36
msgstream.hh
cvdebug
CDebugSink & cvdebug()
Definition msgstream.hh:226