Journey through a two-tone corridor

This article brings a Python code, with and without numpy, to generate a video of a camera moving down an infinite two-tone rectangular corridor, with a checkerboard pattern (with quite large panels).

When I started to code for real, what I did was “demos”. Those small programs that produce cool animations. I never did nothing too fancy, just basic Computer Graphics stuff like translucent rotating cubes. I started with Pascal (yuck!) and then moved to C… Today I use mostly Python, and I work with Computer Vision what is more or less the opposite of Computer Graphics. While CG is concerned with producing images from a given model, CV tackles the inverse problem of figuring out what model would produce that image. Or images, in the case of analyzing a video.

It is useful in CV research to have simulations made with CG so we can test our programs with a known model, to check our answers and debug the systems starting way back at the image processing steps. I have just finished creating a program, and found myself needing some artificial images to make a test. I needed a video of a camera moving along a corridor while its orientation and position changes slightly in time. And the corridor needed to have edges on the directions of the coordinate axes (x, y, z, i.e. vertical, horizontal and “along the corridor” directions).

So I settled to make that in Python… This is simple because it requires nothing fancy such as complicated textures, illumination models, ray tracing… The only complicated thing is to understand all the geometry stuff, know about camera models, perspective, rotation matrices, etc. Here is one frame produced by the program. Of course, you can also watch the final video, on YouTube.

So, without further ado, here is the code. It is commented, but if anyone has any questions, feel free to ask!

#!/usr/bin/python
#####################################################################
## Code to draw an animation of a camera moving down a
## checkerboad-pattern corridor.
## 
## By Nicolau Werneck <nwerneck@gmail.com>
## https://xor0110.wordpress.com
## http://nwerneck.sdf.org
## 

from pylab import *
import numpy

from itertools import product


#####################################################################
## Produces a rotation matrix from the 3 last components of a
## quaternion.
def quaternion_to_matrix(myx):
    xb,xc,xd = myx
    xnormsq = xb*xb+xc*xc+xd*xd
    assert xnormsq < 1
    b,c,d = xb,xc,xd
    a = numpy.sqrt(1-xnormsq)

    assert a >= -1
    assert a <= 1

    return numpy.array([ [(a*a+b*b-c*c-d*d), (2*b*c-2*a*d),     (2*b*d+2*a*c)      ],
                         [(2*b*c+2*a*d),     (a*a-b*b+c*c-d*d), (2*c*d-2*a*b)      ],
                         [(2*b*d-2*a*c),     (2*c*d+2*a*b),     (a*a-b*b-c*c+d*d)] ]  )  \
                         / (a*a+b*b+c*c+d*d)
##
####################################################################

sig = 1./(640.-1.0)
intrinsic = array([ [ sig,    0, 0.],
                    [   0,  sig, 0.],
                    [-0.5, -0.37480438184663539, 1.]  ])

## Simple initial and sluggish implementation.
def draw_stuff(frame, pos, rotM ):
    ## pos is the position of the camera, and rotM the rotation matrix
    ## that models its orientation.

    ## This loop iterates over each pixel of the image
    for j,k in product(range(frame.shape[0]), range(frame.shape[1])):
        ## Converts from image coordinates to camera coordinates using
        ## the "intrinsic parameters"
        vimg = dot(array([k,j,1]), intrinsic)
        ## Converts to world coordinates, using the "extrinsic parameters"
        vspc = dot(vimg, rotM)

        ## Test if ray points to left or right, up or down planes, and
        ## calculates the coordinates of the point where the ray
        ## crosses the plane.
        if vspc[0]>0:
            ptx = pos + vspc * (1-pos[0])/vspc[0]
        else:
            ptx = pos + vspc * (-1-pos[0])/vspc[0]
        if vspc[1]>0:
            pty = pos + vspc * (1-pos[1])/vspc[1]
        else:
            pty = pos + vspc * (-1-pos[1])/vspc[1]
        
        ## Now test if the ray is crossing a vertical or horizontal
        ## plane, and selects the value accordingly.
        if ptx[2]<pty[2]:
            frame[j,k]=int(ptx[2])%2
        else:
            frame[j,k]=(1+int(pty[2]))%2


## Now a faster numpy implementation. First of all, we define this
## list with a vector for each pixel.
thejks = array([(k,j,1) for j,k in product(range(480), range(640))])
def draw_stuff2(pos, rotM ):
    ##Calculates ray direction in camera reference frame
    vimg = dot(thejks, intrinsic)
    ##Calculates ray direction in world reference frame
    vspc = dot(vimg, rotM)

    ## This value is added in the following expresison according to
    ## the ray being pointing to the left or right plane
    incx= (vspc[:,0]>0)*3-1.5
    ## Calculates point where ray crosses the plane
    ptx = pos[2] + vspc[:,2] * (incx-pos[0])/vspc[:,0]
    ## Intensity according to z coordinate of the point
    valx= ((ptx/2)%2) > 1

    incy= (vspc[:,1]>0)*2-1
    pty = pos[2] + vspc[:,2] * (incy-pos[1])/vspc[:,1]
    ## All the same thing here, except we "flip" the intensity on the
    ## horizontal planes
    valy= ((pty/2)%2) < 1

    ## What are the pixels that came from the vertical planes
    sel = ptx<pty
    ## valy will now hold all the pixel values
    valy[ sel ] = valx[ sel ]

    ## Grays out the pixels that are too far away (not working!!!)
    sel = (numpy.min(c_[ptx,pty],1)-pos[2]) > 100.
    valy[sel]=0.1

    ## Now just reshape valy for the image array
    return reshape(valy, (480, 640) )


## No rotation
# rotM = quaternion_to_matrix([0.,0.,0.])

## Loop to generate each frame
for ind in range(500):
    print ind

    ## Calculates new position. z coordinate is proportional to image
    ## index on sequence, creating a forward-moving track. On the
    ## other coordinates we use sinusoidal functions with different
    ## frequencies
    pos = array([0.5*cos(ind*0.05),-0.2*sin(ind*0.15),ind*0.2])

    # Calculates caemra orientation. Torns mostly around x axis, some
    # around y and a little around z (spinning aroudnthe motion
    # direction)
    rotM = quaternion_to_matrix([-0.03 * sin(ind*0.1) ,0.01*cos(ind*0.2),0.008*sin(ind*0.04)])

    ## Calculate simage and stores on frame
    frame = draw_stuff2(pos, rotM)

    ## Saves frame into numbered png image
    imsave('qq%04d.png'%ind, frame, cmap='gray')

One last note, if you are curious about the quaternions thing, I recommend the Wikipedia article Quaternions and spatial rotation

Advertisements

One thought on “Journey through a two-tone corridor

  1. Pingback: Cythonic CG corridor trip « xor

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s