title: C++ - Fastest possible bilinear interpolation?
body: <p>I would like to speed up my c++ bilinear interpolation code.</p>
  
  <p>The setup is as following: <strong>From a grayscale image <em>img</em> I want to extract a rectangular patch <em>pat</em> at position <em>cent</em> with unit spacing and no up-/downsampling.</strong></p>
  
  <p>Since <em>cent</em> will generally not be integer, I have to bilinearly interpolate the extracted patch.</p>
  
  <p>Image <em>img</em>, extracted patch <em>pat</em> and position <em>cent</em> are stored as floats.
  The patch is of size [2*pad+1], with <em>pad</em> being the padding left and right of position <em>cent</em>.</p>
  
  <p>The current solution looks like this:</p>
  
  <pre><code>void function(Eigen::Matrix&lt;float, Eigen::Dynamic, 1&gt;* pat, 
                const float* img, 
                const Eigen::Vector2f* cent)
  {
  
    Eigen::Vector4f we; // bilinear weight vector
    // ... [CROPPED: compute bilinear weights]
  
    float *pat_it = pat-&gt;data();
    for (y=cent[1]-pad; y &lt;= cent[1]+pad; ++y)    
    {
      int postmp_a = y        * image_width;
      int postmp_b = (y-1)    * image_width;
  
      for (x=cent[0]-pad; x &lt;= cent[0]+pad; ++x, ++pat_it)    
      {
  
            (*pat_it)     = we[0] * img[ x    +  postmp_a] +  
                            we[1] * img[x-1   +  postmp_a] +
                            we[2] * img[ x    +  postmp_b] +
                            we[3] * img[x-1   +  postmp_b]; 
      }
    }
  }
  </code></pre>
  
  <ol>
  <li><p><strong>Is it possible to speed this up further?</strong> This function will be called several millions of times in a real-time signal processing pipeline. <strong>There is no memory constraint.</strong></p></li>
  <li><p>Are there maybe specific Eigen functions for that?</p></li>
  <li><p>Since this is the most crucial bottleneck of my code, I would also be willing to consider moving the code to different programming languages/architectures (Assembler, CUDA, etc...). Any thoughts/hints on this?</p></li>
  <li><p>More generally, how would you systematically approach this for profiling?</p></li>
  </ol>
  
  <p>Some more details: The Code is compiled with '-Ofast -std=c++11' and is already running in parallel using OpenMP. The image size is in the order of ~1000x1200 pixels, and <em>pad</em> between 5-10 pixels.</p>
  
  <p><strong>EDIT</strong></p>
  
  <p>I already managed to get a ~ 6 percent speed-up by using pointers to the 4 corresponding image locations directly.</p>
  
  <pre><code>...
  for (x=cent[0]-pad; x &lt;= cent[0]+pad; ++x,++pat_it,
       ++img_a,++img_b,++img_c,++img_d)    
  {
  
        (*pat_it)   = we[0] * (*img_a) +  
                      we[1] * (*img_b) +
                      we[2] * (*img_c) +
                      we[3] * (*img_d); 
  }
  ...
  </code></pre>