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! 🙂
xor 