triquad.py

# ##### BEGIN GPL LICENSE BLOCK #####
#
#  This program 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 2
#  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 this program; if not, write to the Free Software Foundation,
#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# ##### END GPL LICENSE BLOCK #####

# <pep8 compliant>


from . import geom
import math
import random
from math import sqrt

# Points are 3-tuples or 2-tuples of reals: (x,y,z) or (x,y)
# Faces are lists of integers (vertex indices into coord lists)
# After triangulation/quadrangulation, the tris and quads will
# be tuples instead of lists.
# Vmaps are lists taking vertex index -> Point

TOL = 1e-7     # a tolerance for fuzzy equality
GTHRESH = 75   # threshold above which use greedy to _Quandrangulate
ANGFAC = 1.0   # weighting for angles in quad goodness measure
DEGFAC = 10.0  # weighting for degree in quad goodness measure

# Angle kind constants
Ang0 = 1
Angconvex = 2
Angreflex = 3
Angtangential = 4
Ang360 = 5


def TriangulateFace(face, points):
    """Triangulate the given face.

    Uses an easy triangulation first, followed by a constrained delauney
    triangulation to get better shaped triangles.

    Args:
      face: list of int - indices in points, assumed CCW-oriented
      points: geom.Points - holds coordinates for vertices
    Returns:
      list of (int, int, int) - 3-tuples are CCW-oriented vertices of
          triangles making up the triangulation
    """

    if len(face) <= 3:
        return [tuple(face)]
    tris = EarChopTriFace(face, points)
    bord = _BorderEdges([face])
    triscdt = _CDT(tris, bord, points)
    return triscdt


def TriangulateFaceWithHoles(face, holes, points):
    """Like TriangulateFace, but with holes inside the face.

    Works by making one complex polygon that has segments to
    and from the holes ("islands"), and then using the same method
    as TriangulateFace.

    Args:
      face: list of int - indices in points, assumed CCW-oriented
      holes: list of list of int - each sublist is like face
          but CW-oriented and assumed to be inside face
      points: geom.Points - holds coordinates for vertices
    Returns:
      list of (int, int, int) - 3-tuples are CCW-oriented vertices of
          triangles making up the triangulation
    """

    if len(holes) == 0:
        return TriangulateFace(face, points)
    allfaces = [face] + holes
    sholes = [_SortFace(h, points) for h in holes]
    joinedface = _JoinIslands(face, sholes, points)
    tris = EarChopTriFace(joinedface, points)
    bord = _BorderEdges(allfaces)
    triscdt = _CDT(tris, bord, points)
    return triscdt


def QuadrangulateFace(face, points):
    """Quadrangulate the face (subdivide into convex quads and tris).

    Like TriangulateFace, but after triangulating, join as many pairs
    of triangles as possible into convex quadrilaterals.

    Args:
      face: list of int - indices in points, assumed CCW-oriented
      points: geom.Points - holds coordinates for vertices
    Returns:
      list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
          quadrilaterals and triangles making up the quadrangulation.
    """

    if len(face) <= 3:
        return [tuple(face)]
    tris = EarChopTriFace(face, points)
    bord = _BorderEdges([face])
    triscdt = _CDT(tris, bord, points)
    qs = _Quandrangulate(triscdt, bord, points)
    return qs


def QuadrangulateFaceWithHoles(face, holes, points):
    """Like QuadrangulateFace, but with holes inside the faces.

    Args:
      face: list of int - indices in points, assumed CCW-oriented
      holes: list of list of int - each sublist is like face
          but CW-oriented and assumed to be inside face
      points: geom.Points - holds coordinates for vertices
    Returns:
      list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
          quadrilaterals and triangles making up the quadrangulation.
    """

    if len(holes) == 0:
        return QuadrangulateFace(face, points)
    allfaces = [face] + holes
    sholes = [_SortFace(h, points) for h in holes]
    joinedface = _JoinIslands(face, sholes, points)
    tris = EarChopTriFace(joinedface, points)
    bord = _BorderEdges(allfaces)
    triscdt = _CDT(tris, bord, points)
    qs = _Quandrangulate(triscdt, bord, points)
    return qs


def _SortFace(face, points):
    """Rotate face so leftmost vertex is first, where face is
    list of indices in points."""

    n = len(face)
    if n <= 1:
        return face
    lefti = 0
    leftv = face[0]
    for i in range(1, n):
        # following comparison is lexicographic on n-tuple
        # so sorts on x first, using lower y as tie breaker.
        if points.pos[face[i]] < points.pos[leftv]:
            lefti = i
            leftv = face[i]
    return face[lefti:] + face[0:lefti]


def EarChopTriFace(face, points):
    """Triangulate given face, with coords given by indexing into points.
    Return list of faces, each of which will be a triangle.
    Use the ear-chopping method."""

    # start with lowest coord in 2d space to try
    # to get a pleasing uniform triangulation if starting with
    # a regular structure (like a grid)
    start = _GetLeastIndex(face, points)
    ans = []
    incr = 1
    n = len(face)
    while n > 3:
        i = _FindEar(face, n, start, incr, points)
        vm1 = face[(i - 1) % n]
        v0 = face[i]
        v1 = face[(i + 1) % n]
        face = _ChopEar(face, i)
        n = len(face)
        incr = - incr
        if incr == 1:
            start = i % n
        else:
            start = (i - 1) % n
        ans.append((vm1, v0, v1))
    ans.append(tuple(face))
    return ans


def _GetLeastIndex(face, points):
    """Return index of coordinate that is leftmost, lowest in face."""

    bestindex = 0
    bestpos = points.pos[face[0]]
    for i in range(1, len(face)):
        pos = points.pos[face[i]]
        if pos[0] < bestpos[0] or \
                (pos[0] == bestpos[0] and pos[1] < bestpos[1]):
            bestindex = i