Rotating greyscale images

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(ax + by + e, 3), round(cx + dy + 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

1. Comments on your code

-
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.

-
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.)

-
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.

-
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.

-
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.

1. 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 i

StackExchange Code Review Q#41688, answer score: 14