A python implementation of delaunay algorithm

delaunay.py

#!/usr/bin/env python
# -*- coding:UTF-8 -*-
# File Name : face_morphing.py
# Creation Date : 04-04-2018
# Created By : Jeasine Ma [jeasinema[at]gmail[dot]com]

import cv2 
import json 
import numpy as np 

class Delaunay(object):
    def __init__(self):
        pass 

    @staticmethod 
    def load_landmark(fname):
        js = json.load(open(fname))
        ld = [[i['x'], i['y']] for i in list(js['faces'][0]['landmark'].values())]
        return np.array(ld)
    
    @staticmethod 
    def cal_delaunay_triangle(image_size, verticies):
        # Input:
        #   (h, w)
        #   (N, 2)
        h, w = image_size
        origin_length = verticies.shape[0]
        verticies = verticies.copy() 
        triangles, temp_triangles = set(), set()
        indices = np.arange(verticies.shape[0])
        indices = indices[np.argsort(verticies[:, 0])]
        # super triangle 
        a, b, c = [-10, h], [10+w, h], [w/2, -np.tan(np.pi/2 - np.arctan(10/h))*w/2]
        verticies = np.concatenate([verticies, np.array([a,b,c])], axis=0)
        indices = list(indices)
        [indices.append(len(indices)) for i in range(3)]
        triangles.add((indices[-3], indices[-1], indices[-2])) # a,c,b
        temp_triangles.add((indices[-3], indices[-1], indices[-2])) # a,c,b
        # main loop 
        for ind in indices[:-3]:
            edge_buffer = set()
            buf_temp_triangles = temp_triangles.copy()
            for tri_ind in temp_triangles:
                tri = verticies[list(tri_ind)]
                center, r = FaceMorphing.cal_circum_circle(tri)
                res = FaceMorphing.judge_circle_point_relation(center, r, verticies[ind, :])
                if res == 0:
                    triangles.add(tri_ind)
                    buf_temp_triangles.remove(tri_ind)
                elif res == 1:
                    continue 
                else: # a, c, b => (a,c), (a,b), (c,b)
                    if (tri_ind[0], tri_ind[1]) not in edge_buffer:
                        edge_buffer.add((tri_ind[0], tri_ind[1]))
                    else: 
                        edge_buffer.remove((tri_ind[0], tri_ind[1]))
                    if (tri_ind[0], tri_ind[2]) not in edge_buffer:
                        edge_buffer.add((tri_ind[0], tri_ind[2]))
                    else:
                        edge_buffer.remove((tri_ind[0], tri_ind[2]))
                    if (tri_ind[1], tri_ind[2]) not in edge_buffer:
                        edge_buffer.add((tri_ind[1], tri_ind[2]))
                    else:
                        edge_buffer.remove((tri_ind[1], tri_ind[2]))

                    buf_temp_triangles.remove(tri_ind)
            for edge in edge_buffer: 
                tmp_ind = np.array([ind, edge[0], edge[1]])
                points = np.array([verticies[tmp_ind[0]], verticies[tmp_ind[1]], verticies[tmp_ind[2]]])
                tmp_ind = tmp_ind[np.argsort(points[:, 0])]
                buf_temp_triangles.add((tmp_ind[0], tmp_ind[1], tmp_ind[2]))
            temp_triangles = buf_temp_triangles  

        triangles = triangles.union(temp_triangles)
        final_triangles = triangles.copy()
        for tri_ind in triangles:
            if np.sum(np.array(tri_ind) >= origin_length) > 0:
                final_triangles.remove(tri_ind)
        res = np.zeros((len(final_triangles), 3, 2))
        for ind, tri_ind in enumerate(final_triangles):
            res[ind] = verticies[list(tri_ind)]

        return res.astype(np.int32)
        
    @staticmethod 
    def cal_circum_circle(triangle):
        # Input:
        #   (3,2)
        # Output:
        #   (x,y), r
        x1, x2, x3 = triangle[:, 0]
        y1, y2, y3 = triangle[:, 1]
        a = np.sqrt(np.sum((triangle[0, :] - triangle[1, :])**2))
        b = np.sqrt(np.sum((triangle[0, :] - triangle[2, :])**2))
        c = np.sqrt(np.sum((triangle[1, :] - triangle[2, :])**2))
        S = (1/2)*a*b*np.sqrt(1 - ((a**2 + b**2 - c**2)/(2*a*b))**2)
        if S == 0:
            print(triangle)
            r = 0
            x, y = x1, y1
        else:
            r = (a*b*c)/(4*S)
            x = np.linalg.det([[x1**2 + y1**2, y1, 1],[x2**2 + y2**2, y2, 1],[x3**2 + y3**2, y3, 1]])/(2*np.linalg.det([[x1, y1, 1],[x2, y2, 1],[x3, y3, 1]]))
            y = np.linalg.det([[x1, x1**2 + y1**2, 1],[x2, x2**2 + y2**2, 1],[x3, x3**2 + y3**2, 1]])/(2*np.linalg.det([[x1, y1, 1],[x2, y2, 1],[x3, y3, 1]]))

        return (x, y), r

    @staticmethod 
    def judge_circle_point_relation(center, r, point):
        # Input:
        #   (x,y), r, (x,y)
        x1, y1 = center 
        x2, y2 = point 
        if (x2-x1)**2 + (y2-y1)**2 <= r**2:
            return 2
        elif x2 > x1:
            return 0
        else:
            return 1
    
    @staticmethod 
    def draw_delaunay(img, triangles):
        # Input:
        #   (h, w, 3), (N, 3, 2)
        img = img.copy()
        for ind in range(triangles.shape[0]):
            img = cv2.line(img, tuple(triangles[ind, 0, :]), tuple(triangles[ind, 1, :]), (255,0,0), 1)
            img = cv2.line(img, tuple(triangles[ind, 0, :]), tuple(triangles[ind, 2, :]), (255,0,0), 1)
            img = cv2.line(img, tuple(triangles[ind, 1, :]), tuple(triangles[ind, 2, :]), (255,0,0), 1)

        return img 

if __name__ == '__main__':
    alg = Delaunay()
    landmark = facemorph.load_landmark('landmark.txt')
    img = cv2.imread('img.png')

    ret = alg.cal_delaundry_triangle(img.shape[0:2], landmark)
    img = alg.draw_delaunay(img, ret)
    cv2.imwrite('res.png', img)