HW04: Convex polygon with random order

hw04.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from functools import cmp_to_key
from pylab import *
from scipy import stats


def standard_simplex_draw_point(n):
    sample = np.random.random_sample(2 * n)
    for a, b in zip(*[iter(sample)]*2):
        if a > b:
            a, b = b, a
        yield np.array([a, 1 - b])


def transition_matrix(tri):
    (x1, y1), (x2, y2), (x3, y3) = tri
    Q = np.array([
        [x2 - x1, y2 - y1],
        [x3 - x1, y3 - y1]
    ])
    return Q


def generate(points, n, plot_=False):
    center = [np.mean(points[:, 0]), np.mean(points[:, 1])]

    def points_cmp(p1, p2):
        angle1 = np.arctan2(p1[1] - center[1], p1[0] - center[0])
        angle2 = np.arctan2(p2[1] - center[1], p2[0] - center[0])
        return angle1 - angle2

    points = np.array(sorted(points, key=cmp_to_key(points_cmp)))
    if plot_:
        plot(points[:, 0], points[:, 1], '.')

    tris = []
    areas = []
    for i in range(1, points.shape[0] - 1):
        if (np.linalg.norm(points[i + 1] - points[i]) + np.linalg.norm(points[i] - points[0])) == \
                np.linalg.norm(points[i + 1] - points[0]):
            print('Degenerate triangle found')
            continue
        tri = np.array([points[0], points[i], points[i + 1]])
        Q = transition_matrix(tri)
        area = 0.5 * np.linalg.det(Q)
        tris.append(tri)
        areas.append(area)

    total = sum(areas)
    probs = [area / total for area in areas]

    tri_chooser = stats.rv_discrete(values=(np.arange(len(probs)), probs))

    if plot_:
        for tri in tris:
            plot(np.append(tri[:, 0], tri[0, 0]), np.append(tri[:, 1], tri[0, 1]))

    tri_indices = tri_chooser.rvs(size=n)
    standard_points = standard_simplex_draw_point(n)

    result = []
    for tri_index, standard_point in zip(tri_indices, standard_points):
        tri = tris[tri_index]

        Q = transition_matrix(tri)
        point = np.dot(standard_point, Q) + tri[0]

        result.append(point)

        if plot_:
            plot(point[0], point[1], '.', ms=10)

    if plot_:
        show()

    return result


def main():
    # assume that polygon is convex
    points = np.array([
        [-1, 4],
        [4, 0],
        [-0.5, 8],
        [3, 8],
        [1, 10.5],
        [2, 11],
        [1.5, 2]
    ])
    result = generate(points, 1000, True)
    print(result)


if __name__ == '__main__':
    main()