Delaunay tesselation in 2 dimensions.

delaunay.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import logging
import sys
from enum import Enum


logger = logging.getLogger(__name__)


class Orientation(Enum):
    COLLINEAR=0,
    CLOCKWISE=1,
    COUNTER_CLOCKWISE=2


def find_orientation(p1, p2, p3):
    det = (p3[0] - p2[0]) * (p2[1] - p1[1]) - (p3[1] - p2[1]) * (p2[0] - p1[0])
    if det < 0.0:
        return Orientation.COUNTER_CLOCKWISE

    if det > 0.0:
        return Orientation.CLOCKWISE

    return Orientation.COLLINEAR


class AABB(object):
    '''Axis-Aligned Bounding Box.'''
    def __init__(self, lower_bound, upper_bound):
        self._lower_bound = lower_bound
        self._upper_bound = upper_bound

    @property
    def lower_bound(self):
        return self._lower_bound

    @property
    def upper_bound(self):
        return self._upper_bound

    @property
    def width(self):
        return abs(self.upper_bound[0] - self.lower_bound[0])

    @property
    def height(self):
        return abs(self.upper_bound[1] - self.lower_bound[1])

    def calc_circumcircle(self):
        '''Calculates the center and the radius of the circumcircle of a rectangle.'''
        w = self.width
        h = self.height
        r = 0.5 * (w**2.0 + h**2.0)**0.5
    
        return (self.lower_bound[0] + w * 0.5, self.lower_bound[1] + h * 0.5), r

    @classmethod
    def from_points(cls, points):
        min_x = sys.float_info.max
        min_y = sys.float_info.max
        max_x = -sys.float_info.max
        max_y = -sys.float_info.max

        for point in points:
            if point[0] < min_x:
                min_x = point[0]
            if point[1] < min_y:
                min_y = point[1]
            if point[0] > max_x:
                max_x = point[0]
            if point[1] > max_y:
                max_y = point[1]

        return cls(lower_bound=(min_x, min_y), upper_bound=(max_x, max_y))


class Delaunay(object):
    '''Delaunay tesselation in 2 dimensions.'''
    def __init__(self, points):

        # Make an AABB from given points.
        aabb = AABB.from_points(points)

        # Make a super triangle from the AABB.
        center, radius = aabb.calc_circumcircle()
        p1, p2, p3 = self.calc_super_triangle(center=center, radius=radius)

        # Add the super triangle as the first element to the list.
        temporary_points = points + [p1, p2, p3]
        num_points = len(temporary_points)
        super_triangle = (num_points-3, num_points-2, num_points-1)
        self.triangles = [super_triangle, ]

        # Incremental method.
        for i, point in enumerate(points):
            logger.debug('-'*30)
            logger.debug('point: {}'.format(i))
            # A list of triangles newly generated by triangulation.
            new_triangles = list()

            j = 0
            while j < len(self.triangles):
                triangle = self.triangles[j]
                # Calculate the circumcircle of a triangle.
                p1 = temporary_points[triangle[0]]
                p2 = temporary_points[triangle[1]]
                p3 = temporary_points[triangle[2]]
                center, radius = self.calc_circumcircle_of_triangle(p1, p2, p3)

                if self.distance(point, center) < radius:
                    # The circumcircle of a triangle contains a point.
                    # A triangle used for triangulation is removed from the list.
                    logger.debug('*** src triangle: {}'.format(self.triangles[j]))
                    del self.triangles[j]
                    self.triangulation(triangles=new_triangles, triangle=triangle, index=i)
                else:
                    j += 1

            self.triangles.extend(new_triangles)
            logger.debug('triangles: {} \n{}'.format(len(self.triangles), self.triangles))

        # Remove all triangles with vertices of the super triangle.
        i = 0
        while i < len(self.triangles):
            triangle = self.triangles[i]
            if (super_triangle[0] in triangle) or (super_triangle[1] in triangle) or (super_triangle[2] in triangle):
                del self.triangles[i]
            else:
                i += 1

        # Remove all vertices of the super triangle.
        del temporary_points[super_triangle[0]:]

        # Validation.
        self.validate(points=temporary_points, triangles=self.triangles)

        logger.debug('-'*30)
        logger.debug('triangles: {} \n{}'.format(len(self.triangles), self.triangles))
        #self.points = temporary_points

    @classmethod
    def triangulation(cls, triangles, triangle, index):
        for n in range(3):
            tri = (triangle[n], triangle[(n+1)%3], index)
            logger.debug('new triangle: {}'.format(tri))
            i = 0
            while i < len(triangles):
                if set(tri) == set(triangles[i]):
                    # same permutation.
                    logger.debug('same permutation {} {}'.format(tri, triangles[i]))
                    del triangles[i]
                    break
                else:
                    i += 1
            else:
                triangles.append(tri)

    @staticmethod
    def calc_super_triangle(center, radius):
        sqrt3 = 3.0**0.5
        a = (center[0], center[1] + 2.0 * radius)
        b = (center[0] - sqrt3 * radius, center[1] - radius)
        c = (center[0] + sqrt3 * radius, center[1] - radius)

        return a, b, c

    @staticmethod
    def distance(p1, p2):
        return ((p2[0] - p1[0])**2.0 + (p2[1] - p1[1])**2.0)**0.5

    @classmethod
    def calc_circumcircle_of_triangle(cls, p1, p2, p3):
        '''Calculates the center and the radius of the circumcircle of a triangle.'''
        p3_p1 = (p3[0] - p1[0], p3[1] - p1[1])
        p2_p1 = (p2[0] - p1[0], p2[1] - p1[1])
        # When a triangle is degenerate the denominator is 0.
        denominator = 2.0 * (p2_p1[0] * p3_p1[1] - p2_p1[1] * p3_p1[0])

        t1 = p2[0]**2.0 - p1[0]**2.0 + p2[1]**2.0 - p1[1]**2.0
        t2 = p3[0]**2.0 - p1[0]**2.0 + p3[1]**2.0 - p1[1]**2.0

        c = (
            ( p3_p1[1] * t1 - p2_p1[1] * t2) / denominator,
            (-p3_p1[0] * t1 + p2_p1[0] * t2) / denominator
        )
        r = cls.distance(c, p1)

        return c, r

    @classmethod
    def validate(cls, points, triangles):
        #
        for i, c in enumerate(triangles):
            center, radius = cls.calc_circumcircle_of_triangle(
                points[c[0]],
                points[c[1]],
                points[c[2]]
            )
            for j, o in enumerate(triangles):
                if i == j:
                    continue
                if o[0] not in c:
                    assert cls.distance(points[o[0]], center) >= radius
                if o[1] not in c:
                    assert cls.distance(points[o[1]], center) >= radius
                if o[2] not in c:
                    assert cls.distance(points[o[2]], center) >= radius
        
        # Validate the orientation of each triangle.
        for triangle in triangles:
            p1 = points[triangle[0]]
            p2 = points[triangle[1]]
            p3 = points[triangle[2]]
            orientation = find_orientation(p1, p2, p3)
            assert orientation == Orientation.COUNTER_CLOCKWISE

main.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import colorsys
import logging
import random
from PIL import Image, ImageDraw, ImageFont
import delaunay


def make_random_points(num_points, lower_bound, upper_bound):
    points = set()
    while len(points) < num_points:
        x = random.uniform(lower_bound[0], upper_bound[0])
        y = random.uniform(lower_bound[1], upper_bound[1])
        points.add((x, y))

    return list(points)


def gen_distinct_colors(num_colors):
    '''Generate distinct colors.'''
    def hsv_to_rgb(h, s, v):
        (r, g, b) = colorsys.hsv_to_rgb(h, s, v)
        return int(r * 255), int(g * 255), int(b * 255)

    hue_partition = 1.0 / (num_colors + 1)
    return (hsv_to_rgb(value * hue_partition, 1.0, 1.0) for value in range(0, num_colors))


def display_as_image(width, height, points, tri):

    def calc_scale_and_offset(width, height):
        aspect = width / height
        hw = width * 0.5
        hh = height * 0.5
        return (hw / aspect, -hh), (hw, hh)

    aa_width = width * 2
    aa_height = height * 2

    scale, offset = calc_scale_and_offset(width=aa_width, height=aa_height)
    screen_coords = [
        (point[0] * scale[0] + offset[0], point[1] * scale[1] + offset[1])
        for point in points
    ]

    im = Image.new('RGB', (aa_width, aa_height), (128, 128, 128))
    draw = ImageDraw.Draw(im)

    #font = ImageFont.load_default()
    font = ImageFont.truetype('arial.ttf', 32)

    distinct_colors = list(gen_distinct_colors(len(tri.triangles)))

    for i, triangle in enumerate(tri.triangles):
        p1 = screen_coords[triangle[0]]
        p2 = screen_coords[triangle[1]]
        p3 = screen_coords[triangle[2]]
        draw.polygon(
            [p1, p2, p3],
            fill=distinct_colors[i],
            outline=(0, 0, 0)
        )

    for i, triangle in enumerate(tri.triangles):
        p1 = screen_coords[triangle[0]]
        p2 = screen_coords[triangle[1]]
        p3 = screen_coords[triangle[2]]
        c, r = delaunay.Delaunay.calc_circumcircle_of_triangle(p1, p2, p3)
        draw.arc(
            [c[0] - r, c[1] - r, c[0] + r, c[1] + r],
            start=0,
            end=360,
            fill=distinct_colors[i]
        )

    for i, screen_coord in enumerate(screen_coords):
        draw.text((screen_coord[0], screen_coord[1]), str(i), font=font)

    # super triangle.
    aabb = delaunay.AABB.from_points(points)
    center, radius = aabb.calc_circumcircle()
    p1, p2, p3 = delaunay.Delaunay.calc_super_triangle(center=center, radius=radius)
    s1 = (p1[0] * scale[0] + offset[0], p1[1] * scale[1] + offset[1])
    s2 = (p2[0] * scale[0] + offset[0], p2[1] * scale[1] + offset[1])
    s3 = (p3[0] * scale[0] + offset[0], p3[1] * scale[1] + offset[1])
    draw.line((s1, s2, s3, s1), fill=(0, 0, 0), width=1)
    draw.text(s1, str(len(points)+0), fill=(192, 192, 192), font=font)
    draw.text(s2, str(len(points)+1), fill=(192, 192, 192), font=font)
    draw.text(s3, str(len(points)+2), fill=(192, 192, 192), font=font)

    im = im.resize((width, height), Image.LANCZOS) # for anti-alias.
    im.show()


def main():
    points = make_random_points(
        num_points=2,
        lower_bound=(-0.2, -0.6),
        upper_bound=(+0.2, +0.2),
    )
    points = [(-0.3, -0.7), (0.3, -0.7), (0.3, 0.3), (-0.3, 0.3)] + points

    tri = delaunay.Delaunay(points)

    # draw
    display_as_image(width=640, height=480, points=points, tri=tri)


if __name__ == '__main__':
    logging.basicConfig(level=logging.DEBUG)
    main()