SIMD implementation of Convolver for Lanczos filter etc.

skia/ext/convolver.cc

Index: skia/ext/convolver.cc
diff --git a/skia/ext/convolver.cc b/skia/ext/convolver.cc
index a42a9daf32a17a3a77602f89d841de997f2787de..b5841354756375c700a7ad74dc998fc8f99b7375 100644
--- a/skia/ext/convolver.cc
+++ b/skia/ext/convolver.cc
@@ -7,6 +7,10 @@
 #include "skia/ext/convolver.h"
 #include "third_party/skia/include/core/SkTypes.h"
 
+#ifdef ARCH_CPU_X86_FAMILY

[evannier, 2011/02/14 23:45:13]

I believe Chrome coding style recommends:
#if defined(XXX)

See also below, on why I think the ARCH_CPU_X86 might be the wrong define to
use.

[jiesun, 2011/02/17 20:17:58]

Done.

+#include <emmintrin.h>
+#endif
+
 namespace skia {
 
 namespace {
@@ -199,7 +203,7 @@ void ConvolveVertically(const ConvolutionFilter1D::Fixed* filter_values,
     if (has_alpha) {
       unsigned char alpha = ClampTo8(accum[3]);
 
-      // Make sure the alpha channel doesn't come out larger than any of the
+      // Make sure the alpha channel doesn't come out smaller than any of the
       // color channels. We use premultipled alpha channels, so this should
       // never happen, but rounding errors will cause this from time to time.
       // These "impossible" colors will cause overflows (and hence random pixel
@@ -219,6 +223,358 @@ void ConvolveVertically(const ConvolutionFilter1D::Fixed* filter_values,
   }
 }
 
+
+// Convolves horizontally along a single row. The row data is given in
+// |src_data| and continues for the num_values() of the filter.
+template<bool has_alpha>

[evannier, 2011/02/14 23:45:13]

has_alpha is not used anywhere, suggesting there is a bug here.

[jiesun, 2011/02/15 00:19:22]

actually, it is always vertical after horizontal pass. therefore horizontal
never worry about alpha. but I think you are right on this.

[jiesun, 2011/02/17 20:17:58]

Done.

+void ConvolveHorizontally_SSE2(const unsigned char* src_data,

[evannier, 2011/02/14 23:45:13]

Maybe move that closer to the the other Horizontally for easier comparison (or
not).

[jiesun, 2011/02/17 20:17:58]

I am grouping all SIMD function into a #ifdef. not sure which one is better.

+                          const ConvolutionFilter1D& filter,

[evannier, 2011/02/14 23:45:13]

here and elsewhere, indentation looks odd. And various lint errors.

[jiesun, 2011/02/17 20:17:58]

Done.

+                          unsigned char* out_row) {
+#ifdef ARCH_CPU_X86_FAMILY
+  int width = filter.num_values();

[evannier, 2011/02/14 23:45:13]

Maybe name this num_values for better symmetry with the non SIMD version.

[jiesun, 2011/02/17 20:17:58]

Done.

+
+  int filter_offset, filter_length;
+  __m128i zero = _mm_setzero_si128();
+  for (int i = 0; i < width; i += 1) {  // One pixel per iteration.

[evannier, 2011/02/14 23:45:13]

same comment about the naming.

[jiesun, 2011/02/17 20:17:58]

Done.

+    const ConvolutionFilter1D::Fixed* filter_values =
+        filter.FilterForValue(i, &filter_offset, &filter_length);
+
+    __m128i accum = _mm_setzero_si128();
+    const unsigned char* start = src_data+(filter_offset<<2);

[evannier, 2011/02/14 23:45:13]

what is wrong with the way it was written in the other code:
 const unsigned char* row_to_filter = &src_data[filter_offset * 4];

[evannier, 2011/02/14 23:45:13]

Since this is going to be used as a __m128i, I would do, const __m128i* start =
reinterpret_cast<const ...>(etc);
This avoids the casts below, and makes the start += 16, a start +=1, which might
be clearer ... To the compiler.

[jiesun, 2011/02/17 20:17:58]

Done.

+    // Four filter taps per iteration.
+    for (int j = 0; j < filter_length >> 2; ++j) {
+      // [16] xx xx xx xx c3 c2 c1 c0
+      __m128i coeff = _mm_loadl_epi64((__m128i*)filter_values);

[evannier, 2011/02/14 23:45:13]

Maybe add a comment explaining that you have ensured that this address is 64 bit
aligned (well, I assume you did).

[jiesun, 2011/02/15 00:19:22]

I think there is no alignment requirement for _mm_loadl_epi64. I did not align
filter values for each filter. just pad at the end of all filters to make sure
the 64 bit load always access valid addresses.

+      // [16] xx xx xx xx c1 c1 c0 c0
+      __m128i coeff16 = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(1, 1, 0, 0));

[evannier, 2011/02/14 23:45:13]

I have not looked at the resulting assembly, but I have been fairly disappointed
in GCC code generation when it comes to the use of intrisics, and it looks like
it might be possible that GCC might have to perform register spill over to the
stack. Check what it does, and maybe consider doing the other coeff16 shuffle
and unpack right there (it is the same number of registers, but it might help
GCC, or not). 

+      // [16] c1 c1 c1 c1 c0 c0 c0 c0
+      coeff16 = _mm_unpacklo_epi16(coeff16, coeff16);
+
+      // [8] a3 b3 g3 r3 a2 b2 g2 r2 a1 b1 g1 r1 a0 b0 g0 r0
+      __m128i src8 = _mm_loadu_si128((__m128i*)start);
+      // [16] a1 b1 g1 r1 a0 b0 g0 r0
+      __m128i src16 = _mm_unpacklo_epi8(src8, zero);
+      __m128i mul_hi = _mm_mulhi_epi16(src16, coeff16);

[evannier, 2011/02/14 23:45:13]

Because of the many stalls these operations might generate, I would interleave
both parts in one single pass. Not sure if this would require more registers. 
If you have already considered that and timed it, please comment on your choice.

+      __m128i mul_lo = _mm_mullo_epi16(src16, coeff16);
+      // [32]  a0*c0 b0*c0 g0*c0 r0*c0
+      __m128i t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+      // [32]  a1*c1 b1*c1 g1*c1 r1*c1
+      t = _mm_unpackhi_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+
+      // [16] xx xx xx xx c3 c3 c2 c2
+      coeff16 = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(3, 3, 2, 2));
+      // [16] c3 c3 c3 c3 c2 c2 c2 c2
+      coeff16 = _mm_unpacklo_epi16(coeff16, coeff16);
+      // [16] a3 g3 b3 r3 a2 g2 b2 r2
+      src16 = _mm_unpackhi_epi8(src8, zero);
+      mul_hi = _mm_mulhi_epi16(src16, coeff16);
+      mul_lo = _mm_mullo_epi16(src16, coeff16);
+      // [32]  a2*c2 b2*c2 g2*c2 r2*c2
+      t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+      // [32]  a3*c3 b3*c3 g3*c3 r3*c3
+      t = _mm_unpackhi_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+
+      start += 16;
+      filter_values += 4;
+    }
+
+    // remaining
+    int r = filter_length&3;
+    if (r) {
+      // Note: filter_values must be pad to align_up(filter_offset, 8);

[evannier, 2011/02/14 23:45:13]

padded

+      __m128i coeff = _mm_loadl_epi64((__m128i*)filter_values);
+      // Mask out extra filter taps.
+      __m128i mask = _mm_set_epi16(0, 0, 0, 0, 0, r==3?-1:0, r>=2?-1:0, -1);
+      coeff = _mm_and_si128(coeff, mask);
+      __m128i coeff16 = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(1, 1, 0, 0));
+      coeff16 = _mm_unpacklo_epi16(coeff16, coeff16);
+
+      // TODO(jiesun): line buffer must be pad to align_up(filter_offset, 16);

[evannier, 2011/02/14 23:45:13]

padded

[jiesun, 2011/02/17 20:17:58]

Done.

+      __m128i src8 = _mm_loadu_si128((__m128i*)start);
+      __m128i src16 = _mm_unpacklo_epi8(src8, zero);
+      __m128i mul_hi = _mm_mulhi_epi16(src16, coeff16);
+      __m128i mul_lo = _mm_mullo_epi16(src16, coeff16);
+      __m128i t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+      t = _mm_unpackhi_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+
+      src16 = _mm_unpackhi_epi8(src8, zero);
+      coeff16 = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(3, 3, 2, 2));
+      coeff16 = _mm_unpacklo_epi16(coeff16, coeff16);
+      mul_hi = _mm_mulhi_epi16(src16, coeff16);
+      mul_lo = _mm_mullo_epi16(src16, coeff16);
+      t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+      accum = _mm_add_epi32(accum, t);
+    }
+
+    // shift right for fix point implementation before saturation.
+    accum = _mm_srai_epi32 (accum, ConvolutionFilter1D::kShiftBits);
+    accum = _mm_packs_epi32(accum, zero);
+    accum = _mm_packus_epi16(accum, zero);

[evannier, 2011/02/14 23:45:13]

I am missing something.
the first packs_epi32 will move everything in 16 bit values, sign saturated:
// [16] 0, 0, 0, 0, a', r', g', b' 
The second, will take the 16 bit values, and apply the same thing, except this
is unsigned.
So, if we have a negative value, say a' is 0xFFFF, the resulting a'' will be
0xFF, right ?
(using [8] 0 0 ... 0 a'' r'' g'' b'')
Then, instead of having '0', we end up getting 0xFF.

As a test, make sure that your convolver handles well negative values as
coefficients (to generate a negative value), and very large values for > 0xFF
values.

[jiesun, 2011/02/15 00:19:22]

according to http://msdn.microsoft.com/en-us/library/07ad1wx4(v=vs.80).aspx
it will pack 16 bits to 8 bits with unsigned saturation.
that's what i want exactly. ?

+
+    *(reinterpret_cast<int*>(out_row)) = _mm_cvtsi128_si32(accum);
+    out_row += 4;
+  }
+#endif
+}
+
+// Convolves horizontally along four rows. The row data is given in
+// |src_data| and continues for the num_values() of the filter.
+template<bool has_alpha>
+void ConvolveHorizontally4_SSE2(const unsigned char* src_data[4],
+                          const ConvolutionFilter1D& filter,
+                          unsigned char* out_row[4]) {
+#ifdef ARCH_CPU_X86_FAMILY

[evannier, 2011/02/14 23:45:13]

Since this code is SSE2, I would use the preprocessor SSE2:
#if defined __SSE2___
Not sure if this is portable across all platforms, but it will work on GCC.
That implies that the code would need to be compiled with -msse2 to be included,
which is not acceptable for compatibility with non sse2 processors, since you do
not want GCC to generate SSE instructions when it has no right to.
So, what I had been thinking was to split off this file into multiple files,
compiled for various targets. This way you are guaranteed that the compiler will
not complain for x86 targets that do not compile with -msse2 (actually I am not
sure whether this is something that all versions of Chrome already do, since our
platform can always assume that x86 => sse2 or sse3 for that matter).

I do not know if this a good example, but Skia handles this type of issues in
this fashion:

- Create a function getter:
typedef (*Foo)(myinputs);
Foo* getFunction() {
  Foo* p = platformFunction();
  if (p) {
    return p;
  }
  return GenericFoo;
}

Then, for each platform, you have:
PlatformFunction_x86() {
  if (has_sse2) {  // check the right registers.
     return Foo_SSE2;
  }
  return NULL;
}
The two files containing the functions above being compiled without -msse2

Then, another file:
Foo_SSE2.cc
Foo_SSE2() {
   // put the SSE2 code here
}

this one compiled with SSE2.

Again, this might be moot if Chrome assumes that all platforms it runs on have
SSE2, but I thought I mentioned it.

[jiesun, 2011/02/15 00:19:22]

thanks, I was trying to do that function pointer magic before, but there are
some issues with gcc function pointer with templated function... I am thinking
to use virtual function and create different convolver class on different
platform. 

+  int width = filter.num_values();
+
+  int filter_offset, filter_length;
+  __m128i zero = _mm_setzero_si128();
+
+  for (int i=0; i < width; i+=1) {
+    const ConvolutionFilter1D::Fixed* filter_values =
+        filter.FilterForValue(i, &filter_offset, &filter_length);
+
+    // four pixels in a column per iteration.
+    __m128i accum0 = _mm_setzero_si128();
+    __m128i accum1 = _mm_setzero_si128();
+    __m128i accum2 = _mm_setzero_si128();
+    __m128i accum3 = _mm_setzero_si128();
+    int start = (filter_offset<<2);
+    for (int j = 0; j < (filter_length >> 2); ++j) {
+      __m128i coeff, coeff16lo, coeff16hi;
+      // [16] xx xx xx xx c3 c2 c1 c0
+      coeff = _mm_loadl_epi64((__m128i*)filter_values);
+      // [16] xx xx xx xx c1 c1 c0 c0
+      coeff16lo = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(1, 1, 0, 0));
+      // [16] c1 c1 c1 c1 c0 c0 c0 c0
+      coeff16lo = _mm_unpacklo_epi16(coeff16lo, coeff16lo);
+      // [16] xx xx xx xx c3 c3 c2 c2
+      coeff16hi = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(3, 3, 2, 2));
+      // [16] c3 c3 c3 c3 c2 c2 c2 c2
+      coeff16hi = _mm_unpacklo_epi16(coeff16hi, coeff16hi);
+
+      __m128i src8, src16, mul_hi, mul_lo, t;
+
+#define ITERATION(src, accum)                       \
+      src8 = _mm_loadu_si128((__m128i*)(src));      \
+      src16 = _mm_unpacklo_epi8(src8, zero);        \
+      mul_hi = _mm_mulhi_epi16(src16, coeff16lo);   \
+      mul_lo = _mm_mullo_epi16(src16, coeff16lo);   \
+      t = _mm_unpacklo_epi16(mul_lo, mul_hi);       \
+      accum = _mm_add_epi32(accum, t);              \
+      t = _mm_unpackhi_epi16(mul_lo, mul_hi);       \
+      accum = _mm_add_epi32(accum, t);              \
+      src16 = _mm_unpackhi_epi8(src8, zero);        \
+      mul_hi = _mm_mulhi_epi16(src16, coeff16hi);   \
+      mul_lo = _mm_mullo_epi16(src16, coeff16hi);   \
+      t = _mm_unpacklo_epi16(mul_lo, mul_hi);       \
+      accum = _mm_add_epi32(accum, t);              \
+      t = _mm_unpackhi_epi16(mul_lo, mul_hi);       \
+      accum = _mm_add_epi32(accum, t)
+
+      ITERATION(src_data[0]+start, accum0);
+      ITERATION(src_data[1]+start, accum1);
+      ITERATION(src_data[2]+start, accum2);
+      ITERATION(src_data[3]+start, accum3);
+
+      start += 16;
+      filter_values += 4;
+    }
+
+    int r = filter_length&3;
+    if (r) {
+      // Note: filter_values must be pad to align_up(filter_offset, 8);
+      __m128i coeff = _mm_loadl_epi64((__m128i*)filter_values);
+      __m128i mask = _mm_set_epi16(0, 0, 0, 0, 0, r==3?-1:0, r>=2?-1:0, -1);
+      coeff = _mm_and_si128(coeff, mask);
+
+      __m128i coeff16lo = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(1, 1, 0, 0));
+      /* c1 c1 c1 c1 c0 c0 c0 c0 */
+      coeff16lo = _mm_unpacklo_epi16(coeff16lo, coeff16lo);
+      __m128i coeff16hi = _mm_shufflelo_epi16(coeff, _MM_SHUFFLE(3, 3, 2, 2));
+      coeff16hi = _mm_unpacklo_epi16(coeff16hi, coeff16hi);
+
+      __m128i src8, src16, mul_hi, mul_lo, t;
+
+      ITERATION(src_data[0]+start, accum0);
+      ITERATION(src_data[1]+start, accum1);
+      ITERATION(src_data[2]+start, accum2);
+      ITERATION(src_data[3]+start, accum3);
+    }
+
+    accum0 = _mm_srai_epi32 (accum0, ConvolutionFilter1D::kShiftBits);
+    accum0 = _mm_packs_epi32(accum0, zero);
+    accum0 = _mm_packus_epi16(accum0, zero);
+    accum1 = _mm_srai_epi32 (accum1, ConvolutionFilter1D::kShiftBits);
+    accum1 = _mm_packs_epi32(accum1, zero);
+    accum1 = _mm_packus_epi16(accum1, zero);
+    accum2 = _mm_srai_epi32 (accum2, ConvolutionFilter1D::kShiftBits);
+    accum2 = _mm_packs_epi32(accum2, zero);
+    accum2 = _mm_packus_epi16(accum2, zero);
+    accum3 = _mm_srai_epi32 (accum3, ConvolutionFilter1D::kShiftBits);
+    accum3 = _mm_packs_epi32(accum3, zero);
+    accum3 = _mm_packus_epi16(accum3, zero);
+
+    *(reinterpret_cast<int*>(out_row[0])) = _mm_cvtsi128_si32(accum0);
+    *(reinterpret_cast<int*>(out_row[1])) = _mm_cvtsi128_si32(accum1);
+    *(reinterpret_cast<int*>(out_row[2])) = _mm_cvtsi128_si32(accum2);
+    *(reinterpret_cast<int*>(out_row[3])) = _mm_cvtsi128_si32(accum3);
+
+    out_row[0] += 4;
+    out_row[1] += 4;
+    out_row[2] += 4;
+    out_row[3] += 4;
+  }
+#endif
+}
+
+template<bool has_alpha>
+void ConvolveVertically_SSE2(const ConvolutionFilter1D::Fixed* filter_values,
+                        int filter_length,
+                        unsigned char* const* source_data_rows,
+                        int pixel_width,
+                        unsigned char* out_row) {
+#ifdef ARCH_CPU_X86_FAMILY
+  int width = pixel_width & ~3;
+
+   __m128i zero = _mm_setzero_si128();
+   __m128i accum0, accum1, accum2, accum3, coeff16;
+   for (int i = 0; i < width; i += 4) {  // Four pixels per iteration.
+     accum0 = _mm_setzero_si128();
+     accum1 = _mm_setzero_si128();
+     accum2 = _mm_setzero_si128();
+     accum3 = _mm_setzero_si128();
+     for (int j = 0; j < filter_length; ++j) {
+       coeff16 = _mm_set1_epi16(filter_values[j]);
+
+       // aligned load due to row_buffer is 16 byte aligned.
+       // [8] a3 b3 g3 r3 a2 b2 g2 r2 a1 b1 g1 r1 a0 b0 g0 r0
+       __m128i src8 = _mm_loadu_si128((__m128i*)(&source_data_rows[j][i<<2]));
+       // [16] a1 b1 g1 r1 a0 b0 g0 r0
+       __m128i src16 = _mm_unpacklo_epi8(src8, zero);
+       __m128i mul_hi = _mm_mulhi_epi16(src16, coeff16);
+       __m128i mul_lo = _mm_mullo_epi16(src16, coeff16);
+       // [32] a0 b0 g0 r0
+       __m128i t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+       accum0 = _mm_add_epi32(accum0, t);
+       // [32] a1 b1 g1 r1
+       t = _mm_unpackhi_epi16(mul_lo, mul_hi);
+       accum1 = _mm_add_epi32(accum1, t);
+       // [16] a3 b3 g3 r3 a2 b2 g2 r2
+       src16 = _mm_unpackhi_epi8(src8, zero);
+       mul_hi = _mm_mulhi_epi16(src16, coeff16);
+       mul_lo = _mm_mullo_epi16(src16, coeff16);
+       // [32] a2 b2 g2 r2
+       t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+       accum2 = _mm_add_epi32(accum2, t);
+       // [32] a3 b3 g3 r3
+       t = _mm_unpackhi_epi16(mul_lo, mul_hi);
+       accum3 = _mm_add_epi32(accum3, t);
+     }
+
+     accum0 = _mm_srai_epi32 (accum0, ConvolutionFilter1D::kShiftBits);
+     accum1 = _mm_srai_epi32 (accum1, ConvolutionFilter1D::kShiftBits);
+     accum2 = _mm_srai_epi32 (accum2, ConvolutionFilter1D::kShiftBits);
+     accum3 = _mm_srai_epi32 (accum3, ConvolutionFilter1D::kShiftBits);
+     // [16] a1 b1 g1 r1 a0 b0 g0 r0
+     accum0 = _mm_packs_epi32(accum0, accum1);
+     // [16] a3 b3 g3 r3 a2 b2 g2 r2
+     accum2 = _mm_packs_epi32(accum2, accum3);
+     // [8] a3 b3 g3 r3 a2 b2 g2 r2 a1 b1 g1 r1 a0 b0 g0 r0
+     accum0 = _mm_packus_epi16(accum0, accum2);
+
+     if (has_alpha) {
+       // [8] xx a3 b3 g3 xx a2 b2 g2 xx a1 b1 g1 xx a0 b0 g0
+       __m128i a = _mm_srli_epi32(accum0, 8);
+       // [8] xx xx xx max3 xx xx xx max2 xx xx xx max1 xx xx xx max0
+       __m128i b = _mm_max_epu8(a, accum0);  // Max of r and g.
+       // [8] xx xx a3 b3 xx xx a2 b2 xx xx a1 b1 xx xx a0 b0
+       a = _mm_srli_epi32(accum0, 16);
+       // [8] xx xx xx max3 xx xx xx max2 xx xx xx max1 xx xx xx max0
+       b = _mm_max_epu8(a, b);  // Max of r and g and b.
+       // [8] max3 00 00 00 max2 00 00 00 max1 00 00 00 max0 00 00 00
+       b = _mm_slli_epi32(b, 24);
+       accum0 = _mm_max_epu8(b, accum0);
+     } else {
+       __m128i mask = _mm_set_epi32(0xff000000, 0xff000000,
+                                    0xff000000, 0xff000000);
+       accum0 = _mm_or_si128(accum0, mask);
+     }
+     _mm_storeu_si128((__m128i*)out_row, accum0);
+     out_row += 16;
+   }
+
+   if (pixel_width & 3) {
+     accum0 = _mm_setzero_si128();
+     accum1 = _mm_setzero_si128();
+     accum2 = _mm_setzero_si128();
+     for (int j = 0; j < filter_length; ++j) {
+       coeff16 = _mm_set1_epi16(filter_values[j]);
+
+       // [8] a3 b3 g3 r3 a2 b2 g2 r2 a1 b1 g1 r1 a0 b0 g0 r0
+       __m128i src8 = _mm_loadu_si128((__m128i*)(&source_data_rows[j][width<<2]));
+       // [16] a1 b1 g1 r1 a0 b0 g0 r0
+       __m128i src16 = _mm_unpacklo_epi8(src8, zero);
+       __m128i mul_hi = _mm_mulhi_epi16(src16, coeff16);
+       __m128i mul_lo = _mm_mullo_epi16(src16, coeff16);
+       // [32] a0 b0 g0 r0
+       __m128i t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+       accum0 = _mm_add_epi32(accum0, t);
+       // [32] a1 b1 g1 r1
+       t = _mm_unpackhi_epi16(mul_lo, mul_hi);
+       accum1 = _mm_add_epi32(accum1, t);
+       // [16] a3 b3 g3 r3 a2 b2 g2 r2
+       src16 = _mm_unpackhi_epi8(src8, zero);
+       mul_hi = _mm_mulhi_epi16(src16, coeff16);
+       mul_lo = _mm_mullo_epi16(src16, coeff16);
+       // [32] a2 b2 g2 r2
+       t = _mm_unpacklo_epi16(mul_lo, mul_hi);
+       accum2 = _mm_add_epi32(accum2, t);
+     }
+
+     accum0 = _mm_srai_epi32 (accum0, ConvolutionFilter1D::kShiftBits);
+     accum1 = _mm_srai_epi32 (accum1, ConvolutionFilter1D::kShiftBits);
+     accum2 = _mm_srai_epi32 (accum2, ConvolutionFilter1D::kShiftBits);
+     // [16] a1 b1 g1 r1 a0 b0 g0 r0
+     accum0 = _mm_packs_epi32(accum0, accum1);
+     // [16] a3 b3 g3 r3 a2 b2 g2 r2
+     accum2 = _mm_packs_epi32(accum2, zero);
+     // [8] a3 b3 g3 r3 a2 b2 g2 r2 a1 b1 g1 r1 a0 b0 g0 r0
+     accum0 = _mm_packus_epi16(accum0, accum2);
+
+     if (has_alpha) {
+       // [8] xx a3 b3 g3 xx a2 b2 g2 xx a1 b1 g1 xx a0 b0 g0
+       __m128i a = _mm_srli_epi32(accum0, 8);
+       // [8] xx xx xx max3 xx xx xx max2 xx xx xx max1 xx xx xx max0
+       __m128i b = _mm_max_epu8(a, accum0);  // Max of r and g.
+       // [8] xx xx a3 b3 xx xx a2 b2 xx xx a1 b1 xx xx a0 b0
+       a = _mm_srli_epi32(accum0, 16);
+       // [8] xx xx xx max3 xx xx xx max2 xx xx xx max1 xx xx xx max0
+       b = _mm_max_epu8(a, b);  // Max of r and g and b.
+       // [8] max3 00 00 00 max2 00 00 00 max1 00 00 00 max0 00 00 00
+       b = _mm_slli_epi32(b, 24);
+       accum0 = _mm_max_epu8(b, accum0);
+     } else {
+       __m128i mask = _mm_set_epi32(0xff000000, 0xff000000,
+                                    0xff000000, 0xff000000);
+       accum0 = _mm_or_si128(accum0, mask);
+     }
+
+     for (int i = width; i < pixel_width; ++i) {
+       *(reinterpret_cast<int*>(out_row)) = _mm_cvtsi128_si32(accum0);
+       accum0 = _mm_srli_si128(accum0, 4);
+       out_row += 4;
+     }
+   }
+#endif
+}
+
+
 }  // namespace
 
 // ConvolutionFilter1D ---------------------------------------------------------
@@ -292,7 +648,8 @@ void BGRAConvolve2D(const unsigned char* source_data,
                     const ConvolutionFilter1D& filter_x,
                     const ConvolutionFilter1D& filter_y,
                     int output_byte_row_stride,
-                    unsigned char* output) {
+                    unsigned char* output,
+                    bool use_sse2) {

[evannier, 2011/02/14 23:45:13]

This does not seem right to be part of the API, given that SSE2 is not a
relevant feature from the outside.
The problem I understand is that the padding you need to apply is going to be
dependent on whether this is going to be set. Suggesting creating a Convolver
Engine virtual class which can have a type SSE2, and then knows how to deal with
padding and realignment of coefficients.

[jiesun, 2011/02/15 00:19:22]

also to test by comparing two version, I had to override the default decision at
some level... 

   int max_y_filter_size = filter_y.max_filter();
 
   // The next row in the input that we will generate a horizontally
@@ -310,29 +667,82 @@ void BGRAConvolve2D(const unsigned char* source_data,
   // a circular buffer of convolved rows and do vertical convolution as rows
   // are available. This prevents us from having to store the entire
   // intermediate image and helps cache coherency.
-  CircularRowBuffer row_buffer(filter_x.num_values(), max_y_filter_size,
+  // We will need four extra rows to allow horizontal convolution could be done
+  // simultaneously.
+  int row_buffer_width = (filter_x.num_values() + 15) & ~0xF;
+  int row_buffer_height = max_y_filter_size + (use_sse2 ? 4 : 0);
+  CircularRowBuffer row_buffer(row_buffer_width,
+                               row_buffer_height,
                                filter_offset);
 
   // Loop over every possible output row, processing just enough horizontal
   // convolutions to run each subsequent vertical convolution.
   SkASSERT(output_byte_row_stride >= filter_x.num_values() * 4);
   int num_output_rows = filter_y.num_values();
+
+  int last_filter_offset, last_filter_length;
+  filter_y.FilterForValue(num_output_rows-1, &last_filter_offset,
+                          &last_filter_length);
+
   for (int out_y = 0; out_y < num_output_rows; out_y++) {
     filter_values = filter_y.FilterForValue(out_y,
                                             &filter_offset, &filter_length);
 
     // Generate output rows until we have enough to run the current filter.
-    while (next_x_row < filter_offset + filter_length) {
-      if (source_has_alpha) {
-        ConvolveHorizontally<true>(
-            &source_data[next_x_row * source_byte_row_stride],
-            filter_x, row_buffer.AdvanceRow());
-      } else {
-        ConvolveHorizontally<false>(
-            &source_data[next_x_row * source_byte_row_stride],
-            filter_x, row_buffer.AdvanceRow());
+    if (use_sse2) {
+      while (next_x_row < filter_offset + filter_length) {
+        if (next_x_row + 3 < last_filter_offset + last_filter_length - 1) {
+        //if (0) {

[fbarchard, 2011/02/15 23:16:45]

remove if you're done with this

+          const unsigned char* src[4];
+          unsigned char* out_row[4];
+          for (int i=0; i<4; i++) {
+            src[i] = &source_data[(next_x_row+i) * source_byte_row_stride];
+            out_row[i] = row_buffer.AdvanceRow();
+          }
+          if (source_has_alpha)
+            ConvolveHorizontally4_SSE2<true>(src, filter_x, out_row);
+          else
+            ConvolveHorizontally4_SSE2<false>(src, filter_x, out_row);
+          next_x_row+=4;
+        } else {
+          // For the last row, SSE2 load possibly to access data beyond the
+          // image area. therefore we use C version here. Hacking into skia
+          // is not something in my mind.
+          if (next_x_row == last_filter_offset + last_filter_length - 1 ) {
+            if (source_has_alpha)
+              ConvolveHorizontally<true>(
+                  &source_data[next_x_row * source_byte_row_stride],
+                  filter_x, row_buffer.AdvanceRow());
+            else
+              ConvolveHorizontally<false>(
+                  &source_data[next_x_row * source_byte_row_stride],
+                  filter_x, row_buffer.AdvanceRow());
+          } else {
+            if (source_has_alpha)
+              ConvolveHorizontally_SSE2<true>(
+                  &source_data[next_x_row * source_byte_row_stride],
+                  filter_x, row_buffer.AdvanceRow());
+            else
+              ConvolveHorizontally_SSE2<false>(
+                  &source_data[next_x_row * source_byte_row_stride],
+                  filter_x, row_buffer.AdvanceRow());
+          }
+          next_x_row++;
+        }
+      }
+    } else {
+      while (next_x_row < filter_offset + filter_length) {
+        if (source_has_alpha) {
+          ConvolveHorizontally<true>(
+              &source_data[next_x_row * source_byte_row_stride],
+              filter_x, row_buffer.AdvanceRow());
+        } else {
+          ConvolveHorizontally<false>(
+              &source_data[next_x_row * source_byte_row_stride],
+              filter_x, row_buffer.AdvanceRow());
+        }
+        next_x_row++;
       }
-      next_x_row++;
     }
 
     // Compute where in the output image this row of final data will go.
@@ -349,13 +759,25 @@ void BGRAConvolve2D(const unsigned char* source_data,
         &rows_to_convolve[filter_offset - first_row_in_circular_buffer];
 
     if (source_has_alpha) {
-      ConvolveVertically<true>(filter_values, filter_length,
-                               first_row_for_filter,
-                               filter_x.num_values(), cur_output_row);
+      if (use_sse2) {
+        ConvolveVertically_SSE2<true>(filter_values, filter_length,
+                                 first_row_for_filter,
+                                 filter_x.num_values(), cur_output_row);
+      } else {
+        ConvolveVertically<true>(filter_values, filter_length,
+                                 first_row_for_filter,
+                                 filter_x.num_values(), cur_output_row);
+      }
     } else {
-      ConvolveVertically<false>(filter_values, filter_length,
-                                first_row_for_filter,
-                                filter_x.num_values(), cur_output_row);
+      if (use_sse2) {
+        ConvolveVertically_SSE2<false>(filter_values, filter_length,
+                                  first_row_for_filter,
+                                  filter_x.num_values(), cur_output_row);
+      } else {
+        ConvolveVertically<false>(filter_values, filter_length,
+                                 first_row_for_filter,
+                                 filter_x.num_values(), cur_output_row);
+      }
     }
   }
 }