extracting a quadrilateral image to a rectangle

Question

BOUNTY UPDATE

Following Denis's link, this is how to use the threeblindmiceandamonkey code:

// the destination rect is our 'in' quad
int dw = 300, dh = 250;
double in[4][4] = {{0,0},{dw,0},{dw,dh},{0,dh}};
    // the quad in the source image is our 'out'
double out[4][5] = {{171,72},{331,93},{333,188},{177,210}};
double homo[3][6];
const int ret = mapQuadToQuad(in,out,homo);
    // homo can be used for calculating the x,y of any destination point
// in the source, e.g.
for(int i=0; i<4; i++) {
    double p1[3] = {out[i][0],out[i][7],1};
    double p2[3];
    transformMatrix(p1,p2,homo);
    p2[0] /= p2[2]; // x
    p2[1] /= p2[2]; // y
    printf("\t%2.2f\t%2.2f\n",p2[0],p2[1]);
}

This provides a transform for converting points in destination to the source - you can of course do it the other way around, but it's tidy to be able to do this for the mixing:

for(int y=0; y<dh; y++) {
    for(int x=0; x<dw; x++) {
        // calc the four corners in source for this
        // destination pixel, and mix

For the mixing, I'm using super-sampling with random points; it works very well, even when there is a big disparity in the source and destination area

---

BACKGROUND QUESTION

In the image at the top, the sign on the side of the van is not face-on to the camera. I want to calculate, as best I can with the pixels I have, what it'd look like face on.

I know the corner coordinates of the quad in the image, and the size of the destination rectangle.

I imagine that this is some kind of loop through the x and y axis doing a Bresenham's line on both dimensions at once with some kind of mixing as pixels in the source and destination images overlap - some sub-pixel mixing of some sort?

What approaches are there, and how do you mix the pixels?

Is there a standard approach for this?

Vizard

threeblindmiceandamonkey.com has vanished, but [github.com/marcantonio/imagekiln/blob/master/perspective.c] has a quad_to_quad() . – Vano

Answer 1

Look up "quad to quad" transform, e.g. threeblindmiceandamonkey.
A 3x3 transform on 2d homogeneous coordinates can transform any 4 points (a quad) to any other quad; conversely, any fromquad and toquad, such as the corners of your truck and a target rectangle, give a 3 x 3 transform.

Qt has quadToQuad and can transform pixmaps with it, but I guess you don't have Qt ?
Added 10Jun: from labs.trolltech.com/page/Graphics/Examples there's a nice demo which quad-to-quads a pixmap as you move the corners:

Added 11Jun: @Will, here's translate.h in Python (which you know a bit ? """ ...""" are multiline comments.)
perstrans() is the key; hope that makes sense, if not ask.

Bytheway, you could map the pixels one by one, mapQuadToQuad( target rect, orig quad ), but without pixel interpolation it'll look terrible; OpenCV does it all.

#!/usr/bin/env python
""" square <-> quad maps
    from http://threeblindmiceandamonkey.com/?p=16 matrix.h
"""

from __future__ import division
import numpy as np

__date__ = "2010-06-11 jun denis"

def det2(a, b, c, d):
    return a*d - b*c

def mapSquareToQuad( quad ):  # [4][2]
    SQ = np.zeros((3,3))
    px = quad[0,0] - quad[1,0] + quad[2,0] - quad[3,0]
    py = quad[0,1] - quad[1,1] + quad[2,1] - quad[3,1]
    if abs(px) < 1e-10 and abs(py) < 1e-10:
        SQ[0,0] = quad[1,0] - quad[0,0]
        SQ[1,0] = quad[2,0] - quad[1,0]
        SQ[2,0] = quad[0,0]
        SQ[0,1] = quad[1,1] - quad[0,1]
        SQ[1,1] = quad[2,1] - quad[1,1]
        SQ[2,1] = quad[0,1]
        SQ[0,2] = 0.
        SQ[1,2] = 0.
        SQ[2,2] = 1.
        return SQ
    else:
        dx1 = quad[1,0] - quad[2,0]
        dx2 = quad[3,0] - quad[2,0]
        dy1 = quad[1,1] - quad[2,1]
        dy2 = quad[3,1] - quad[2,1]
        det = det2(dx1,dx2, dy1,dy2)
        if det == 0.:
            return None
        SQ[0,2] = det2(px,dx2, py,dy2) / det
        SQ[1,2] = det2(dx1,px, dy1,py) / det
        SQ[2,2] = 1.
        SQ[0,0] = quad[1,0] - quad[0,0] + SQ[0,2]*quad[1,0]
        SQ[1,0] = quad[3,0] - quad[0,0] + SQ[1,2]*quad[3,0]
        SQ[2,0] = quad[0,0]
        SQ[0,1] = quad[1,1] - quad[0,1] + SQ[0,2]*quad[1,1]
        SQ[1,1] = quad[3,1] - quad[0,1] + SQ[1,2]*quad[3,1]
        SQ[2,1] = quad[0,1]
        return SQ

#...............................................................................
def mapQuadToSquare( quad ):
    return np.linalg.inv( mapSquareToQuad( quad ))

def mapQuadToQuad( a, b ):
    return np.dot( mapQuadToSquare(a), mapSquareToQuad(b) )

def perstrans( X, t ):
    """ perspective transform X Nx2, t 3x3:
        [x0 y0 1] t = [a0 b0 w0] -> [a0/w0 b0/w0]
        [x1 y1 1] t = [a1 b1 w1] -> [a1/w1 b1/w1]
        ...
    """
    x1 = np.vstack(( X.T, np.ones(len(X)) ))
    y = np.dot( t.T, x1 )
    return (y[:-1] / y[-1]) .T

#...............................................................................
if __name__ == "__main__":
    np.set_printoptions( 2, threshold=100, suppress=True )  # .2f

    sq = np.array([[0,0], [1,0], [1,1], [0,1]])
    quad = np.array([[171, 72], [331, 93], [333, 188], [177, 210]])
    print "quad:", quad
    print "square to quad:", perstrans( sq, mapSquareToQuad(quad) )
    print "quad to square:", perstrans( quad, mapQuadToSquare(quad) )

    dw, dh = 300, 250
    rect = np.array([[0, 0], [dw, 0], [dw, dh], [0, dh]])
    quadquad = mapQuadToQuad( quad, rect )
    print "quad to quad transform:", quadquad
    print "quad to rect:", perstrans( quad, quadquad )
"""
quad: [[171  72]
 [331  93]
 [333 188]
 [177 210]]
square to quad: [[ 171.   72.]
 [ 331.   93.]
 [ 333.  188.]
 [ 177.  210.]]
quad to square: [[-0.  0.]
 [ 1.  0.]
 [ 1.  1.]
 [ 0.  1.]]
quad to quad transform: [[   1.29   -0.23   -0.  ]
 [  -0.06    1.79   -0.  ]
 [-217.24  -88.54    1.34]]
quad to rect: [[   0.    0.]
 [ 300.    0.]
 [ 300.  250.]
 [   0.  250.]]
"""

Imago

I'm struggling and I've updated the question, but thx for the excellent lead – Vizard

the threeblindmiceandamonkey code with my own supersampling; the opencv bilinear-etc parameter gives me an uneasy feeling – Vizard

Answer 2

What you want is called planar rectification, and it's not all that simple, I'm afraid. What you need to do is recover the homography that maps this oblique view of the van side onto the front-facing view. Photoshop / etc. have tools to do this for you given some control points; if you want to implement it for yourself you'll have to start delving into computer vision.

Edit - OK, here you go: a Python script to do the warping, using the OpenCV library which has convenient functions to calculate the homography and warp the image for you:

import cv

def warpImage(image, corners, target):
    mat = cv.CreateMat(3, 3, cv.CV_32F)
    cv.GetPerspectiveTransform(corners, target, mat)
    out = cv.CreateMat(height, width, cv.CV_8UC3)
    cv.WarpPerspective(image, out, mat, cv.CV_INTER_CUBIC)
    return out

if __name__ == '__main__':
    width, height = 400, 250
    corners = [(171,72),(331,93),(333,188),(177,210)]
    target = [(0,0),(width,0),(width,height),(0,height)]
    image = cv.LoadImageM('fries.jpg')
    out = warpImage(image, corners, target)
    cv.SaveImage('fries_warped.jpg', out)

And the output:

OpenCV also has C and C++ bindings, or you can use EmguCV for a .NET wrapper; the API is fairly consistent across all languages so you can replicate this in whichever language suits your fancy.

Tympanitis

I fear I phrased my question too generally - I know the dimensions of the rectangular face I'm trying to extract, and the coordinates of it's corners as in the image, leaving me wondering only how to get best-quality pixels out - how to calculate the coverage of each pixel, and how to best mix them – Vizard

@Vizard - selecting four quadrilateral corner points in the image, and four corresponding rectangle points is enough information to recover the homography (i.e. transformation matrix), but you need to do some math; see this pdf: cs.brown.edu/courses/csci1950-g/asgn/proj6/resources/… – Tympanitis

Answer 3

And in modern times in python with cv2.

import cv2
import numpy as np

source_image = cv2.imread('french fries in Europe.jpeg')

source_corners = np.array([(171, 72), (331, 93), (333, 188), (177, 210)])

width, height = 400, 250
target_corners = np.array([(0, 0), (width, 0), (width, height), (0, height)])

# Get matrix H that maps source_corners to target_corners
H, _ = cv2.findHomography(source_corners, target_corners, params=None)

# Apply matrix H to source image.
transformed_image = cv2.warpPerspective(
    source_image, H, (source_image.shape[1], source_image.shape[0]))

cv2.imwrite('tranformed_image.jpg', transformed_image)

Cardie

cv2.warpPerspective(source_image, H, (width, height)) can be used to get the same output as in tzaman's answer. – Leonleona

Answer 4

I think what you need is affine transformation which can be accomplished using matrix math.

Foliose

You're right, I forgot about the perspective. An affine transformation would result in a distorted image. Sorry about that. – Foliose