title: Fast SSE threshold algorithm
body: <p>I'm trying to come up with a very fast threshold algorithm using SSE to replace this:</p>
  
  <pre><code>uint8_t *pSrc, *pDst;
  
  // Assume pSrc and pDst point to valid data
  
  // Handle left edge
  *pDst++ = *pSrc++;
  
  // Likeness filter
  for (uint32_t k = 2; k &lt; width; k++, pSrc++, pDst++)
      if ((*pDst - *pSrc) * (*pDst - *pSrc) &gt; 100 /*THRESHOLD_SQUARED*/) {
          *pDst = *pSrc;
      }
  }
  
  // Handle right edge
  *pDst++ = *pSrc++;
  </code></pre>
  
  <p>So far I have this:</p>
  
  <pre><code>const uint8_t THRESHOLD = 10;
  
  __attribute__((aligned (16))) static const uint8_t mask[16] = {
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD,
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD,
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD,
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD
  };
  
  __m128i xmm1, xmm3, xmm4, xmm5, xmm6, xmm7, xmm8, xmm9;
  xmm1 = _mm_load_si128((__m128i const *)mask);
  xmm6 = _mm_setzero_si128();
  
  uint8_t *pSrc, *pDst;
  
  // Assume pSrc and pDst point to valid data
  
  // I have other code with another mask for the first 16 entries
  
  for (uint32_t k = 16; k &lt; (width - 16); k += 16, pSrc += 16, pDst += 16) {
      xmm3 = _mm_load_si128((__m128i const *)pDst);
      xmm4 = _mm_load_si128((__m128i const *)pSrc);
      xmm5 = _mm_unpacklo_epi8(xmm3, xmm6);
      xmm7 = _mm_unpackhi_epi8(xmm3, xmm6);
      xmm8 = _mm_unpacklo_epi8(xmm4, xmm6);
      xmm9 = _mm_unpackhi_epi8(xmm4, xmm6);
      xmm5 = _mm_sub_epi16(xmm5, xmm8);
      xmm7 = _mm_sub_epi16(xmm7, xmm9);
      xmm5 = _mm_abs_epi16(xmm5);
      xmm7 = _mm_abs_epi16(xmm7);
      xmm5 = _mm_packs_epi16(xmm5, xmm7);
      xmm5 = _mm_cmpgt_epi8(xmm5, xmm1);
      xmm3 = _mm_blendv_epi8(xmm3, xmm4, xmm5);
      _mm_store_si128((__m128i *)pDst, xmm3);
  }
  
  // I have other code with another mask for the last 16 entries
  </code></pre>
  
  <p>I have idea using another type of algorithm to handle the absolute value of the difference of two values (mainly to stay in U8 (uchar) space):</p>
  
  <pre><code>a' = a &gt;&gt; 1;
  b' = b &gt;&gt; 1;
  diff = (abs(sub(a' - b')) &lt;&lt; 1) + ((a ^ b) &amp; 1);
  </code></pre>
  
  <p>This would take 8 SSE instructions instead of the above 9 (not including any extra register moves the compiler generates) but I'm not sure if it is faster due to dependency latencies.</p>
  
  <p>Does any other SSE expert have any better suggestions (using up to SSE 4.2)?</p>
  
  <p><strong>Update 1 - Thanks to Yves's suggestion!</strong></p>
  
  <pre><code>const uint8_t THRESHOLD = 10;
  
  __attribute__((aligned (16))) static const uint8_t mask[16] = {
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD,
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD,
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD,
      THRESHOLD, THRESHOLD, THRESHOLD, THRESHOLD
  };
  
  __m128i xmm1, xmm3, xmm4, xmm5, xmm6, xmm7;
  xmm1 = _mm_load_si128((__m128i const *)mask);
  xmm6 = _mm_setzero_si128();
  
  uint8_t *pSrc, *pDst;
  
  // Assume pSrc and pDst point to valid data
  
  // I have other code with another mask for the first 16 entries
  
  for (uint32_t k = 16; k &lt; (width - 16); k += 16, pSrc += 16, pDst += 16) {
      xmm3 = _mm_load_si128((__m128i const *)pDst);
      xmm4 = _mm_load_si128((__m128i const *)pSrc);
      xmm5 = _mm_subs_epu8(xmm3, xmm4);
      xmm7 = _mm_subs_epu8(xmm4, xmm3);
      xmm5 = _mm_adds_epu8(xmm5, xmm7);
      xmm5 = _mm_subs_epu8(xmm5, xmm1);
      xmm5 = _mm_cmpeq_epi8(xmm5, xmm6);
      xmm4 = _mm_blendv_epi8(xmm4, xmm3, xmm5);
      _mm_store_si128((__m128i *)pDst, xmm4);
  }
  
  // I have other code with another mask for the last 16 entries
  </code></pre>