# -*- coding: utf-8 -*-

#  Copyright © 2013-2017  B. Clausius <barcc@gmx.de>
#
#  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 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 this program.  If not, see <http://www.gnu.org/licenses/>.

from collections import deque
from math import atan2, sin, cos, pi, sqrt
import weakref
import itertools

import pybiklib.utils
from pybiklib.utils import epsdigits, epsilon


sid_cache = weakref.WeakKeyDictionary()

def roundeps(val):
    return round(val, epsdigits+1)
    
    
class Vector (pybiklib.utils.Vector):
    def rounded(self):
        return self.__class__(roundeps(s) for s in self)
        
    def scaled(self, vector):
        return self.__class__(s*v for s,v in zip(self, vector))
        
    def equalfuzzy(self, other):
        return all((-epsilon<s<epsilon) for s in self-other)
        
    def inversfuzzy(self, other): return (self+other).length_squared() < epsilon
        
    @classmethod
    def polygon(cls, point, n):
        x0, y0 = point
        yield x0, y0
        r = cls(point).length()
        angle0 = atan2(y0, x0)
        angled = 2 * pi / n
        
        for k in range(1, n):
            angle = angle0 + k * angled
            yield r * cos(angle), r * sin(angle)
        
        
def distance_axis_vert(axis, vert):
    ''' axis:  vector that defines a line through (0,0,0)
        return: the shortest distance between axis and vert
    '''
    axis = Vector(axis)
    base_point = axis * vert.point.dot(axis) / axis.dot(axis)
    return (vert.point - base_point).length()
    
def distance_axis_edge(axis, verts):
    ''' axis:  vector that defines a line through (0,0,0)
        verts: two points that define a line segment
        return: the shortest distance between axis and edge
    '''
    vec_axis = axis             # S1.P1 - S1.P0, S1.P0 == (0,0,0)
    vec_edge = verts[1].point - verts[0].point  # S2.P1 - S2.P0
    vec_vert0 = -verts[0].point  # S1.P0 - S2.P0
    a = vec_axis.dot(vec_axis)  # always >= 0
    b = vec_axis.dot(vec_edge)
    c = vec_edge.dot(vec_edge)  # always >= 0
    d = vec_axis.dot(vec_vert0)
    e = vec_edge.dot(vec_vert0)
    D = a*c - b*b   # always >= 0
    sD = D      # sc = sN / sD, default sD = D >= 0
    tD = D      # tc = tN / tD, default tD = D >= 0
    
    # compute the line parameters of the two closest points
    if D < epsilon: # the lines are almost parallel
        sN = 0.0    # force using point P0 on segment S1
        sD = 1.0    # to prevent possible division by 0.0 later
        tN = e
        tD = c
    else:   # get the closest points on the infinite lines
        sN = b*e - c*d
        tN = a*e - b*d
        
    if tN < 0.0:  # tc < 0 => the t=0 edge is visible
        tN = 0.0
        # recompute sc for this edge
        sN = -d
        sD = a
    elif tN > tD:   # tc > 1  => the t=1 edge is visible
        tN = tD
        # recompute sc for this edge
        sN = (-d + b)
        sD = a
    # finally do the division to get sc and tc
    sc = 0.0 if abs(sN) < epsilon else sN / sD
    tc = 0.0 if abs(tN) < epsilon else tN / tD
    
    # get the difference of the two closest points
    dP = vec_vert0 + (vec_axis * sc) - (vec_edge * tc)    # =  S1(sc) - S2(tc)
    
    return dP.length()     # return the closest distance
    
    
class Plane:
    def __init__(self, normal):
        self.normal = Vector(normal)
        self.distance = 0
        
    def __repr__(self):
        return '{0}({1.normal}, {1.distance})'.format(self.__class__.__name__, self)
        
    def signed_distance(self, point):
        return self.normal.dot(point) - self.distance
        
    def intersect_segment(self, p, q):
        pvalue = self.signed_distance(p.point)
        qvalue = self.signed_distance(q.point)
        if abs(pvalue - qvalue) < epsilon:
            return -1, None
        if abs(pvalue) < epsilon:
            return 0, p
        if abs(qvalue) < epsilon:
            return 0, q
        if pvalue < 0 < qvalue or qvalue < 0 < pvalue:
            d = pvalue / (pvalue - qvalue)
            return 1, Vert(p.point + (q.point - p.point) * d)
        else:
            return -1, None
            
        
class GeomBase:
    def __euler(self):
        chars = 'vefc456'
        halfsets, fullsets = self.get_all_ldims()
        part = lambda c,h,f: ('{}={}'.format(c,f) if h==f else '{}={}/{}'.format(c,f,h))
        return ' '.join(part(c,len(hs),len(fs)) for c,hs,fs in zip(chars,halfsets,fullsets))
        
    def __str__(self):
        return '{}({})'.format(self.__class__.__name__, self.__euler())
        
    def __repr__(self):
        sid = sid_cache.setdefault(self, len(sid_cache))
        return '{}({}, {})'.format(self.__class__.__name__, sid, self.__euler())
        
    def clone(self, clone_data):
        if id(self) in clone_data:
            return
        obj = clone_data[id(self)] = self.__class__()
        vals = [(attr, getattr(self, attr)) for attr in reversed(self.fields_geom)]
        def func(clone_data):
            for attr in self.fields_kwarg:
                if attr not in self.fields_geom:
                    setattr(obj, attr, getattr(self, attr))
            for attr, val in vals:
                if isinstance(val, list):
                    newval = [clone_data[id(o)] for o in val]
                elif val is not None:
                    newval = clone_data[id(val)]
                else:
                    newval = None
                setattr(obj, attr, newval)
        for attr, val in vals:
            if isinstance(val, list):
                for o in val:
                    yield o.clone
            elif val is not None:
                yield val.clone
        yield func
        
    def center(self):
        c = Vector()
        for i, v in enumerate(self.verts):
            c += v.point
        return c / (i+1)
        
        
class HalfGeom (GeomBase):
    fields_geom = '_ldims', '_hdim', '_fullgeom'
    
    def __init__(self, ldims, fullgeom):
        self._ldims = []
        for l in ldims:
            l.hdim = self
        #assert all(l.dim+1==self.dim for l in self._ldims)
        self._hdim = None
        self._fullgeom = None
        self.fullgeom = fullgeom
        
    def assert_valid(self):
        assert all(l.dim+1==self.dim for l in self._ldims)
        assert all(l._hdim is self for l in self._ldims)
        assert self._hdim is None or self in self._hdim._ldims
        if self.fullgeom is not None:
            assert self.fullgeom.assert_valid()
        return True
        
    @property
    def fullgeom(self):
        return self._fullgeom
    @fullgeom.setter
    def fullgeom(self, fg):
        if self._fullgeom is not None:
            self._fullgeom._halfgeoms.remove(self)
        self._fullgeom = fg
        if fg is not None:
            assert self not in fg._halfgeoms
            fg._halfgeoms.append(self)
            
    @property
    def ldims(self):
        return self._ldims
        
    @property
    def hdim(self):
        return self._hdim
    @hdim.setter
    def hdim(self, hf):
        if self._hdim is not None:
            self._hdim._ldims.remove(self)
        self._hdim = hf
        if hf is not None:
            assert self not in hf._ldims
            hf._ldims.append(self)
            
    def insert_before(self, hg):
        #TODO: try to remove this function, the order shouldnt matter,
        #      but currently needed to preserve content of model data files
        if self._hdim is not None:
            self._hdim._ldims.remove(self)
        self._hdim = hg._hdim
        assert self not in hg._hdim._ldims
        i = hg._hdim._ldims.index(hg) or len(hg._hdim._ldims)
        hg._hdim._ldims.insert(i, self)
        
    def insert_after(self, hg):
        #TODO: try to remove this function, the order shouldnt matter,
        #      but currently needed to preserve content of model data files
        if self._hdim is not None:
            self._hdim._ldims.remove(self)
        self._hdim = hg._hdim
        assert self not in hg._hdim._ldims
        i = hg._hdim._ldims.index(hg)
        hg._hdim._ldims.insert(i+1, self)
        
    def insert_after_first(self, hg):
        #TODO: try to remove this function, the order shouldnt matter,
        #      but currently needed to preserve content of model data files
        if hg._ldims:
            self.insert_after(hg._ldims[0])
        else:
            self._hdim = hg
            
    def get_all_ldims(self, halfsets=None, fullsets=None):
        if halfsets is None:
            halfsets = [set() for i in range(self.dim+1)]
        if fullsets is None:
            fullsets = [set() for i in range(self.dim+1)]
        for ldim in self._ldims:
            ldim.get_all_ldims(halfsets, fullsets)
        halfsets[self.dim].add(self)
        fullsets[self.dim].add(self._fullgeom)
        return halfsets, fullsets
        
    def rotate_ldims_to(self, hg):
        i = self._ldims.index(hg)
        if i:
            self._ldims = self._ldims[i:] + self._ldims[:i]
        return self
            
    @property
    def other(self):
        '''other halfgeom of fullgeom in the same hdim'''
        for hg in self.fullgeom.halfgeoms:
            if hg is not self and hg.hdim.hdim is self.hdim.hdim:
                return hg
        return None  # hole in surface of self.hdim.hdim
        
        
class FullGeom (GeomBase):
    fields_geom = '_halfgeoms',
    
    def __init__(self):
        self._halfgeoms = []
        
    def assert_valid(self):
        assert all(hg.fullgeom is self for hg in self.halfgeoms)
        return True
        
    @property
    def halfgeoms(self):
        return self._halfgeoms
        
    def get_all_ldims(self):
        halfsets = fullsets = None
        for hg in self._halfgeoms:
            halfsets, fullsets = hg.get_all_ldims(halfsets, fullsets)
        return halfsets, fullsets
        
        
class HalfVert (HalfGeom):
    dim = 0
    fields_arg = ()
    fields_kwarg = ()
    
    def __init__(self, *, vert=None):
        super().__init__([], vert)
        
    vert = HalfGeom.fullgeom
    halfedge = HalfGeom.hdim
    
    def assert_valid(self):
        super().assert_valid()
        return True
        
        
class Vert (FullGeom):
    fields_arg = ()
    fields_kwarg = 'point', 'id',
    
    def __init__(self, point=None, *, id=None):
        super().__init__()
        self.point = None if point is None else Vector(point)
        self.id = id
        
    def assert_valid(self):
        super().assert_valid()
        import numbers
        assert all(isinstance(i, numbers.Number) for i in self.point)
        return True
        
    halfverts = FullGeom.halfgeoms
    
        
class HalfEdge (HalfGeom):
    dim = 1
    fields_arg = ()
    fields_kwarg = ()
    
    def __init__(self, verts=None, edge=None):
        #XXX: different from other Half* classes where Half*-objects are used for initialisation
        verts = list(verts or [])
        assert all(isinstance(v, FullGeom) for v in verts)
        halfverts = [HalfVert(vert=v) for v in verts]
        super().__init__(halfverts, edge)
        
    halfverts = HalfGeom.ldims
    edge = HalfGeom.fullgeom
    halfface = HalfGeom.hdim
            
    @property
    def verts(self):
        return tuple(hv.vert for hv in self.halfverts)
        
    def create_other(self):
        return self.__class__(reversed(self.verts), self.edge)
        
    def assert_valid(self):
        assert super().assert_valid()
        assert all(isinstance(hv, HalfVert) for hv in self.halfverts)
        assert isinstance(self.edge, Edge)
        assert isinstance(self.halfface, HalfFace)
        assert all(hv.assert_valid() for hv in self.halfverts)
        assert len(self.halfverts) == 2
        assert self.edge.assert_valid()
        assert self.halfverts[1].vert is self.nextatface.halfverts[0].vert
        assert self.other is None or isinstance(self.other, self.__class__), self.other.__class__
        assert self.other is None or self.other.other is self
        return True
        
    @property
    def nextatface(self):
        ldims = self.hdim.ldims
        return ldims[(ldims.index(self)+1) % len(ldims)]
        
    @classmethod
    def polygon(cls, points, ids=None):
        vert0 = Vert(next(points))
        vertp = vert0
        ids = itertools.repeat(None) if ids is None else iter(ids)
        for point, eid in zip(points, ids):
            vertc = Vert(point)
            yield cls((vertp, vertc), Edge(id=eid))
            vertp = vertc
        yield cls((vertp, vert0), Edge(id=next(ids)))
        
    def close_hole(self, nedges, id=None):
        def next_at_hole(v):
            for hv in v.halfverts:
                he = hv.halfedge
                if he.halfface is not hf and he.other is None:
                    # only works if at the vert is only one hole
                    return he
            return None
        assert self.other is None
        hf = HalfFace(face=Face(type='face', id=id))
        hf.halfcell = self.halfface.halfcell
        he = self
        v0 = he.halfverts[1].vert
        v = he.halfverts[0].vert
        is_attach_prev = len(v0.halfverts) > 2
        is_attach_next = len(v.halfverts) > 2
        he.create_other().halfface = hf
        # attach hf to existing edges
        while v is not v0 and is_attach_next:
            he = next_at_hole(v)
            v = he.halfverts[0].vert
            is_attach_next = len(v.halfverts) > 2
            he.create_other().halfface = hf
        # attach hf to existing edges in the other direction
        halfedges = []
        while v is not v0 and is_attach_prev:
            he = next_at_hole(v0)
            halfedges.append(he)
            v0 = he.halfverts[1].vert
            is_attach_prev = len(v0.halfverts) > 2
        # if the hole needs more faces, create some edges
        for i in range(len(hf.halfedges)+len(halfedges), nedges):
            vn = Vert() if i < nedges-1 else v0
            HalfEdge([v, vn], edge=Edge()).halfface = hf
            v = vn
        for he in halfedges:
            he.create_other().halfface = hf
        return hf
        
        
class Edge (FullGeom):
    fields_arg = ()
    fields_kwarg = 'id',
    
    def __init__(self, *, id=None):
        super().__init__()
        self.id = id
        
    def assert_valid(self):
        super().assert_valid()
        assert all(isinstance(he, HalfEdge) for he in self.halfedges)
        assert len(self.halfedges) > 0
        assert all(isinstance(he, HalfEdge) for he in self.halfedges)
        assert all(he.edge is self for he in self.halfedges)
        return True
        
    halfedges = FullGeom.halfgeoms
    
    def split(self, vert):
        verts = self.halfedges[0].verts
        edge2 = self.__class__()
        for he1 in self.halfedges:
            if he1.verts[0] == verts[0]:
                he2 = HalfEdge([vert, he1.verts[1]], edge2)
                he1.halfverts[1].vert = vert
                he2.insert_after(he1)
            else:
                assert he1.verts[1] == verts[0], (he1.verts, verts)
                he2 = HalfEdge([he1.verts[0], vert], edge2)
                he1.halfverts[0].vert = vert
                he2.insert_before(he1)
                
        
class HalfFace (HalfGeom):
    dim = 2
    fields_arg = 'halfedges',
    fields_kwarg = ()
    
    def __init__(self, *, halfedges=None, face=None):
        super().__init__(halfedges or [], face)
        
    def assert_valid(self):
        assert super().assert_valid()
        assert set(self.halfedges) == set(self.ldims), (set(self.halfedges), self.ldims)
        assert isinstance(self.face, Face)
        assert self.face.assert_valid()
        assert all(he.halfface is self for he in self.ldims)
        assert all(isinstance(he, HalfEdge) for he in self.halfedges)
        assert all(he.assert_valid() for he in self.halfedges)
        return True
        
    halfedges = HalfGeom.ldims
    face = HalfGeom.fullgeom
    halfcell = HalfGeom.hdim
        
    @property
    def edges(self):
        for he in self.halfedges:
            yield he.edge
            
    @property
    def verts(self):
        for he in self.halfedges:
            yield he.verts[0]
            
    def normal(self):
        verts = self.verts
        try:
            v1 = next(verts).point
            v2 = next(verts).point
            v3 = next(verts).point
        except StopIteration:
            raise ValueError('Too few vertices')
        return (v2-v1).cross(v3-v1).normalised()
        
    def sort(self):
        try:
            he = self.halfedges[0]
        except IndexError:
            pass
        else:
            hv_last = he.halfverts[0].other
            while he.halfverts[1] is not hv_last:
                he = he.halfverts[1].other.halfedge
                he.hdim = self
        return self
        
    @classmethod
    def polygon(cls, point, n, ids=None):
        points = Vector.polygon(point, n)
        halfedges = HalfEdge.polygon(points, ids)
        halfface = cls(face=Face())
        
        for he in halfedges:
            he.halfface = halfface
        return halfface
        
    def extrude_y(self, upy, downy, ids=None):
        cell = Cell()
        self.halfcell = cell
        # down face
        face_down = Face()
        if ids is not None:
            self.face.type = face_down.type = 'face'
            self.face.id, face_down.id = ids
        halfface_down = HalfFace(face=face_down)
        
        def extrude_vert(vert_up, upy, downy):
            x1, y1 = vert_up.point
            vert_up.point = Vector((x1, upy, -y1))
            vert_down = Vert((x1, downy, -y1))
            return vert_up, vert_down
            
        vp_down = None
        for he_up in self.halfedges:
            # vertices
            vc_up, vc_down = extrude_vert(he_up.verts[1], upy, downy)
            # edges
            e_down = Edge()
            ec_side = Edge()
            if vp_down is None:
                v1_up, v1_down = vc_up, vc_down
                e0_down = e_down
                e1_side = ec_side
            else:
                # down halfedge
                HalfEdge([vc_down, vp_down], e_down).halfface = halfface_down
                # side face
                halfface_side = HalfFace(face=Face(type='face', id=he_up.edge.id))
                HalfEdge([vp_down, vc_down], e_down).halfface = halfface_side
                HalfEdge([vc_down, vc_up], ec_side).halfface = halfface_side
                HalfEdge([vc_up, vp_up], he_up.edge).halfface = halfface_side
                HalfEdge([vp_up, vp_down], ep_side).halfface = halfface_side
                halfface_side.halfcell = cell
            # next
            vp_up, vp_down = vc_up, vc_down
            ep_side = ec_side
        he_up = self.ldims[0]
        # joining down halfedge
        he_down = HalfEdge([v1_down, vp_down], e0_down)
        he_down.halfface = halfface_down
        halfface_down.rotate_ldims_to(he_down).sort()
        # joining side face
        halfface_side = HalfFace(face=Face(type='face', id=he_up.edge.id))
        HalfEdge([vp_down, v1_down], e0_down).halfface = halfface_side
        HalfEdge([v1_down, v1_up], e1_side).halfface = halfface_side
        HalfEdge([v1_up, vp_up], he_up.edge).halfface = halfface_side
        HalfEdge([vp_up, vp_down], ep_side).halfface = halfface_side
        halfface_side.halfcell = cell
        
        halfface_down.halfcell = cell
        return cell
        
    def pyramid(self, upy, downy, id=None):
        cell = Cell()
        self.halfcell = cell
        vert_up = Vert((0., upy, 0.))
        
        def change_vert(vert_down, downy):
            x1, y1 = vert_down.point
            vert_down.point = Vector((-x1, downy, -y1))
            return vert_down
            
        vp_down = None
        for he in self.halfedges:
            # vertices
            vc_down = change_vert(he.verts[1], downy)
            # edges
            ec_side = Edge()
            if vp_down is None:
                v1_down = vc_down
                e1_side = ec_side
            else:
                # side face
                halfface_side = HalfFace(face=Face(type='face', id=he.edge.id))
                HalfEdge([vc_down, vp_down], he.edge).halfface = halfface_side
                HalfEdge([vp_down, vert_up], ep_side).halfface = halfface_side
                HalfEdge([vert_up, vc_down], ec_side).halfface = halfface_side
                halfface_side.halfcell = cell
            # next
            vp_down = vc_down
            ep_side = ec_side
        he = self.ldims[0]
        # joining side face
        halfface_side = HalfFace(face=Face(type='face', id=he.edge.id))
        HalfEdge([v1_down, vp_down], he.edge).halfface = halfface_side
        HalfEdge([vp_down, vert_up], ep_side).halfface = halfface_side
        HalfEdge([vert_up, v1_down], e1_side).halfface = halfface_side
        halfface_side.halfcell = cell
        # down face
        if id is not None:
            self.face.type = 'face'
            self.face.id = id
        return cell
        
        
class Face (FullGeom):
    fields_arg = ()
    fields_kwarg = 'type', 'id'
    
    def __init__(self, *, type=None, id=None):
        super().__init__()
        self.type = type
        self.id = id
        
    def assert_valid(self):
        super().assert_valid()
        assert all(isinstance(f, HalfFace) for f in self.halffaces)
        assert all(hf.face is self for hf in self.halffaces)
        assert 1 <= len(self.halffaces) <= 2
        lenhf0 = len(list(self.halffaces[0].halfedges))
        assert all(lenhf0 == len(list(hf.halfedges)) for hf in self.halffaces)
        return True
        
    halffaces = FullGeom.halfgeoms
    
    @property
    def edges(self):
        return self.halffaces[0].edges
        
    @property
    def verts(self):
        return self.halffaces[0].verts
        
    def split(self, split_verts):
        assert len(split_verts) == 2, len(split_verts)
        edge = Edge()
        newface = Face(type=self.type, id=self.id)
        edgef1 = None
        for halfface in self.halffaces:
            he = halfface.ldims[0]
            # make sure to start on the same edge for all halffaces
            if edgef1 is None:
                edgef1 = he.edge
            else:
                he = next(_he for _he in edgef1.halfedges if _he.halfface is halfface)
                halfface.rotate_ldims_to(he)
            hedge2f1 = None
            hedge1f2 = None
            for he in halfface.halfedges[:]:
                if hedge2f1 is None:
                    if he.verts[1] in split_verts:
                        # end of the old halfface
                        hedge2f1 = he
                elif hedge1f2 is None:
                    # start of the new halfface
                    newhalfface = HalfFace(face=newface)
                    newhalfface.insert_after(halfface)
                    he.halfface = newhalfface
                    hedge1f2 = he
                elif he.verts[1] not in split_verts:
                    # proceed to the end of the new halfface
                    he.halfface = newhalfface
                else:
                    # end of the new halfface
                    he.halfface = newhalfface
                    newhalfface.rotate_ldims_to(hedge1f2)
                    hedge2f2 = he
                    break
            # insert a new halfedge between the splitted halfface
            hedge1 = HalfEdge(split_verts, edge)
            hedge2 = HalfEdge([split_verts[1], split_verts[0]], edge)
            if hedge2f1.verts[1] is split_verts[1]:
                hedge1, hedge2 = hedge2, hedge1
            hedge1.insert_after(hedge2f1)
            hedge2.insert_after(hedge2f2)
        return edge
        
        
class HalfCell (HalfGeom):
    dim = 3
    fields_arg = 'halffaces',
    fields_kwarg = 'indices',
    
    def __init__(self, *, halffaces=None):
        super().__init__(halffaces or [], None)
        self.indices = None
        
    def assert_valid(self):
        assert super().assert_valid()
        assert list(self.halffaces) == list(self.ldims), (list(self.halffaces), list(self.ldims))
        assert all(isinstance(f, HalfFace) for f in self.halffaces)
        assert all(hf.assert_valid() for hf in self.halffaces)
        assert len(list(self.edges)) > 3
        assert len(list(self.verts)) > 3
        assert all(he.hdim.hdim is not None for e in self.edges for he in e.halfedges)
        return True
        
    halffaces = HalfGeom.ldims
                
    @property
    def faces(self):
        for hf in self.halffaces:
            yield hf.face
            
    @property
    def edges(self):
        cache = []
        for f in self.faces:
            for e in f.edges:
                if e not in cache:
                    cache.append(e)
                    yield e
                    
    @property
    def verts(self):
        cache = []
        for hf in self.halffaces:
            for v in hf.verts:
                if v not in cache:
                    cache.append(v)
                    yield v
                    
    def split(self, split_edges, id=None):
        halffaces, halffaces1, halffaces2 = [self.ldims[0]], [self.ldims[0]], []
        # new face
        face = Face(type='cut', id=id)
        halfface1 = HalfFace(face=face)
        halfface2 = HalfFace(face=face)
        
        while halffaces:
            for he in halffaces.pop().halfedges:
                hf = he.other.halfface
                if he.edge in split_edges:
                    he.create_other().halfface = halfface1
                    if hf not in halffaces2:
                        halffaces2.append(hf)
                else:
                    if hf not in halffaces1:
                        halffaces.append(hf)
                        halffaces1.append(hf)
        halffaces = halffaces2[:]
        while halffaces:
            for he in halffaces.pop().halfedges:
                if he.edge in split_edges:
                    he.create_other().halfface = halfface2
                else:
                    hf = he.other.halfface
                    if hf not in halffaces2:
                        halffaces.append(hf)
                        halffaces2.append(hf)
        # update cell
        halfface1.sort()
        phf = halfface1
        for hf in halffaces1:
            phf = hf
            #XXX: just to reorder ldims to the order of halfface1, try to remove this line:
            hf.hdim = self
        halfface1.hdim = self
        
        # new cell
        halfface2.sort()
        halffaces2.append(halfface2)
        return Cell(halffaces=halffaces2)
            
    @classmethod
    def dodecahedron(cls, ids=None):
        #rfi = sqrt(25 + 10*sqrt(5)) / 10
        rfu = sqrt(50 + 10*sqrt(5)) / 10
        ri = sqrt((25 + 11*sqrt(5))/10) / 2
        #ru = sqrt(3) * (1 + sqrt(5)) / 4
        
        #1: h01² + (r1u-rfu)² = 1²
        #2: (ri-h01)² + r1u² = ru²
        # h01 = rfu
        # r1u² = sqrt(45 + 20*sqrt(5)) / 5
        r1u = sqrt(sqrt(45 + 20*sqrt(5)) / 5)
        
        ids = itertools.repeat(None) if ids is None else iter(ids)
        cell = cls()
        
        hf_down = HalfFace.polygon((0.,-rfu), 5)
        hf_down.halfcell = cell
        hf_down.face.type = 'face'
        hf_down.face.id = next(ids)
        for v in hf_down.verts:
            x,z = v.point
            v.point = Vector((x, -ri, z))
        hfs1 = [he.close_hole(5, id=fid) for he, fid in zip(hf_down.halfedges, ids)]
        hfs2 = [hf.halfedges[2].close_hole(5, id=fid) for hf, fid in zip(hfs1, ids)]
        hf_up = hfs2[0].halfedges[3].close_hole(5, id=next(ids))
        
        hes01 = [he for e in cell.edges for he in e.halfedges if he.verts[0].point is not None and he.verts[1].point is None]
        for he in hes01:
            v0 = he.verts[0]
            v1 = he.verts[1]
            x,y,z = v0.point
            v1.point = Vector((x/rfu*r1u, rfu-ri, z/rfu*r1u))
        hes12 = [he for e in cell.edges for he in e.halfedges if he.verts[0].point is not None and he.verts[1].point is None]
        for he in hes12:
            he.verts[1].point = True  #XXX: marker, any value not None would do it
        hes23 = [he for e in cell.edges for he in e.halfedges if he.verts[0].point is not None and he.verts[1].point is None]
        for i, he in enumerate(hes23):
            he.verts[0].point = -hes01[(i+2)%5].verts[1].point
            he.verts[1].point = -hes01[(i+2)%5].verts[0].point
            
        return cell
        
Cell = HalfCell


class Polyhedron:
    def __init__(self):
        self._cells = []
        
    def clone(self):
        clone_data = {}
        obj = self.__class__()
        funcs = deque(f for o in self._cells for f in o.clone(clone_data))
        obj._cells = [clone_data[id(o)] for o in self._cells]
        while funcs:
            f = funcs.popleft()
            r = f(clone_data)
            if r is not None:
                funcs.extendleft(reversed(list(r)))
        return obj
        
    def assert_valid(self):
        import itertools
        assert all(isinstance(c, Cell) and c.assert_valid() for c in self.cells)
        assert all(isinstance(f, Face) and f.assert_valid() for f in self.faces)
        assert all(isinstance(e, Edge) and e.assert_valid() for e in self.edges)
        assert all(isinstance(v, Vert) and v.assert_valid() for v in self.verts)
        assert all(not v1.point.equalfuzzy(v2.point) for v1, v2 in itertools.combinations(self.verts, 2)), self.verts
        return True
        
    @property
    def cells(self):
        return self._cells
    @cells.setter
    def cells(self, cells):
        self._cells = cells
        
    @property
    def faces(self):
        cache = []
        for c in self.cells:
            for f in c.faces:
                if f not in cache:
                    cache.append(f)
        return cache
        
    @property
    def edges(self):
        cache = []
        for f in self.faces:
            for e in f.edges:
                if e not in cache:
                    cache.append(e)
        return cache
        
    @property
    def verts(self):
        cache = []
        for f in self.faces:
            for v in f.verts:
                if v not in cache:
                    cache.append(v)
        return cache
        
    def set_verts(self, *args):
        self._polytopes = [Vert(arg) for arg in args]
            
    def set_edges(self, *args):
        halfedges = []
        for arg in args:
            verts = [self._polytopes[i] for i in arg]
            halfedges.append(HalfEdge(verts=verts, edge=Edge()))
        self._polytopes = halfedges
        
    def set_faces(self, *args, ids=None):
        halffaces = []
        if ids is None:
            faces = [Face() for unused in args]
        else:
            assert len(args) == len(ids)
            faces = [Face(type='face', id=fid) for fid in ids]
        for arg, face in zip(args, faces):
            halfface = HalfFace(face=face)
            for i in arg:
                # this only works if args describe one closed cell
                if i > 0:
                    self._polytopes[i-1].halfface = halfface
                elif i < 0:
                    self._polytopes[-i-1].create_other().halfface = halfface
                else:
                    raise ValueError()
            halffaces.append(halfface)
        del self._polytopes
        self._cells = [Cell(halffaces=halffaces)]
            
    def split_plane(self, plane, indexpos=None, split_id=None):
        # split intersected edges
        new_verts = []
        for edge in self.edges[:]:
            pos, vert = plane.intersect_segment(*edge.halfedges[0].verts)
            if pos > 0:
                edge.split(vert)
                new_verts.append(vert)
            elif pos == 0 and vert not in new_verts:
                new_verts.append(vert)
        # split intersected faces
        new_edges = []
        for face in self.faces[:]:
            split_verts = [v for v in face.verts for nv in new_verts if nv is v]
            assert 0 <= len(split_verts) <= 2, split_verts
            if len(split_verts) == 2:
                for he in face.halffaces[0].halfedges:
                    if set(he.verts) == set(split_verts):
                        new_edges.append(he.edge)
                        break
                else:
                    new_edge = face.split(split_verts)
                    new_edges.append(new_edge)
        # split cells
        for cell in self.cells[:]:
            split_edges = [e for e in cell.edges if e in new_edges]
            if len(split_edges) > 2:
                new_cell = cell.split(split_edges, split_id)
                self.cells.append(new_cell)
            else:
                new_cell = None
            if indexpos is not None:
                for hf in cell.halffaces:
                    for he in hf.halfedges:
                        if he.edge not in split_edges:
                            sd1 = -plane.signed_distance(he.verts[0].point)
                            # the verts of the edge e cannot be both in range epsilon to the plane,
                            # otherwise e would be in split_edges, so only test one vert for epsilon
                            if sd1 > epsilon or sd1 >= -epsilon and plane.signed_distance(he.verts[1].point) < 0:
                                if new_cell is not None:
                                    new_cell.indices = cell.indices
                                new_cell = cell
                            elif new_cell is None:
                                break
                            indices = list(cell.indices)
                            indices[indexpos] += 1
                            new_cell.indices = tuple(indices)
                            break
                    else:
                        continue
                    break
            
    def distances(self, plane):
        ddict = {}
        for v in self.verts:
            dist = plane.signed_distance(v.point)
            ddict.setdefault(roundeps(dist), []).append(dist)
        return sorted(sum(dists)/len(dists) for dists in ddict.values())
        
    def scale(self, value):
        for v in self.verts:
            v.point = v.point.scaled(value)
            
    def remove_cells(self, func):
        self._cells = [c for c in self._cells if not func(c)]