In a previous post I 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 ?
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 in 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 moving points. In the top image, 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 shown 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?
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. You get a Delaunay triangulation if you connect the points in neighboring Voronoi regions!
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])
Create an instance of Subdiv2D with the rectangle obtained in the previous step.
C++
Subdiv2D subdiv(rect);
Python
subdiv = cv2.Subdiv2D(rect);
Insert the points into subdiv using subdiv.insert(point). The video above shows an animation of triangulation as we add points to subdiv.Use subdiv.getTriangleList to get a list of Delaunay triangles.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 with OpenCV and simplified and modified it to suit our purpose. The python example that comes with OpenCV uses the old (and ugly) interface, 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
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)
Image Credits
- The image of President Obama was is licensed under CC By 3.0 [ link ]
- The images used in Figure 2. are two frames of a gif licensed under CC By-SA 3.0 [ link ]
- The image in Figure 3. is licensed under CC BY-SA 4.0 [ link ]
