title: Faster Method for Multiple Bilinear Interpolation?
body: <p>I am writing a program in C++ to reconstruct a 3D object from a set of projected 2D images, the most computation-intensive part of which involves magnifying and shifting each image via bilinear interpolation. I currently have a pair of functions for this task; "blnSetup" defines a handful of parameters outside the loop, then "bilinear" applies the interpolation point-by-point within the loop:</p>
  
  <p>(NOTE: 'I' is a 1D array containing ordered rows of image data)</p>
  
  <pre><code>//Pre-definition structure (in header)
  struct blnData{
      float* X;
      float* Y;
      int* I;
      float X0;
      float Y0;
      float delX;
      float delY;
  };
  
  //Pre-definition function (outside the FOR loop)
  extern inline blnData blnSetup(float* X, float* Y, int* I)
  {
      blnData bln;
      //Create pointers to X, Y, and I vectors
      bln.X = X;
      bln.Y = Y;
      bln.I = I;
  
      //Store offset and step values for X and Y
      bln.X0 = X[0];
      bln.delX = X[1] - X[0];
      bln.Y0 = Y[0];
      bln.delY = Y[1] - Y[0];
  
      return bln;
  }
  
  //Main interpolation function (inside the FOR loop)
  extern inline float bilinear(float x, float y, blnData bln)
  {
      float Ixy;
  
      //Return -1 if the target point is outside the image matrix
      if (x &lt; bln.X[0] || x &gt; bln.X[-1] || y &lt; bln.Y[0] || y &gt; bln.Y[-1])
          Ixy = 0;
      //Otherwise, apply bilinear interpolation
      else
      {
          //Define known image width
          int W = 200;
  
          //Find nearest indices for interpolation
          int i = floor((x - bln.X0) / bln.delX);
          int j = floor((y - bln.Y0) / bln.delY);
  
          //Interpolate I at (xi, yj)
          Ixy = 1 / ((bln.X[i + 1] - bln.X[i])*(bln.Y[j + 1] - bln.Y[j])) *
              (
              bln.I[W*j + i] * (bln.X[i + 1] - x) * (bln.Y[j + 1] - y) +
              bln.I[W*j + i + 1] * (x - bln.X[i]) * (bln.Y[j + 1] - y) +
              bln.I[W*(j + 1) + i] * (bln.X[i + 1] - x) * (y - bln.Y[j]) +
              bln.I[W*(j + 1) + i + 1] * (x - bln.X[i]) * (y - bln.Y[j])
              );
      }
  
      return Ixy;
  }
  </code></pre>
  
  <p>EDIT: The function calls are below. 'flat.imgdata' is a std::vector containing the input image data and 'proj.imgdata' is a std::vector containing the transformed image.</p>
  
  <pre><code>int Xs = flat.dim[0];
  int Ys = flat.dim[1];
  
  int* Iarr = flat.imgdata.data();
  float II, x, y;
  
  bln = blnSetup(X, Y, Iarr);
  
  for (int j = 0; j &lt; flat.imgdata.size(); j++)
  {
      x = 1.2*X[j % Xs];
      y = 1.2*Y[j / Xs];
      II = bilinear(x, y, bln);
      proj.imgdata[j] = (int)II;
  }
  </code></pre>
  
  <p>Since I started optimizing, I have been able to reduce computation time by ~50x (!) by switching from std::vectors to C arrays within the interpolation function, and another 2x or so by cleaning up redundant computations/typecasting/etc, but assuming O(n) with n being the total number of processed pixels, the full reconstruction (~7e10 pixels) should still take 40min or so--about an order of magnitude longer than my goal of &lt;5min.</p>
  
  <p>According to Visual Studio's performance profiler, the interpolation function call ("II = bilinear(x, y, bln);") is unsurprisingly still the majority of my computation load. I haven't been able to find any linear algebraic methods for fast multiple interpolation, so my question is: is this basically as fast as my code will get, short of applying more or faster CPUs to the task? Or is there a different approach that might speed things up?</p>
  
  <p>P.S. I've also only been coding in C++ for about a month now, so feel free to point out any beginner mistakes I might be making.</p>