Several Ways to Express a Convolution in Halide

There are at least two ways to express a convolution operation in Halide; more if the kernel is separable.  I'll review these in this post.

Let's examine a simple Gaussian 3x3 lowpass (smoothing) filter (also known as a Gaussian Blur):

The straight-forward method sums up the pixel neighborhood, using the weights in the convolution kernel.

Halide::Func gaussian_3x3_1(Halide::Func input) {

    Halide::Func k, gaussian;

    Halide::Var x,y,c;

   

    gaussian(x,y,c) = input(x-1, y-1, c) * 1   + input(x, y-1, c) * 2  + input(x+1, y-1, c) * 1 +

                                 input(x-1, y, c) * 2      + input(x, y, c) * 4     + input(x+1, y, c) * 2 +

                                 input(x-1, y+1, c) * 1  + input(x, y+1, c) * 2 + input(x+1, y+1, c) * 1;

    gaussian(x,y,c) /= 16;

    return gaussian;

}

We have to watch out not to overflow the summation of the pixel values.  In the gaussian_3x3_1 function below, the type of the summation (gaussian(x,y,c)) is deduced from the type of input(x,y,c) and if this is an 8-bit type for example, then it will most likely overflow and the output will be wrong without emitting any errors.  One way to handle this is to explicitly set the kernel weights to floats, but this will most likely hurt performance because it will require a cast and arithmetic operations on the float type.

gaussian(x,y,c) = input(x-1, y-1, c) * 1.f + input(x, y-1, c) * 2.f + input(x+1, y-1, c) * 1.f +

                  input(x-1, y, c) * 2.f   + input(x, y, c) * 4.f   + input(x+1, y, c) * 2.f +

                  input(x-1, y+1, c) * 1.f + input(x, y+1, c) * 2.f + input(x+1, y+1, c) * 1.f;

I prefer to cast the input type so that the caller of the function has the control on when to cast and when not to cast.  Here's an example which loads a 24-bit RGB image (8 bit per color channel), clamps the image values and converts from uint8_t to int32_t.

Halide::Image input = load("images/rgb.png");

Halide::Func padded, padded32;

padded(x,y,c) = input(clamp(x, 0, input.width()-1), clamp(y, 0, input.height()-1), c);

padded32(x,y,c) = Halide::cast(padded(x,y,c));

Halide::Func gaussian_3x3_fn = gaussian_3x3_1(padded32);

Another method to perform that convolution uses a 2-dimensional reduction domain for the convolution kernel.  We define a 3x3 RDom which spans from -1 to +1 in both width and height. When the RDom in the code below is evaluated, r.x takes the values {-1, 0, 1} and r.y similarly takes the values {-1, 0, 1}. Therefore, x+r.x takes the values {x-1, x, x+1}.  

Halide::Func gaussian_3x3_2(Halide::Func input) {

    Halide::Func k, gaussian;

    Halide::RDom r(-1,3,-1,3);

    Halide::Var x,y,c;

    

    k(x,y) = 0;

    k(-1,-1) = 1;    k(0,-1) = 2;    k(1,-1) = 1;

    k(-1, 0) = 2;    k(0, 0) = 4;    k(1, 0) = 2;

    k(-1, 1) = 1;    k(0, 1) = 2;    k(1, 1) = 1;

    gaussian(x,y,c) = sum(input(x+r.x, y+r.y, c) * k(r.x, r.y));

    gaussian(x,y,c) /= 16;

    return gaussian;

}

Because a Gaussian kernel is separable (that is, it can be expressed as the outer product of two vectors), we can express it in yet another way:

Halide::Func gaussian_3x3_3(Halide::Func input) {

    Halide::Func gaussian_x, gaussian_y;

    Halide::Var x,y,c;

    

    gaussian_x(x,y,c) = (input(x-1,y,c) + input(x,y,c) * 2 + input(x+1,y,c)) / 4;

    gaussian_y(x,y,c) = (gaussian_x(x,y-1,c)  + gaussian_x(x,y,c) * 2 + gaussian_x(x,y+1,c) ) / 4;

    return gaussian_y;

}

Of course, we can also use a reduction domain here.  In this case we need a 1-dimensional RDom spanning {-1, 0, 1}:

Halide::Func gaussian_3x3_4(Halide::Func input) {

    Halide::Func k, gaussian_x, gaussian_y;

    Halide::Var x,y,c;

    Halide::RDom r(-1,3);

    k(x) = 0;

    k(-1) = 1;    k(0) = 2;    k(1) = 1;

    gaussian_x(x,y,c) = sum(input(x+r.x, y, c) * k(r)) / 4;

    gaussian_y(x,y,c) = sum(gaussian_x(x, y+r, c) * k(r)) / 4;

    return gaussian_y;

}

And we can also play a bit with the syntax, replacing the Halide::Sum function with explicit summation over the reduction domain:

Halide::Func gaussian_3x3_5(Halide::Func input) {

    Halide::Func k, gaussian_x, gaussian_y;

    Halide::Var x,y,c;

    Halide::RDom r(-1,3);

    k(x) = 0;

    k(-1) = 1;    k(0) = 2;    k(1) = 1;

    gaussian_x(x,y,c) += input(x+r.x, y, c) * k(r) / 4;

    gaussian_y(x,y,c) += gaussian_x(x, y+r, c) * k(r) / 4;

    return gaussian_y;

}

So does it matter how we specify the algorithm?  The premise of Halide says 'no': write the algorithm once, and then experiment with different schedules until you get the best performance.  Intuitively, gaussian_3x3_2 is better than gaussian_3x3_1 because the Halide::RDom should have been optimized by Halide's compiler.  And gaussian_3x3_3 should perform better than gaussian_3x3_2 because it provides another degree of freedom when scheduling.  However, this is intuition and what we care about is actual performance measurements. 

I haven't measured this yet, so I owe you the results soon...  ;-)

Halide Excursions

I finally took the time to start a Halide-based github project, called halide-excursions.
I'm new to computer imaging and vision, and the algorithms and applications of this domain are a new frontier for my curiosity.  The halide-excursions project is an attempt to create a large, open source repository of Halide computer-vision, computational-photography, and image processing functions.  Anyone and everyone is more than welcome to join.

Halide is a new language for image processing and computer vision.  It was developed by several PhD students at MIT and it is actively maintained. There's a developer mailing list with less than a handful of messages a day (so it is easy to lazily eavesdrop on the conversation) where you can communicate directly with the Halide developers.  The response time is very short and the guys on the other end are very nice and eager to help.

Halide's succinct, functional syntax is very appealing and perhaps this is what drew me in.  If, like me, you are starting with little knowledge in the domain, and the word 'kernel' makes you think about the Linux kernel; and 'convolution' is a long forgotten concept from the university, then it might take a little more energy to get to smooth sailing.  But that's part of the fun.  Moving from Wikipedia to implementation is always a nice feeling and I find Halide a great platform to do this.  Here's the gradient magnitude of the Sobel operator (source code in halide-excursions):

The source code comes with a bunch of tutorials, sample "applications" and a large set of unit tests that can serve as a jumping board to start learning and messing around.  There's no user guide or official language specification, but there is doxygen-generated documentation.  All in all, I think there's plenty of resources to get you started.

Halide is supported on several OSs and cores (incl. GPUs) and promises the same performance as hand-optimized native code, with less code lines,less mess and with portability.  Hand optimized code - using vectorization and parallelization intrinsics and a sort of other tricks - is hard to read, hard to maintain, hard/impossible to port and makes exploring the scheduling optimization space very slow.  Halide's ability to separate the algorithm from the scheduling policy is very appealing and works well, most of the time.  For example, when implementing the Viola-Jones face detection algorithm, I found that implementing the classifier cascade phase in Halide cannot be done optimally because of Halide's poor support for control code. In a future post I'll provide more examples showing where Halide shines, and where a hybrid native-Halide solution works better.
Until then, I hope you join the project.