Delaunay Triangulation and Voronoi Diagram using OpenCV ( C++ / Python )

By | November 5, 2015
opencv delaunay vornoi subdiv example

Figure 1. Left : Image of President Obama with landmarks detected using dlib. Center : Delaunay triangulation of the landmarks. Right : Corresponding Voronoi Diagram.

In a previous post I had discussed two libraries for facial landmark detection, and had pointed to several interesting applications like Face Morphing, Face Replacement etc. that use facial landmarks. In many such applications a triangulation of facial landmarks is first found (See Figure 1), and these triangles are warped to do something interesting. This post will help us understand Delaunay triangulation and Voronoi diagrams ( a.k.a Voronoi tesselation, Voronoi decomposition, Voronoi partition, and Dirichlet tessellation ), and help us uncover a barely documented function in OpenCV.

Boris Nikolaevich Delaunay, a Russian mathematician, spelled his last name in two different ways — Delaunay, for French and German publications, and Delone elsewhere. Few people have algorithms or concepts named after them. But Delaunay has one mathematical concept named after each spelling of his last name! — Delaunay Triangulation and Delone Sets. As if all this was not enough, by 1913, he became one of the top three Russian mountain climbers!

Delaunay’s Ph.D. advisor was Georgy Voronoy after whom the Voronoi Diagrams are named. And here is another interesting bit of trivia — Voronoy’s Ph.D. advisor was Andrey Markov (yes, yes Markov of Markov chains, and Markov processes fame).

What is Delaunay Triangulation ?

Delaunay triangulation favors small angles

Figure 2 : Delaunay triangulation favors small angles

Given a set of points in a plane, a triangulation refers to the subdivision of the plane into triangles, with the points as vertices. In Figure 1, we see a set of landmarks on the left image, and the triangulation in the middle image. A set of points can have many possible triangulations, but Delaunay triangulation stands out because it has some nice properties. In a Delaunay triangulation, triangles are chosen such that no point is inside the circumcircle of any triangle. Figure 2. shows Delaunay triangulation of 4 points A, B, C and D.  In the top image, for the triangulation to be a valid Delaunay triangulation, point C should be outside the circumcircle of triangle ABD, and point A should be outside the circumcircle of triangle BCD.

An interesting property of Delaunay triangulation is that it does not favor “skinny” triangles ( i.e. triangles with one large angle ).

Figure 2 shows how the triangulation changes to pick “fat” triangles when the points are moved. In the top image, the points B and D have their x-coordinates at x = 1.5, and in the bottom image they are moved to the right to x = 1.75. In the top image angles ABC and ABD are large, and Delaunay triangulation creates an edge between B and D splitting the two large angles into smaller angles ABD, ADB, CDB, and CBD. On the other hand in the bottom image, the angle BCD is too large, and Delaunay triangulation creates an edge AC to divide the large angle.

There are many algorithms to find the Delaunay triangulation of a set of points. The most obvious ( but not the most efficient ) one is to start with any triangulation, and check if the circumcircle of any triangle contains another point. If it does, flip the edges ( as show in Figure 2. ) and continue until there are no triangles whose circumcircle contains a point.

Any discussion on Delaunay triangulation has to include Voronoi diagrams because the Voronoi diagram of a set of points is mathematical dual to its Delaunay triangulation.

What is a Voronoi Diagram ?

Voronoi-Diagram

Figure 3. Voronoi Diagram

Given a set of points in a plane, a Voronoi diagram partitions the space such that the boundary lines are equidistant from neighboring points. Figure 3. shows an example of a Voronoi diagram calculated from the points shown as black dots. You will notice that every boundary line passes through the center of two points. If you connect the points in neighboring Voronoi regions, you get a Delaunay triangulation!

Delaunay triangulation and Voronoi diagram are related in more ways than one.  Georgy Voronoy, the mathematician after which Voronoi diagram is named, was Boris Delaunay’s Ph.D. advisor.

Delaunay Triangulation & Voronoi Diagram in OpenCV

Given a set of points, you can calculate the Delaunay Triangulation or Voronoi Diagram using the class Subdiv2D. Here are the steps. A complete working example is shown in the next section.

  1. Collect all the points in a vector.
    C++

    vector<Point2f> points;
    // This is how you can add one point. 
    points.push_back(Point2f(x,y));
    

    Python

    points = []
    # This is how you can add one point. 
    points.append((x, y))
    
  2. Define the space you want to partition using a rectangle (rect). If the points you have defined in the previous step are defined on an image, this rectangle can be ( 0, 0, width, height ). Otherwise, you can choose a rectangle that bounds all the points.

    C++

     
    Mat img = imread("image.jpg");
    Size size = img.size();
    Rect rect(0, 0, size.width, size.height);
    

    Python

    img = cv2.imread("image.jpg");
    size = img.shape
    rect = (0, 0, size[1], size[0])
    
  3. Create an instance of Subdiv2D with the rectangle obtained in the previous step.
    C++

     
    Subdiv2D subdiv(rect); 
    

    Python

     
    subdiv  = cv2.Subdiv2D(rect); 
    
  4. Insert the points into subdiv using subdiv.insert(point). The video above shows an animation of triangulation as we add points to subdiv.
  5. Use subdiv.getTriangleList to get a list of Delaunay triangles.
  6. Use subdiv.getVoronoiFacetList to get the list of Voronoi facets.

OpenCV Example for Delaunay Triangulation & Voronoi Diagram

Here is a complete working example. I have copied some of this code from examples that come with OpenCV and simplified and modified it to suit our purpose. The python example that comes with OpenCV uses the old (and ugly) interface, and so I wrote it from scratch. This code assumes an image is stored in image.jpg and the points are stored in points.txt. Each row of points.txt contains the x and y coordinates of a point separated by a space. E.g.

207 242
210 269
214 297
220 322
229 349

You can also download this code along with the associated data ( obama example in this post ) by subscribing to our newsletter here.

C++ Example

#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <iostream>
#include <fstream>

using namespace cv;
using namespace std;

// Draw a single point
static void draw_point( Mat& img, Point2f fp, Scalar color )
{
    circle( img, fp, 2, color, CV_FILLED, CV_AA, 0 );
}

// Draw delaunay triangles
static void draw_delaunay( Mat& img, Subdiv2D& subdiv, Scalar delaunay_color )
{

    vector<Vec6f> triangleList;
    subdiv.getTriangleList(triangleList);
    vector<Point> pt(3);
    Size size = img.size();
    Rect rect(0,0, size.width, size.height);

    for( size_t i = 0; i < triangleList.size(); i++ )
    {
        Vec6f t = triangleList[i];
        pt[0] = Point(cvRound(t[0]), cvRound(t[1]));
        pt[1] = Point(cvRound(t[2]), cvRound(t[3]));
        pt[2] = Point(cvRound(t[4]), cvRound(t[5]));
        
        // Draw rectangles completely inside the image.
        if ( rect.contains(pt[0]) && rect.contains(pt[1]) && rect.contains(pt[2]))
        {
            line(img, pt[0], pt[1], delaunay_color, 1, CV_AA, 0);
            line(img, pt[1], pt[2], delaunay_color, 1, CV_AA, 0);
            line(img, pt[2], pt[0], delaunay_color, 1, CV_AA, 0);
        }
    }
}

//Draw voronoi diagram
static void draw_voronoi( Mat& img, Subdiv2D& subdiv )
{
    vector<vector<Point2f> > facets;
    vector<Point2f> centers;
    subdiv.getVoronoiFacetList(vector<int>(), facets, centers);

    vector<Point> ifacet;
    vector<vector<Point> > ifacets(1);

    for( size_t i = 0; i < facets.size(); i++ )
    {
        ifacet.resize(facets[i].size());
        for( size_t j = 0; j < facets[i].size(); j++ )
            ifacet[j] = facets[i][j];

        Scalar color;
        color[0] = rand() & 255;
        color[1] = rand() & 255;
        color[2] = rand() & 255;
        fillConvexPoly(img, ifacet, color, 8, 0);

        ifacets[0] = ifacet;
        polylines(img, ifacets, true, Scalar(), 1, CV_AA, 0);
        circle(img, centers[i], 3, Scalar(), CV_FILLED, CV_AA, 0);
    }
}


int main( int argc, char** argv)
{

    // Define window names
    string win_delaunay = "Delaunay Triangulation";
    string win_voronoi = "Voronoi Diagram";
    
    // Turn on animation while drawing triangles
    bool animate = true;

    // Define colors for drawing.
    Scalar delaunay_color(255,255,255), points_color(0, 0, 255);
    
    // Read in the image.
    Mat img = imread("image.jpg");
    
    // Keep a copy around
    Mat img_orig = img.clone();
    
    // Rectangle to be used with Subdiv2D
    Size size = img.size();
    Rect rect(0, 0, size.width, size.height);

    // Create an instance of Subdiv2D
    Subdiv2D subdiv(rect);
    
    // Create a vector of points.
    vector<Point2f> points;

    // Read in the points from a text file
    ifstream ifs("points.txt");
    int x, y;
    while(ifs >> x >> y)
    {
        points.push_back(Point2f(x,y));
    }
    
    // Insert points into subdiv
    for( vector<Point2f>::iterator it = points.begin(); it != points.end(); it++)
    {
        subdiv.insert(*it);
        // Show animation
        if (animate)
        {
            Mat img_copy = img_orig.clone();
            // Draw delaunay triangles
            draw_delaunay( img_copy, subdiv, delaunay_color );
            imshow(win_delaunay, img_copy);
            waitKey(100);
        }
        
    }
    
    // Draw delaunay triangles
    draw_delaunay( img, subdiv, delaunay_color );

    // Draw points
    for( vector<Point2f>::iterator it = points.begin(); it != points.end(); it++)
    {
        draw_point(img, *it, points_color);
    }
    
    // Allocate space for Voronoi Diagram
    Mat img_voronoi = Mat::zeros(img.rows, img.cols, CV_8UC3);
    
    // Draw Voronoi diagram
    draw_voronoi( img_voronoi, subdiv );

    // Show results.
    imshow( win_delaunay, img);
    imshow( win_voronoi, img_voronoi);
    waitKey(0);

    return 0;
}

Python Example

 
#!/usr/bin/python

import cv2
import numpy as np
import random

# Check if a point is inside a rectangle
def rect_contains(rect, point) :
    if point[0] < rect[0] :
        return False
    elif point[1] < rect[1] :
        return False
    elif point[0] > rect[2] :
        return False
    elif point[1] > rect[3] :
        return False
    return True

# Draw a point
def draw_point(img, p, color ) :
    cv2.circle( img, p, 2, color, cv2.cv.CV_FILLED, cv2.CV_AA, 0 )


# Draw delaunay triangles
def draw_delaunay(img, subdiv, delaunay_color ) :

    triangleList = subdiv.getTriangleList();
    size = img.shape
    r = (0, 0, size[1], size[0])

    for t in triangleList :
        
        pt1 = (t[0], t[1])
        pt2 = (t[2], t[3])
        pt3 = (t[4], t[5])
        
        if rect_contains(r, pt1) and rect_contains(r, pt2) and rect_contains(r, pt3) :
        
            cv2.line(img, pt1, pt2, delaunay_color, 1, cv2.CV_AA, 0)
            cv2.line(img, pt2, pt3, delaunay_color, 1, cv2.CV_AA, 0)
            cv2.line(img, pt3, pt1, delaunay_color, 1, cv2.CV_AA, 0)


# Draw voronoi diagram
def draw_voronoi(img, subdiv) :

    ( facets, centers) = subdiv.getVoronoiFacetList([])

    for i in xrange(0,len(facets)) :
        ifacet_arr = []
        for f in facets[i] :
            ifacet_arr.append(f)
        
        ifacet = np.array(ifacet_arr, np.int)
        color = (random.randint(0, 255), random.randint(0, 255), random.randint(0, 255))

        cv2.fillConvexPoly(img, ifacet, color, cv2.CV_AA, 0);
        ifacets = np.array([ifacet])
        cv2.polylines(img, ifacets, True, (0, 0, 0), 1, cv2.CV_AA, 0)
        cv2.circle(img, (centers[i][0], centers[i][1]), 3, (0, 0, 0), cv2.cv.CV_FILLED, cv2.CV_AA, 0)


if __name__ == '__main__':

    # Define window names
    win_delaunay = "Delaunay Triangulation"
    win_voronoi = "Voronoi Diagram"

    # Turn on animation while drawing triangles
    animate = True
    
    # Define colors for drawing.
    delaunay_color = (255,255,255)
    points_color = (0, 0, 255)

    # Read in the image.
    img = cv2.imread("image.jpg");
    
    # Keep a copy around
    img_orig = img.copy();
    
    # Rectangle to be used with Subdiv2D
    size = img.shape
    rect = (0, 0, size[1], size[0])
    
    # Create an instance of Subdiv2D
    subdiv = cv2.Subdiv2D(rect);

    # Create an array of points.
    points = [];
    
    # Read in the points from a text file
    with open("points.txt") as file :
        for line in file :
            x, y = line.split()
            points.append((int(x), int(y)))

    # Insert points into subdiv
    for p in points :
        subdiv.insert(p)
        
        # Show animation
        if animate :
            img_copy = img_orig.copy()
            # Draw delaunay triangles
            draw_delaunay( img_copy, subdiv, (255, 255, 255) );
            cv2.imshow(win_delaunay, img_copy)
            cv2.waitKey(100)

    # Draw delaunay triangles
    draw_delaunay( img, subdiv, (255, 255, 255) );

    # Draw points
    for p in points :
        draw_point(img, p, (0,0,255))

    # Allocate space for Voronoi Diagram
    img_voronoi = np.zeros(img.shape, dtype = img.dtype)

    # Draw Voronoi diagram
    draw_voronoi(img_voronoi,subdiv)

    # Show results
    cv2.imshow(win_delaunay,img)
    cv2.imshow(win_voronoi,img_voronoi)
    cv2.waitKey(0)

Subscribe & Download Code

If you liked this article and would like to download code and example images used in this post, please subscribe to our newsletter. You will also receive a free Computer Vision Resource guide. In our newsletter we share OpenCV tutorials and examples written in C++/Python, and Computer Vision and Machine Learning algorithms and news.

Subscribe Now

Image Credits

  1. The image of President Obama was is licensed under CC By 3.0 [ link ]
  2. The images used in Figure 2. are two frames of a gif licensed under CC By-SA 3.0 [ link ]
  3. The image in Figure 3. is licensed under CC BY-SA 4.0 [ link ]
Category: Face how-to Tags: , , ,

About Satya Mallick

I am an entrepreneur with a love for Computer Vision and Machine Learning with a dozen years of experience (and a Ph.D.) in the field. In 2007, right after finishing my Ph.D., I co-founded TAAZ Inc. with my advisor Dr. David Kriegman and Kevin Barnes. The scalability, and robustness of our computer vision and machine learning algorithms have been put to rigorous test by more than 100M users who have tried our products.

  • Sun Steven

    Great! Expect your next post!

  • Sathish Anand

    Great Work.

  • Ruchika

    Nice Article. But how did you get the x, y points in obama.txt?

  • Shubham Jaiswal

    A very useful article for beginners.I wanted to know that can we extract the image region based on denaulay triangulation

    • Thanks Shubham. Let’s say you want to extract a single triangle from the traingulation. You can use cvFillConvexPoly to fill the triangle with white, and leave the rest of the image black. This serves as a mask to extract the triangle.

  • Hamid Kavianathar

    thanks for your post. excuse me, when I run your code, I received this error:
    AttributeError: ‘module’ object has no attribute ‘CV_AA’

    what should I do? thank you.

    • alpha

      replace it to cv2.LINE_AA but have many wrong module 😀

    • Try 16.

  • neenu

    its great………….

    • Thanks Neenu.

      • neenu

        Dear…. i am working with snake identification from it’s scales. so that i want to make Voronoi diagram with respect to the scale. That’s help me to identify each species…pls help me for making that………………….

  • Mario Lukas

    Great Blog, thank you for those interesting tutorials.Very helpful! Keep on posting.

  • neenu

    Dear…. i am working with snake identification from it’s scales. so that i want to make Voronoi diagram with respect to the scale. That’s help me to identify each species…pls help me for making that………………

    • Sorry for the late reply. I don’t think you actually need or should use voronoi diagrams for this problem. I would love to see the data. If you have the time, please send me an email at spmallick@learnopencv.com and I will suggest something.

  • Alf

    Dear Mr. Mallick.
    First of all, thank you for this amazing blog!

    I have a question about .txt-files which take place in projects such as FaceMorph, Delaunay etc. How to get necessary points from a picture and then put them in a .txt-file?
    Thank you.

    • Thanks for the kind words.

      I used dlib to get the landmark points. You can see my post for more details.

      http://www.learnopencv.com/facial-landmark-detection

      I am in the process of adding a few new files to dlib that will generate this file more easily. You will receive an email in a few days about that.

      • Alf

        Dear Mr. Mallick.
        Thank you so much for prompt response and interesting information!
        Thank you.

  • Mahmoud Mohamed Ouf

    Dear Dr.Mallick,
    I have already subscribed but I still can’t download the text file.
    Thanks a lot for this great website .

    • Hi Mahmoud,

      Not sure why you did not receive a link to the code. Please send an email to spmallick@learnopencv.com and I will respond with a link to the code. I am assuming you are talking about obama.txt ?

      • Mahmoud Mohamed Ouf

        Thanks a lot, I have already received an e-mail with the link.
        thanks again .

  • IBTISSEM

    WHY I FIND THIS ERROR `Subdiv2D’ has not been declared ‘

  • Abhijit Bhattacharyya

    Hello! Your post illustrates the point in simple way. I understand that use of opencv is the theme. Just I like to ask if you have any reference for voronoi and delaunay using simple C++ code only i.e. not using opencv classes. Thanks

  • pranjal

    How did you generate the points for an image?

  • Valent

    Hi Mr. Satya Mallick, thank you very much for this contribution.

    I have a question for you. Is there a way I can choose a certain edge on a triangle and swap it with another edge, here’s an example: http://i.imgur.com/RqAQKN5.png

    There’s a member function of a subdiv2D called swapedges but i’m not sure how it’s supposed to be used.

    If you could help me, I would be very thankful.

    • Delaunay triangulation will always choose the triangles shown in “before” because it is designed to choose fat triangles. But if you need the triangles to be like in “After” image, you can simply go through the list of triangles and remove the two triangles in “before” and add the two triangles in “after” to the list.

      • Valent

        Thank you, kind sir. I’ve already found another solution for it 🙂

  • Mohanraj V

    I want to know detail about how the facial landmarks points are detected.

  • Yarif

    Hi, I would appreciate any help about how to extract features from delaunay triangulation. I was trying to get a area from each triangle. However, I got a distinct total number of triangles and not in order. So, how can I deal with this? Thank you.

    • Hi Yarif,

      I am sorry I don’t understand the question. The triangles are in no specific order.

      • Yarif

        Hi Satya, thank you for your response. So, how can I use the delaunay triangulation to perform object recognition? Thank you.

  • Nibin Jose

    Hi, I would appreciate any help about triangulate one small object in a 2d image and how can I get the area of that triangle? How can I achieve that? Thank You.

  • Sıradan Kullanıcı

    First of all thank you very much for your work.
    How do I use these codes in the android app?
    Thanks in advance…

  • Sıradan Kullanıcı

    Error:
    Traceback (most recent call last):
    img_orig = img.copy();
    AttributeError: ‘NoneType’ object has no attribute ‘copy’

  • Sıradan Kullanıcı
  • Thái Nguyễn

    Hi, This post is very useful. But I face up to some problems. Can you help me?
    I used Dlib for ficial landmark detection and I can get 68 points as the photo. https://uploads.disquscdn.com/images/1d6b0e7996ff25031438caafa25ef26c3cd4a4316ae25c384dbce0f176afbb3d.png
    When I use these points to make Delaunay Triangulation, I can get list of 828 items => 138 triangulations. I don’t know why I get some wrong point as the photo I show you. Can you tell me the reason. Many thanks.

  • soyoung shin

    Hello. I want to download the ‘text.txt’ file.
    I click the subscribe button but I still can not get the file.
    Please help..

    • When you subscribe and confirm you email, we send you a link to all code and other files using in the post.

  • Diego Victor

    Hello, Satya.

    In this section:
    ifacets[0] = ifacet;
    polylines(img, ifacets, true, Scalar(), 1, CV_AA, 0);

    Are you using ifacets, which is a vector of a vector of points, because you need to pass it to polylines as the accepted argument?