Voronoi.FCMacro

"""
This macro is inteneded to be used in the Sketcher Workbench.

The idea is to place down some points, which are the center points of the Voronoi cells.
The macro will compute the convex hull of the points, a Voronoi tesselation
and additionally offset all Voronoi cells in such a way, that they can be used
as pockets.

An additional switch exists to draw the convex hull as line segments too, which
can be used to create a extrudeable sketch.

Edit the settings in the main() function for the creation of the Voronoi tesselation
and run the macro inside an active sketch.
"""
__Name__ = 'PDVoronoiFace'
__Comment__ = 'Create Voronoi Pattern on Face'
__Author__ = 'Sebastian Bachmann'
__Version__ = '0.0.1'
__Date__ = '2020-02-28'
__License__ = 'MIT'
__Web__ = 'https://github.com/reox/FreeCAD_Macros'
__Wiki__ = ''
__Icon__ = 'PDVoronoi.svg'
__Help__ = 'Select a face'
__Status__ = 'Beta'
__Requires__ = '>=0.19.18234; py3 only'
__Communication__ = 'https://github.com/reox/FreeCAD_Macros/issues'
"""
Copyright (c) 2019, Sebastian Bachmann <hello@reox.at>

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
import Part
import numpy as np
from numpy.linalg import det, norm
from scipy.spatial import Voronoi, Delaunay, ConvexHull

from PySide.QtGui import (
        QMainWindow,
        QDialog,
        QLineEdit,
        QPushButton,
        QGridLayout,
        QLabel,
        QGroupBox,
        QSpacerItem,
        QMessageBox,
        )
from PySide import QtCore


def msg(title, message):
    diag = QMessageBox(QMessageBox.Warning, title, message)
    diag.setWindowModality(QtCore.Qt.ApplicationModal)
    diag.exec_()


def sum_conserving_round(inp):
    """
    Round a list but conserve the overall sum
    by minimizing the roundoff error
    """
    frac, res = np.modf(inp)
    items = np.argsort(frac)[::-1]
    for i in range(int(round(np.sum(inp) - np.sum(res)))):
        res[items[i]] += 1.0
    return res


def triangle_area(v0, v1, v2):
    """
    Returns the surface area of a triangle with given points

    :param np.array v0: point 1
    :param np.array v1: point 2
    :param np.array v2: point 3
    """
    return norm(np.cross(v1 - v0, v2 - v0)) * 0.5


def Det(u, v):
    """
    Calculate the determinate of two vectors concatenated to matrix
    """
    return det(np.vstack((u, v)))


def point_in_triangle(p, v0, v1, v2):
    """
    Checks if the point p is part of the triangle formed by v0, v1, v2
    v0 is a point and v1 and v2 are vectors from v0 to the other two vertices

    From http://mathworld.wolfram.com/TriangleInterior.html
    """

    a = (Det(p, v2) - Det(v0, v2)) / Det(v1, v2)
    b = -(Det(p, v1) - Det(v0, v1)) / Det(v1, v2)

    return a > 0 and b > 0 and (a + b) < 1


def mirror_point_at_point(p, s):
    """
    Mirrors the point p at the point s
    """
    return p - 2*(p - s)


def fill_triangle(v0, v1, v2, n=100):
    """
    Fills a triangle with points.
    The points stem from a uniform distribution.

    From: http://mathworld.wolfram.com/TrianglePointPicking.html
    """
    a1 = np.random.random(n)[:, np.newaxis]
    a2 = np.random.random(n)[:, np.newaxis]

    # symmetry point at the parallelogram's edge
    s = v1 + ((v2 - v1) * 0.5)

    # Get vectors from v0 to v1 and v2
    v01 = (v1 - v0)
    v02 = (v2 - v0)

    x = v0 + a1 * v01 + a2 * v02
    x = np.array([p if point_in_triangle(p, v0, v01, v02) else mirror_point_at_point(p, s) for p in x])
    return x.T

"""
Function by @pv (Pauli Virtanen)
https://gist.github.com/pv/8036995
"""
def voronoi_finite_polygons_2d(vor, radius=None):
    """
    Reconstruct infinite voronoi regions in a 2D diagram to finite
    regions.
    Parameters
    ----------
    vor : Voronoi
        Input diagram
    radius : float, optional
        Distance to 'points at infinity'.
    Returns
    -------
    regions : list of tuples
        Indices of vertices in each revised Voronoi regions.
    vertices : list of tuples
        Coordinates for revised Voronoi vertices. Same as coordinates
        of input vertices, with 'points at infinity' appended to the
        end.
    """

    if vor.points.shape[1] != 2:
        raise ValueError("Requires 2D input")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max()*2

    # Construct a map containing all ridges for a given point
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    # Reconstruct infinite regions
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]

        if all(v >= 0 for v in vertices):
            # finite region
            new_regions.append(vertices)
            continue

        # reconstruct a non-finite region
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # finite ridge: already in the region
                continue

            # Compute the missing endpoint of an infinite ridge

            t = vor.points[p2] - vor.points[p1] # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius

            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())

        # sort region counterclockwise
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        # finish
        new_regions.append(new_region.tolist())
        
    return new_regions, np.asarray(new_vertices)


def fill_voronoi_offset(sketch, points, origin_face, offset):
    """
    Create Voronoi cells from a list of points given by XY coordinates
    Z will always be set to zero.
    These are sketch coordinates!

    .. todo::
        The method seems to miss out some points, which means that the surface might not completely
        filled with voronoi regions.

    :param Sketcher.SketchObject sketch: The sketch object to place the cells into
    :param List[Tuple[float]] points: A list of 2D points to form the centers of the cells
    :param Part.Face origin_face: Originating face in sketch coordinates
    :param float offset: the offset in Units to offset all lines. Positive offsets make the cells smaller
    """
    vor = Voronoi(points)
    
    regions, vertices = voronoi_finite_polygons_2d(vor)

    # Debug: show all points
    #for p in points:
    #    sketch.addGeometry(Part.Point(App.Vector(*p, 0)))

    for r in regions:
        # Simply remove points which are outside of the face...
        #while -1 in r:
        #   r.remove(-1)

        if r == []:
            continue

        # Make a wire and calculate the offset
        wire = Part.makePolygon([App.Vector(*vertices[x], 0) for x in r]+ [App.Vector(*vertices[r[0]], 0)])
        face = Part.makeFace(wire, 'Part::FaceMakerSimple')

        section = face.common(origin_face)

        # If the face has holes, there might be a different number of wires
        # than one...
        for new_wire in section.Wires:
            try:
                newshape = new_wire.makeOffset2D(-offset, openResult=False, intersection=True)
            except:
                print("Shape was probly too small... ignoring")
            else:
                for e in newshape.Edges:
                    a, b = e.Vertexes
                    sketch.addGeometry(Part.LineSegment(a.Point, b.Point))


class VoronoiForm(QDialog):
    def __init__(self, parent=None):
        super().__init__(parent)
        self.setWindowTitle('Voronoi Pattern on Face')

        layout = QGridLayout()

        layout.addWidget(QLabel('Points:', self), 0, 0)
        layout.addWidget(QLabel('Offset:', self), 1, 0)
        layout.addWidget(QLabel('Tessellation Tol:', self), 2, 0)
        layout.addWidget(QLabel('Seed:', self), 3, 0)

        newseed = QPushButton("Random")
        newseed.clicked.connect(self.setnewseed)
        layout.addWidget(newseed, 3, 2)

        self.points = QLineEdit('100')
        self.offset = QLineEdit('0.3')
        self.tol = QLineEdit('0.1')
        self.seed = QLineEdit('1337')

        layout.addWidget(self.points, 0, 1)
        layout.addWidget(self.offset, 1, 1)
        layout.addWidget(self.tol, 2, 1)
        layout.addWidget(self.seed, 3, 1)

        ok = QPushButton("OK")
        ok.clicked.connect(self.do)
        abort = QPushButton("Abort")
        abort.clicked.connect(self.exit)

        layout.addWidget(ok, 4, 0)
        layout.addWidget(abort, 4, 1)

        self.setLayout(layout)

        self.show()

    def setnewseed(self):
        self.seed.setText(str(np.random.randint(0, np.iinfo(np.int32).max)))

    @staticmethod
    def get_field_as(field, fun):
        err = False
        n = None
        try:
            n = fun(field.text())
        except ValueError:
            msg('Error', 'Invalid Value "{}", must be of type {}!'.format(field.text(), fun))
            err = True
        return n, err

    def do(self):
        err = False
        n, e1 = self.get_field_as(self.points, int)
        err |= e1
        offset, e1 = self.get_field_as(self.offset, float)
        err |= e1
        tol, e1 = self.get_field_as(self.tol, float)
        err |= e1
        seed, e1 = self.get_field_as(self.seed, int)

        if not err:
            np.random.seed(seed)
            voronoi_on_face(n, offset, tol)
            self.close()

    def exit(self, event):
        self.close()


def voronoi_on_face(n, offset, tol):
    """
    Create a Voronoi Pattern on a Face by creating random points.

    A sketch is created which contains the Voronoi cell geometry.

    The face is tessellated first, then the triangles are filled with random points
    choosen from an uniform distribution.
    Afterwards, each Voronoi Cell is offsetted in such a way, that ridges form between
    the cells. The width of each ridge is offset/2.

    Note, that the tessellation of the face might produce unwanted artefacts!
    Voronoi cells might intersect with features (holes) of the face.
    This can be explained by the tessellation: Consider a perfect circle, which is then
    tessellated. Now a finite number of triangles will be used, and the form
    of the circle will be approximated with a n-gon.
    During the calculation these n-gons are used! As the edges of the n-gon
    are always inside the circle, some Voronoi cells might intersect the circle.

    The face will in most cases also not filled completely with cells.
    The choise of points influences the end result. Unfortunately, I'm not entirely sure
    how to control the cells in such a way that cells are created on all the face.

    Right now, there is also no way to control the minimum size of each cell.
    One method might be to distribute the points more evenly - which also depends
    on the tessellation!

    :param int n: Number of points to generate
    :param float offset: Offset for voronoi cells, larger than 0
    :param float tol: Tessellation tolerance
    """
    sel = Gui.Selection.getSelectionEx()

    if len(sel) != 1:
        msg('Error', 'Please select a single face first!')
        return

    so = sel[0].SubObjects
    if len(so) != 1:
        msg('Error', 'Please select a single face first!')
        return

    face = so[0]
    print(face.Wires)

    poly_points, simplices = face.tessellate(tol)
    # We need lists instead if tuples for numpy to be able to return items from
    # an array
    simplices = [list(x) for x in simplices]
    total_area = face.Area

    # Get the object which is selected
    feature = Gui.Selection.getSelection()[0]
    face_name = sel[0].SubElementNames[0]

    body = Gui.ActiveDocument.ActiveView.getActiveObject('pdbody')
    if body is None:
        msg('Error', 'Need an active body!')
        return

    sketch = App.ActiveDocument.addObject('Sketcher::SketchObject', 'Voronoi_Sketch')
    body.addObject(sketch)
    sketch.MapMode = 'FlatFace'
    sketch.Support = (feature, face_name)

    # We need the placement to calculate the points in the sketch...
    # This transform will convert all points on the surface into 2D sketch
    # coordinates
    transform = sketch.Placement.toMatrix().inverse()

    # NOTE: The transformation should always make the point have zero at the z
    # coordinate. However, sometimes the z coordinate is <<1 but not 0.
    # But we are confident that this should work for all faces ;)
    poly = np.array([[p.x, p.y] for p in [transform * p for p in poly_points]])

    n_points = sum_conserving_round(np.array([triangle_area(*poly[s]) / total_area * n for s in simplices])).astype(np.uint)
    if np.sum(n_points) != n:
        msg('Warning', 'Point numbers do not match after distribution step!')

    points = np.empty((2, 0))
    for i, s in enumerate(simplices):
        # If the tesselation produces too small triangles, they might not get
        # filled with points
        if n_points[i] == 0:
            continue
        points = np.hstack((points, fill_triangle(*poly[s], n=n_points[i])))
    points = points.T

    fill_voronoi_offset(sketch, points, face.transformed(transform), offset)

    App.ActiveDocument.recompute()



if __name__ == "__main__":
    form = VoronoiForm()