Dijkstra on boost grid graph

There is a whole science of finding tubular structures in
3D-images. It started sometime in the late nineties. One of the most
cited papers¹ in this area shows how to use eigenvalues of the
local second order derivatives to identify tubular structures
heuristically. It also shows how to combine the results on multiple scales.

Raw Data

After Pre-processing

That is quite neat, but the approach has a major problem: It does not
assume any global structure. Therefore, when the data does not have
evidence for tubular structures, you still do not see them, and
furthermore, you have to have a method that allows you to say "points
x and y lie on the same tube". This, and determining properties of the
tubes, are the problems which most of the followupr research strives
to solve.

I wanted to have a quick-and-dirty approach to the problem. For a
given seed point, I would like to know the connectivity of the tubular
structure to this seed point. I pre-processed the raw image using the
method of Sato et al.¹, determined the seed point, the graph is
defined implicitly over the 3D-image using voxel-neighborhoods, so we
could in theory just run Dijkstra on it. I was lazy, however and did
not want to implement it all.

BGL Grid Graph

A quick google search turned up boost::grid_graph, a relatively new
member of the boost graph library (BGL). It defines a graph on a grid
without representing everything. Neat. The documentation is still a
bit shortish and sadly, only a 6-neighborhood (in 3D) is supported,
but who'd expect a free lunch?


Usability is (at first) easy as pie:

typedef grid_graph<3> graph_t;
typedef graph_traits<graph_t> Traits;
typedef Traits::edge_descriptor edge_descriptor;
typedef Traits::vertex_descriptor vertex_descriptor;

// ...

boost::array<vidx_t, 3> lengths = { { 256, 256, 256 } };
graph_t graph(lengths, false); // do not wrap any dimensions

Defining properties of edges in a grid graph

Contrary to graphs with internal properties, we have to define
external properties for graph_t. This is currently not documented
well, and it took me some time to figure out:

• In most BGL examples, nodes are indexed using an integer type. The
  grid_graph, however, uses n-tuples of type boost::array<n>.
  This means, most examples cannot by applied to grid_graph in a
  straight-forward way: You cannot index an array by a tuple without
  doing some tricks first.
  
  
• We would like to have an arbitrary function of the nodes as the
  weight. Most BGL examples simply use an array for whichever
  property they put in their graph.
  

The first problem can be solved using shared array property maps
(which do not seem to be well-documented either):

shared_array_property_map<vertex_descriptor,
  property_map<graph_t, vertex_index_t>::const_type>
  p_map(num_vertices(graph), get(vertex_index, graph));
shared_array_property_map<double,
  property_map<graph_t, vertex_index_t>::const_type>
  d_map(num_vertices(graph), get(vertex_index, graph));

We will use p_map as our predecessor map, which maps a vertex index (integer)
to a vertex descriptor (n-tuple).

The second map d_map is a mapping from vertex index to double, which
will be our distance map.


Finally, we need to create a edge_weight for graph_t, and it
should be a function, not an array. For this, I peaked W.P. McNeills
Boost-Implicit-Graph Example on GitHub.

struct edge_weight_map;
namespace boost {
  template<>
  struct property_map< graph_t, edge_weight_t > {
    typedef edge_weight_map type;
    typedef edge_weight_map const_type;
  };
}

/*
   Map from edges to weight values
*/
struct edge_weight_map {
  typedef double value_type;
  typedef value_type reference;
  typedef edge_descriptor key_type;
  typedef boost::readable_property_map_tag category;
  const graph_t& m_graph;
  edge_weight_map(const graph_t& g)
  :m_graph(g) { }
  reference operator[](key_type e) const; // implemented below
};

typedef boost::property_map<graph_t, boost::edge_weight_t>::const_type
        const_edge_weight_map;
typedef boost::property_traits<const_edge_weight_map>::reference
        edge_weight_map_value_type;
typedef boost::property_traits<const_edge_weight_map>::key_type
        edge_weight_map_key;

namespace boost{
    // PropertyMap valid expressions
    edge_weight_map_value_type get(const_edge_weight_map pmap, edge_weight_map_key e) {
        return pmap[e]; }
    // ReadablePropertyGraph valid expressions
    const_edge_weight_map get(boost::edge_weight_t, const graph_t&g) {
        return const_edge_weight_map(g); }
    edge_weight_map_value_type get(boost::edge_weight_t tag, const graph_t& g, edge_weight_map_key e) {
        return get(tag, g)[e]; }
}

// whoa, lots of typedefs, but now we can write the function that we
// wanted to write: Map edges to weights!
edge_weight_map::reference
edge_weight_map::operator[](key_type e) const {
    vertex_descriptor t = target(e,m_graph);
    vertex_descriptor s = source(e,m_graph);
    // f = f(t,s)
    return f;
}

That was more nasty than I expected, but it is still not that much
code to write.

Representing the data in a Boost Multi-Array

Finally, we need to represent 3D-data. I hate dealing with raw memory,
so I always wrap it in some convenient abstracting structure first
chance I get. 3D-arrays are definitely a case for Boost Multi-Array.
Its not a big deal using it, but I want to advertise its simplicity here:

typedef boost::multi_array_ref<unsigned char, 3> array_type;

std::ifstream dat("raw.dat", std::ios::in | std::ios::binary);
unsigned char* data = new unsigned char[256*256*256];
dat.read((char*)data,256*256*256);
array_type A(data,boost::extents[256][256][256]);
dat.close();

By the way, you can also change the base of this array as well as the
memory order, so if you need to port Matlab (TM) code to a programming
language (no adjective here on purpose), you can't get much simpler.

Finish line: Running Dijkstra

This is trivial now and straight from the book:

boost::array<vidx_t, 3> start  = { { 109, 129,  24 } };
dijkstra_shortest_paths(graph, start
        ,predecessor_map(p_map)
        .distance_map(d_map)
        );

Note the "." before distance_map, this is a named argument called
distance_map, which is set to d_map. The same goes for p_map.

Thresholded Distances

We can now trace paths starting at points which have a high value in
the raw image:

Traced paths starting at thresholded raw data

Footnotes:

¹ "3D multi-scale line filter for segmentation and visualization
of curvilinear structures in medical images" by Sato et al.,
CVRMed-MRCAS, 1997 (PDF)

Nearest Neighbor Classifier in CUV

After reading this blog post by Alex Smola I implemented a
one-nearest-neighbor classifier in CUV. Apart from a bug that popped
up, I was happy to see that no changes to the library where needed,
and the code is very short, so I'll share the implementation here.


There is a small annoyance which I'll fix later:

• there is no argmin_col, so we instead resort to multiplying by
  -1 and then use argmax_col for now. I'll clean this mess up
  later, probably by adding a argFunc_col functor.

import cuv_python as cp
import pyublas
import numpy as np

class KNN:
    """ calculates the labels of a test set by determining the closest
    instance in a training set and returning the corresponding label."""
    def __init__(self, data, data_l):
        """
        data:   a (number_instances x dimensionality) numpy matrix
        data_l: a number_instances numpy vector containing the labels
        """
        self.data   = cp.push(data)
        self.data_l = data_l
        self.dsq    = cp.dev_matrix_cmf(self.data.h,1)
        cp.reduce_to_col(self.dsq.vec,self.data,cp.reduce_functor.ADD_SQUARED)
    def __get_distance_matrix(self, test):
        """
        test: a (number_test_instances x dimensionality) numpy matrix
        returns: a (number_instances x number_test_instances) CUV distance matrix
        """
        t   = cp.push(test)
        assert t.w == self.data.w
        tsq = cp.dev_matrix_cmf(t.h, 1)
        cp.reduce_to_col(tsq.vec,t,cp.reduce_functor.ADD_SQUARED)
        p   = cp.dev_matrix_cmf(self.data.h, t.h)
        cp.prod(p, self.data, t, 'n','t',-2, 0)
        cp.matrix_plus_col(p,self.dsq.vec)
        cp.matrix_plus_row(p,tsq.vec)
        return p
    def run(self,test):
        """
        test:    a (number_test_instances x dimensionality) numpy matrix
        returns: the estimated labels of the test instances
        """
        p = self.__get_distance_matrix(test)
        p *= -1.                # no argmin supported yet, sorry
        idx = cp.dev_matrix_cmi(test.shape[0],1)
        cp.argmax_to_row(idx.vec, p)
        hidx  = idx.np.reshape(idx.h)
        return self.data_l.reshape(self.data.h)[hidx]

As suggested in the original post, this is a very fast method to
calculate the nearest neighbor. Instead of looping over all possible
pairs $x∈ X$ and $z∈ Z$ calculating the distance

\[D_{ij} = \|x_i-z_j\|^2,\]


we first precalculate $\|x_i\|^2$ and $\|z_i\|^2$ and then use a fast
matrix multiplication and two additions to determine the distance
above (derived by expanding the second binomial formula)


\[D_{ij} = \|x_i\|^2 + \|z_i\|^2 - 2 XZ^T.\]


We can measure the speed of the implementation by running the above
program using GPU or CPU. The CPU implementation uses cBLAS, I'm using
a 3.2 GHz CPU and a GTX580 for comparison.

Method	Time (s)	Relative speed
CUV CPU (cBLAS)	58.35	34
CUV GPU (cuBLAS)	1.71	1

To test this on MNIST, save the above as knn.py and run the following source:

import os
from knn import KNN
import numpy as np

class MNIST:
  """ This simply reads the MNIST files to main memory and converts them to float """
  def __init__(self,dir):
      from scipy.io.numpyio import fread
      fd = open(dir+'/train-labels.idx1-ubyte')
      fread(fd,8,'c')
      self.data_labels = np.fromfile(file=fd, dtype=np.uint8).reshape( 60000 )
      fd.close()

      fd = open(dir+'/train-images.idx3-ubyte')
      fread(fd,16,'c')
      self.data = np.fromfile(file=fd, dtype=np.uint8).reshape( (60000,784) )
      fd.close()

      fd = open(dir+'/t10k-images.idx3-ubyte')
      fread(fd,16,'c')
      self.test = np.fromfile(file=fd, dtype=np.uint8).reshape( (10000,784) )
      fd.close()

      fd = open(dir+'/t10k-labels.idx1-ubyte')
      fread(fd,8,'c')
      self.test_labels = np.fromfile(file=fd, dtype=np.uint8).reshape( 10000 )
      fd.close()
  def get_test(self):
      v = self.test.astype('float32').T.copy("F")
      t = self.test_labels
      return v,t
  def get(self):
      v = self.data.astype('float32').T.copy("F")
      t = self.data_labels
      return v,t

pg = MNIST("/home/local/datasets/MNIST");

data, data_l  = pg.get()
test, test_l = pg.get_test()

knn = KNN(data.T.copy("F"),data_l)

off, err_cnt = 5000, 0
for i in xrange(0,10000,off):
    pred = knn.run(test[:,i:i+off].T.copy("F"))
    err_cnt += (pred!=test_l[i:i+off]).sum()

print "Errors: ", err_cnt

Hurricane in the Server Room

Our new toy arrived, with twelve cores, four GTX580 GPUs, four
terabyte hard disk and 24 Gigabytes RAM. The thing looks quite
impressive from the inside:

Note the four GPUs on the left hand side and the vertical stack of
fans in the center. Together with the CPU fans above the GPUs it is
impossible to sit next to this machine and think clearly. It blows
your head off, or at least it sounds as if it would. Even when almost
unused, this machine produced so much heat our server room could not
discharge the heat anymore. Amazing.

If you ever wondered: A current Ubuntu Server version seems to have no
problems with this specification.

Benchmarking GPUs using the CUV library

Since I started programming GPUs, a few generations of these cards
have been released, most recently the GTX580. I always wondered how
good these actually are and how programs written for one GPU scale to
the next generation.


We have a test suite available here, which is at least of importance
to us: The CUV library. Apart from unit-tests checking correctness of
the implementation, the library also has a few "tests" which measure
execution speed on GPU and CPU for comparison.


Instead of inventing new tests for the benchmark, we simply reuse the
speed tests which come with CUV, as they are probably relevant
use-cases that the programmers optimized for, anyway. A small perl
script that now resides in the scripts/ directory of CUV now
identifies speed tests by their name (*_speed), runs them and parses
their output. We simply collect all values in a defined order and save
them to a file. A larger number of these files can then be analyzed
using a python script included in the same directory, which uses the
superb matplotlib library to draw a bar chart. We use a reference GPU
to compare relative timings.


To cut a long story short, here are the results comparing GTX285,
GTX295, GX2-9800, GTX480 and GTX580:

As expected, Fermi generation cards perform a lot better than the
older generation, the GTX580 also improves on GTX480. Some operations,
which are apparently not implemented well, perform worse with newer
generation cards. The real work horses of our library have definitely
improved a lot over generations, even though we did not spend time on
optimizing them for the later cards.


We're not the first ones to compare these cards of course. Legit Reviews compares the frame rate rendered by the GTX480 and GTX580,
finding only marginal improvements. Brightsideofnews runs many
rendering related benchmarks as well (3DMark Vantage, Unigine Heaven,
Pripyat), with mixed results. However, we're more interested in
general purpose computing (GPGPU) and in writing algorithms for GPU
hardware which then scale with the new hardware, as it becomes
available.

Directory Color Embedding

Do you also save all your experiment data in folder names which
represent the settings? I do. If you have many experiments, however,
this strategy will defy comparison of the experiments, as it becomes
hard to put everything in one plot. The default settings in plotting
programs such as gnuplot or matplotlib do not have a large enough
range of colors, and if they do, it becomes hard to assign the colors
in a meaningful way for exploratory data analysis. The script here
might come to rescue.

Reading pre-processing directory names

We first transform directory names such that test-paramA_0.1 becomes
test-paramA_0000.1, allowing easier string comparison.

import Levenshtein as L
import numpy as np
from colormath.color_objects import LabColor
import os, re

def fillfunc(mo):
    """ returns the matched number, with padded zeros to the front """
    s = mo.group(1)
    while(len(s)<6): s = '0'+s
    return s
def get_unified_names(fns):
    n = len(fns)
    strs = [x for x in fns]
    for i in xrange(len(strs)):
        strs[i] = re.sub(r'(?<!\d)([\d.]+)(?!\d)',fillfunc,strs[i])
    return strs

Creating the distance matrix based on Levenstein string distance

The Levenshtein string distance is a measure of how many
edits/insertions/deletions one needs to transform one string into the
other. With this we automatically derive a dissimilarity measure for
the unified directory strings:

mat = np.zeros((n,n))
for i in xrange(n):
    for j in xrange(n):
        if j<=i: continue
        mat[i,j] = L.distance(strs[i],strs[j])
        mat[j,i] = mat[i,j]
return np.exp(-mat)

Embedding the dissimilarity matrix using Multidimensional Scaling (MDS) in R

The distance matrix does not necessarily directly map to coordinates
in 3D, so we need an embedding algorithm which deals with distance
matrices only.


For simplicity, we call R from python with a saved dissimilarity matrix.

def emb(mat):
    np.savetxt("file.dat", mat)
    os.system("./analyse/emb.R")
    res = np.loadtxt("mds.dat")
    return res

Now to the part written in R, which only does the embedding and saves
the result in a text file:

#!/usr/bin/Rscript
library(MASS)
library(vegan)
tab <- read.table("file.dat", header = FALSE, sep=" ")
data.m    <- as.matrix(tab)
data.mds <- vegan::metaMDS(data.m, k=3, trymax=50)
write.table(data.mds$points, file="mds.dat", quote=F, row.names=F, col.names=F)

Transforming embedded coordinates into colors

The coordinates are in some – it seems arbitrary – range and
therefore need to be mapped to colors. We want similar colors to have
similar edit distance, we therefore map our coordinates directly to
Lab Color Space and then transform the Lab color coordinates to RGB:

def getcolor(embedded,i):
    #lab = LabColor(0.8,*embedded[i,:])
    lab = LabColor(*embedded[i,:])
    rgb = lab.convert_to('RGB', debug=False).get_rgb_hex()
    return rgb

The resulting string can be used directly in the color specification
of a matplotlib plot and is then best combined with picking to find
out more about a particularly interesting plotline.

Example Embedding of directory names