Piecewise affine warping using SciPy

As we saw in the previous post, three point correspondences uniquely defines an affine transformation and therefore affine warping of triangles can be done exactly. Let's look at the most common form of warping between a set of corresponding points, piecewise affine warping. Given any image with landmark points we can warp that image to corresponding landmarks in another image by triangulating the points in a mesh and warping each triangle with an affine transform. This is standard stuff for any graphics and image processing library. It is also possible to do it using SciPy and matplotlib.

To triangulate points, Delaunay triangulation is often used. This triangulation comes included in Matplotlib (but outside the PyLab part) and can be used like this.

import matplotlib.delaunay as md 
from pylab import *
from numpy import *

x,y =  array(random.standard_normal((2,100)))
centers,edges,tri,neighbors = md.delaunay(x,y)

figure()
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] #add first point to end
    plot(x[t_ext],y[t_ext],'r')

plot(x,y,'*')
axis('off')
show()

The image below shows some example points and the resulting triangulation. Delaunay triangulation chooses the triangles so that the edges are the dual graph of a Voronoi diagram. There are four outputs of delaunay() of which we only need the list of triangles. Create a function in warp.py for the triangulation.

import matplotlib.delaunay as md
    
def triangulate_points(x,y):
    """delaunay triangulation of 2D points"""

    centers,edges,tri,neighbors = md.delaunay(x,y)
    return tri

The output is an array with each row containing the indices for the three points of each triangle.



Let's look at an example of warping an image to a non-flat object in another image using 30 control points in a 5 by 6 grid. The example below shows an image to be warped to the facade of the "turning torso". The target points were manually selected using ginput().

First we need a general warp function for piecewise affine image warping. The code below does the trick, where we also take the opportunity to learn how to warp color images (you simply warp each color channel).

def pw_affine(fromim,toim,fp,tp,tri):
    """ warp triangular patches from an image.
        fromim = image to warp 
        toim = destination image
        fp = from points in hom. coordinates
        tp = to points in hom.  coordinates
        tri = triangulation"""

    im = toim.copy()

    #check if image is grayscale or color
    is_color = len(fromim.shape) == 3

    #create image to warp to (needed if iterate colors)
    im_t = zeros(im.shape, 'uint8') 

    for t in tri:
        #compute affine transformation
        H = homography.Haffine_from_points(tp[:,t],fp[:,t])

        if is_color:
            for col in range(3):
                im_t[:,:,col] = ndimage.affine_transform(fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        else:
            im_t = ndimage.affine_transform(fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])

        #alpha for triangle
        alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
        #add triangle to image
        if is_color:
            for col in range(3):
                im[:,:,col] = (1-alpha)*im[:,:,col] + alpha*im_t[:,:,col]
        else:
            im = (1-alpha)*im + alpha*im_t

    return im

Add this function to warp.py (here homography.Haffine_from_points() refers to the function in the previous post). You also need:

def alpha_for_triangle(points,m,n):
    """ creates alpha map of size (m,n) 
        for a triangle with corners defined by points
        (given in normalized homogeneous coordinates)."""

    alpha = zeros((m,n))
    for i in range(min(points[0]),max(points[0])):
        for j in range(min(points[1]),max(points[1])):
            x = linalg.solve(points,[i,j,1])
            if min(x) > 0: #all coefficients positive
                alpha[i,j] = 1
                
    return alpha

To use this on the current example, just do

from PIL import Image
from pylab import *
from numpy import *
import homography
import warp

#open image to warp
fromim = array(Image.open('jes.jpg'))
x,y = meshgrid(range(5),range(6))
x = (fromim.shape[1]/4) * x.flatten()
y = (fromim.shape[0]/5) * y.flatten()

#triangulate
tri = warp.triangulate_points(x,y)

#open image and destination points
im = array(Image.open('turningtorso1.jpg'))
tp = loadtxt('turningtorso1_points.txt') #destination points

#convert points to hom. coordinates
fp = vstack((y,x,ones((1,len(x))))) 
tp = vstack((tp[:,1],tp[:,0],ones((1,len(tp)))))

#warp triangles
im = warp.pw_affine(fromim,im,fp,tp,tri)

The triangles can be plotted, as in the simple example above, with the following code (add this to warp.py).

def plot_mesh(x,y,tri):
    """ plot triangles""" 

    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]] #add first point to end
        plot(x[t_ext],y[t_ext],'r')

Example of placing an image on the "turning torso" using a 5 by 6 grid of landmark points. Image credits: rutgerblom from Flickr (original, license).

If you want to try this yourself, here is the target image and control points (you can use any source image).