Rotating greyscale images
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:

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


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

&gt;&gt;&gt; v = np.arange(1000)
&gt;&gt;&gt; w = np.arange(1000)
&gt;&gt;&gt; timeit(lambda:dot(v, w), number=1000)
1.4182810601778328
&gt;&gt;&gt; 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:

&gt;&gt;&gt; grid = np.meshgrid(np.arange(1000), np.arange(1000))
&gt;&gt;&gt; coords = np.vstack(grid).reshape(2, -1).T
&gt;&gt;&gt; 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.

&gt;&gt;&gt; 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 &lt;= sx) &amp; (sx &lt; sw) &amp; (0 &lt;= sy) &amp; (sy &lt; 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:

&gt;&gt;&gt; from scipy.misc import imread, imsave
&gt;&gt;&gt; img = imread('fabio.jpg')
&gt;&gt;&gt; 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 1/60 of a second:

&gt;&gt;&gt; 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'.)