// Fun with 3x3 box filters on SSE2.
// Written by Nils Liaaen Corneliusen 2012.
// License: CC0 1.0 Universal (CC0 1.0) Public Domain Dedication license
// https://www.ignorantus.com

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

/* "Clean C", according to those guys
void blur(const Image &in, Image &blurred) {
 Image tmp(in.width(), in.height());

 for (int y = 0; y < in.height(); y++)
  for (int x = 0; x < in.width(); x++)
   tmp(x, y) = (in(x-1, y) + in(x, y) + in(x+1, y))/3;

 for (int y = 0; y < in.height(); y++)
  for (int x = 0; x < in.width(); x++)
   blurred(x, y) = (tmp(x, y-1) + tmp(x, y) + tmp(x, y+1))/3;
}
*/

/* "Fast C++" by the same guys. I'm not making this up. It actually looked like that.
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);
}}}}}
*/

// First rewrite by me. It's amazing: Spaces can be used to improve readability!

void fast_blur_halide( const uint16_t *in, uint16_t *blurred, int width, int height )
{
    int x, y;
    int xTile, yTile;
    __m128i one_third = _mm_set1_epi16(21846);

    for( yTile = 0; yTile < height; yTile += 32 ) {

        __m128i a, b, c;
        __m128i sum, avg;
        __m128i tmp[(256/8)*(32+2)];

        for( xTile = 0; xTile < width; xTile += 256 ) {

            __m128i *tmpPtr = tmp;

            for( y = -1; y < 32+1; y++ ) {

                const uint16_t *inPtr = &in[(yTile+y)*width + xTile];

                for (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( y = 0; y < 32; y++ ) {

                __m128i *outPtr = (__m128i *)&blurred[(yTile+y)*width + xTile];

                for( 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 );
                }

            }
        }
    }
}

// But assuming you're not using a 15 year old computer, this is much simpler. Also faster.

void fast_blur_horiz2d( const uint16_t *in, uint16_t *blurred, int width, int height )
{
    int x, y;
    __m128i one_third;
    __m128i *dst0, *dst1;
    
    one_third = _mm_set1_epi16(21846);
    dst0 = (__m128i *)blurred;
    dst1 = (__m128i *)(blurred+width);

    for( y = 0; y < height; y += 2 ) {
        const uint16_t *row0, *row1, *row2, *row3;

        row1 = in + y*width;
        row0 = row1 - width;
        row2 = row1 + width;
        row3 = row2 + width;

        for( x = 0; x < width; x += 8 ) {
            __m128i s00 ,s01 ,s02;
            __m128i r00, r01, r02;
            
            s00 = _mm_loadu_si128( (__m128i*)(row0-1) );
            s01 = _mm_loadu_si128( (__m128i*)(row0+1) );
            s02 = _mm_load_si128(  (__m128i*)(row0) );
            r00 = _mm_mulhi_epi16( _mm_add_epi16( _mm_add_epi16( s00, s01 ), s02 ), one_third );

            s00 = _mm_loadu_si128( (__m128i*)(row1-1) );
            s01 = _mm_loadu_si128( (__m128i*)(row1+1) );
            s02 = _mm_load_si128(  (__m128i*)(row1) );
            r01 = _mm_mulhi_epi16( _mm_add_epi16( _mm_add_epi16( s00, s01 ), s02 ), one_third );

            s00 = _mm_loadu_si128( (__m128i*)(row2-1) );
            s01 = _mm_loadu_si128( (__m128i*)(row2+1) );
            s02 = _mm_load_si128(  (__m128i*)(row2) );
            r02 = _mm_mulhi_epi16( _mm_add_epi16( _mm_add_epi16( s00, s01 ), s02 ), one_third );

            _mm_store_si128( dst0++, _mm_mulhi_epi16( _mm_add_epi16( _mm_add_epi16( r00, r01 ), r02 ), one_third ) );

            s00 = _mm_loadu_si128( (__m128i*)(row3-1) );
            s01 = _mm_loadu_si128( (__m128i*)(row3+1) );
            s02 = _mm_load_si128(  (__m128i*)(row3) );
            r00 = _mm_mulhi_epi16( _mm_add_epi16( _mm_add_epi16( s00, s01 ), s02 ), one_third );

            _mm_store_si128( dst1++, _mm_mulhi_epi16( _mm_add_epi16( _mm_add_epi16( r00, r01 ), r02 ), one_third ) );

            row0 += 8;
            row1 += 8;
            row2 += 8;
            row3 += 8;
        }

        dst0 += width>>3;
        dst1 += width>>3;

    }

}