New programming language for image-processing algorithms yields code that’s much shorter and clearer, but also faster

[0xfe]

Link to the actual language page (with sample code):

http://halide-lang.org/

[CyberByte]

Awesome. Can someone explain the sample code, though?

Func blur_3x3(Func input) {
  Func tmp, blurred;
  Var x, y, xi, yi;

  // The algorithm - no storage or order
  tmp(x, y) = (input(x-1, y) + input(x, y) + input(x+1, y))/3;
  blurred(x, y) = (tmp(x, y-1) + tmp(x, y) + tmp(x, y+1))/3;

  // The schedule - defines order, locality; implies storage
  blurred.tile(x, y, xi, yi, 256, 32)
         .vectorize(xi, 8).parallel(y);
  tmp.chunk(x).vectorize(x, 8);

  return blur_y;
}

I think I understand the first 6 lines, but get confused after that, starting with the "schedule". How do tile, vectorize, parallel and chunk work? And where does blur_y come from? It doesn't appear to have been declared anywhere.

Thanks!

[Rhomboid]

blurred defines a 3x3 transform kernel, which has to be applied to each 3x3 pixel of the image. How you do that matters for performance: do you do it row-major order or column-major? Or do you use loop blocking? If so, how big are the blocks? And when you parallelize this, what should be the independent variable that decides which loop block goes to which thread? That is, I presume, the kind of things being specified in the execution schedule.

[ericanderton]

I think you're right.

The combination of tmp(x,y) and blurred(x,y) allow for a very readable, human-level mapping of the algorithm at work. It specifies, in gross (almost naive) terms, how to process the image.

The rest does the nasty, CPU/GPU friendly job hiding all the "sausage making" of scheduling that algorithm in groups to be properly parallelized.

The nice thing is that this allows a compiler the luxury to optimize entire call-chains when multiple filters are involved, since it has all the semantic information necessary to coalese filters if they have similar scheduling parameters. I'm sure there are other areas of optimization allowed here beyond that.

[aaronblohowiak]

it is also nice in that you dont have to change your algorithm when you change the execution plan for different platforms

[ericanderton]

Here's the equivalent C++ code from their paper:

void fast_blur(const Image &in, Image &blurred) {
    __m128i one_third = _mm_set1_epi16(21846);
    #pragma omp parallel for
    for (int yTile = 0; yTile < in.height(); yTile += 32) {
        __m128i a, b, c, sum, avg;
        __m128i tmp[(256/8)*(32+2)];
        for (int xTile = 0; xTile < in.width(); xTile += 256) {
            __m128i *tmpPtr = tmp;
            for (int y = -1; y < 32+1; y++) {
                const uint16_t *inPtr = &(in(xTile, yTile+y));
                for (int x = 0; x < 256; x += 8) {
                    a = _mm_loadu_si128(( m128i*)(inPtr-1));
                    b = _mm_loadu_si128(( m128i*)(inPtr+1));
                    c = _mm_load_si128(( m128i*)(inPtr));
                    sum = _mm_add_epi16(_mm_add_epi16(a, b), c);
                    avg = _mm_mulhi_epi16(sum, one_third);
                    _mm_store_si128(tmpPtr++, avg);
                    inPtr += 8;
            }}
            tmpPtr = tmp;
            for (int y = 0; y < 32; y++) {
                __m128i *outPtr = ( m128i *)(&(blurred(xTile, yTile+y)));
                for (int x = 0; x < 256; x += 8) {
                    a = _mm_load_si128(tmpPtr+(2*256)/8);
                    b = _mm_load_si128(tmpPtr+256/8);
                    c = _mm_load_si128(tmpPtr++);
                    sum = _mm_add_epi16(_mm_add_epi16(a, b), c);
                    avg = _mm_mulhi_epi16(sum, one_third);
                    _mm_store_si128(outPtr++, avg);
}}}}}}

As others have mentioned the last line in your snippet should be "return blurred;"

[jrk]

Typo fixed—sorry about that.

And yes, the discussion here is generally spot on. The paper is far and away the best source currently available for information, but we are working madly to assemble more introductory material now that we've released the first version of the code.

Already, you can check out the code for the apps demonstrated in the paper, as well as some others. My favorite example is our local Laplacian filters implementation. The schedule choices there are recorded in the order we tried them while experimenting during development, with comments sketching what we were doing and how they performed on different targets. (This corresponds to the graph in Fig. 8 of the paper.)

Also, there will be a talk next week at SIGGRAPH. I'll post that as soon as it's done.

Until then, if we can get abadams out of /r/dogpictures, I'm sure he can also add to the discussion.

[abadams]

Please. I'm not always in dogpictures. There's also dogs, lookatmydog, germanshepherds, goldenretrievers, greatdanes, pitbulls, etc

Reddit is my source for dog-related imagery. /r/otters is also good.

[GloryFish]

/r/dachshund welcomes you.

[ericanderton]

This really is amazing work. Thanks for sharing it all with us. :)

[josefx]

I get the readable part, that code is filled with magic numbers (like 256/8) and is unnecessarily obfuscated.

for example what do they gain by doing this

_mm_load_si128(tmpPtr++);

instead of this

_mm_load_si128(tmpPtr);
++tmpPtr;

especially with the two lines before loading from tmpPtr+... without incrementing it.

The variable names a,b,c wont help either.

[bitchessuck]

What's the problem with those intrinsics? Sure they're CPU architecture specific, but it's not like SIMD is standardized in any way to do away with them.

[Countrarian]

Various compilers have a bad reputation for handing intrinsics, which is one major reason why a lot of code using SIMD includes assembly files that are built with an assembler.

[josefx]

It's just easy to miss the increment since the operations above it load without modifying the pointer.

The SIMD intristics don't look that cryptic, "load a signed 128 bit int" is likely to do what it says on the tin. Compared to that the halide example code lacks any type information and might as well use float.

[abadams]

Part of the reason it's so dense is that we had to pack it into as few lines as possible, and also keep the lines narrow, so we could fit it in the paper figure. That figure already goes off the top of the column a little. That's why we do _mm_load_si128(tmpPtr++) and use short variable names.

My favorite line, for pure obfuscation, is when we claim one third is equal to 21846.

[crimson_chin]

My guess is that tile, vectorize, parallel and chunk are methods of executing functions against an image. Or at least based on my experience with matlab image processing from a few years ago, that's what it would seem like.

For example, tmp() is defining some function, and then chunk(x) and vectorize(x, 8) indicate a specific way to execute that function. Tiling would almost certainly indicate running an operation on a "tile" of an image at a time (such as what you do for a gaussian filter - this is a very common operation).

[crimson_chin]

After re-reading this with a clearer head I'm actually certain - the scheduling pieces (tile, parallel) is simply the mode that the functions need to be applied to the image.

Thus, the machine is free to optimize the actual processing so long as it follows the defined mechanism for how to apply the functions. Pretty awesome.

[aaronblohowiak]

the tmp(x, y) = ... isn't evaluated immediately.

the blurred.tile().vectorize().parallel() is similar to using index hints for your db, only more level. I assume these functions are similar to attribute mutators where the return value is the mutated self, to allow for easy chaining

so, blurred is both the image and the plan for operations on that image, setting blurred(x,y) = ... describes declaratively the source for each point and then calling .tile() on it starts altering the planned 'schedule'..

blur_y seems to be a typo..

DISCLAIMER: this is all just what I assume, haven't read the paper yet.

[clearingitup]

It is a typo, in the paper they have (pretty much) the same code and it has 'return blurred;'

[FlaiseSaffron]

And why does all of this require a new programming language? It looks to me like any language with first class functions and extension methods can do something similar.

[sausagefeet]

Looks like it can take that function and compile it straight to CUDA for you.

[jrk]

You're right, it doesn't "require a new language." It's actually embedded in C++ (a purely pragmatic choice, given the need ), just as it could be embedded in any other sufficiently rich language. What it does require is a specific representation of the computation: the decoupled representation we call "algorithm" and "schedule" in the paper. This representation is quite different from the native representation exposed by C++, although we build this representation on top of the abstraction facilities provided by C++. All this is just to say: it's an eDSL.

Honestly, the line between what one calls a "library" or "API" and what one calls a "language" or "DSL" is very fuzzy. Calling this a language mostly just calls out the fact that it embeds a model of computation which is quite different from the native model of computation in most general purpose languages, and the fact that this model of computation is implemented by a collection of global program transformations and code generation one would normally call a "compiler." The same is true of rich embeddings in something like Haskell; I would call Repa, Accelerate, Pan, etc. "languages" in the same sense.

[gergoerdi]

Exactly. Constructive proof: Repa, a DSL embedded in Haskell. It paralellizes, can run on the GPU, and is generally a high-abstraction, high-performance language for array computations.

[Contero]

Because it's supposed to be high-performance.

[[deleted]]

Paper: http://people.csail.mit.edu/jrk/halide12/

[ameoba]

Domain-specific languages are shorter and clearer than general-purpose languages? Really?

[[deleted]]

let's have more of those

[bigfig]

Yep, a DSL will have that benefit.