14
\$\begingroup\$

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:

enter image description here

Using Bilinear:

enter image description here

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

\$\endgroup\$
12
+200
\$\begingroup\$

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

Fabio Lanzoni

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

Fabio rotated by 45 degrees Fabio rotated by 75 degrees Fabio rotated by 90 degrees Fabio rotated by 120 degrees

The rotation of this 200×200 image takes about 1/60 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'.)

\$\endgroup\$
  • \$\begingroup\$ Perfect. I did find it odd that numpy was taking longer than my function. I had never done any image manipulation so the only thing I could think of was to do it point by point. As for the comments, global variables, etc I didn't bother since it was just an exercise, I know it's messy. Thanks! (I can't award the bounty yet because of the min time limit, I'll do it as soon as possible). \$\endgroup\$ – Alex Feb 17 '14 at 21:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.