model.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>

"""Manipulations of Models.
"""

__author__ = "howard.trickey@gmail.com"

from . import geom
from . import triquad
from . import offset
import math


def PolyAreasToModel(polyareas, bevel_amount, bevel_pitch, quadrangulate):
    """Convert a PolyAreas into a Model object.

    Assumes polyareas are in xy plane.

    Args:
      polyareas: geom.PolyAreas
      bevel_amount: float - if > 0, amount of bevel
      bevel_pitch: float - if > 0, angle in radians of bevel
      quadrangulate: bool - should n-gons be quadrangulated?
    Returns:
      geom.Model
    """

    m = geom.Model()
    if not polyareas:
        return m
    polyareas.points.AddZCoord(0.0)
    m.points = polyareas.points
    for pa in polyareas.polyareas:
        PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate)
    return m


def PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate):
    if bevel_amount > 0.0:
        BevelPolyAreaInModel(m, pa, bevel_amount, bevel_pitch, quadrangulate,
            False)
    elif quadrangulate:
        if len(pa.poly) == 0:
            return
        qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points)
        m.faces.extend(qpa)
        m.face_data.extend([pa.data] * len(qpa))
    else:
        m.faces.append(pa.poly)
        # TODO: just the first part of QuadrangulateFaceWithHoles, to join
        # holes to outer poly
        m.face_data.append(pa.data)


def ExtrudePolyAreasInModel(mdl, polyareas, depth, cap_back):
    """Extrude the boundaries given by polyareas by -depth in z.

    Assumes polyareas are in xy plane.

    Arguments:
      mdl: geom.Model - where to do extrusion
      polyareas: geom.Polyareas
      depth: float
      cap_back: bool - if True, cap off the back
    Side Effects:
      For all edges in polys in polyareas, make quads in Model
      extending those edges by depth in the negative z direction.
      The application data will be the data of the face that the edge
      is part of.
    """

    for pa in polyareas.polyareas:
        back_poly = _ExtrudePoly(mdl, pa.poly, depth, pa.data, True)
        back_holes = []
        for p in pa.holes:
            back_holes.append(_ExtrudePoly(mdl, p, depth, pa.data, False))
        if cap_back:
            qpa = triquad.QuadrangulateFaceWithHoles(back_poly, back_holes,
              polyareas.points)
            # need to reverse each poly to get normals pointing down
            for i, p in enumerate(qpa):
                t = list(p)
                t.reverse()
                qpa[i] = tuple(t)
            mdl.faces.extend(qpa)
            mdl.face_data.extend([pa.data] * len(qpa))


def _ExtrudePoly(mdl, poly, depth, data, isccw):
    """Extrude the poly by -depth in z

    Arguments:
      mdl: geom.Model - where to do extrusion
      poly: list of vertex indices
      depth: float
      data: application data
      isccw: True if counter-clockwise
    Side Effects
      For all edges in poly, make quads in Model
      extending those edges by depth in the negative z direction.
      The application data will be the data of the face that the edge
      is part of.
    Returns:
      list of int - vertices for extruded poly
    """

    if len(poly) < 2:
        return
    extruded_poly = []
    points = mdl.points
    if isccw:
        incr = 1
    else:
        incr = -1
    for i, v in enumerate(poly):
        vnext = poly[(i + incr) % len(poly)]
        (x0, y0, z0) = points.pos[v]
        (x1, y1, z1) = points.pos[vnext]
        vextrude = points.AddPoint((x0, y0, z0 - depth))
        vnextextrude = points.AddPoint((x1, y1, z1 - depth))
        if isccw:
            sideface = [v, vextrude, vnextextrude, vnext]
        else:
            sideface = [v, vnext, vnextextrude, vextrude]
        mdl.faces.append(sideface)
        mdl.face_data.append(data)
        extruded_poly.append(vextrude)
    return extruded_poly


def BevelPolyAreaInModel(mdl, polyarea,
    bevel_amount, bevel_pitch, quadrangulate, as_percent):
    """Bevel the interior of polyarea in model.

    This does smart beveling: advancing edges are merged
    rather than doing an 'overlap'.  Advancing edges that
    hit an opposite edge result in a split into two beveled areas.

    If the polyarea is not in the xy plane, do the work in a
    transformed model, and then transfer the changes back.

    Arguments:
      mdl: geom.Model - where to do bevel
      polyarea geom.PolyArea - area to bevel into
      bevel_amount: float - if > 0, amount of bevel
      bevel_pitch: float - if > 0, angle in radians of bevel
      quadrangulate: bool - should n-gons be quadrangulated?
      as_percent: bool - if True, interpret amount as percent of max
    Side Effects:
      Faces and points are added to model to model the
      bevel and the interior of the polyareas.
    """

    pa_norm = polyarea.Normal()
    if pa_norm == (0.0, 0.0, 1.0):
        m = mdl
        pa_rot = polyarea
    else:
        (pa_rot, inv_rot, inv_map) = _RotatedPolyAreaToXY(polyarea, pa_norm)
        # don't have to add the original faces into model, just their points.
        m = geom.Model()
        m.points = pa_rot.points
    vspeed = math.tan(bevel_pitch)
    off = offset.Offset(pa_rot, 0.0, vspeed)
    if as_percent:
        bevel_amount = bevel_amount * off.MaxAmount() / 100.0
    off.Build(bevel_amount)
    inner_pas = AddOffsetFacesToModel(m, off, polyarea.data)
    for pa in inner_pas.polyareas:
        if quadrangulate:
            if len(pa.poly) == 0:
                continue
            qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes,
                pa.points)
            m.faces.extend(qpa)
            m.face_data.extend([pa.data] * len(qpa))
        else:
            m.faces.append(pa.poly)
            m.face_data.append(pa.data)
    if m != mdl:
        _AddTransformedPolysToModel(mdl, m.faces, m.points, m.face_data,
             inv_rot, inv_map)


def AddOffsetFacesToModel(mdl, off, data=None):
    """Add the faces due to an offset into model.

    Returns the remaining interiors of the offset as a PolyAreas.

    Args:
      mdl: geom.Model - model to add offset faces into
      off: offset.Offset
      data: any - application data to be copied to the faces
    Returns:
      geom.PolyAreas
    """

    mdl.points = off.polyarea.points
    assert(len(mdl.points.pos) == 0 or len(mdl.points.pos[0]) == 3)
    o = off
    ostack = []
    while o:
        if o.endtime != 0.0:
            for face in o.facespokes:
                n = len(face)
                for i, spoke in enumerate(face):
                    nextspoke = face[(i + 1) % n]
                    v0 = spoke.origin
                    v1 = nextspoke.origin
                    v2 = nextspoke.dest
                    v3 = spoke.dest
                    if v2 == v3:
                        mface = [v0, v1, v2]
                    else:
                        mface = [v0, v1, v2, v3]
                    mdl.faces.append(mface)
                    mdl.face_data.append(data)
        ostack.extend(o.inneroffsets)
        if ostack:
            o = ostack.pop()
        else:
            o = None
    return off.InnerPolyAreas()


def BevelSelectionInModel(mdl, bevel_amount, bevel_pitch, quadrangulate,
        as_region, as_percent):
    """Bevel all the faces in the model, perhaps as one region.

    If as_region is False, each face is beveled individually,
    otherwise regions of contiguous faces are merged into
    PolyAreas and beveled as a whole.

    TODO: something if extracted PolyAreas are not approximately
    planar.

    Args:
      mdl: geom.Model
      bevel_amount: float - amount to inset
      bevel_pitch: float - angle of bevel side
      quadrangulate: bool - should insides be quadrangulated?
      as_region: bool - should faces be merged into regions?
      as_percent: bool - should amount be interpreted as a percent
          of the maximum amount (if True) or an absolute amount?
    Side effect:
      Beveling faces will be added to the model
    """

    pas = []
    if as_region:
        pas = RegionToPolyAreas(mdl.faces, mdl.points, mdl.face_data)
    else:
        for f, face in enumerate(mdl.faces):
            pas.append(geom.PolyArea(mdl.points, face, [],
                mdl.face_data[f]))
    for pa in pas:
        BevelPolyAreaInModel(mdl, pa,
            bevel_amount, bevel_pitch, quadrangulate, as_percent)


def RegionToPolyAreas(faces, points, data):
    """Find polygonal outlines induced by union of faces.

    Finds the polygons formed by boundary edges (those not
    sharing an edge with another face in region_faces), and
    turns those into PolyAreas.
    In the general case, there will be holes inside.
    We want to associate data with the region PolyAreas.
    Just choose a representative element of data[] when
    more than one face is combined into a PolyArea.

    Args:
      faces: list of list of int - each sublist is a face (indices into points)
      points: geom.Points - gives coordinates for vertices
      data: list of any - parallel to faces, app data to put in PolyAreas
    Returns:
      list of geom.PolyArea
    """

    ans = []
    (edges, vtoe) = _GetEdgeData(faces)
    (face_adj, is_interior_edge) = _GetFaceGraph(faces, edges, vtoe, points)
    (components, ftoc) = _FindFaceGraphComponents(faces, face_adj)
    for c in range(len(components)):
        boundary_edges = set()
        betodata = dict()
        vstobe = dict()
        for e, ((vs, ve), f) in enumerate(edges):
            if ftoc[f] != c or is_interior_edge[e]:
                continue
            boundary_edges.add(e)
            # vstobe[v] is set of edges leaving v
            # (could be more than one if boundary touches itself at a vertex)
            if vs in vstobe:
                vstobe[vs].append(e)
            else:
                vstobe[vs] = [e]
            betodata[(vs, ve)] = data[f]
        polys = []
        poly_data = []
        while boundary_edges:
            e = boundary_edges.pop()
            ((vstart, ve), face_i) = edges[e]
            poly = [vstart, ve]
            datum = betodata[(vstart, ve)]
            while ve != vstart:
                if ve not in vstobe:
                    print("whoops, couldn't close boundary")
                    break
                nextes = vstobe[ve]
                if len(nextes) == 1:
                    nexte = nextes[0]
                else:
                    # find a next edge with face index face_i
                    # TODO: this is not guaranteed to work,
                    # as continuation edge may have been for a different
                    # face that is now combined with face_i by erasing
                    # interior edges. Find a better algorithm here.
                    nexte = -1
                    for ne_cand in nextes:
                        if edges[ne_cand][1] == face_i:
                            nexte = ne_cand
                            break
                    if nexte == -1:
                        # case mentioned in TODO may have happened;
                        # just choose any nexte - may mess things up
                        nexte = nextes[0]
                ((_, ve), face_i) = edges[nexte]
                if nexte not in boundary_edges:
                    print("whoops, nexte not a boundary edge", nexte)
                    break
                boundary_edges.remove(nexte)
                if ve != vstart:
                    poly.append(ve)
            polys.append(poly)
            poly_data.append(datum)
        if len(polys) == 0:
            # can happen if an entire closed polytope is given
            # at least until we do an edge check
            return []
        elif len(polys) == 1:
            ans.append(geom.PolyArea(points, polys[0], [], poly_data[0]))
        else:
            outerf = _FindOuterPoly(polys, points, faces)
            pa = geom.PolyArea(points, polys[outerf], [], poly_data[outerf])
            pa.holes = [polys[i] for i in range(len(polys)) if i != outerf]
            ans.append(pa)
    return ans


def _GetEdgeData(faces):
    """Find edges from faces, and some lookup dictionaries.

    Args:
      faces: list of list of int - each a closed CCW polygon of vertex indices
    Returns:
      (list of ((int, int), int), dict{ int->list of int}) -
        list elements are ((startv, endv), face index)
        dict maps vertices to edge indices
    """

    edges = []
    vtoe = dict()
    for findex, f in enumerate(faces):
        nf = len(f)
        for i, v in enumerate(f):
            endv = f[(i + 1) % nf]
            edges.append(((v, endv), findex))
            eindex = len(edges) - 1
            if v in vtoe:
                vtoe[v].append(eindex)
            else:
                vtoe[v] = [eindex]
    return (edges, vtoe)


def _GetFaceGraph(faces, edges, vtoe, points):
    """Find the face adjacency graph.

    Faces are adjacent if they share an edge,
    and the shared edge goes in the reverse direction,
    and if the angle between them isn't too large.

    Args:
      faces: list of list of int
      edges: list of ((int, int), int) - see _GetEdgeData
      vtoe: dict{ int->list of int } - see _GetEdgeData
      points: geom.Points
    Returns:
      (list of  list of int, list of bool) -
        first list: each sublist is adjacent face indices for each face
        second list: maps edge index to True if it separates adjacent faces
    """

    face_adj = [[] for i in range(len(faces))]
    is_interior_edge = [False] * len(edges)
    for e, ((vs, ve), f) in enumerate(edges):
        for othere in vtoe[ve]:
            ((_, we), g) = edges[othere]
            if we == vs:
                # face g is adjacent to face f
                # TODO: angle check
                if g not in face_adj[f]:
                    face_adj[f].append(g)
                    is_interior_edge[e] = True
                # Don't bother with mirror relations, will catch later
    return (face_adj, is_interior_edge)


def _FindFaceGraphComponents(faces, face_adj):
    """Partition faces into connected components.

    Args:
      faces: list of list of int
      face_adj: list of list of int - see _GetFaceGraph
    Returns:
      (list of list of int, list of int) -
        first list partitions face indices into separate lists,
            each a component
        second list maps face indices into their component index
    """

    if not faces:
        return ([], [])
    components = []
    ftoc = [-1] * len(faces)
    for i in range(len(faces)):
        if ftoc[i] == -1:
            compi = len(components)
            comp = []
            _FFGCSearch(i, faces, face_adj, ftoc, compi, comp)
            components.append(comp)
    return (components, ftoc)


def _FFGCSearch(findex, faces, face_adj, ftoc, compi, comp):
    """Depth first search helper function for _FindFaceGraphComponents

    Searches recursively through all faces connected to findex, adding
    each face found to comp and setting ftoc for that face to compi.
    """

    comp.append(findex)
    ftoc[findex] = compi
    for otherf in face_adj[findex]:
        if ftoc[otherf] == -1:
            _FFGCSearch(otherf, faces, face_adj, ftoc, compi, comp)


def _FindOuterPoly(polys, points, faces):
    """Assuming polys has one CCW-oriented face when looking
    down average normal of faces, return that one.

    Only one of the faces should have a normal whose dot product
    with the average normal of faces is positive.

    Args:
      polys: list of list of int - list of polys given by vertex indices
      points: geom.Points
      faces: list of list of int - original selected region, used to find
          average normal
    Returns:
      int - the index in polys of the outermost one
    """

    if len(polys) < 2:
        return 0
    fnorm = (0.0, 0.0, 0.0)
    for face in faces:
        if len(face) > 2:
            fnorm = geom.VecAdd(fnorm, geom.Newell(face, points))
    if fnorm == (0.0, 0.0, 0.0):
        return 0
    # fnorm is really a multiple of the normal, but fine for test below
    for i, poly in enumerate(polys):
        if len(poly) > 2:
            pnorm = geom.Newell(poly, points)
            if geom.VecDot(fnorm, pnorm) > 0:
                return i
    print("whoops, couldn't find an outermost poly")
    return 0


def _RotatedPolyAreaToXY(polyarea, norm):
    """Return a  PolyArea rotated to xy plane.

    Only the points in polyarea will be transferred.

    Args:
      polyarea: geom.PolyArea
      norm: the normal for polyarea
    Returns:
      (geom.PolyArea, (float, ..., float), dict{ int -> int }) - new PolyArea,
          4x3 inverse transform, dict mapping new verts to old ones
    """

    # find rotation matrix that takes norm to (0,0,1)
    (nx, ny, nz) = norm
    if abs(nx) < abs(ny) and abs(nx) < abs(nz):
        v = (vx, vy, vz) = geom.Norm3(0.0, nz, - ny)
    elif abs(ny) < abs(nz):
        v = (vx, vy, vz) = geom.Norm3(nz, 0.0, - nx)
    else:
        v = (vx, vy, vz) = geom.Norm3(ny, - nx, 0.0)
    (ux, uy, uz) = geom.Cross3(v, norm)
    rotmat = [ux, vx, nx, uy, vy, ny, uz, vz, nz, 0.0, 0.0, 0.0]
    # rotation matrices are orthogonal, so inverse is transpose
    invrotmat = [ux, uy, uz, vx, vy, vz, nx, ny, nz, 0.0, 0.0, 0.0]
    pointmap = dict()
    invpointmap = dict()
    newpoints = geom.Points()
    for poly in [polyarea.poly] + polyarea.holes:
        for v in poly:
            vcoords = polyarea.points.pos[v]
            newvcoords = geom.MulPoint3(vcoords, rotmat)
            newv = newpoints.AddPoint(newvcoords)
            pointmap[v] = newv
            invpointmap[newv] = v
    pa = geom.PolyArea(newpoints)
    pa.poly = [pointmap[v] for v in polyarea.poly]
    pa.holes = [[pointmap[v] for v in hole] for hole in polyarea.holes]
    pa.data = polyarea.data
    return (pa, invrotmat, invpointmap)


def _AddTransformedPolysToModel(mdl, polys, points, poly_data,
        transform, pointmap):
    """Add (transformed) the points and faces to a model.

    Add polys to mdl.  The polys have coordinates given by indices
    into points.pos; those need to be transformed by multiplying by
    the transform matrix.
    The vertices may already exist in mdl.  Rather than relying on
    AddPoint to detect the duplicate (transform rounding error makes
    that dicey), the pointmap dictionar is used to map vertex indices
    in polys into those in mdl - if they exist already.

    Args:
      mdl: geom.Model - where to put new vertices, faces
      polys: list of list of int - each sublist a poly
      points: geom.Points - coords for vertices in polys
      poly_data: list of any - parallel to polys
      transform: (float, ..., float) - 12-tuple, a 4x3 transform matrix
      pointmap: dict { int -> int } - maps new vertex indices to old ones
    Side Effects:
      The model gets new faces and vertices, based on those in polys.
      We are allowed to modify pointmap, as it will be discarded after call.
    """

    for i, coords in enumerate(points.pos):
        if i not in pointmap:
            p = geom.MulPoint3(coords, transform)
            pointmap[i] = mdl.points.AddPoint(p)
    for i, poly in enumerate(polys):
        mpoly = [pointmap[v] for v in poly]
        mdl.faces.append(mpoly)
        mdl.face_data.append(poly_data[i])