Getting there

After fixing some inlining problems and comparing the Python and C++ versions in a fairer way (read: with equal amounts of makeUnion calls), the C++ implementation is now blazingly fast. Compared to raw pointer arithmetics, it takes just 1.5 times as long, and the numpy version is much slower when the arrays are non-trivial. This means I actually got time to write lazyflow operators now!


More Iterations

I experimented a bit with different iteration schemes in vigra, and ended up in stealing the one used in labelVolume: Some nested functions that split the multi array into a list of iterators, shapes and accessors. I am not quite sure why this works so much better than the CoupledIterator approach (about a factor of 10 I’d say), but I will stick to that. The original algorithm from which I’m stealing is regrettably just for 3-dimensional arrays, but I wanted to be a bit more generic.

The new and improved mergeLabels function calls itself recursively for every slice of the outermost dimension. This takes about 1.5 times as long, but we don’t have to keep several versions of the merging function for different dimensions lying around. What also helped in improving speed is to bind singleton dimensions — this could be an indicator that the recursive function calls are not inline’d and therefore causing the .5 overhead.

I keep this improvements separate from the current developments in lazyflow (aka OpLazyConnectedComponents) because I wanted to just mess around with this. For the next week I will return to making the operator work smoothly and for general dimensions — I got the feeling that I spend to much time looking at template error messages.




Numpy is pretty impressive

As it turns out, my C++ functions suck. Really bad. Since my python functions turned out to be alright, I started developing their vigra counterparts. After a long night I finally got the python interface to the UnionFind data structure somewhat right, and I turned to the implementation of the mergeLabels function.

This function should simply inspect two adjacent chunks, find labeled regions that cross the boundary and merge the corresponding labels in the global UnionFind structure. A fairly simple task, as you can see from my python implementation:

def mergeLabels(hyperplane_a, hyperplane_b,
                label_hyperplane_a, label_hyperplane_b,
                mapping_a, mapping_b, GUF):
    # the indices where objects are adjacent
    idx = hyperplane_a == hyperplane_b
    # merge each pair of labels
    for label_a, label_b in zip(label_hyperplane_a[idx],
        GUF.makeUnion(mapping_a[label_a], mapping_b[label_b])

But converting this into a vigra function turned out to be not that trivial. Consider this naive port (skip to line 24 if you are not familiar with vigra) :

template void mergeLabels(MultiArrayView const & left,
                          MultiArrayView const & right,
                          MultiArrayView const & leftLabels,
                          MultiArrayView const & rightLabels,
                          MultiArrayView<1, LabelType> const & leftMap,
                          MultiArrayView<1, LabelType> const & rightMap,
                          detail::UnionFindArray & unionFind
                          ) {   
    vigra_precondition(left.shape() == right.shape(), "Shape mismatch");
    vigra_precondition(leftLabels.shape() == rightLabels.shape(), "Shape mismatch");
    vigra_precondition(leftLabels.shape() == left.shape(), "Shape mismatch");

    const MultiArrayIndex LEFT_PIXEL = 1;
    const MultiArrayIndex RIGHT_PIXEL = 2;
    const MultiArrayIndex LEFT_LABEL = 3;
    const MultiArrayIndex RIGHT_LABEL = 4;

    typedef typename CoupledIteratorType::type Iterator;

    Iterator start = createCoupledIterator(left, right, leftLabels, rightLabels);
    Iterator end = start.getEndIterator();

    for (Iterator it = start; it < end; it++)
        if (it.get() == it.get())

This does actually do what we want, and we are using a C-loop over the blocks so we should be fine performance-wise, shouldn’t we? Guess again. A little benchmark showed these results:

pyMergeLabels() for shape (10, 100, 100):    0.014scMergeLabels() for shape (10, 100, 100):    0.293s

This is just sad. Numpy is actually about 20 times faster than my C++ function, even with a favourable axis order. I guess I’m gonna stick to my python function until I figure out why my array traversal is so slow.


Region Growing

Yay, I got accepted into the Google Summer of Code project ‘Lazy Connected Components’! The official start of development will be on May 17, but I have lectures to visit and excercises to do, so I decided to extend the community bonding phase and do some real development. In fact I feel quite bonded to the HCI already, so I guess this is fine.

For my first steps I created a repo outside ilastik / vigra so I can do some experimenting. In the latter stages of the project my code will have to be ported from there to the respective project’s git repository though.

Here’s a little sketch of what my plans are for tackling the lazy connected components problem. I am going to explain things in 2d, but the extension to arbitrary dimensions should be clear.

Suppose we got an input image X, and a user wants a ROI of that image labeled. As I have explained in the initial post, the labels should be consistent if the user requests some other ROI in the future. This means we have to make sure that all objects (connected components) which are visible in the current ROI have to be processed in this step. Consider some horse shoe shaped object: If I label the one end red, the other end independently blue, and later I find out that red and blue are actually the same. I cannot tell the user ‘Hey, I got the labeling wrong, just think of blue as red’.

We dubbed the process of extending the region such that all objects are within the extension as ‘region growing’. The algorithm would look something like this:

Given a labeled chunk and a set of labels to finalize (initially this would be all labels in the chunk)

  • for every adjacent chunk

    • label classically
    • find objects with labels to be finalized that are on the chunk boundary
    • merge the global labels for these objects
    • repeat from the beginning with the neighbour chunk and the labels that where merged

Someone could state that this looks pretty much like some classical connected components algorithm, and indeed the only part that differs is the ‘keep track of labels’ part. If we look at the problem from a graph theory point of view (that’s where the term ‘connected components’ originated, I guess), we can interpret the image as a grid graph and compute spanning trees. One could argue that the problem is essentially the same for chunks. We can take each chunk as a node and compute spanning trees for them as well. The difference is that connectedness in terms of chunks is not transitive: If ROI A is connected to ROI B via object 1 and ROI B is connected to ROI C via object 2, A and C are in general not connected. That’s why we need to keep a list of labels that we are actually interested in, or we would most likely label the whole image every time.

There is already a good amount of python code and tests that implements this strategy in my repo. The probem is split into these parts:

  • classical labeling of chunks
  • merging of adjacent chunks
  • region growing

The original idea was that only the latter should be implemented in python, but as it turns out the merge function is a three-liner which can be handled with numpy quite efficiently. Which means that the lazyflow operator will be all python (except the UnionFind structure from vigra). For the later stages of the project, Ulli and I agreed on implementing a vigra version of the lazy connected components algorithm which will make use of the new ChunkedArrays. And it will be implemented in terms of general graphs, which is a prospect that I really enjoy!