GraphTech Expert Blog

Today’s main stream CPU’s have SIMD (Single Instruction Multiple Data) compute capabilities. Intel has SSE, with AMD its 3DNow and ARM has NEON.

However, it seems, from my very narrow point of view, that these instruction set extensions do not nearly get the attention that they deserve in terms of opportunities to accelerate computational code. This post is about exploring and demonstrating the use of these capabilities using Intel’s SSE instruction set.

As a learning exercise, we chose to use a simple pixel alpha composition kernel. The nice thing about this kernel is that it is very simple on one hand, and has visual representation on the other. Our test function takes a RGBA ’source’ and a RGB ‘destination’ images and blends the source image into the destination image. It uses the following expression to calculate the composed pixel color:

Pdst = (P_src * α_src) + P_dst * (1 – α_src)

Where:

P_src – Pixel color of the source pixel

α_src – Alpha value of the source pixel

P_dst – Pixel color of the destination pixel

Data layout:

The first step with designing SIMD code, is to get your data arranged in a way that is SIMD friendly. In our case this is fairly straight forward – we represent our images one channel after the other.  An RGBA image looks like:

Figure 1 Image data layout

The component data element, in our case, is unsigned 8 bit.

Baseline Implementation:

As a base line, we used a simple composition loop, leaving the compiler to do it’s tricks. The code (for a single line) looks like:

for (unsigned int x = X0; x < X1; x++) {
	short diff;
	short tmp;

	diff = *pSrc[0] - *pDst[0];
	tmp = short(*pSrc[3] * diff) >> 8;
	*pDst[0] = tmp + *pDst[0];

	diff = *pSrc[1] - *pDst[1];
	tmp = short(*pSrc[3] * diff) >> 8;
	*pDst[1] = tmp + *pDst[1];

	diff = *pSrc[2] - *pDst[2];
	tmp = short(*pSrc[3] * diff) >> 8;
	*pDst[2] = tmp + *pDst[2];

	pSrc[0] += 1;
	pSrc[1] += 1;
	pSrc[2] += 1;
	pSrc[3] += 1;

	pDst[0] += 1;
	pDst[1] += 1;
	pDst[2] += 1;
}

The observant reader would notice that this code actually implements an alternative variant of the blending equation presented above. This modified form is:

Pdst = α_{src *  (}P_src - P_dst) + P_dst

The advantage of this form is that it uses a single multiply operation. However, I never really verified that this version is more efficient than the original form.

Running this function on my laptop which has:
Intel Core 2 Duo SP9300 @ 2.26 Mhz, 2 * 32KB L1, 6 MB L2

The compiler used is Microsoft Visual Studio 2008

We get the following results:

The column that is interesting is the pixel time, which is the pixel processing time in nanoseconds.

The ‘line overhead’ is calculated using the following expression:

L = W1 * W2 * (R1 – R2) / (W2 – W1)

Where:

L – Line overhead time

Wn – line width in test case n

Rn – pixel_time for test case n

It is left for the interested (and bored) reader to figure out why this expression represents the per-line calculation time.

SSE Implementation using Microsoft’s Visual Studio 2008 SSE Intrinsics

The next step in our journey was to implement our blending kernel  using Visual Studio SSE intrinsics. In this case, the blending function processes 16 pixels per iteration using the processor’s 128 bit vector processing capabilities. The result code is:

for (unsigned x = X0; x < X1; x += 16) {
	register __m128i s0, s1, d0, d1, a0, a1, r0, r1, zero;
	register __m128i diff0, tmp0, diff1, tmp1, t;
	zero = _mm_setzero_si128();
	// load alpha
	t = _mm_loadu_si128((__m128i *) pSrc[3]);
	a0 = _mm_unpacklo_epi8(t, zero);
	a1 = _mm_unpackhi_epi8(t, zero);

	EMMX_BLEND(0);
	EMMX_BLEND(1);
	EMMX_BLEND(2);

	pSrc[0] += 16;
	pSrc[1] += 16;
	pSrc[2] += 16;
	pSrc[3] += 16;

	pDst[0] += 16;
	pDst[1] += 16;
	pDst[2] += 16;
}

Where EMMX_BLEND is a macro:

#define EMMX_BLEND(comp)
	t = _mm_loadu_si128((__m128i *) pDst[(comp)]);
	d0 = _mm_unpacklo_epi8(t, zero);
	d1 = _mm_unpackhi_epi8(t, zero);
	t = _mm_loadu_si128((__m128i *) pSrc[(comp)]);
	s0 = _mm_unpacklo_epi8(t, zero);
	s1 = _mm_unpackhi_epi8(t, zero);
	/* A * S */
	tmp0 = _mm_mullo_epi16(s0, a0);
	tmp1 = _mm_mullo_epi16(s1, a1);
	/* 255 - A    */
	diff0 = _mm_sub_epi16(ff, a0);
	diff1 = _mm_sub_epi16(ff, a1);
	/* (255 - A) * D */
	diff0 = _mm_mullo_epi16(diff0, d0);
	diff1 = _mm_mullo_epi16(diff1, d1);
	/* r = A * S + (255 - A) * D */
	r0 = _mm_add_epi16(tmp0, diff0);
	r1 = _mm_add_epi16(tmp1, diff1);
	/* shift; */
	r0 = _mm_srli_epi16(r0, 8);
	r1 = _mm_srli_epi16(r1, 8);
	/* and pack */
	t = _mm_packus_epi16(r0, r1);
	_mm_storeu_si128((__m128i *) pDst[(comp)], t)

The results for this test case are:

As can be observed, the results are somewhat disappointing. A speed up of 160%-170% is not much, given that you have 16 times more execution units. A close look at the assembler code that is generated by Visual Studio reveals the culprit. It seems that the Visual studio compiler uses only 2 registers our of the 8 SSE registers that are available in the processor. After scratching my head for while, I thought that I should be more explicit in terms of how one calculation result should be re-used for the next calculation. The downside, of course, is that it requires writing the code in a totally unreadable form (or, also known as, write-only-code). The result EMMX_BLEND macro became:

#define EMMX_BLEND(comp)
	t = _mm_loadu_si128((__m128i *) pDst[(comp)]);
	d0 = _mm_unpacklo_epi8(t, zero);
	d1 = _mm_unpackhi_epi8(t, zero);
	t = _mm_loadu_si128((__m128i *) pSrc[(comp)]);
	_mm_storeu_si128((__m128i *) pDst[(comp)],
            _mm_packus_epi16(_mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(_mm_unpacklo_epi8(t, zero), a0),
            _mm_mullo_epi16(_mm_sub_epi16(ff, a0), d0)), 8),
            _mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(_mm_unpackhi_epi8(t, zero), a1),
            _mm_mullo_epi16(_mm_sub_epi16(ff, a1), d1)), 8)))

And the results:

An extra ~30% over the readable version which translates into a total of up to 200% speedup over the baseline. Not quite there yet!

So, if the brute force method doesn’t work, perhaps we should be even more brutal – time for assembler. At first, I was a little hesitant about writing the code in assembly, mainly when it comes to debugging the code using the Visual Studio debugger. However, as it turned out, it wasn’t that hard and perhaps even a lesson to bare in mind. Besides the ‘mechanical’ translation of the intrinsics into their assembly syntax, all I was left to deal with was the register allocation. And the end result code looks like this:

__asm {
	// initalization & load Alpha
	pxor xmm0, xmm0 // xmm0 <- 0
		mov eax, dword ptr [pSrc + 12]
		movdqu xmm1, [eax]; xmm1 <- *pSrc[3]
		movdqa xmm2, xmm1;
	punpcklbw xmm2, xmm0; // xmm2 <- a0, 16bit
	movdqa xmm3, xmm1;
	punpckhbw xmm3, xmm0; // xxm3 <- a1, 16bit

	// blending the red;
	mov eax, dword ptr[pDst + 0];
	movdqu xmm1, [eax]; // xmm1 = pDst[0]
	movdqa xmm6, xmm1;
	punpcklbw xmm6, xmm0; // xmm6 <- pDst[0] low 16bit
	movdqa xmm7, xmm1;
	punpckhbw xmm7, xmm0; // xmm7 <- pDst[0] high, 16 bit
	// load the ff constant
	movdqu xmm4, [ffconst]; // xmm4 <- ff
	movdqa xmm5, xmm4;
	psubw  xmm5, xmm2; // xmm5 = ff - a0
	pmullw xmm6, xmm5; // xmm6 = (ff - a0) * d0;
	// now for the upper bits
	movdqa xmm5, xmm4;
	psubw  xmm5, xmm3; // xmm5 = ff - a1
	pmullw xmm7, xmm5; // xmm7 = (ff - a1) * d1;
	// load the source;
	mov eax, dword ptr[pSrc + 0];
	movdqu xmm1, [eax]; // xmm1 = pSrc[0]
	// low bits of pSrc[0]
	movdqa xmm5, xmm1;
	punpcklbw xmm5, xmm0; // xmm5 = pSrc[0], low, 16 bit;
	pmullw xmm5, xmm2; // xmm5 = s0 * a0;
	paddw xmm6, xmm5; // xmm6 = s0 * a0 + (ff - a0) * d0;
	// high bits of pSrc[0]
	movdqa xmm5, xmm1;
	punpckhbw xmm5, xmm0;
	pmullw xmm5, xmm3; // xmm5 = s1 * a1
	paddw xmm7, xmm5; // xmm7 = s1 * a1 + (ff - a1) * d1;
	// shift the results;
	psrlw xmm6, 8;
	psrlw xmm7, 8;
	// pack back
	packuswb xmm6, xmm7; // xmm6 <- xmm6{}xmm7 low bits;
	mov eax, dword ptr [pDst + 0];
	movdqu [eax], xmm6; // done for this component;

	// Similar code goes for the green and blue channles
	// see the code for details

};

As for the results:

Finally it seems that we are getting close. The improvement over the base line is 7.3 X. This is not the desired 16X yet, but given the 4 evenings that I have put into this effort, this is probably at the point where I should call it version 1 and post something!

The Test program & the code

The test program (img_ops.exe) that used for this project is available *here*. It has two running options:

Benchmark mode – Is used to measure the performance with various image sizes.

Syntax: img_ops.exe benchmark

Bubble Tank demo – Demonstrate the performance difference between the serial and the vectorized code. In this mode the blending of multiple instances of the overlay image are blended on the background images.

Syntax: img_ops.exe bubble_tank <background image> <overlay image>

Disclaimers, next steps and conclusions

Ok, there are many of them so lets get started.

First, you might have noticed that the SIMD implementation assumes that source image line width is a multiply of 16 bit. Obviously this restriction is not realistic and should be relieved by adding wind up/down code to handle the head and the tail of a line.

Next, As can be observed from the results above, the ‘per-line’ overhead is quite considerable. With the fastest version, it seems that we add overhead of ~80 pixels per line. Obviously this should be reviewed for potential optimization. In general, both the serial and the SIMD codes have not been reviewed and there is probably big potential for shaving a few cycles off the inner loops.

In terms of next steps (once I get to finish up this one) is to test this code on other architectures. Especially interesting one is ARM/Neon. The reason is that while there is a growing market desire for richer user interfaces, there are many ARM based embedded systems that can’t afford to include a GPU due to area/costs considerations. Such implementations could use the SIMD extensions to deliver a rich and compelling user experience.

To conclude this phase of the experiment, there are few takeaways:

First, as expected, the SIMD instruction set indeed offers great optimization potential. We have presented 7.3 X performance boost, with strong feeling that there is plenty of room for more.

Second, SIMD code, especially written in assembly language and using Visual Studio, is not debugging friendly. So, you better have a solid, debugged serial or vectorized implementation before you start your assembly programming.

This post will discuss some of the theory of image shift detection. What’s shift detection for images and what is it good for?

I think it’s best to explain by example, which will be our reference for the rest of the post. We’ll refer to a typical server-client distributed visualisation environment, in which the server runs the applications, the image is grabbed from the server, processed and encoded if needed, and then send to the client, where it’s decoded and displayed. The time consuming parts are server & client processing and network limitations. See this post for more details.

Now lets consider a scenario where the server runs an application like a browser or image gallery, and the user scrolls down the screen, so the entire screen is updated, but actually only small part of it is getting changed, most of it is just being shifted upwards. Another example could be a map application with a small blinking location mark on it, and the application updates the entire map every time the small mark blinks or is changed. In these examples it seems intuitive to save some processing time and bandwidth by “reusing” the images from previous frames. This is where we could use shift detection.
Seeing as some parts of the frame are changed while others aren’t, we need to divide the full frame to tiles. Tiling usually causes some overhead so the size of the tiles should be optimised to the environment. If for example we’re using libjpeg to encode frames, there is overhead for each tile, therefore we don’t want to make them too small. The essence of shift detection is trying to find tiles from the current frame in the previous frames (possibly with some shift) and then reusing them to save time. In order to do this, we need copies of the previous and current frames on the server and client. We need to allocate memory for these and maintain them correctly and consistently (the server and client copies should be identical).
The shift detection process has 3 main tasks:
1. Create an anchor – a reference point for matching.
2. Find shift by matching the anchor.
3. Match tiles, based on the found shift.

Note: the shift detection discussed here only handles vertical and horizontal shifting. Random direction shift detection is much more expensive, therefore not always useful.

1. Creating an anchor

We need an efficient way to detect potential shift. Comparing the entire tile is very slow, so we need to define some unique fragments which will give us a good reference point. For example, we could use “dirty” stripes, stripes with lots of pixels changes in them. An anchor could be defined by N dirty stripes of length L. The anchor shouldn’t be too simple – it should represent a unique area in the frame, otherwise we’re prone to mismatches – where the anchor is matched, but tiles actually aren’t. We also want the anchor to be somewhere close to the center of the frame, to insure, in good probability, it’s not scrolled out. Creating an anchor could be integrated within the encoding, in order to save time. For example, for every tile we encode we could try create a potential anchor, by counting different pixels and creating dirty stripes from them (we don’t need to save stripes data, only the coordinates). If we managed to create a potential anchor in the encoding process, we can compare it to the current anchor, if exists, and replace it if it’s closer to center.
2. Shift matching
When process a new frame, first of all we want to see if it’s possible to reuse data from previous frames. If we have a valid anchor from previous frames, it might be possible. We want to see where the anchor from the previous frame is located in our current frame, thus calculate the shift from the previous frame. As mentioned, we’re only handling vertical & horizontal shifts.

Note: Zero shift, is also a shift – therefore “diff” between two frames is also covered by shift detection.

So if our anchor is N stripes, we can start by matching the first stripe, vertically or horizontally (we could save previous shift direction as hint). If we matched the first stripe, we can then use the found shift to match the rest of the anchor’s stripes. If we matched the entire anchor, then we can declare we’ve found a shift, otherwise we have no match, and we encode the frame normally.

3. Tile matching
Now we’ve detected a potential shift, we still need to make sure the shifted tiles in the current frame are matched with the previous frame, because even if we matched the anchor, the tiles still can have some changes between them (like mouse cursor, or a letter in an editor application). So now, for every tile in the current frame, we try to match an area in the previous frame, based on the shift we’ve detected. We have to go over all the pixels and compare them to make sure we can match the entire tile. If we have a match, then instead of sending tile data to the client, we can send only the shift coordinates, and the client can reproduce the tile from the previous frame copy. If shift is zero, then we don’t have to send anything, because the area is unchanged. If we don’t have a match, then we process the tile normally and send the new image data to the client.

Performance wise, the shift detection takes a least one full pass over all pixels per tile ( the shift matching is cheap, because of the vertical & horizontal restriction ), but it’s still significantly less than most of the image compression algorithms, and if our target applications matches the “shift” scenario, it could be a very effective strategy. We can also save a lot of bandwidth, by sending only shift coordinates instead of full image data, even in compressed form.