Active Contours are mathematical models that enable detection of objects within images, and are extensively used in Computer Vision as self-adapting frameworks for the delineation and tracking of objects.

To demonstrate Jacket's cross-platform versatility, we implemented the Chan Vese contour tracking app on Android. The video can be viewed <a title="Jacket Mobile" href="http://www.accelereyes.com/products/mobile">here</a>.

Today, however, we'd like to use a MATLAB implementation of active contours as an example of how to take a large project, and with minimal changes, achieve speedups with Jacket.

We'll dangle the proverbial carrot first: the GPU Chan-Vese implementation contains only <em>three kinds of changes overall,</em> and the computational code is <em>exactly the same</em> for both CPU and GPU versions. Plus, take a look at the speed-ups below! How did we accomplish that?
<p style="text-align: center;"><a href="http://blog.accelereyes.com/blog/wp-content/uploads/2011/06/chanvese-jacket-1.8-benchmark-chart.png"><img class="aligncenter size-full wp-image-1394" title="Runtimes for Chan-Vese Contours " src="http://blog.accelereyes.com/blog/wp-content/uploads/2011/06/chanvese-jacket-1.8-benchmark-chart.png" alt="" width="700" height="512" /></a></p>
Before we jump into the code, you might want to read up on Jacket's <a href="http://wiki.accelereyes.com/wiki/index.php/CLASS">CLASS</a> functionality and a <a href="http://blog.accelereyes.com/blog/2011/03/17/write_once_run_everywhere/">blog post</a> that discusses how to exploit it.

<span style="text-decoration: underline;"><em>0. Identify the bottleneck</em></span>

As with any project, the biggest question is where to begin. Fortunately, MATLAB's Profiler will give you a good sense of where the bottlenecks are:

[caption id="attachment_1371" align="aligncenter" width="668" caption="MATLAB Profile for CPU Chan-Vese function"]<a href="http://blog.accelereyes.com/blog/wp-content/uploads/2011/06/profile-overall.png"><img class="size-full wp-image-1371" title="MATLAB Profile for CPU Chan-Vese function" src="http://blog.accelereyes.com/blog/wp-content/uploads/2011/06/profile-overall.png" alt="Shows kappa.m to be the bottleneck" width="668" height="86" /></a>[/caption]

From the profile, we have our target: the function kappa() that consumes about 60% of the time.

We now formulate a set of objectives to ensure proper Jacketization:
<ul>
	<li>Change CPU code as necessary to ensure optimized GPU performance. Target bottleneck first.</li>
	<li>Ensure that all <em>large</em> intermediate variables are generated on the GPU</li>
	<li>Minimize operations between mixed CPU and GPU types (scalars and small-sized variables are okay)</li>
	<li>Keep the computational chunk of code exactly the same between CPU and GPU runs.</li>
</ul>
<span style="text-decoration: underline;"><em>1. Change CPU code</em></span>

Some investigation reveals that in kappa.m:
<ul>
	<li>The same matrix, P, is indexed in many different patterns (an operation called "<a href="http://wiki.accelereyes.com/wiki/index.php?title=SUBSREF"><span style="font-family: monospace;">subsref</span></a>")</li>
	<li>Some index patterns repeat themselves.</li>
</ul>
While straight-up subsrefs are generally fast, the GPU starts to slow-down when it's asked to make "uncoalesced" memory accesses. This most often happens during an indexing operation like in kappa.m. Other examples include conditional indexing, such as <span style="font-family: monospace;">A(A&lt;1)</span>.

Jacket has many built-in workarounds for commonly-used subsrefs where the GPU is slow, but a piece of code with lesser indexing operations might win over code that has more.

Here is the original kappa.m:
<pre>fy = P(3:end,2:n+1)-P(1:m,2:n+1);
fx = P(2:m+1,3:end)-P(2:m+1,1:n);
fyy = P(3:end,2:n+1)+P(1:m,2:n+1)-2*I;
fxx = P(2:m+1,3:end)+P(2:m+1,1:n)-2*I;
fxy = 0.25.*(P(3:end,3:end)-P(1:m,3:end)+P(3:end,1:n)-P(1:m,1:n));</pre>
And here's the modified version:
<pre>subP1 = P(3:end,2:n+1); subP2 = P(1:m,2:n+1);
subP3 = P(2:m+1,3:end); subP4 = P(2:m+1,1:n);
subP5 = P(3:end,3:end); subP6 = P(1:m,3:end);
subP7 = P(3:end,1:n); subP8 = P(1:m,1:n);
fy = subP1-subP2; fx = subP3-subP4;
fyy = subP1+subP2-2*I; fxx = subP3+subP4-2*I;
fxy = 0.25.*(subP5-subP6+subP7-subP8);</pre>
<em><span style="text-decoration: underline;">2. Generate intermediate variables on GPU</span>
</em>

Having changed our M code, we now need to ensure that all <em>large</em> helper variables are generated on the GPU. This is where the CLASS construct comes in handy.

When using CLASS,  the focus is on keeping the functional part of the code identical for both CPU and GPU. To accomplish that, we need to follow the control flow through the code and change variables as necessary.

At a very high level, the flow proceeds as shown:
<ol>
	<li>Generate inputs</li>
	<li>Call chenvese() on those inputs</li>
	<li>Generate masks if necessary (maskcircle2)</li>
	<li>Call heaviside(), kappa()</li>
	<li>Plot images</li>
</ol>
Some careful inspection of the code will show that the helper functions (<span style="font-family: monospace;">maskcircle2, heaviside, kappa</span>) all contain helper variables that are hard-coded to be CPU variables. You can make them change their data-types according to the primary input, as follows:

In function heaviside():
<pre>H=zeros(size(z,1),size(z,2));</pre>
is modified to follow the data type of that function's primary input, z:
<pre>H=zeros(size(z,1),size(z,2),class(z));</pre>
<span style="text-decoration: underline;"><em>3. Minimize operations between mixed types</em></span>

The original M-code contained lots of hard-coded casts to <span style="font-family: monospace;">uint8</span>. For example, in function chenvese():
<pre>P = rgb2gray(uint8(I));</pre>
Here, I is the input image. You need P to be of type <span style="font-family: monospace;">uint8</span> or <span style="font-family: monospace;">guint8</span>, depending on whether I is CPU or GPU respectively. How do you make this statement follow the data-type of I?

This is accomplished using a simple function, smartcast(), with the syntax:
<pre>op = smartcast(ip, type)</pre>
For example,
<pre>op = smartcast(ip, 'uint8')</pre>
will mean that OP is a <span style="font-family: monospace;">guint8</span> or <span style="font-family: monospace;">uint8</span>, depending on whether IP is GPU or CPU, respectively.
You can download <span style="font-family: monospace;">smartcast.m</span> <a href="http://blog.accelereyes.com/blog/wp-content/uploads/2011/06/smartcast.zip">here</a>.

And that's it! Lather, rinse, repeat.

<span style="text-decoration: underline;"><em>4. Keep CPU and GPU code same</em></span>

Our fourth objective is automatically achieved if we follow the steps outlined above!

This code was run on two different CPUs with a Tesla C2070 and Jacket 1.8rc (downloadable <a href="http://www.accelereyes.com/nightly">here</a>). It gave a <strong>12x</strong> speedup on a 1500x1200 image.

You can download the Jacketized versions <a href="http://blog.accelereyes.com/blog/wp-content/uploads/2011/06/chanvese.zip">here</a>. Run <span style="font-family: monospace;">demochanvesecpu</span> for the CPU code, and <span style="font-family: monospace;">demochanvesegpu</span> for the GPU version (they both call the same functions). You can also search the M-files for the word "accelereyes" for the points where the original code was modified.

The original M-code implementation of Chan-Vese active contours was written by Yue Wu at Tufts University. Yue maintains a <a href="http://sites.google.com/site/rexstribeofimageprocessing/">Google page</a> on Image Processing research, where the original, unmodified code can be downloaded. We would like to thank Yue for his efforts and wish him all the best in his research endeavors.

<em>References:</em>

<em>Chan, T. F., &amp; Vese, L. A. Active contours without edges. IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 10, NO. 2, FEBRUARY 2001</em>