title: Transpose SSE2 Vectors
body: <p>I try to convolve an image for wavelet decomposition using SSE2 and C. This image has 4 channels (Lab + alpha) stored contiguously in memory : [LabA][LabA][LabA]… The alpha channel is irrelevant for what I do here.</p>
  
  <p>Accessing a pixel is then straightforward by loading the content of a pointer incremented by 4 interatively:</p>
  
  <pre><code>static void eaw_decompose_sse2(float *const out, 
                                 const float *const in, 
                                 float *const detail, 
                                 const int scale,
                                 const float sharpen, 
                                 const size_t width, 
                                 const size_t height)
  {
  /* Convolve rows */
  #ifdef _OPENMP
  #pragma omp parallel for collapse(2)
  #endif
    for(size_t j = 0; j &lt; height; j++)
    {
      for(size_t i = 0; i &lt; width; i++)
      {
        const size_t inc = (j * width + i) * 4;
        float *pdetail = detail + inc;
        float *pcoarse = tmp + inc;
  
        // pixel to be convolved
        const __m128 pin0 = _mm_load_ps(in + inc);
        const __m128 w_0 = _mm_set1_ps(filter[2]);
  
        // neighbours
        const __m128 pin1 = _mm_load_ps(in + ASAN_ROW(i, j, -2, mult, max_height_i, width));
        const __m128 pin2 = _mm_load_ps(in + ASAN_ROW(i, j, -1, mult, max_height_i, width));
        const __m128 pin3 = _mm_load_ps(in + ASAN_ROW(i, j, 1, mult, max_height_i, width));
        const __m128 pin4 = _mm_load_ps(in + ASAN_ROW(i, j, 2, mult, max_height_i, width));
  
        // neighbours contribution
        const __m128 w_1 = _mm_set1_ps(filter[0]) * weight_sse2(pin0, pin1, sharpen);
        const __m128 w_2 = _mm_set1_ps(filter[1]) * weight_sse2(pin0, pin2, sharpen);
        const __m128 w_3 = _mm_set1_ps(filter[3]) * weight_sse2(pin0, pin2, sharpen);
        const __m128 w_4 = _mm_set1_ps(filter[4]) * weight_sse2(pin0, pin3, sharpen);
  
        // Filter computation
        const __m128 wgt = w_1 + w_2 + w_3 + w_4 + w_0;
        const __m128 sum = (w_1 * pin1 + w_2 * pin2 + w_3 * pin3 + w_4 * pin4 + w_0 * pin0) * _mm_rcp_ps(wgt);
  
        // High frequency layer
        _mm_stream_ps(pdetail, pin0 - sum);
  
        // Low frequency layer
        _mm_stream_ps(pcoarse, sum);
      }
    }
  }
  </code></pre>
  
  <p>The function <code>ASAN_ROW</code> slides the pointer along the rows ensuring we stay in the bounds, if not, it takes the nearest neighbour. <code>weight_sse2</code> is a gaussian weight that does complicated bit-shifts because L and a/b have different weightings.</p>
  
  <p>So, instead of operating on 4 Lab SSE vectors, with the last element lost, I feel it would be faster to operate on 3 SSE vecors, each vector being a Lab channel, of which each element is a neighbouring pixel. So that would become:</p>
  
  <pre><code>static void eaw_decompose_sse2(float *const out, 
                                 const float *const in, 
                                 float *const detail, 
                                 const int scale,
                                 const float sharpen, 
                                 const size_t width, 
                                 const size_t height)
  {
  /* Convolve rows */
  #ifdef _OPENMP
  #pragma omp parallel for collapse(2)
  #endif
    for(size_t j = 0; j &lt; height; j++)
    {
      for(size_t i = 0; i &lt; width; i++)
      {
        const size_t inc = (j * width + i) * 4;
        float *pdetail = detail + inc;
        float *pcoarse = tmp + inc;
  
        // pixel to be convolved
        const __m128 pin0 = _mm_load_ps(in + inc);
        const __m128 w_0 = _mm_set1_ps(filter[2]);
  
        // neighbours
        const __m128 pin1 = _mm_load_ps(in + ASAN_ROW(i, j, -2, mult, max_height_i, width));
        const __m128 pin2 = _mm_load_ps(in + ASAN_ROW(i, j, -1, mult, max_height_i, width));
        const __m128 pin3 = _mm_load_ps(in + ASAN_ROW(i, j, 1, mult, max_height_i, width));
        const __m128 pin4 = _mm_load_ps(in + ASAN_ROW(i, j, 2, mult, max_height_i, width));
  
        // Lab extraction - pixel to be convolved
        __m128 L_0 =  _mm_set1_ps( pin0[0] ); // ?
        __m128 a_0 =  _mm_set1_ps( pin0[1] ); // ?
        __m128 b_0 =  _mm_set1_ps( pin0[2] ); // ?
  
        // Lab extraction - neighbours
        __m128 L_f = _mm_set_ps ({ pin1[0], pin2[0], pin3[0], pin4[0] }); // ?
        __m128 a_f = _mm_set_ps ({ pin1[1], pin2[1], pin3[1], pin4[1] }); // ?
        __m128 b_f = _mm_set_ps ({ pin1[2], pin2[2], pin3[2], pin4[2] }); // ?
  
        // neighbours contribution
        const __m128 filter = _mm_load_ps(filter_coeff);
        const __m128 w_L = filter * weight_sse(L_0, L_f, sharpen);
        const __m128 w_c = filter * weight_sse(a_0 + b_0, a_f + b_f, sharpen);
  
        // Filter computation
        const __m128 wgt = _mm_set_ps( { sum_of_elts_sse(w_L), 
                                         sum_of_elts_sse(w_c), 
                                         sum_of_elts_sse(w_c), 
                                         0.0f } );
        const __m128 w1 = _mm_set_ps ({ w_L[0], w_c[0], w_c[0], 0.0f }); // ?
        const __m128 w2 = _mm_set_ps ({ w_L[1], w_c[1], w_c[1], 0.0f }); // ?
        const __m128 w3 = _mm_set_ps ({ w_L[2], w_c[2], w_c[2], 0.0f }); // ?
        const __m128 w4 = _mm_set_ps ({ w_L[3], w_c[3], w_c[3], 0.0f }); // ?
        const __m128 sum = (w_1 * pin1 + w_2 * pin2 + w_3 * pin3 + w_4 * pin4 + w_0 * pin0) * _mm_rcp_ps(wgt);
  
        // High frequency layer
        _mm_stream_ps(pdetail, pin0 - sum);
  
        // Low frequency layer
        _mm_stream_ps(pcoarse, sum);
      }
    }
  }
  </code></pre>
  
  <p>What is the most cache-efficient way to switch from the channel-based vectors (pixels vectors <code>pin0</code> to <code>pin4</code>) to the neighbour-based vectors (<code>L_0</code>, <code>L_f</code>), and the other way around (<code>w_L</code> and <code>w_c</code> to <code>w_1</code>-<code>w_4</code>) ? Would the second version be faster ?</p>