Convolutions

[mg262]

So... I have some SIMD code to perform general 4x4 separable convolutions as part of another filter, and I was thinking of pulling it out into a separate plug-in. A couple of questions:

Does anyone actually have any use for 4x4 (as opposed e.g. 3x3)?
What edge behaviour is wanted (shrink frame, mirror, assume 0)?
What chroma behaviour is wanted (discard, copy, convolve)?

These are particularly targeted at those writing large scripts...

[Manao]

How's the 4x4 kernel centered on a pixel ? Basically, it amounts to a 5x5 kernel with constraints on the kernel coefficients, doesn't it ?

I guess for script gurus, and for chroma processing, the more versatile, the better ( hence, discard, copy & convolve )

For the edge behavior, i'd do padding, perhaps mirror, and definitely not 0. What is shrink frame ?

[mg262]

The way I phrased things was probably coloured by the fact that I'm currently using it for subpixel translations in a stabiliser... the notion of "centered on a pixel" doesn't come in in this situation; you run the 4x4 window across the frame, with each window position giving an output pixel. Currently the output frame is smaller than the input frame by three pixels each way (which is what I meant by shrink frame) Edit: so border issues don't come in at all. I rather doubted 4x4 was useful for anything else, but I thought I'd check, since the code is already there...

By padding, do you mean repeating the border value outwards?

I can think of one more edge method, which is to truncate the kernel and scale the remaining values up to normalise. (Obviously this is only appropriate for some kernels.)

Multiple border treatments will probably require a separate version for each one (unless I wait for the AVISynth padding support, but initialising the padding can have a substantial time penalty for reasons you pointed out) -- so I wanted to check what might actually be used.

[IanB]

@mg262,

Paste the guts of the code here, those with mild interest (like me) can browse the algo, those more keen can wrap it appropriate to their need.

IanB

[mg262]

With assembly, I would generally like to squeeze out as much speed as I can before posting -- and then when someone finds a way to speed it up, I learn something. In this case I may not have time to work on it for a while... so here it is. Please bear in mind that this is not a finished product. (BTW, this kind of thing: movq mm0, [esi]/movq2dq xmm0, mm0 really should be replaced by a single instruction... I just haven't had time to go through and do it.)

start:
//eightpixelshorizontal:
pxor xmm7, xmm7
movdqa xmm5, short_xkernel

//first load 11 consecutive pixels of NEW data from the source (esi), needed for eight (overlapping) 4x4 convolutions
movq mm0, [esi]
movq mm1, [esi+1]
movq mm2, [esi+2]
movq mm3, [esi+3]
movq2dq xmm0, mm0
movq2dq xmm1, mm1
movq2dq xmm2, mm2
movq2dq xmm3, mm3
punpcklbw xmm0, xmm7
punpcklbw xmm1, xmm7
punpcklbw xmm2, xmm7
punpcklbw xmm3, xmm7

//now apply the horizontal convolutions
pmaddwd xmm0, xmm5
pmaddwd xmm1, xmm5
pmaddwd xmm2, xmm5
pmaddwd xmm3, xmm5

//now xmm0 holds the results for pixels 0 and 4
// xmm1 holds the results for pixels 1 and 5, etc

//remaining code interleaved not to save registers but to match the execution units/subunits
//empirically much faster than the alternative

//xmm0,1: combine the two DWORD outputs from pmaddwd
pshufd xmm4, xmm0, 0xb1
pshufd xmm6, xmm1, 0xb1
add esi, 8
add edx, 64//lea edx, [ edx+64]
paddd xmm0, xmm4
paddd xmm1, xmm6
pshufd xmm4, xmm2, 0xb1
pshufd xmm6, xmm3, 0xb1

//xmm0 pixel results bytewise: ....4444/....0000, etc. (scaled by -32768)
psllq xmm0, 48-9
psllq xmm1, 48-9
//xmm0 pixel results bytewise: 44....../00......, etc. (scaled by -64)
//ignore rounding at this stage; empirically, luma error about 0.1
paddd xmm2, xmm4
paddd xmm3, xmm6

//store data in the 73 62 51 40 format

movdqa xmm7, mask //YY....../YY......

//load the convolution results for the three previous rows, b(-1),c(-2),d
movdqa xmm4, [edx] //bbccdd../bbccdd..
movdqa xmm5, [edx+16]
psllq xmm2, 48-9
psllq xmm3, 48-9

//mask data for this row to prevent overwriting data for previous rows
//n.b. this could be done with pshufw but that uses the wrong execution subunit
pand xmm0,xmm7 //now 44....../00...... for row a
pand xmm1,xmm7 //now 55....../11...... for row a

//shift the three previous rows' data backwards
psrlq xmm4, 16 //now ..bbccdd/..bbccdd
psrlq xmm5, 16

//combine the data for this row
paddd xmm0, xmm4
paddd xmm1, xmm5
movdqa xmm4, [edx+32]
movdqa xmm5, [edx+48]
//and put the data back
movdqa [edx], xmm0
movdqa [edx+16],xmm1

psrlq xmm4, 16
psrlq xmm5, 16
pand xmm2, xmm7
pand xmm3, xmm7
paddd xmm2,xmm4
paddd xmm3,xmm5
movdqa [edx+32],xmm2
movdqa [edx+48],xmm3

//we have now loaded eight pixels worth of intermediate (horizontally convolved) data
//
//xmm0 contains data for pixels 0 and 4, 44444444/00000000
//xmm1 contains data for pixels 1 and 5, 55555555/11111111
//all four registers are in the format aabbccdd/aabbccdd, where a is the LOWEST row
//interleaving code following this point with earlier code has no empirical benefit
__asm movdqa xmm5, short_ykernel
__asm movdqa xmm6, rounder


//convolve, add rounder and combine the two DWORD outputs of pmaddwd
pmaddwd xmm0, xmm5
pmaddwd xmm1, xmm5
pmaddwd xmm2, xmm5
pmaddwd xmm3, xmm5
paddd xmm0, xmm6
paddd xmm1, xmm6
paddd xmm2, xmm6
paddd xmm3, xmm6
pshufd xmm4, xmm0, 0xb1
pshufd xmm5, xmm1, 0xb1
pshufd xmm6, xmm2, 0xb1
pshufd xmm7, xmm3, 0xb1
paddd xmm0, xmm4
paddd xmm1, xmm5
paddd xmm2, xmm6
paddd xmm3, xmm7
//xmm0 pixel results bytewise: ....4444/....0000, etc. (scaled by -32768*-64)
//xmm0 pixel results bytewise: ....44../....00.., etc. (scaled by 32)

movdqa xmm5, [lowcap] //lowcap; set to zero for all words except ......xx/......xx
movdqa xmm6, [highcap] //highcap; set to zero for all words except ....xx../....xx..
pxor xmm7, xmm7

psrad xmm2, 5+16
psrad xmm3, 5
psrad xmm0, 5+16
psrad xmm1, 5
//xmm0 pixel results bytewise: ......44/......00,
//xmm1 pixel results bytewise: ....55../....11..,
//xmm2 pixel results bytewise: ......66/......22,
//xmm3 pixel results bytewise: ....77../....33..

pmaxsw xmm2, xmm7
pmaxsw xmm3, xmm7
pmaxsw xmm0, xmm7
pmaxsw xmm1, xmm7
pminsw xmm2, xmm5
pminsw xmm3, xmm6
pminsw xmm0, xmm5
pminsw xmm1, xmm6
//each result is now at most 255, so lies in a single byte

paddb xmm2, xmm3
paddb xmm0, xmm1
movq mm6, [edi]
psllq xmm2,32
movq2dq xmm6, mm6 //xmm6 = xxxxxxxx/76543210 of compare
paddb xmm0, xmm2
pshufd xmm5, xmm6, 0xe4 //copy

//xmm0 pixel results bytewise: 77665544/33221100
packuswb xmm0, xmm0
//xmm0 pixel results bytewise: xxxxxxxx/76543210

pminub xmm6, xmm0
pmaxub xmm0, xmm5
add edi, 8
psubusb xmm0, xmm6
//now the differences are stored in the lower eight bytes of xmm0
//this is the place to AND with the main mask
//the metric lookup happens at the beginning of the loop
cmp esi, [sourceend]//moving this up causes slowdown (!)
jne eightpixels
//END OF INNER LOOP

I personally find that a lot less enlightening than an English description... but there you go ;)