Cythonic CG corridor trip

This is a sequel to my previous post, Journey through a two-tone corridor. I just enhanced the code a bit, implementing the main procedure in Cython. I also implemented an anti-aliasing technique that is only used when necessary.

I will first just spill out the code, to leave as reference to anyone trying to do something similar. First, as usual when working with Cython, it is useful to have a small setup script to help you compile the code and generate the python module you will sue later. Here is what I used… The file xadrez_aux-setup.py is simply

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_modules = [Extension("xadrez_aux", ["xadrez_aux.pyx"], include_dirs=['/usr/share/pyshared/numpy/core/include/','/usr/local/lib/python2.6/site-packages/numpy/core/include/numpy/'], extra_objects=[],libraries=['m'])]

setup(
  name = 'xadrez aux functions',
  cmdclass = {'build_ext': build_ext},
  ext_modules = ext_modules
)

The module code goes in the file xadrez_aux.pyx, and to compile it you just have to run

python2.6 xadrez_aux-setup.py build_ext --inplace --force

This uses C to compile the xadrez_aux.so module that you can use from inside Python. The code at xadrez_aux.pyx is:

from itertools import product
import numpy as np
cimport numpy as np
cimport cython


@cython.boundscheck(False)
@cython.wraparound(False)
def draw_stuff(width,height,
               np.ndarray[np.int_t, ndim=2, mode="c"] frame,
               np.ndarray[np.float32_t, ndim=3, mode="c"] frameaux,
               np.ndarray[np.float32_t, ndim=1, mode="c"] pos,
               np.ndarray[np.float32_t, ndim=2, mode="c"] fullmat):
    cdef float level
    cdef int j, k
    cdef float jj, kk
    cdef float incx, incy, ptx,pty

    cdef float * fm = <float*>fullmat.data

    cdef float * pt = <float*>frameaux.data

    cdef float vspc[3]

    cdef int valx=0, valy=0


    cdef float amin, amax
    cdef float a1,a2,a3,a4
    cdef int v1,v2,v3,v4


    ## First loop calcualtes the frameaux array. it stores the ptx and
    ## pty values of the z coordinates of the point where the camera
    ## ray crosses the plane.
    for j,k in product(range(0,height+2,1), range(0,width+2,1)):
        ##Calculates ray direction in world reference frame
        vspc[0] = (k-1.0)*fm[0] + (j-1.0)*fm[3] + fm[6]
        vspc[1] = (k-1.0)*fm[1] + (j-1.0)*fm[4] + fm[7]
        vspc[2] = (k-1.0)*fm[2] + (j-1.0)*fm[5] + fm[8]

        ## Test againt zero division, and calculates the z coordinate
        ## of the point where the ray cast form the camera crosses the
        ## left or right plane, depending on the direction of the x or
        ## y component.
        if vspc[0] != 0:
            incx= 4.0 if (vspc[0]>0) else -4.0
            ptx = pos[2] + vspc[2] * (incx-pos[0])/vspc[0]

        if vspc[1] != 0:
            incy= 1.5 if (vspc[1]>0) else -15
            # incy= 4.0 if (vspc[1]>0) else -4.0
            pty = pos[2] + vspc[2] * (incy-pos[1])/vspc[1]

        # Store component for vertical and horizontal planes.
        frameaux[j,k,0] = ptx
        frameaux[j,k,1] = pty


    ## The next loop calculates final pixel value according to
    ## frameaux. First it tests the neighboring pixels, and check if
    ## they belong to the same plane. If they do, there is a second
    ## check to see if they have the same color, and if the space
    ## range of each pixel is not too large. Everythin OK, the color
    ## is calculated frmo the z coordinate, with color expression
    ## depending if it's a verticla or horizontal plane. If not we use
    ## the aliasing procedure.
    for j,k in product(range(0,height,1), range(0,width,1)):
        if frameaux[j+1,k+1,0] < frameaux[j  ,k+1,1] and \
           frameaux[j+1,k+1,0] < frameaux[j+2,k+1,1] and \
           frameaux[j+1,k+1,0] < frameaux[j+1,k  ,1] and \
           frameaux[j+1,k+1,0] < frameaux[j+1,k+2,1] :
            
            a1,a2,a3,a4 = frameaux[j+1,k,0],frameaux[j+1,k+2,0],frameaux[j,k+1,0],frameaux[j+2,k+1,0]

            amin =   a1 if a1<a2 else a2
            amin = amin if amin<a3 else a3
            amin = amin if amin<a4 else a4
            amax =   a1 if a1>a2 else a2
            amax = amax if amax>a3 else a3
            amax = amax if amax>a4 else a4

            v1=(a1%4)>2
            v2=(a2%4)>2
            v3=(a3%4)>2
            v4=(a4%4)>2


            if v1==v2 and v1==v3 and v1==v4 and (amax-amin)<2.0:
                frame[j,k] = 255*v1
                continue

        if frameaux[j+1,k+1,1] < frameaux[j  ,k+1,0] and \
           frameaux[j+1,k+1,1] < frameaux[j+2,k+1,0] and \
           frameaux[j+1,k+1,1] < frameaux[j+1,k  ,0] and \
           frameaux[j+1,k+1,1] < frameaux[j+1,k+2,0] :
            
            ## Is an horizontal plane. So
            a1,a2,a3,a4 = frameaux[j+1,k,1],frameaux[j+1,k+2,1],frameaux[j,k+1,1],frameaux[j+2,k+1,1]

            amin =   a1 if a1<a2 else a2
            amin = amin if amin<a3 else a3
            amin = amin if amin<a4 else a4
            amax =   a1 if a1>a2 else a2
            amax = amax if amax>a3 else a3
            amax = amax if amax>a4 else a4

            v1=(a1%4)<=2
            v2=(a2%4)<=2
            v3=(a3%4)<=2
            v4=(a4%4)<=2


            if v1==v2 and v1==v3 and v1==v4 and (amax-amin)<2.0:
                frame[j,k] = 255*v1
                continue

        level=0
        for jj,kk in product( (np.arange(16)-7.5)/16., (np.arange(16)-7.5)/16.):
            vspc[0] = (kk+float(k))*fm[0] + (jj+float(j))*fm[3] + fm[6]
            vspc[1] = (kk+float(k))*fm[1] + (jj+float(j))*fm[4] + fm[7]
            vspc[2] = (kk+float(k))*fm[2] + (jj+float(j))*fm[5] + fm[8]

            if vspc[0] != 0:
                incx= 4.0 if (vspc[0]>0) else -4.0
                ptx = pos[2] + vspc[2] * (incx-pos[0])/vspc[0]
            else:
                ptx=0

            if vspc[1] != 0:
                incy= 1.5 if (vspc[1]>0) else -15
                # incy= 4.0 if (vspc[1]>0) else -4.0
                pty = pos[2] + vspc[2] * (incy-pos[1])/vspc[1]
            else:
                pty=0
            

            if ptx<pty:
                level+= 1 if ((ptx%4)>2) else 0
            else:
                level+= 1 if ((pty%4)<=2) else 0

        frame[j,k]= round( 255 * (level/256.) )

And finally, there is a scrit called geraxader.py that just sets things up and repeatedly calls this Cython procedure to cenerate each frame from the film. It is where we calculate the camera position at each time.

#!/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
import sys
from itertools import product

from xadrez_aux import draw_stuff

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


###########################
## Image formats
## 1280 x 720
## 640 x 480
###########################



# height,width = 720,
height,width = 480,640
# height,width = 905,640


sig = 1./(width-1.0)
intrinsic = array([ [ sig,    0, 0.],
                    [   0,  sig, 0.],
                    [-0.5, -0.5*sig*(height-1.0), 1.]  ], dtype=numpy.float32)



print "alocando imagem"
frame = zeros( (height,width), dtype=numpy.int  )
frameaux = zeros( (height+2,width+2,2), dtype=numpy.float32  )
print "imagem alocada"
## Loop to generate each frame
for ind in range(int(sys.argv[1]), 100+int(sys.argv[1])):
    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.6*cos(ind*0.05),-0.3*sin(ind*0.15),ind*0.2], dtype=numpy.float32)
    # pos = array([0.2+1.6*cos(pi*ind*0.05),0.1-0.5*sin(pi*ind*0.05), ind*0.4 +0.02*cos(ind*0.10)], dtype=numpy.float32)
    # pos=array([-0.1,-0.2,0.1],dtype=numpy.float32)
    # pos=array([0.,0.,0.],dtype=numpy.float32)
    pos=array([1.5+1.5*cos(pi*ind*0.042),  0.6*cos(pi*ind*0.031) ,  ind*0.3],dtype=numpy.float32)
    
    # 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.04 * sin(ind*0.1) ,0.01*cos(ind*0.25),0.008*sin(ind*0.06)])
    rotM = quaternion_to_matrix([-0.03, 0.03*sin(pi*ind*0.042)  ,0.02*sin(pi*ind*0.03)])

    # dd=sin(  0.5*pi*((ind*0.0025)%1)   )
    # dd = 0
    # dd=sin(  0.5*pi*0.25   )
    # rotM = quaternion_to_matrix([0,0,dd])

    rotM = dot(intrinsic, rotM)
    ## Calculate simage and stores on frame
    # frame = draw_stuff2(pos, rotM)
    draw_stuff(width,height,frame, frameaux, pos, rotM.astype(numpy.float32) ) 

    ## Saves frame into numebred png image
    imsave('/home/nlw/ciencia/DADOS/xadrez/qw%04d.png'%ind, frame, cmap='gray')

This is a resulting frame. This time I am simulating a very tall corridor. You can check the full video in YouTube. The came moves in a sinusoidal track, with little oscillations to all directions.

The code can probably still be enhanced in may ways… Using SSE instructions, for example. But that is good enough for me. Good luck if you are pursuing a similar project! 🙂

Advertisements

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