Rotating greyscale images

Question

For educational purposes I wrote a little piece of code to rotate greyscale images as "low level" as possible, that is, not using any rotate() function, but doing the math. I was wondering if it could be improved in any way, specially in order to achieve better performance. I'm not concerned with the libraries or the programming language chosen to perform this task, I am aware that they aren't the best option if what I want is performance. I'm only concerned with the code itself.

It basically maps the pixels on the new image to the ones in the original image and interpolates using nearest neighbor or bilenear.

def affine_t(x, y, a, b, c, d, e, f):
    # A lot faster than numpy.dot([[a,b],[c,d]],[x,y])+[e,f]
    return round(a*x + b*y + e, 3), round(c*x + d*y + f, 3)

def mrotate(r, cx=0, cy=0):
    # Returns the coefficients to an affine transformation which rotates
    # a point clockwise by r radians in respect to a central point. For counter 
    # clockwise, just pass r as -r.
    return (math.cos(r), -math.sin(r), math.sin(r), math.cos(r), cx, cy)

def lin_interp(x, x0, x1, y0, y1):
    # Faster than numpy.interp()
    return y0 + (y1 - y0)*((x-x0)/(x1-x0))

def bilinear(bmp, ox, oy):
    # Try to interpolate. If the pixel falls on the image boundary, use the 
    # nearest pixel value (nearest neighbor).
    try:
        x0, x1, y0, y1 = int(ox), int(ox)+1, int(oy), int(oy)+1
        a = lin_interp(ox, x0, x1, bmp[y0][x0], bmp[y0][x1])
        b = lin_interp(ox, x0, x1, bmp[y1][x0], bmp[y1][x1])
    except IndexError:
        return nn(bmp, ox, oy)
    return int(lin_interp(oy, y0, y1, a, b))

def nn(bmp, ox, oy):
    return bmp[oy][ox] # bmp is numpy.array, it casts float indexes to int

def rotate(bmp, r, mx=0, my=0, filename=None, interpol=None):
    """Rotate bitmap bmp r radians clockwise from the center. Move it mx, my."""

    # Get the bitmap's original dimensions and calculate the new ones
    oh, ow = len(bmp), len(bmp[0])
    nwl = ow * math.cos(r)
    nwr = oh * math.sin(r)
    nhl = ow * math.sin(r)
    nhu = oh * math.cos(r)
    nw, nh = int(math.ceil(nwl + nwr)), int(math.ceil(nhl+nhu))
    cx, cy = ow/2.0 - 0.5, oh/2.0 - 0.5 # The center of the image
    # Some rotations yield pixels offscren. They will be mapped anyway, so if 
    # the user moves the image he gets what was offscreen. 
    xoffset, yoffset = int(math.ceil((ow-nw)/2.0)), int(math.ceil((oh-nh)/2.0))
    for x in xrange(xoffset,nw):
        for y in xrange(yoffset,nh):
            ox, oy = affine_t(x-cx, y-cy, *mrotate(-r, cx, cy))
            if ox > -1 and ox < ow and oy > -1 and oy < oh:
                pt = bilinear(bmp, ox, oy) if interpol else nn(bmp, ox, oy)
                draw.point([(x+mx,y+my),],fill=pt)
    if filename is not None:
        im.save(filename)

The complete code, used to generate the examples below:

http://pastebin.com/2j7BQyzi

Some example rotations, first image is the original, to the left is one rotated 90 degrees and the one below is rotated 45 degrees.

Using Nearest Neighbors:

Using Bilinear:

Note: The original "Lenna.png" image was first resized to 0.5 the original size using bilinear interpolation and flatten to a greyscale.

Answer 1

1. Comments on your code

1. For most of your functions, you've written a comment describing what it does. It's usual in Python to put this in a docstring, so that a user can get at it from the interactive interpreter using the help function.

2. The function rotate relies on a global variable draw. This makes it hard to reuse and test. (And you couldn't use it in a multi-threaded program.)

3. The function rotate carries out multiple tasks: it rotates the bitmap; translates it; draws the result to an Image object; and saves it to a file. It would make the code easier to use and more general if you separated these functions.

4. The comments suggest that you have not understood how to use NumPy. For example, you write that your function lin_interp is "Faster than numpy.interp()". Which no doubt it is for the case of a single point. But when properly vectorized, NumPy will beat pure Python.

   For a very simple example, consider a pure-Python implementation of the dot product operation:

   def dot(v, w):
       return sum(x * y for x, y in zip(v, w))
   

   Let's compare the speed with numpy.dot. The pure-Python version is substantially faster for very small arrays:

   >>> import numpy as np
   >>> from timeit import timeit
   >>> timeit(lambda:dot([1, 2], [3, 4]))
   2.1183970817364752
   >>> timeit(lambda:np.dot([1, 2], [3, 4]))
   11.150392828974873
   

   but the NumPy version is hundreds of times faster for large arrays:

   >>> v = np.arange(1000)
   >>> w = np.arange(1000)
   >>> timeit(lambda:dot(v, w), number=1000)
   1.4182810601778328
   >>> timeit(lambda:np.dot(v, w), number=1000)
   0.0053491611033678055
   

   So in order to get a performance benefit out of NumPy, you have to vectorize your algorithm so that each step of the algorithm operates on many elements simultaneously.

5. Here's another example. In this extract you iterate over a grid of coordinates (x, y) and rotate each coordinate by −r radians about the point (cx, cy):

   for x in xrange(xoffset,nw):
       for y in xrange(yoffset,nh):
           ox, oy = affine_t(x-cx, y-cy, *mrotate(-r, cx, cy))
   

   Let's put that in a function for timing purposes:

   def rotate1(a, theta, ox, oy):
       """Rotate array of 2-dimensional coordinates by theta radians about
       the point (ox, oy).
   
       """
       for x, y in a:
           affine_t(x-ox, y-oy, *mrotate(theta, ox, oy))
   

   To be as generous as possible with the timing results, I haven't even attempted to store or return a result. Let's try that out on an array of a million coordinates:

   >>> grid = np.meshgrid(np.arange(1000), np.arange(1000))
   >>> coords = np.vstack(grid).reshape(2, -1).T
   >>> timeit(lambda:rotate1(coords, np.pi / 4, 500, 500), number=1)
   94.44321689195931
   

   Here's how the same operation looks when vectorized:

   def rotate2(x, y, theta, ox, oy):
       """Rotate arrays of coordinates x and y by theta radians about the
       point (ox, oy).
   
       """
       s, c = np.sin(theta), np.cos(theta)
       x, y = x - ox, y - oy
       return x * c - y * s + ox, x * s + y * c + oy
   

   Instead of looping over the coordinates and applying all the stages of the transformation to one coordinate at a time, we apply each stage of the transformation to all the coordinates simultaneously.

   >>> timeit(lambda:rotate2(grid[0], grid[1], np.pi / 4, 500, 500), number=1)
   0.06842728098854423
   

   So although your comment says that affine_t is "A lot faster than numpy", in fact, when properly vectorized, NumPy runs more than 1300 times faster.

2. Revised code

If I was doing this for real, then obviously I'd use scipy.misc.imrotate or skimage.transform.rotate or whatever. But considered only as an exercise, here's how I'd program this (in the nearest neighbour case, which is quite enough for one answer).

import numpy as np

def rotate_coords(x, y, theta, ox, oy):
    """Rotate arrays of coordinates x and y by theta radians about the
    point (ox, oy).

    """
    s, c = np.sin(theta), np.cos(theta)
    x, y = np.asarray(x) - ox, np.asarray(y) - oy
    return x * c - y * s + ox, x * s + y * c + oy

def rotate_image(src, theta, ox, oy, fill=255):
    """Rotate the image src by theta radians about (ox, oy).
    Pixels in the result that don't correspond to pixels in src are
    replaced by the value fill.

    """
    # Images have origin at the top left, so negate the angle.
    theta = -theta

    # Dimensions of source image. Note that scipy.misc.imread loads
    # images in row-major order, so src.shape gives (height, width).
    sh, sw = src.shape

    # Rotated positions of the corners of the source image.
    cx, cy = rotate_coords([0, sw, sw, 0], [0, 0, sh, sh], theta, ox, oy)

    # Determine dimensions of destination image.
    dw, dh = (int(np.ceil(c.max() - c.min())) for c in (cx, cy))

    # Coordinates of pixels in destination image.
    dx, dy = np.meshgrid(np.arange(dw), np.arange(dh))

    # Corresponding coordinates in source image. Since we are
    # transforming dest-to-src here, the rotation is negated.
    sx, sy = rotate_coords(dx + cx.min(), dy + cy.min(), -theta, ox, oy)

    # Select nearest neighbour.
    sx, sy = sx.round().astype(int), sy.round().astype(int)

    # Mask for valid coordinates.
    mask = (0 <= sx) & (sx < sw) & (0 <= sy) & (sy < sh)

    # Create destination image.
    dest = np.empty(shape=(dh, dw), dtype=src.dtype)

    # Copy valid coordinates from source image.
    dest[dy[mask], dx[mask]] = src[sy[mask], sx[mask]]

    # Fill invalid coordinates.
    dest[dy[~mask], dx[~mask]] = fill

    return dest

Here's fabio.jpg, a 200×200 image (courtesy of Oasisisgood at en.wikipedia, licensed under CC-BY-SA):

and here are some rotated versions:

>>> from scipy.misc import imread, imsave
>>> img = imread('fabio.jpg')
>>> for angle in (45, 75, 90, 120):
...     imsave('fabio-{}.jpg'.format(angle),
...            rotate_image(img, angle * np.pi / 180, 100, 100))

The rotation of this 200×200 image takes about ¹/₆₀ of a second:

>>> timeit(lambda:rotate_image(img, np.pi / 4, 100, 100), number=1)
0.01648401003330946

(scipy.misc.imrotate is about 4 times faster than this when using interp='nearest'.)