solem's vision blog

Polar panoramas with Python

2012-02-11T22:59:00.000-08:00

Today I saw (again) Filipe Calegario's nice photos and instructions for making polar panoramas. To show how to warp images with Python I wrote a Python version of the polar warp. Here's a complete code example that uses scipy.ndimage (except for the final hole filling/interpolation).


import Image
from numpy import *

from scipy.ndimage import geometric_transform
from scipy.misc import imsave

def polar_warp(im,sz,r):
 """ Warp an image to a square polar version. """
 
 h,w = im.shape[:2]
  
 def polar_2_rect(xy):
  """ Convert polar coordinates to coordinates in 
   the original image for geometric_transform. """
   
  x,y = float(xy[1]),float(xy[0])
  theta = 0.5*arctan2(x-sz/2,y-sz/2)
  R = sqrt( (x-sz/2)*(x-sz/2) + (y-sz/2)*(y-sz/2) ) - r
  
  return (2*h*R/(sz-2*r),w/2+theta*w/pi+pi/2)
 
 # if color image, warp each channel, otherwise there's just one
 if len(im.shape)==3: 
  warped = zeros((sz,sz,3))
  for i in range(3):
   warped[:,:,i] = geometric_transform(im[:,:,i],polar_2_rect,output_shape=(sz,sz),mode='nearest')
 else:
  warped = geometric_transform(im,polar_2_rect,output_shape=(sz,sz),mode='nearest')
 
 return warped


if __name__ == "__main__":
 
 # load image from http://www.filipecalegario.com/comppho/3assignment/
 im = array(Image.open('autostitch-stitch.jpg'))
 h,w = im.shape[:2]
 
 # flip vertical to make math easier for warp
 im = flipud(im)
 
 # warp image
 wim = uint8(polar_warp(im,h,5))
 
 # save
 imsave('polar.png',wim)

This reads a regular 360 panorama (in this example the "autostitch-stitch.jpg" file), flips the image and warps so that the sky faces out. Save this as a file and run from the command line and you should get an image like the one above.

Here's the file if you don't want to copy/paste the above.

Installing the OpenCV Python Interface on Lion

2011-12-20T16:51:00.000-08:00

On some machines I ran into "Segmentation fault 11" when I use the OpenCV Python interface on Mac OS X 10.7 Lion. These problem machines were all upgraded from Snow Leopard which had Python 2.6 (Lion has 2.7). Anyway, the simplest solution for me was to do a clean Lion install and then install OpenCV. Here's the whole process, as a big fat "note to self" and hopefully a help for others.

1. Make a clean install of Lion by copying the install .dmg to a DVD or USB stick. See process here.

2. Install Xcode through the Mac App Store.

3. Get the Scipy/Matplotlib/iPython/... stack from the Scipy Superpack.

4. Install Cmake. Download the .dmg from here and run.

4. Download OpenCV, the latest version is 2.3.1.

5. Unpack, open Terminal and go to the folder OpenCV-2.3.1. Run the following commands:


mkdir build
cd build
cmake -G "Unix Makefiles" ..
make -j8
sudo make install

6. Add the build/lib folder that contains cv2.so to your PYTHONPATH.

I have tried this on a number of computers and it works great for me. If you are having OpenCV-Python problems and want to use MacPorts or other variants, the clean Lion install might help as well.

Good luck.

Book draft

2011-08-23T09:19:00.000-07:00

The other night I put up the first early draft of my upcoming book "Programming Computer Vision with Python".

The book is meant as an entry point to hands-on computer vision for students, researchers and enthusiasts using python. I'm thinking introductory courses in image analysis and computer vision. Also OpenCV python uses who need to do something outside what's in OpenCV will find the book useful. There is a growing python community around computer vision with projects like pythonvision and SimpleCV. The book can serve as a computer vision introduction for these users as well. And let's not forget all the hackers and enthusiasts out there.

It is not meant as a textbook, it is a manual for getting started with python and/or computer vision. I'm putting everything online in order to get early feedback and help hammering out errors and improve the text and code.

Draft versions together with code and data are available from the book webpage.

I'll take all the feedback and comments I can get! Preferably via email.

UPDATE (Sept 25): Thanks so much for all the feedback and comments! I'm now on the third iteration of the book pdf since putting it online a month ago. Check the book page for updates and keep sending me your thoughts. /JES

Another Python Interface for SIFT

2011-06-08T20:46:00.000-07:00

I previously wrote about using David Lowe's SIFT code from Python here. There is a good open source alternative for those on Windows or Mac OS (Lowe's binaries are Linux only) called VLFeat. The library is written in c but has a command line interface that we can use much in the same way.

To install VLFeat, download and unpack the latest binary package from the download page (currently the latest version is 0.9.9). Add the paths to your environment or copy the binaries to a directory in your path. The binaries are in the bin/ directory, just pick the sub-directory for your platform. The use of the VLFeat command line binaries is described in the src/ sub-directory. Alternatively you can find the documentation online.

I wrote a similar Python interface as the one I had for Lowe's code. You can download it here (vlfeat.py).

Try it out like this:


from PIL import Image
from pylab import *

process_image('box.pgm','tmp.sift')
l,d = read_features_from_file('tmp.sift')

im = array(Image.open('box.pgm'))
figure()
plot_features(im,l,True)

process_image('scene.pgm','tmp2.sift')
l2,d2 = read_features_from_file('tmp2.sift')
im2 = array(Image.open('scene.pgm')) 

m = match_twosided(d,d2)
figure()
plot_matches(im,im2,l,l2,m)

show()

The result looks like this.

For those who want more control, there is also the option of using the Python wrapper I wrote about before.

Read/Write World

2011-05-31T05:12:00.000-07:00

Last week at ARE2011 Blaise Aguera y Arcas and Avi Bar-Zeev gave two talks that got my attention. The Read/Write World initiative from Bing. This is interesting from a number of angles:

RML - Reality Markup Language

RML is a markup language for describing the shape, location, and content of the world (and images of it). All the details are not posted yet but supposedly supports transforms between views, images, video, panoramas, geo etc. Worth checking in to see where this goes in the near future.

Open Source Viewers

The demo apps show some interesting ways of viewing content. The project will contain HTML5 based viewers available under open source licenses. Again, worth checking in on later.

Services

Hosted by Bing, services will include indexing, matching between images (e.g. homographies), 20GB of storage and possibly other image and video services.

All in all it looks ambitious and very interesting. Not much to see right now but definitely a project to keep an eye on and check out in a few months time.

Adjacency matrix for image pixel graphs

2011-04-08T18:52:00.000-07:00

One of my favorite NumPy functions is roll(). Here's an example of using this function to get neighborhood indices to create adjacency matrices for images.

An adjacency matrix is a way of representing a graph and shows which nodes of the graph are adjacent to which other nodes. For n nodes is is an n*n matrix A where a_ij is the number of edges from vertex i to vertex j. The number of edges can also be replaced with edge weights depending on the application.

For images, a pixel neighborhood defines a graph which is sparse since each pixel is only connected to its neighbors (usually 4 or 8). A sparse representation, for example using dictionaries, is preferable if only the edges are needed. For clustering using spectral graph theory, a full matrix is needed. The following function creates an adjacency matrix in sparse or full matrix form given an image:


from numpy import *

def adjacency_matrix(im,connectivity=4,sparse=False):
 """ Create a pixel connectivity matrix with
  4 or 8 neighborhood. If sparse then a dict is
  returned, otherwise a full array. """
  
 m,n = im.shape[:2]

 # center pixel index
 c = arange(m*n)
 ndx = c.reshape(m,n)

 if sparse:
  A = {}
 else:
  A = zeros((m*n,m*n))

 if connectivity==4: # 4 neighborhood
  hor = roll(ndx,1,axis=1).flatten() # horizontal
  ver = roll(ndx,1,axis=0).flatten() # vertical
  for i in range(m*n):
   A[c[i],hor[i]] = A[hor[i],c[i]] = 1
   A[c[i],ver[i]] = A[ver[i],c[i]] = 1

 elif connectivity==8: # 8 neighborhood
  hor = roll(ndx,1,axis=1).flatten() # horizontal
  ver = roll(ndx,1,axis=0).flatten() # vertical
  diag1 = roll(roll(ndx,1,axis=0),1,axis=1).flatten() # diagonal
  diag2 = roll(roll(ndx,1,axis=0),1,axis=0).flatten() # diagonal
  for i in range(m*n):
   A[c[i],hor[i]] = A[hor[i],c[i]] = 1
   A[c[i],ver[i]] = A[ver[i],c[i]] = 1
   A[c[i],diag1[i]] = A[diag1[i],c[i]] = 1
   A[c[i],diag2[i]] = A[diag2[i],c[i]] = 1

 return A

Use it like this:


from pylab import *

im = random.random((10,10))
A = adjacency_matrix(im,8,False)

figure()
imshow(1-A)
gray()
show()

Which should give you an image like the one below.

RQ Factorization of Camera Matrices

2011-03-22T22:58:00.000-07:00

A common operation on camera matrices is to use RQ factorization to obtain the camera calibration matrix given a 3*4 projection matrix. The simple example below shows how to do this using the scipy.linalg module. Assuming the camera is modeled as P = K [R | t], the goal is to recover K and R by factoring the first 3*3 part of P.

The scipy.linalg module actually contains RQ factorization although this is not always clear from the documentation (here is a page that shows it though). To use this version, import rq like this:


from scipy.linalg import rq

Alternatively, you can use the more common QR factorization and with some modifications write your own RQ function.


from scipy.linalg import qr

def rq(A): 
 Q,R = qr(flipud(A).T)
 R = flipud(R.T)
 Q = Q.T 
 return R[:,::-1],Q[::-1,:]

RQ factorization is not unique. The sign of the diagonal elements can vary. In computer vision we need them to be positive to correspond to focal length and other positive parameters. To get a consistent result with positive diagonal you can apply a transform that changes the sign. Try this on a camera matrix like this:


# factor first 3*3 part of P
K,R = rq(P[:,:3])

# make diagonal of K positive
T = diag(sign(diag(K)))

K = dot(K,T)
R = dot(T,R) #T is its own inverse

Multidimensional meshgrid with Python

2011-03-09T22:50:00.000-08:00

Looking for an easy way to plot cubes in 3D I found NumPy's mgrid that can do meshgrid in any number of dimensions. The syntax is a bit confusing but once you understand the trick, it works.

Here's a simple example that creates points for a unit cube and plots the points.


from pylab import *
from numpy import *
from mpl_toolkits.mplot3d import axes3d

# create 3D points
x,y,z = mgrid[0:2,0:2,0:2]
xx = x.flatten()
yy = y.flatten()
zz = z.flatten()

# plot 3D points
fig = figure()
ax = fig.gca(projection='3d')
ax.plot(xx,yy,zz,'o')

show()

3D plotting with Matplotlib

2011-02-08T08:59:00.000-08:00

Earlier, Matplotlib (and Python) lacked proper support for 3D plotting. This is no longer the case as the mplot3d toolkit provides 3D plotting of points, lines, contours, surfaces and all other basic components as well as 3D rotation and scaling.
Making a plot 3D is done by adding the projection=‘3d’ keyword to the axes object like this:


from pylab import *
from numpy import *
from mpl_toolkits.mplot3d import axes3d

fig = figure()
ax = fig.gca(projection='3d')

# plot points in 3D
class1 = 0.6 * random.standard_normal((200,3))
ax.plot(class1[:,0],class1[:,1],class1[:,2],'o')
class2 = 1.2 * random.standard_normal((200,3)) + array([5,4,0])
ax.plot(class2[:,0],class2[:,1],class2[:,2],'o')
class3 = 0.3 * random.standard_normal((200,3)) + array([0,3,2])
ax.plot(class3[:,0],class3[:,1],class3[:,2],'o')

3D scatter plots of three distributions.

Plotting surfaces is just as easy. The following example uses get_test_data() to create sample data.


# generate 3D sample data
X,Y,Z = axes3d.get_test_data(0.05)

fig = figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y,Z,alpha=0.5)

Surface plot with transparency.

There are many options for formatting the appearance of the plots, including color and transparency of surfaces. The mplot3d website has the details.

Installing LibSVM

2011-01-27T12:41:00.000-08:00

LibSVM is a good Support Vector Machine (SVM) implementation which works well with Python. The current release of LibSVM is version 3.0 (released September 2010). Download the zip file from the LibSVM website. Unzip the file (a directory "libsvm-3.0" will be created). In a terminal window go to this directory and type "make".


$ cd libsvm-3.0
$ make

Then go to the "python" directory and do the same:


$ cd python/
$ make

This should be all you need to do. To test your installation, start python from the command line and try


>>> import svm

(You might have to change the line libsvm = CDLL('../libsvm.so.2') in svm.py to contain the location of libsvm if you have problems importing svm from other directories on your computer.)

The authors of LibSVM wrote a practical guide for using LivSVM. This is a good starting point if you are new to SVM.

Converting video to NumPy arrays

2011-01-10T22:55:00.000-08:00

Using OpenCV it is possible to read video frames from a file and convert them to NumPy arrays. Frames returned by QueryFrame() are in OpenCV's IplImage format where the pixels are stored as arrays of color in BGR order. Here's an example of converting the 100 first frames of a video to standard RGB and storing them in a NumPy array.


from numpy import *
import cv

capture = cv.CaptureFromFile('skyline.mov')
frames = []
for i in range(100):
    img = cv.QueryFrame(capture)
    tmp = cv.CreateImage(cv.GetSize(img),8,3)
    cv.CvtColor(img,tmp,cv.CV_BGR2RGB)
    frames.append(asarray(cv.GetMat(tmp))) 
frames = array(frames)

Here the CvtColor() function converts from BGR to RGB and NumPy's asarray() does the final conversion. The array will have dimensions (100,width,height,3) in this case.

There is no good way to determine the end of the file, QueryFrame() will just loop around to the beginning. If you want to you can add a check to break when the first frame appears a second time.

Here's the video I used:

OpenCV with Python on Mac OS X

2010-12-14T17:39:00.000-08:00

Open source computer vision library OpenCV is gaining a lot of traction and now comes with what looks like a pretty good Python interface. There are installers for Windows and Linux but Mac OS X support has been lacking. There are several ways to install from source as described here. MacPorts is one option that I usually like but since I'm not using the MacPort version of the Python/Numpy/Scipy/Matplotlib stack I decided to try to build using cmake. I followed the instructions to the point:


svn co https://code.ros.org/svn/opencv/trunk/opencv 
cd opencv
sudo cmake -G "Unix Makefiles" .
sudo make -j8
sudo make install

Now everything builds and installs properly. When testing the python interface I got:


>>> import cv
Traceback (most recent call last):
  File "", line 1, in 
ImportError: No module named cv

This was solved by adding the directory containing cv.so to my PYTHONPATH. In my case:


export PYTHONPATH=/usr/local/lib/python2.6/site-packages/

Then it all works fine. I can start python and run the cookbook examples without any problems. That's it. Hope these instructions can be of use.

Check out the online OpenCV Python reference guide for more examples and details.

Transparent figures

2010-12-07T18:52:00.000-08:00

Matplotlib (Pylab) has an interesting option in the savefig() command; the figure background can be set to be transparent. Very useful if you want to put your figures on any non-white background (e.g. presentation slides). You can also set degrees of transparency for the different figure elements using an alpha value. See the documentation for more details.

Here's an example:


from pylab import *

# plot sin(x) for some interval
x = linspace(-2*pi,2*pi,200)
plot(x,sin(x))

# plot marker for every 4th point
plot(x[::4] ,sin(x[::4] ),'r*')
title('Function sin(x) and some points plotted')

savefig('test.pdf', transparent=True)

And here's the figure on a sample presentation slide, pretty nice.

Plotting Geo Data using Basemap

2010-09-01T07:28:00.000-07:00

Matplotlib has an interesting toolkit called Basemap for plotting geo data and maps. Basemap handles all the possible map projections and uses Matplotlib to do the plotting. Basemap also includes data like countries, coastlines and background satellite images from NASA Blue Marble.

Plotting geo data using Basemap is actually very simple. Here is an example with geo data from Twitter using the script in my previous post.


from mpl_toolkits.basemap import Basemap
from pylab import *
import twitter_geo

# setup Lambert Conformal basemap.
m = Basemap(width=400000,height=400000,projection='lcc',resolution='f',lat_0=55.5,lon_0=13)
m.bluemarble() #background
m.drawcoastlines() #coastlines

# get some geo data from Twitter
geo = '55.583333,13.033333,100km' #Malmo, Sweden
tweets = twitter_geo.parse(twitter_geo.search(geo=geo))

for t in tweets:
    coords = t[2]['coordinates']
    x,y = m(coords[1],coords[0]) #project to map coords
    plot(x,y,'ro')
    
show() #you might not need this

This example takes tweets in a 100km radius around Malmö, Sweden and plots the location with red dots on top of a NASA Blue Marble image with coastlines drawn. The result looks like this:

Pretty simple. If you had some dense data like weather or topographic data, that could be placed on the map (e.g. using Matplotlib's contour()) and opens up interesting possibilities for data visualization.

Geo Data from Twitter using Python

2010-08-30T07:46:00.000-07:00

I wanted to write an example of plotting geo data with maps but I didn't have any good example data. Having played around with Twitter's API before, I decided to use that. There are several Python libraries for interfacing with Twitter but to use the Twitter Search API you actually don't need any of them. Since the result is returned in JSON you just need something like simplejson (or the built in JSON modules in later versions of Python).

Here is a simple example (twitter_geo.py) that searches for tweets with a query string and/or a geographic location. Pretty simple.


import urllib2
import simplejson as json

def search(q='',geo=''):
    """ use search api to find tweets matching a query string and/or 
        location as string of "lat,lon,radius", see api documentation"""
    query = 'http://search.twitter.com/search?q=%s&format=json&rpp=100&result_type=recent&geocode=%s' % (q,geo)
    f = urllib2.urlopen(query)
    r = json.loads(f.read())
    return r["results"]

def parse(res):
    """ take the relevant parts of the result"""
    #list of tuples (user,time,geo)
    return [(r['from_user'],r['created_at'],r['geo']) 
            for r in res if r['geo'] != None]


if __name__ == "__main__":
    tag = '%23fail' 
    geo = '55.583333,13.033333,100km' #Malmo, Sweden
    
    #query Search API
    print parse(search(q=tag,geo=geo))

Installing Numpy and Matplotlib on Mac OS X 10.6 Snow Leopard

2010-08-22T11:33:00.000-07:00

This is partly an update to my previous post on installing for Leopard, partly a note to myself. After upgrading to Snow Leopard (and 64bit) I could still use my MacPort installs but not upgrade or install new ones (since all previous builds were for 32bit). So, time for another fresh install.

With MacPorts:


sudo port install python26
sudo port install py26-numpy
sudo port install py26-matplotlib

After this I can start python and load pylab but I get problems with the plotting window totally missing. My colleague Jerome Piovano found a fix over at stackoverflow. Replace the last line with


sudo port install py26-matplotlib +cairo+gtk2

Then edit ~/.matplotlib/matplotlibrc (it should be in your home directory, if no such file exists just make a new file) by adding:


backend: MacOSX

Now everything works as it should.

Computational Mathematics with Python

2010-06-04T13:24:00.000-07:00

Posting has been a little slow the last months. I have been busy at Polar Rose and in parallel I've been teaching a course in computational mathematics with python. I hope to increase posting soon as I have a couple of unfinished posts and examples I'd like to write about.

If you are interested in numerical analysis or how to program mathematics using python, check out the course page. Lecture slides and course material is available here.

Computer Vision and the New iPhone SDK

2010-04-19T13:28:00.000-07:00

A look at the newly released iPhone SDK for iPhone OS 4 shows some very interesting new features for people building computer vision and augmented reality applications.

Video playback & Capture
"You now have full programmatic control over video playback and capture, using new APIs in the AV Foundation framework." This is the one AR people have been waiting for. No more hacks to get the video feed. Should give more video based apps in the open and better video performance.

Photo Library Access
"Applications now have direct access to user photos and videos with the Media Library APIs." Also very useful. Build auto-tagging tools or photo browsers for your on-phone photos! (Cooliris, are you watching this? You finally have access.)

And finally...

Accelerate
"Accelerate provides hundreds of mathematical functions optimized for iPhone and iPod touch, including signal-processing routines, fast Fourier transforms, basic vector and matrix operations, and industry-standard functions for factoring matrices and solving systems of linear equations." Will need to try it to find out how well it works but this is potentially extremely useful. Some basic optimized linear algebra will take computer vision or AR far on this device! Is it using OpenCL? I don't know. Do you?

sparray - sparse n-dimensional arrays in Python

2010-02-09T14:09:00.000-08:00

Yesterday I wrote about sparse arrays for Python using dictionaries. I couldn't get it out of my head so tonight I actually made a complete module containing a class for sparse n-dimensional arrays. The module supports any number of dimensions and any size. Sparray has all basic operations like adding, subtracting, multiplication, division, mod, pow, printing, sum() and can easily be converted to full NumPy arrays.

The sparray module can be downloaded from here. Running sparray.py will show a series of examples and print the results in the console.

Enjoy. Bug reports and feedback welcome!

Below is some output from the example:


adding...
[[  3.   3.   3.]
 [  3.  13.   3.]
 [  3.   3.   3.]]
subtracting...
[[-3. -3. -3.]
 [-3.  7. -3.]
 [-3. -3. -3.]]
multiplication...
[[  0.   0.   0.]
 [  0.  30.   0.]
 [  0.   0.   0.]]
[[ 9.  9.  9.]
 [ 9.  9.  9.]
 [ 9.  9.  9.]]
division...
[[ 0.  0.  0.]
 [ 0.  3.  0.]
 [ 0.  0.  0.]]
mod...
[[ 0.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  0.]]
power...
[[    0.     0.     0.]
 [    0.  1000.     0.]
 [    0.     0.     0.]]
iadd...
[[  3.   3.   3.]
 [  3.  13.   3.]
 [  3.   3.   3.]]
sum of elements...
74
mix with NumPy arrays...
[[  6.   6.   6.]
 [  6.  26.   6.]
 [  6.   6.   6.]]
Frobenius norm...
601.0

Sparse arrays in Python using dictionaries

2010-02-08T15:35:00.000-08:00

UPDATE: I made a complete module for sparse n-dimensional array. See the next post.

A great application of Python's dictionary type is to use it for representing sparse arrays and matrices (in arbitrary dimensions). Consider the following code example:


a = {}
a[0,3] = 12
a[2,1] = 5
a[3,3] = 7

Now, this dictionary looks like this:


{(0, 3): 12, (3, 3): 7, (2, 1): 5}

It is easy to define operations on these dictionaries. For example, the Euclidean (L2) distance between two arrays (Frobenius norm in the case of matrices) is simply:


def L2(a1,a2):
    ndx = set(a1.keys()+a2.keys())
    return sum([ (a1.get(i,0)-a2.get(i,0))**2 for i in ndx])

Here we used set() to remove duplicates after concatenating the keys and the dictionary get() method that returns a default value (in this case zero) if a key is missing.

With another example array like this:


b = {}
b[0,3] = 10
b[1,1] = 15
b[1,3] = 4

the value of L2(a,b) is 319.

There are some interesting efforts for sparse multidimensional arrays in Python out there. Ndsparse is a new one, which uses lists instead of dictionaries and looks very nice and clean. Unfortunately I found lots of basic things broken when testing today. Hope to see a new release with bugs fixed soon! There is also the SciPy sparse module with data structures for 2D sparse matrices but that looks a little messy (and is only 2D).

CVXOPT - convex optimization with Python

2010-01-21T01:04:00.000-08:00

The other day I was looking around for a good optimization package for Python and found something that looks great.

CVXOPT is a free software package (GPL license) for convex optimization based on Python and developed by Joachim Dahl and Lieven Vandenberghe. CVXOPT is NumPy compatible which means that you can use your arrays as usual.

The CVXOPT website contains lots of examples from the book Convex Optimization by Boyd and Vandenberghe. For example this total variation reconstruction which is quite neat.

After playing around with it and trying some examples I really like it and it looks solid. Give it a try next time you have a linear program or something in need of optimization.

Looking for the best of the best

2010-01-01T10:49:00.000-08:00

(image CC by TinyFizzyPop)

If you are interested (or know anyone who might be) in pursuing a PhD in computer vision there are now several open positions at the Mathematical Imaging Group at Lund University. PhD advisers will be either myself or the amazing Fredrik Kahl.

We are looking for strong applicants interested in one or several of the following topics: optimization, medical imaging, image retrieval, recognition and machine learning. Personally my thoughts and ideas revolve a lot around large scale retrieval and graph based methods these days. Experience with applied mathematics, computer vision, optimization, graph algorithms is a plus!

If you are interested, send me an email or a tweet. The official posting is here.

Deadline for applying is January 25, 2010.

VLFeat - now also with Python interfaces

2009-10-21T15:00:00.000-07:00

VLFeat is an open source library that implements some popular computer vision algorithms including SIFT, MSER, k-means, hierarchical k-means and more. It is written in C (by Andrea Vedaldi and Brian Fulkerson) and comes with interfaces in Matlab. The license is GNU GPL v2.

If you are using local descriptors like SIFT or MSER or if you need fast clustering code you should definitely try out VLFeat. The easiest is to download the version with both source and binaries from here. There is a Git repository with the most current development code.

The good news is that there is now a fairly complete Python wrapper, open sourced and developed by my Polar Rose colleague Mikael Rousson. The wrapper is currently in a Git fork, get it here. Try it out! (You need Boost:Python installed).

Is your NumPy using the right ATLAS?

2009-10-08T11:10:00.000-07:00

Here's a tip I got from my colleague Martin. Making sure that your NumPy installation uses the right back-end can drastically improve performance. ATLAS is the software library used by NumPy (and Matlab, Octave etc.) to do linear algebra.

Try running the following 1000*1000 matrix multiplication test:


>>> from numpy import *
>>> import time
>>> A = random.random((1000,1000))
>>> B = random.random((1000,1000))
>>> t = time.time(); dot(A,B); print time.time()-t

Disappointed about the time this takes? Take a minute to check that you are running the right ATLAS for your computer. On my Mac this test was fast (0.2-0.3s) but on my desktop running Ubuntu very slow (5-6s). In Ubuntu, open the Package Manager and search for "atlas". If you find something like "libatlas3gf-base" you are using a generic version that works on most hardware. Unless your computer is antique, you can probably do better. Try selecting "libatlas3gf-sse2" (or even libatlas3gf-sse1 if this doesn't work) instead and run the test again. (on other linux platforms the process should be similar)

Tests I did on a couple of machines show speedup of a factor 4-10. Well worth the minute it takes to check and fix. Thanks Martin!

Exporting from NumPy to Weka

2009-09-16T12:19:00.000-07:00

Weka is a very useful tool for experimenting with different classifiers, pre-processing of data and much more. Weka is an open source (GPL) collection of machine learning algorithms for data mining written in Java. Not only algorithms, Weka also has a GUI for exploring data, running classification experiments and visualization without having to write any code. Some screen shots:

Loading and exploring some data (in this case the iris sample data set).

Running cross-correlation with different classifiers.

Visualization.

If you are used to working in NumPy, the following function will write NumPy arrays to Weka .arff data files that can be loaded directly from the GUI.


def write_to_weka(filename,relationname,attributenames,attributes,comment=''):
    """ writes NumPy arrays with data to WEKA format .arff files
    
        input: relationname (string with a description), attributenames (list 
        of the names of different attributes), attributes (array of attributes, 
        one row for each attribute, WEKA treats last row as classlabels by 
        default), comment (short description of the content)."""
    
    nbrattributes = len(attributenames)
    if attributes.shape[1] != nbrattributes:
        raise Exception('Number of attribute names is not equal to length of attributes')
    
    f = open(filename, 'w')
    f.write('% '+comment+'\n')
    f.write('\n')
    f.write('@RELATION '+relationname+'\n')
    
    for a in attributenames:
        f.write('@ATTRIBUTE '+str(a)+' NUMERIC\n') #assume values are numeric
    
    f.write('\n')    
    f.write('@DATA\n') #write the data, one attribute vector per line
    for i in range(attributes.shape[0]):
        for j in range(nbrattributes):
            f.write(str(attributes[i,j]))
            if j < nbrattributes-1:
                f.write(', ')
        f.write('\n') 
    f.close()

Weka is one of those tools I always forget I have, I hope this post will help me remember and hopefully introduce Weka to new people.

Tip: If you run out of memory, start Weka with "java -Xmx512m weka.jar weka.gui.GUIChooser" and you will have 512MB... or whatever you choose.

"Cover Flow" with PyGame and OpenGL

2009-08-31T12:34:00.000-07:00

The other day I was curious at PyOpenGL and PyGame and decided to try to implement a crude "cover flow"-like image browser in Python. PyGame is a set of Python modules designed for writing games. PyOpenGL is a free Python binding to OpenGL. The NEHE OpenGL demos (python version available at here) and the examples found at this page contains lots of useful examples to learn from.

The result is here. Very crude and unpolished but actually something that works. The demo loads all .jpg images in the current folder and shows them with "periodic boundary conditions". (there are some limitations: it doesn't work with too big images but you could just add a check and resize, images are assumed to have 2:3 aspect ratio, all images are loaded in memory + many more I'm sure)

Anyway, it was a fun exercise.

Report from CIVR 2009

2009-08-20T00:10:00.000-07:00

I'm slowly getting into posting after a stint of traveling and vacation. As a warmup, I'll give a short report (way overdue!) on CIVR 2009 where I gave a talk early July. The conference was held in Thira on the island of Santorini, Greece.

The conference twitter hashtag gives a nice record of what people were tweeting about at the time. The conference program is here.

One of the highlights was Luis von Ahn's keynote "Human Computation". The keynote was more or less the same as his talk at Library of Congress mixed with some from his Google tech talk at the end. Both presentations are well worth seeing if you are not familiar with Luis' work.

Some interesting papers:
* Jasper Uijlings, Arnold Smeulders and Remko Scha, Real-time Bag of Words, Approximately (Speedups for bag of word pipeline, fast computation of histograms using matrix multiplication. Use RBF kernel, random forest etc, giver higher accuracy and much faster.) Awarded best paper of CIVR 2009.
* Rainer Lienhart, Stefan Romberg and Eva Hörster, Multilayer pLSA for Multimodal Image Retrieval (Combine visual words and tags a la Flickr)
* Ville Viitaniemi and Jorma Laaksonen. Spatial Extensions to Bag of Visual Words ("soft" assignment)
* Matthijs Douze, Hervé Jégou, Harsimrat Singh, Laurent Amsaleg and Cordelia Schmid. Evaluation of GIST descriptors for web-scale image search (use of GIST as more compact representation to scale to 100M images)
In general there is a clear trend towards using Flickr tags combined with SIFT descriptors.

My presentation at the "practitioner day" was about Polar Rose and lessons learned when deploying large scale computer vision for the web. Slides are here.

A nice conference, smaller than the ICCV/ECCV/CVPRs I try to go to but focused and full of interesting presentation. The next conference will be in Xi'an, China, July 5-7, 2010.

Picture of the speakers at the practitioner day outside the conference building with the volcano in the back.

Robust homography estimation using RANSAC

2009-07-15T00:09:00.000-07:00

I'm just back from lovely Santorini and the ACM CIVR conference. State-of-the-art image retrieval seems to use a bag of feature representation (a topic itself worthy of a few future posts) with re-ranking of the top results based on geometric consistency of local image features. This consistency check is done using homography estimation. As I mentioned in my previous post about RANSAC, this is a useful algorithm for among other things homography estimation. Here's a simple way of using the Python RANSAC implementation to do exactly this.

We can use ransac.py for any model. All that is needed is a class with fit() and get_error() methods.

To fit a homography using RANSAC we first need the following model class (you can add it to the file homography.py from another post which already contains some useful things).


class ransac_model:
    """ class for testing homography fit with ransac.py from
        http://www.scipy.org/Cookbook/RANSAC"""

    def __init__(self,debug=False):
        self.debug = debug

    def fit(self, data):
        """ fit homography to four selected correspondences"""

        #transpose to fit H_from_points()
        data = data.T

        #from points
        fp = data[:3,:4]
        #target points
        tp = data[3:,:4]

        #fit homography, H
        H = H_from_points(fp,tp)

        return H

    def get_error( self, data, H):
        """ apply homography to all correspondences, 
            return error for each transformed point"""

        data = data.T

        #from points
        fp = data[:3]
        #target points
        tp = data[3:]

        #transform fp
        fp_transformed = dot(H,fp)

        #normalize hom. coordinates
        for i in range(3):
            fp_transformed[i] /= fp_transformed[2]

        err_per_point = sqrt( sum((tp-fp_transformed)**2,axis=0) )

        return err_per_point
        
        
def H_from_points(fp,tp):
    """ find homography H, such that fp is mapped to tp
        using the linear DLT method. Points are conditioned
        automatically."""

    if fp.shape != tp.shape:
        raise RuntimeError, 'number of points do not match'


    #condition points (important for numerical reasons)
    #--from points--
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1))
    C1 = diag([1/maxstd, 1/maxstd, 1]) 
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp = dot(C1,fp)

    #--to points--
    m = mean(tp[:2], axis=1)
    #C2 = C1.copy() #must use same scaling for both point sets
    maxstd = max(std(tp[:2], axis=1))
    C2 = diag([1/maxstd, 1/maxstd, 1])
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp = dot(C2,tp)


    #create matrix for linear method, 2 rows for each correspondence pair
    nbr_correspondences = fp.shape[1]
    A = zeros((2*nbr_correspondences,9))
    for i in range(nbr_correspondences):        
        A[2*i] = [-fp[0][i],-fp[1][i],-1,0,0,0,tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]
        A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]

    U,S,V = linalg.svd(A)

    H = V[8].reshape((3,3))    

    #decondition
    H = dot(linalg.inv(C2),dot(H,C1))

    #normalize and return
    return H / H[2][2]

As you can see, this class contains a fit() method which just takes the four correspondences selected by ransac.py (they are the first four in data) and fits a homography. Remember, four points are the minimal number to compute a homography. The method get_error() applies the homography and returns the distance for each correspondence pair so that RANSAC can chose which points to keep as inliers and outliers. This is done with a threshold on this distance.

We also need the function H_from_points() to estimate the homography in each RANSAC iteration. Here we use the DLT algorithm.

For ease of use, add the following function to homography.py.


def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):
    """ robust estimation of homography H from point 
        correspondences using RANSAC (ransac.py from
        http://www.scipy.org/Cookbook/RANSAC).

        input: fp,tp (3*n arrays) points in hom. coordinates"""

    import ransac

    #use ransac class
    model = ransac_model()
    #group corresponding points
    data = vstack((fp,tp))

    H = ransac.ransac(data.T,model,4,maxiter,match_theshold,10)

    return H

The function also lets you supply the threshold and the minimum number of points desired. The most important parameter is the maximum number of iterations, exiting too early might give a worse solution, too many iterations will take more time.

Try this on some sample data to see that it works. Using a set with 20 inliers and 10 outliers you can get a figure like this:

Pretty neat!

CVPR 2009 Word Cloud Summary

2009-07-06T12:19:00.000-07:00

I seem to have missed an exciting CVPR this year. For those of you who also missed out, I have created a word cloud (using Wordle) from the full program page of CVPR 2009. Here it is:

I removed "University", "Time" and "Location" for obvious reasons. Give Wordle a try, it is a nice tool.

RANSAC using Python

2009-06-30T12:21:00.000-07:00

RANSAC, short for "RANdom SAmple Consensus", is an iterative method to fit models to data that can contain outliers. Given a model, e.g. a homography between points, the basic idea is that the data contains inliers, the data points that can be described by the model, and outliers, those that do not fit the model.

The standard example is the case of fitting a line to a set of points that contains outliers. Simple least squares fitting will fail but RANSAC can hopefully single out the inliers and obtain the correct fit. Let's look at using ransac.py which contains this particular example as test case. The image below shows an example of running ransac.test(). As you can see, the algorithm selects only points consistent with a line model and correctly finds the right solution.

RANSAC is a very useful algorithm and is frequently used on computer vision to make homography estimation and structure from motion robust to noise and false image correspondences.

Piecewise affine warping using SciPy

2009-06-27T14:31:00.000-07:00

As we saw in the previous post, three point correspondences uniquely defines an affine transformation and therefore affine warping of triangles can be done exactly. Let's look at the most common form of warping between a set of corresponding points, piecewise affine warping. Given any image with landmark points we can warp that image to corresponding landmarks in another image by triangulating the points in a mesh and warping each triangle with an affine transform. This is standard stuff for any graphics and image processing library. It is also possible to do it using SciPy and matplotlib.

To triangulate points, Delaunay triangulation is often used. This triangulation comes included in Matplotlib (but outside the PyLab part) and can be used like this.


import matplotlib.delaunay as md 
from pylab import *
from numpy import *

x,y =  array(random.standard_normal((2,100)))
centers,edges,tri,neighbors = md.delaunay(x,y)

figure()
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] #add first point to end
    plot(x[t_ext],y[t_ext],'r')

plot(x,y,'*')
axis('off')
show()

The image below shows some example points and the resulting triangulation. Delaunay triangulation chooses the triangles so that the edges are the dual graph of a Voronoi diagram. There are four outputs of delaunay() of which we only need the list of triangles. Create a function in warp.py for the triangulation.


import matplotlib.delaunay as md
    
def triangulate_points(x,y):
    """delaunay triangulation of 2D points"""

    centers,edges,tri,neighbors = md.delaunay(x,y)
    return tri

The output is an array with each row containing the indices for the three points of each triangle.

Let's look at an example of warping an image to a non-flat object in another image using 30 control points in a 5 by 6 grid. The example below shows an image to be warped to the facade of the "turning torso". The target points were manually selected using ginput().

First we need a general warp function for piecewise affine image warping. The code below does the trick, where we also take the opportunity to learn how to warp color images (you simply warp each color channel).


def pw_affine(fromim,toim,fp,tp,tri):
    """ warp triangular patches from an image.
        fromim = image to warp 
        toim = destination image
        fp = from points in hom. coordinates
        tp = to points in hom.  coordinates
        tri = triangulation"""

    im = toim.copy()

    #check if image is grayscale or color
    is_color = len(fromim.shape) == 3

    #create image to warp to (needed if iterate colors)
    im_t = zeros(im.shape, 'uint8') 

    for t in tri:
        #compute affine transformation
        H = homography.Haffine_from_points(tp[:,t],fp[:,t])

        if is_color:
            for col in range(3):
                im_t[:,:,col] = ndimage.affine_transform(fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        else:
            im_t = ndimage.affine_transform(fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])

        #alpha for triangle
        alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
        #add triangle to image
        if is_color:
            for col in range(3):
                im[:,:,col] = (1-alpha)*im[:,:,col] + alpha*im_t[:,:,col]
        else:
            im = (1-alpha)*im + alpha*im_t

    return im

Add this function to warp.py (here homography.Haffine_from_points() refers to the function in the previous post). You also need:


def alpha_for_triangle(points,m,n):
    """ creates alpha map of size (m,n) 
        for a triangle with corners defined by points
        (given in normalized homogeneous coordinates)."""

    alpha = zeros((m,n))
    for i in range(min(points[0]),max(points[0])):
        for j in range(min(points[1]),max(points[1])):
            x = linalg.solve(points,[i,j,1])
            if min(x) > 0: #all coefficients positive
                alpha[i,j] = 1
                
    return alpha

To use this on the current example, just do


from PIL import Image
from pylab import *
from numpy import *
import homography
import warp

#open image to warp
fromim = array(Image.open('jes.jpg'))
x,y = meshgrid(range(5),range(6))
x = (fromim.shape[1]/4) * x.flatten()
y = (fromim.shape[0]/5) * y.flatten()

#triangulate
tri = warp.triangulate_points(x,y)

#open image and destination points
im = array(Image.open('turningtorso1.jpg'))
tp = loadtxt('turningtorso1_points.txt') #destination points

#convert points to hom. coordinates
fp = vstack((y,x,ones((1,len(x))))) 
tp = vstack((tp[:,1],tp[:,0],ones((1,len(tp)))))

#warp triangles
im = warp.pw_affine(fromim,im,fp,tp,tri)

The triangles can be plotted, as in the simple example above, with the following code (add this to warp.py).


def plot_mesh(x,y,tri):
    """ plot triangles""" 

    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]] #add first point to end
        plot(x[t_ext],y[t_ext],'r')

Example of placing an image on the "turning torso" using a 5 by 6 grid of landmark points. Image credits: rutgerblom from Flickr (original, license).

If you want to try this yourself, here is the target image and control points (you can use any source image).

Affine transformations and warping

2009-06-16T00:33:00.000-07:00

This is the first of two posts on warping with Python using SciPy/NumPy. In the last couple of days warping have come up on several occasions when talking to my colleagues at work. This inspired me to write some example code over the weekend. Hope this can be of use for someone. The next post will be about piecewise affine warping. Here goes.

An affine transformation is a linear transformation followed by a translation. For details see e.g. Wikipedia. Affine transformations are a good approximation of image to image mappings if the transformed region is small or if the image is captured with a large focal length.

Affine transforms can be estimated using the following algorithm, described in detail in "Multiple View Geometry" Hartley & Zisserman.

Add the following function to a file homography.py, which computes the affine transformation matrix from point correspondences.


def Haffine_from_points(fp,tp):
    """ find H, affine transformation, such that 
        tp is affine transf of fp"""

    if fp.shape != tp.shape:
        raise RuntimeError, 'number of points do not match'

    #condition points
    #-from points-
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1))
    C1 = diag([1/maxstd, 1/maxstd, 1]) 
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp_cond = dot(C1,fp)

    #-to points-
    m = mean(tp[:2], axis=1)
    C2 = C1.copy() #must use same scaling for both point sets
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp_cond = dot(C2,tp)

    #conditioned points have mean zero, so translation is zero
    A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
    U,S,V = linalg.svd(A.T)

    #create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.
    tmp = V[:2].T
    B = tmp[:2]
    C = tmp[2:4]

    tmp2 = concatenate((dot(C,linalg.pinv(B)),zeros((2,1))), axis=1) 
    H = vstack((tmp2,[0,0,1]))

    #decondition
    H = dot(linalg.inv(C2),dot(H,C1))

    return H / H[2][2]

Applying an affine transformation matrix H on image patches is called warping (or affine warping) and is frequently used in computer graphics but also in several computer vision algorithms. A warp can easily be performed with SciPy using the ndimage package. The command


transformed_im = ndimage.affine_transform(im,A,b,size)

transforms the image patch im with A a linear transformation and b a translation vector as above. The optional argument size can be used to specify the size of the output image. The default is an image with the same size as the original. To see how this works, try running the following commands:


>>>from PIL import Image
>>>from scipy import ndimage
>>>from pylab import *

>>>im = array(Image.open('empire.jpg').convert('L'))
>>>H = array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])
>>>im2 = ndimage.affine_transform(im,H[:2,:2],(H[0,2],H[1,2]))

>>>figure()
>>>gray()
>>>imshow(im2)
>>>show()

This gives you a result like the image below. As you can see, missing pixel values in the result image are filled with zeros.

An example application of affine warping is to place images, or parts of images, inside another image.

Add the function image_in_image() to a file warp.py. This function takes two images and the corner coordinates of where to put the first image in the second.


def image_in_image(im1,im2,tp):
    """ put im1 in im2 with an affine transformation
        such that corners are as close to tp as possible.
        tp are homogeneous and counter-clockwise from top left.""" 

    #points to warp from
    m,n = im1.shape[:2]
    fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])

    #compute affine transform and apply
    H = Haffine_from_points(tp,fp)
    im1_t = ndimage.affine_transform(im1,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])
    alpha = (im1_t > 0)

    return (1-alpha)*im2 + alpha*im1_t

As you can see, there is not much needed to do this. When blending together the warped image and the second image we create an alpha map which defines how much of each pixel to take from each image. Here we use the fact that the warped image is filled with zeros outside the borders of the warped area to create a binary alpha map. (This is not the proper way to do it but works for our purposes here.)

To try this, let's insert an image on a billboard in another image.


from numpy import *
from PIL import Image
from pylab import *
import warp

#example of affine warp of im1 onto im2
im1 = array(Image.open('beatles.jpg').convert('L'))
im2 = array(Image.open('billboard_for_rent.jpg').convert('L'))

#set from points to corners of im
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])

#set to points
tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])

im3 = warp.image_in_image(im1,im2,tp)

figure()
gray()
imshow(im3)
axis('equal')
axis('off')
show()

This puts the image on the upper part of the billboard. Note that the landmark coordinates (tp and fp) are in homogeneous coordinates. Changing the coordinates to


tp = array([[675,826,826,677],[55,52,281,277],[1,1,1,1]])

will put the image on the lower left "for rent" part.

Image credits (CC): billboard_for_rent.jpg by striatic @ Flickr, http://flickr.com/photos/striatic/21671910/, beatles.jpg by oddsock @ Flickr, http://flickr.com/photos/oddsock/82535061/

The function Haffine_from_points() gives the best affine transform for the given point correspondences, in the example above between the image corners and the corners of the billboard. If the perspective effects are small, this will give good results.

That's all for now. More to come soon.

Histogram equalization with Python and NumPy

2009-06-07T11:54:00.000-07:00

I have been looking for a nice and clean implementation of histogram equalization with no luck so far. So I decided to write a short script and share it. You will need to have NumPy installed.

A very useful example of a graylevel transform is histogram equalization. This transform flattens the graylevel histogram of an image so that all intensities are as equally common as possible. This is often a good way to normalize image intensity before further processing and also a way to increase image contrast.

The transform function is in this case a cumulative distribution function (cdf) of the pixel values in the image (normalized to map the range of pixel values to itself).

Here's how to do it.


def histeq(im,nbr_bins=256):

   #get image histogram
   imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
   cdf = imhist.cumsum() #cumulative distribution function
   cdf = 255 * cdf / cdf[-1] #normalize

   #use linear interpolation of cdf to find new pixel values
   im2 = interp(im.flatten(),bins[:-1],cdf)

   return im2.reshape(im.shape), cdf

The function takes a grayscale image and the number of bins to use in the histogram as input and returns an image with equalized histogram together with the cumulative distribution function used to do the mapping of pixel values. To try this on an image,


from PIL import Image
from numpy import *

im = array(Image.open('AquaTermi_lowcontrast.jpg').convert('L'))
im2,cdf = imtools.histeq(im)

The picture below shows an example of histogram equalization (click to enlarge). The top row shows the graylevel histogram before and after equalization together with the cdf mapping. As you can see, the contrast increases and the details of the dark regions appear clearly.

SciPy tip: Saving arrays as images

2009-05-27T14:57:00.000-07:00

A quick tip on a useful function in the scipy.misc module. When manipulating images and doing computations using array objects it is useful to be able to save them directly as image ﬁles. The imsave() function is available through the scipy.misc module. To save an array im to ﬁle just do:


import scipy.misc
scipy.misc.imsave(’test.jpg’,im)

There is also the option of converting the array to a PIL Image object and save it from there but it is tricky and if you don't have the exact right format settings when doing the conversion your image will be scrambled. I definitely prefer the SciPy option.

More datasets - ImageNet

2009-05-13T12:05:00.000-07:00

From Princeton University comes a new extremely exciting dataset called ImageNet. It is a large hierarchical image database with images labeled according to a hierarchical structure provided by WordNet. Right now there are 3,247,902 images but this is expected to grow close to 50 million as more WordNet categories are added.

The dataset is explained in more detail in this paper (to appear in CVPR2009) where it is also compared to many alternative datasets like Caltech, PASCAL, ESP dataset, LabelMe. An interesting detail is that the images are verified by humans using Amazon's Mechanical Turk. These days there is certainly no lack of data for those interested in large-scale object recognition!

First impressions of Google Labs "similar images"

2009-05-07T02:40:00.000-07:00

If you haven't seen it already, the "similar images" feature from Google Labs is pretty cool. The idea is to refine a search with visual similarity. Feedback seems to be generally good and I have to agree. For things that have a distinct visual appearance, the refinement works really well.

I tried a search for "malmö" and then refined with images of turning torso and the öresund bridge, both give nice refinements. Something harder like "fixie bike" also gives good results.

Not all images have the "similar images" link which I take as not all images have been indexed (yet). Overall I like what I see and hope this makes it out of the Lab.

Reading and writing .mat files with Python

2009-05-05T02:22:00.000-07:00

If you have some old data, or find some interesting data set online, stored in Matlab's .mat file format it is possible to read this using the scipy.io module in SciPy. This is how to do it.


data = scipy.io.loadmat('test.mat')

"data" now contains a dictionary with keys corresponding to the variable names saved in the original .mat file. The variables are in array format. Saving to .mat files is equally simple. Just create a dictionary with all variables you want to save and use savemat().


data = {}
data['x'] = x
scipy.io.savemat('test.mat',data)

This saves the array x so that it has the name "x" when read into Matlab. More information on scipy.io can be found in the online documentation.

Hierarchical Clustering in Python

2009-04-30T02:20:00.000-07:00

Continuing on the topic of clustering algorithms, here is another popular choice with a Python implementation. No need to copy paste, the all the code can be found here.

Hierarchical clustering is another simple but powerful clustering algorithm. The idea is to build a similarity tree based on pairwise distances. The algorithm starts with grouping the two closest objects (based on the distance between feature vectors) and creates an "average" node in a tree with the two objects as children. Then the next closest pair is found among the remaining objets but also including any average nodes, and so on. At each node the distance between the two children is also stored. Clusters can then be extracted by traversing this tree and stopping nodes with distance smaller some threshold.

Hierarchical clustering has several benefits. For example the tree structure can be used to visualize relationships and how clusters are related. A good feature vector will give a nice separation in the tree. Another benefit is that the tree can be used with different cluster thresholds without recomputing the tree. The drawback, however, is that one needs to choose a threshold if the actual clusters are needed.

Let's see what this looks like in code. Create a file hcluster.py and add the following code (partially taken and modified from the example in "Programming Collective Intelligence" by Toby Segaran (O'Reilly Media 2007, page 33).


from numpy import *

class cluster_node:
    def __init__(self,vec,left=None,right=None,distance=0.0,id=None,count=1):
        self.left=left
        self.right=right
        self.vec=vec
        self.id=id
        self.distance=distance
        self.count=count #only used for weighted average 

def L2dist(v1,v2):
    return sqrt(sum((v1-v2)**2))

def hcluster(features,distance=L2dist):
    #cluster the rows of the "features" matrix
    distances={}
    currentclustid=-1

    # clusters are initially just the individual rows
    clust=[cluster_node(array(features[i]),id=i) for i in range(len(features))]

    while len(clust)>1:
        lowestpair=(0,1)
        closest=distance(clust[0].vec,clust[1].vec)

        # loop through every pair looking for the smallest distance
        for i in range(len(clust)):
            for j in range(i+1,len(clust)):
                # distances is the cache of distance calculations
                if (clust[i].id,clust[j].id) not in distances: 
                    distances[(clust[i].id,clust[j].id)]=distance(clust[i].vec,clust[j].vec)

                d=distances[(clust[i].id,clust[j].id)]

                if d < closest:
                    closest=d
                    lowestpair=(i,j)

        # calculate the average of the two clusters
        mergevec=[(clust[lowestpair[0]].vec[i]+clust[lowestpair[1]].vec[i])/2.0 \
            for i in range(len(clust[0].vec))]

        # create the new cluster
        newcluster=cluster_node(array(mergevec),left=clust[lowestpair[0]],
                             right=clust[lowestpair[1]],
                             distance=closest,id=currentclustid)

        # cluster ids that weren't in the original set are negative
        currentclustid-=1
        del clust[lowestpair[1]]
        del clust[lowestpair[0]]
        clust.append(newcluster)

    return clust[0]

Running hcluster() on a matrix with feature vectors as rows will create and return the cluster tree. The choice of distance measure depends on the actual feature vectors, here we used the Euclidean distance but you can create any function and use that as parameter.

To extract the clusters from the tree you need to traverse the tree from the top until a node with distance value smaller than some threshold is found. This is easiest done recursively.


def extract_clusters(clust,dist):
    # extract list of sub-tree clusters from hcluster tree with distance < dist
    clusters = {}
    if clust.distance < dist:
        # we have found a cluster subtree
        return [clust]
    else:
        # check the right and left branches
        cl = []
        cr = []
        if clust.left!=None: 
            cl = extract_clusters(clust.left,dist=dist)
        if clust.right!=None: 
            cr = extract_clusters(clust.right,dist=dist)
        return cl+cr

This function will return a list of sub-trees containing the clusters. To get the leaf nodes that contain the object ids, traverse each sub-tree and return a list of leaves using:


def get_cluster_elements(clust):
    # return ids for elements in a cluster sub-tree
    if clust.id>0:
        # positive id means that this is a leaf
        return [clust.id]
    else:
        # check the right and left branches
        cl = []
        cr = []
        if clust.left!=None: 
            cl = get_cluster_elements(clust.left)
        if clust.right!=None: 
            cr = get_cluster_elements(clust.right)
        return cl+cr

Let's try this on a simple example to see it all in action. The file sunsets.zip contains 100 images downloaded from Flickr using the tag "sunset" or "sunsets". For this example we will use the average color of each image as feature vector. This is crude and simple but good enough for illustrating what hierarchical clustering does. Try running the following code in a folder containing the sunset images.


import os
from PIL import Image
from numpy import *
import hcluster


#create a list of images
imlist = []
for filename in os.listdir('./'):
    if os.path.splitext(filename)[1] == '.jpg':
        imlist.append(filename)
n = len(imlist)

#extract feature vector for each image
features = zeros((n,3))
for i in range(n):
    im = array(Image.open(imlist[i]))
    R = mean(im[:,:,0].flatten())
    G = mean(im[:,:,1].flatten())
    B = mean(im[:,:,2].flatten())
    features[i] = array([R,G,B])

tree = hcluster.hcluster(features)

To visualize the cluster tree, one can draw a dendrogram. This often gives useful information on how good a given descriptor vector is and what is considered similar in a particular case. Add the following code (also adapted from Segaran's book).


from PIL import Image,ImageDraw

def drawdendrogram(clust,imlist,jpeg='clusters.jpg'):
    # height and width
    h=getheight(clust)*20
    w=1200
    depth=getdepth(clust)

    # width is fixed, so scale distances accordingly
    scaling=float(w-150)/depth

    # Create a new image with a white background
    img=Image.new('RGB',(w,h),(255,255,255))
    draw=ImageDraw.Draw(img)

    draw.line((0,h/2,10,h/2),fill=(255,0,0))    

    # Draw the first node
    drawnode(draw,clust,10,(h/2),scaling,imlist,img)
    img.save(jpeg)

def drawnode(draw,clust,x,y,scaling,imlist,img):
    if clust.id<0:
        h1=getheight(clust.left)*20
        h2=getheight(clust.right)*20
        top=y-(h1+h2)/2
        bottom=y+(h1+h2)/2
        # Line length
        ll=clust.distance*scaling
        # Vertical line from this cluster to children    
        draw.line((x,top+h1/2,x,bottom-h2/2),fill=(255,0,0))    

        # Horizontal line to left item
        draw.line((x,top+h1/2,x+ll,top+h1/2),fill=(255,0,0))    

        # Horizontal line to right item
        draw.line((x,bottom-h2/2,x+ll,bottom-h2/2),fill=(255,0,0))        

        # Call the function to draw the left and right nodes    
        drawnode(draw,clust.left,x+ll,top+h1/2,scaling,imlist,img)
        drawnode(draw,clust.right,x+ll,bottom-h2/2,scaling,imlist,img)
    else:   
        # If this is an endpoint, draw a thumbnail image
        nodeim = Image.open(imlist[clust.id])
        nodeim.thumbnail((20,20))
        ns = nodeim.size
        img.paste(nodeim,(x,y-ns[1]//2,x+ns[0],y+ns[1]-ns[1]//2))

The dendrogram is then drawn like this.


hcluster.drawdendrogram(tree,imlist,jpeg='sunset.jpg')

This should give an image like the one below.

Hope this was useful. I have one final small thing about clustering I want to write, something for next week maybe.

Things iPhoto Thinks are Faces

2009-04-06T04:35:00.000-07:00

More on the "Flickr as database for computer vision" topic. On Flickr there is now a group pool where people collect pictures with false face detections from Apple's iPhoto'09. Could turn into an interesting set to use for research and training/improving face detection.

I'm planning to write more about faces soon, I just need some time.

(thanks Thijs for the tip)

Clustering using SciPy's k-means

2009-04-03T00:16:00.000-07:00

K-means is a very simple clustering algorithm that tries to partition the input data in k clusters. The main drawback of this algorithm is that the number of clusters needs to be decided beforehand and an inappropriate choice will give poor clustering results. K-means works by iteratively refining an initial guess of class centroids as follows:

Initialize centroids, randomly or with some guess.
Assign each data point to the class of its nearest centroid.
Update the centroids as the average of all data points within that class.

K-means tries to minimize the total within-class variance. The algorithm above is a heuristic refinement algorithm that works fine for most cases but does not guarantee that the best solution is found. To avoid the effects of choosing a bad centroid initialization, the algorithm is often run several times with different initialization. Then the solution with lowest variance is selected.

Although simple to implement, there is no need to. The SciPy vector quantization package comes with a k-means implementation. Here's how to use it.

Let's start with creating some sample data.
>>> from numpy import *
>>> from scipy.cluster.vq import *
>>> class1 = array(random.standard_normal((100,2))) + array([5,5])
>>> class2 = 1.5 * array(random.standard_normal((100,2)))
>>> features = vstack((class1,class2))
This generates two normally distributed classes in two dimensions. To try and cluster the points, run k-means with k=2 like this.
>>> centroids,variance = kmeans(features,2)
The variance is returned but we don't really need it since the SciPy implementation computes several runs (default is 20) and selects the one with smallest variance for us. Now you can check where each data point is assigned using the vector quantization function in the SciPy package.
>>> code,distance = vq(features,centroids)
By checking the value of code we can see if there are any incorrect assignments. To visualize, we can plot the points and the final centroids.
>>> import pylab
>>> pylab.plot([p[0] for p in class1],[p[1] for p in class1],'*')
>>> pylab.plot([p[0] for p in class2],[p[1] for p in class2],'r*')
>>> pylab.plot([p[0] for p in centroids],[p[1] for p in centroids],'go')
>>> pylab.show()
This should give a plot like the one below (click to enlarge).

This is a very useful tool but be aware that the python implementation is not as fast as c/c++ implementations so this can take some time if you are clustering very large data sets. For many cases though, this is good enough.

ROF de-noising in Python

2009-03-24T15:03:00.000-07:00

These days it seems there are lots of applications where the Rudin-Osher-Fatemi (ROF) de-noising model can be applied (more on that later). So when my friend Niels sent me a Matlab implementation the other day, I couldn't resist liberating it from the evil claws of Matlab. Below is a python implementation that uses NumPy. I copied shamelessly from Niels' original m-file.

Anyway, here's the code.


from numpy import *

def denoise(im, U_init, tolerance=0.1, tau=0.125, tv_weight=100):
    """ An implementation of the Rudin-Osher-Fatemi (ROF) denoising model
        using the numerical procedure presented in Eq. (11) of A. Chambolle
        (2005). Implemented using periodic boundary conditions 
        (essentially turning the rectangular image domain into a torus!).
    
        Input:
        im - noisy input image (grayscale)
        U_init - initial guess for U
        tv_weight - weight of the TV-regularizing term
        tau - steplength in the Chambolle algorithm
        tolerance - tolerance for determining the stop criterion
    
        Output:
        U - denoised and detextured image (also the primal variable)
        T - texture residual"""
    
    #---Initialization
    m,n = im.shape #size of noisy image

    U = U_init
    Px = im #x-component to the dual field
    Py = im #y-component of the dual field
    error = 1 
    iteration = 0

    #---Main iteration
    while (error > tolerance):
        Uold = U

        #Gradient of primal variable
        LyU = vstack((U[1:,:],U[0,:])) #Left translation w.r.t. the y-direction
        LxU = hstack((U[:,1:],U.take([0],axis=1))) #Left translation w.r.t. the x-direction

        GradUx = LxU-U #x-component of U's gradient
        GradUy = LyU-U #y-component of U's gradient

        #First we update the dual varible
        PxNew = Px + (tau/tv_weight)*GradUx #Non-normalized update of x-component (dual)
        PyNew = Py + (tau/tv_weight)*GradUy #Non-normalized update of y-component (dual)
        NormNew = maximum(1,sqrt(PxNew**2+PyNew**2))

        Px = PxNew/NormNew #Update of x-component (dual)
        Py = PyNew/NormNew #Update of y-component (dual)

        #Then we update the primal variable
        RxPx =hstack((Px.take([-1],axis=1),Px[:,0:-1])) #Right x-translation of x-component
        RyPy = vstack((Py[-1,:],Py[0:-1,:])) #Right y-translation of y-component
        DivP = (Px-RxPx)+(Py-RyPy) #Divergence of the dual field.
        U = im + tv_weight*DivP #Update of the primal variable

        #Update of error-measure
        error = linalg.norm(U-Uold)/sqrt(n*m);
        iteration += 1;

        print iteration, error

    #The texture residual
    T = im - U
    print 'Number of ROF iterations: ', iteration
    
    return U,T

Save that in a file "rof.py" (copy available here) and you can use it like this:


from PIL import Image
import pylab
import rof

im = array(Image.open('empire.jpg').convert("L"))
U,T = rof.denoise(im,im)

pylab.figure()
pylab.gray()
pylab.imshow(U)
pylab.axis('equal')
pylab.axis('off')
pylab.show()

Here I used empire.jpg as example. The result should look something like this.

This implementation uses an algorithm proposed by A. Chambolle in 2005. The details and what they mean for other problems in vision is a topic for future posts.

Highlights from SSBA2009

2009-03-22T13:01:00.000-07:00

The annual Swedish image analysis conference SSBA was held in Halmstad Thursday-Friday this week. Always good to catch up to meet and see what the Swedish groups are working on at the moment. The highlights were (according to my own personal bias):

Strandmark, Kahl, Overgaard, Optimal Levels for the Two-phase, Piecewise Constant Mumford-Shah Functional. ROF and branch and bound for solving piecewise constant Mumford-Shah.
Läthen, Andersson, Lenz, Borga, Level set based segmentation using gradient descent with momentum. Adding momentum to curve evolution to bridge gaps in segmentation.
Olsson, Byröd, Overgaard, Kahl, Continuous alpha-Expansion Moves. Continuous version of the alpha-expansion moves for multi-class graph-cut segmentation.
Enqvist, Josephson, Kahl, Correspondence from Pairwise Constraints. Graph based 3D registration using vertex cover problem formulation.

(No online papers to link to at the moment, sorry.)

Jacket - GPU Toolbox for Matlab

2009-03-18T13:26:00.000-07:00

Although I'm trying to stay away from Matlab these days, I saw a nice demo of Jacket that made an impression the other day. Jacket is a GPU Toolbox for MATLAB that lets you compile Matlab code (with very small modifications) to code that will run on CUDA-enabled GPUs (this means NVIDIA only). It is incredibly easy to use and could be an interesting option if you are a Matlab-er with a good (NVIDIA) graphics card.

Jacket is made by AccelerEyes and cost a few hundred USD for student and academic licenses and around 2kUSD for a commercial license.

SixthSense at TED

2009-03-12T13:51:00.000-07:00

Pattie Maes from the MIT Media Lab gives an incredibly cool presentation of SixthSense (developed by Pranav Mistry) at TED. The details are on the project webpage and here:

"SixthSense is a wearable gestural interface that augments the physical world around us with digital information and lets us use natural hand gestures to interact with that information."

Extremely cool and very inspiring!

Event Video Categorization added to Multimedia Grand Challenge

2009-03-02T10:59:00.000-08:00

There is a new challenge added to the Multimedia Grand Challenge 2009, and it is a pretty interesting one. Below are parts of the description quoted.

Accenture Challenge: Analysis of Video Footage Captured in Uncontrolled Environments

"Often, it is necessary to analyze this corpus after (for) an event. An event might involve one or more objects (e.g. people, cars, etc.) and the objects’ interaction with each other. We might then want to search the corpus for similar events or objects that were part of the event."

"What are some categories of objects and a set of higher-level events that the objects could help identify?"
"Can the system identify key objects if given footage of an event?"

Great stuff! Will be interesting to follow the developments.

Using Python to download images from Flickr

2009-02-22T10:10:00.000-08:00

The immensely popular photo sharing site Flickr is a gold mine for computer vision researchers and hobbyists. With hundreds of millions of images, many of them tagged by users, it is a great resource to get training data or for doing experiments on real data. Flickr has an API for interfacing with the service that makes it possible to upload, download and annotate images (and much more). A full description of the API is available here and there are kits for many programming languages, including Python.

Let's look at using a library called flickrpy available freely at http://code.google.com/p/flickrpy/. Download the file flickr.py. You will need an API Key from Flickr to get this to work. Keys are free for non-commercial use. Just click the link "Apply for a new API Key" on the Flickr API page and follow the instructions.

Once you have an API key, open flickr.py and replace the empty string on the line
API_KEY = ''
with your key. It should look something like this:
API_KEY = '123fbbb81441231123cgg5b123d92123'

Let's create a simple command line tool that downloads images tagged with a particular tag. Add the following code to a new file called tagdownload.py.


import flickr
import urllib, urlparse
import os
import sys

if len(sys.argv)>1:
    tag = sys.argv[1]
else:
    print 'no tag specified'

# downloading image data
f = flickr.photos_search(tags=tag)
urllist = [] #store a list of what was downloaded

# downloading images
for k in f:
    url = k.getURL(size='Medium', urlType='source')
    urllist.append(url) 
    image = urllib.URLopener()
    image.retrieve(url, os.path.basename(urlparse.urlparse(url).path)) 
    print 'downloading:', url

If you also want to write the list of urls to a text file, add the following lines at the end.


# write the list of urls to file       
fl = open('urllist.txt', 'w')
for url in urllist:
    fl.write(url+'\n')
fl.close()

From the command line, just type

$ python tagdownload.py goldengatebridge

and you will get the 100 latest images tagged with "goldengatebridge". As you can see, we chose to take the "Medium" size. If you want thumbnails or full size originals or something else, there are many sizes available, check the documentation on the Flickr website.

Here we were just interested in downloading images, for API calls that require authentication the process is slightly more complicated. See the API documentation for more information on how to set up authenticated sessions.

Good luck!

Astronomy on Flickr

2009-02-20T11:58:00.000-08:00

Flickr has quickly come to be one of the best databases for computer vision researchers. Now there is also automatic analysis of images. Images in the Astrometry group are automatically matched and annotated by the "blind astrometry solver" (see example here). Besides computing the coordinates of the part of the sky pictured and lots of other data, structures that are recognized are also annotated as photo notes on the images. Very cool!

See the Flickr blog post for more background information on the algorithms used (geometric hashing) and the people behind the project.

(thanks Brian for the tip)

SIFT Python Implementation

2009-02-12T04:42:00.000-08:00

I'd like to share a Python interface I wrote for David Lowe's Scale Invariant Feature Transform (SIFT) implementation. David, the inventor of SIFT, has since several years generously shared binaries with a Matlab interface on his website. Inspired by the Matlab files for reading keypoint descriptor files and for matching between images, I decided to write a Python version. Here it is: sift.py.

The file itself should be self-explanatory, especially together with the documentation that comes with Lowe's zip-file. Anyway, here's an example:

>>>from PIL import Image
>>>from numpy import *
>>>import sift

>>>sift.process_image('basmati.pgm', 'basmati.key')
>>>l1,d1 = sift.read_features_from_file('basmati.key')
>>>im = array(Image.open('basmati.pgm'))
>>>sift.plot_features(im,l1)

Update: I took the corrections from the comments and updated the file. Thanks all. I also wrote a new post on using VLFeat which is worth checking out for those who want binaries that work for Windows and Mac OS.

Multimedia Grand Challenge 2009

2009-02-07T10:34:00.000-08:00

My friend and ex-Yahoo!er Mor is co-chairing the Multimedia Grand Challenge for ACM Multimedia 2009. Industry partners like Yahoo!, Nokia and Google submitted in total eight challenges. Very interesting! My three favorites:

Nokia Challenge: Where was this Photo Taken, and How?
"This challenge focuses on capture device location and orientation, one dimension of content metadata. The problem can be stated simply: try to derive exact camera poses (location and orientation) of given photos that are lacking location annotation."

Yahoo! Challenge: Robust Clustering Guided by User Intent in Image Searc
"The challenge to researchers in the multi-media community is to 1) develop a robust way of understanding user intent and 2) generate highly relevant clusters for the given intent and query."

Google Challenge: Robust, As-Accurate-As-Human Genre Classification for Video
"The goal of this task would be to take user generated videos (along with their sparse and noisy metadata) and automatically classify them into genres."

Harris Corner Detector in Python

2009-01-26T13:46:00.000-08:00

Thought I'd share a simple Python implementation of the Harris corner detector. I have seen people looking for a python implementation for a range of applications so I'm hoping someone finds this useful.

The Harris (or Harris & Stephens) corner detection algorithm is one of the simplest corner indicators available. The general idea is to locate points where the surrounding neighborhood shows edges in more than one direction, these are then corners or interest points. The algorithm is explained here. In short, a matrix W is created from the outer product of the image gradient, this matrix is averaged over a region and then a corner response function is defined as the ratio of the determinant to the trace of W.

Let's see what this looks like in code. First we need to be able to do convolutions of 2D signals. For this NumPy is not enough and we need to use the signal module in SciPy. Create a file filtertools.py and add the following functions needed to create Gaussian derivative kernels and apply them to the image.


from scipy import *
from scipy import signal


def gauss_derivative_kernels(size, sizey=None):
    """ returns x and y derivatives of a 2D 
        gauss kernel array for convolutions """
    size = int(size)
    if not sizey:
        sizey = size
    else:
        sizey = int(sizey)
    y, x = mgrid[-size:size+1, -sizey:sizey+1]

    #x and y derivatives of a 2D gaussian with standard dev half of size
    # (ignore scale factor)
    gx = - x * exp(-(x**2/float((0.5*size)**2)+y**2/float((0.5*sizey)**2))) 
    gy = - y * exp(-(x**2/float((0.5*size)**2)+y**2/float((0.5*sizey)**2))) 

    return gx,gy

def gauss_derivatives(im, n, ny=None):
    """ returns x and y derivatives of an image using gaussian 
        derivative filters of size n. The optional argument 
        ny allows for a different size in the y direction."""

    gx,gy = gauss_derivative_kernels(n, sizey=ny)

    imx = signal.convolve(im,gx, mode='same')
    imy = signal.convolve(im,gy, mode='same')

    return imx,imy

The point of using Gaussian derivative filters is that this computes a smoothing of the image, to a scale defined by the size of the filter, and the derivatives at the same time. The derivatives are less noisy than if computed with a simple difference filter on the original image.

First add the corner response function to a file harris.py which will make use of the Gaussian derivatives above.


def compute_harris_response(image):
    """ compute the Harris corner detector response function 
        for each pixel in the image"""

    #derivatives
    imx,imy = filtertools.gauss_derivatives(image, 3)

    #kernel for blurring
    gauss = filtertools.gauss_kernel(3)

    #compute components of the structure tensor
    Wxx = signal.convolve(imx*imx,gauss, mode='same')
    Wxy = signal.convolve(imx*imy,gauss, mode='same')
    Wyy = signal.convolve(imy*imy,gauss, mode='same')

    #determinant and trace
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    return Wdet / Wtr

This gives an image with each pixel containing the value of the Harris response function. Now it is just a matter of picking out the information needed from this image. Picking all values above a threshold with the additional constraint that corners must be separated with a minimum distance is an approach that often gives good results. To do this, take all candidate pixels, sort them in descending order of corner response values and mark off regions too close to positions already marked as corners. Add this function to harris.py.


def get_harris_points(harrisim, min_distance=10, threshold=0.1):
    """ return corners from a Harris response image
        min_distance is the minimum nbr of pixels separating 
        corners and image boundary"""

    #find top corner candidates above a threshold
    corner_threshold = max(harrisim.ravel()) * threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    #get coordinates of candidates
    candidates = harrisim_t.nonzero()
    coords = [ (candidates[0][c],candidates[1][c]) for c in range(len(candidates[0]))]
    #...and their values
    candidate_values = [harrisim[c[0]][c[1]] for c in coords]

    #sort candidates
    index = argsort(candidate_values)

    #store allowed point locations in array
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_distance:-min_distance,min_distance:-min_distance] = 1

    #select the best points taking min_distance into account
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i][0]][coords[i][1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i][0]-min_distance):(coords[i][0]+min_distance),(coords[i][1]-min_distance):(coords[i][1]+min_distance)] = 0

    return filtered_coords

Now you have all you need to detect corner points in images. To make it easier to show the corner points in the image you can add a plotting function using matplotlib (PyLab) as follows.


def plot_harris_points(image, filtered_coords):
    """ plots corners found in image"""
    from pylab import *

    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    axis('off')
    show()

Try running the following commands on an example image:


>>> im = array(Image.open('empire.jpg').convert("L"))
>>> harrisim = harris.compute_harris_response(im)
>>> filtered_coords = harris.get_harris_points(harrisim,6)
>>> harris.plot_harris_points(im, filtered_coords)

The image is opened and converted to grayscale. Then the response function is computed and points selected based on the response values. Finally, the points are plotted overlaid on the original image. This should give you a plot like this.

An example of corner detection with the Harris corner detector.

An overview of different approaches to corner detection, including improvements on the Harris detector and further developments, see e.g. http://en.wikipedia.org/wiki/Corner_detection.

Mahout - something to watch in 2009

2009-01-21T10:44:00.000-08:00

Tomorrow is the first birthday of the Apache Mahout project which officially launched Jan 22, 2008. The goal of Mahout is "to build scalable, Apache licensed machine learning libraries" on top of Hadoop, the Apache open-source project for distributed computing with MapReduce.

So far it seems to have been mostly collaborative filtering but Mahout will implement several clustering and classification methods. I'm curious to see what comes out of 2009. Will there be a python interface perhaps?

Happy birthday and best of luck for 2009!

Installing NumPy and SciPy on Mac OS X 10.5 Leopard

2009-01-18T04:32:00.000-08:00

I have had problems getting SciPy to work on my MacBook with Mac OS X 10.5. I decided to try a clean install from scratch and see what works. I have previously used NumPy/matplotlib/PIL etc from MacPorts with no problems, except that SciPy didn't build correctly, so I decided to try MacPorts first.

There are some instructions available for python 2.5 but those didn't work for me and I wanted to try python 2.6. Here is a step-by-step recap of what I did:

* Delete opt/local folder.

* Get MacPorts from http://www.macports.org/install.php

Install Python2.6:

$ sudo port install python26

Then some dependencies that in retrospect, I'm not sure I needed to add separately (you could try skipping these).

$ sudo port install fftw-3 gmp libiconv mpfr

$ sudo port install zlib freetype fontconfig jpeg libpng gd2 pdflib ncursesw ncurses readline gnuplot

[note to self:
Gnuplot failed with the following error.
"Image error: /Applications/MacPorts/AquaTerm.app/Contents/Info.plist already exists and does not belong to a registered port. Unable to activate port aquaterm."
Turns out AquaTerm puts a file AquaTerm.app in that folder and I didn't delete the old one. I deleted it and then activated the new AquaTerm with:
$ sudo port install aquaterm ]

The python 2.5 instructions recommends using easy_install to install NumPy and SciPy after this step but that didn't work for me. The installation finishes but it doesn't work - import numpy or import scipy fails. So, skipping easy_install altogether and putting my trust in MacPorts:

$ sudo port install py26-numpy

Works fine. Now for the real test.

$ sudo port install py26-scipy

SciPy for python 2.6 depends on gcc 4.3 and the gcc43 package takes several hours to build, have patience. After gcc43 and some other dependencies, SciPy installs fine!

After testing some code with NumPy and SciPy I found out that matplotlib does not support python 2.6. No packages exist for 2.6 but after some investigations I decided to try the current source from the repository.

$ svn co https://matplotlib.svn.sourceforge.net/svnroot/matplotlib/trunk/matplotlib matplotlib

In the matplotlib direcory, run

$ python setup.py build

and then

$ sudo python setup.py install

It works! Plotting does not need to open X11 like in earlier versions and the plots are faster and look a little nicer.

Finally, I added PIL for image handling:

$ sudo port install py26-pil

That's everything. What a wonderful way to spend a Saturday afternoon & evening.

Thanks Geoff!

2009-01-16T13:10:00.000-08:00

My friend Geoff Parker helped me with styling for my personal page at Lund University. I'm very happy about the way it looks, thanks Geoff!

PCA for Images using Python

2009-01-12T15:49:00.000-08:00

I've looked around for a nice and clean implementation of principal component analysis in python but couldn't really find anything that fits. PCA Module looks useful but seems like a lot of stuff for something this simple. So anyway, I decided to spend part of my Monday evening writing a simple implementation. Here it is.

Principal component analysis (PCA) is a useful technique for dimensionality reduction and is optimal in the sense that it represents the variability if the training data with as few dimensions as possible. The resulting projection matrix can be seen as a change of coordinates to a coordinate system where the coordinates are in descending order of importance.

To do PCA on image data, the images need to be converted to a one-dimensional vector representation, e.g. using flatten() in numpy. A tiny 100*100 pixel grayscale image is then a point in a 10,000 dimensional space, a megapixel image has dimensions in the millions, so dimensionality reduction comes handy in many applications.

The flattened images are collected in a single matrix by stacking them, one row for each image. The rows are then centered around the mean image before the computation of the dominant directions. To find the principal components, singular value decomposition (SVD) is usually used but if the dimensionality is high, there is a useful trick that can be used instead since the
SVD computation will be very slow in that case. Here is what it looks like in python code.


from PIL import Image
from numpy import *

def pca(X):
  # Principal Component Analysis
  # input: X, matrix with training data as flattened arrays in rows
  # return: projection matrix (with important dimensions first),
  # variance and mean

  #get dimensions
  num_data,dim = X.shape

  #center data
  mean_X = X.mean(axis=0)
  for i in range(num_data):
      X[i] -= mean_X

  if dim>100:
      print 'PCA - compact trick used'
      M = dot(X,X.T) #covariance matrix
      e,EV = linalg.eigh(M) #eigenvalues and eigenvectors
      tmp = dot(X.T,EV).T #this is the compact trick
      V = tmp[::-1] #reverse since last eigenvectors are the ones we want
      S = sqrt(e)[::-1] #reverse since eigenvalues are in increasing order
  else:
      print 'PCA - SVD used'
      U,S,V = linalg.svd(X)
      V = V[:num_data] #only makes sense to return the first num_data

  #return the projection matrix, the variance and the mean
  return V,S,mean_X

The threshold for when to switch from SVD to use the trick with computing eigenvectors of the (smaller) covariance matrix X*X^T is a bit arbitrary in this case but that doesn't really matter. The rows of the matrix V are orthogonal and contains the coordinate directions in order of descending variance of the training data.

Let's try this on an example of face images downloaded from the "Labeled Faces in the Wild" data set available at http://vis-www.cs.umass.edu/lfw/index.html.

The file gwb_cropped.zip contains 20 selected images of George W Bush, cropped slightly tighter to the face than the originals. Assuming that the filenames of these cropped face images are stored in a list, imlist, the principal components can be computed and shown like this.


from PIL import Image
import numpy
import pylab

im = numpy.array(Image.open(imlist[0])) #open one image to get the size
m,n = im.shape[0:2] #get the size of the images
imnbr = len(imlist) #get the number of images

#create matrix to store all flattened images
immatrix = numpy.array([numpy.array(Image.open(imlist[i])).flatten() for i in range(imnbr)],'f')

#perform PCA
V,S,immean = pca(immatrix)

#mean image and first mode of variation
immean = immean.reshape(m,n)
mode = V[0].reshape(m,n)

#show the images
pylab.figure()
pylab.gray()
pylab.imshow(immean)

pylab.figure()
pylab.gray()
pylab.imshow(mode)

pylab.show()

Note that the images need to be converted back from the one-dimensional representation using reshape(). Running the example should give two plots like these.

The mean image of the 20 sample images and the first mode, i.e. the direction with most variation.

This example is of course a bit silly since, to have something useful, the images should be properly aligned and cropped or masked to remove the influence of unimportant pixels.

Anyway, hope this is useful to someone. I'll write an update with some better examples (I have some ideas for cool data sets) when there is time.

Face Recognition in iPhoto

2009-01-08T11:40:00.000-08:00

The Apple keynote at Macworld yesterday seems to have gotten a lukewarm response in the press. There was one very interesting announcement though; face recognition for tagging and organizing photos in iPhoto. Similar to Picasa's "name tags" (but oh so much prettier, Apple is the god of user interface) and the Polar Rose partner solution (example here). If you haven't seen the iPhoto demo, it is here.

Associate Professor - Starting Today

2009-01-07T14:23:00.000-08:00

Today I started my new gig as part-time associate professor at the Mathematical Imaging Group at Lund University. Not so much action on this first day but I did get keys, a place to sit and a chance to talk research with some old friends and colleagues.

I'm planning to write here regularly about things that interest me, mainly about computer vision (research as well as applications), about start-ups and occasionally about private/random things.