title: Fast computation of warp matrices
body: <p>For a fixed and given tform, the imwarp command in the Image Processing Toolbox</p>
  
  <pre><code>  B = imwarp(A,tform)
  </code></pre>
  
  <p>is linear with respect to A, meaning there exists some sparse matrix W, <strong>depending on tform but independent of A</strong>, such that the above can be equivalently implemented</p>
  
  <pre><code>  B(:)=W*A(:)
  </code></pre>
  
  <p>for all A of fixed known dimensions [n,n]. My question is whether there are fast/efficient options for computing W. The matrix form is necessary when I need the transpose operation W.'*B(:), or if I need to do W\B(:) or similar linear algebraic things which I can't do directly through imwarp alone.</p>
  
  <p>I know that it is possible to compute W column-by-column by doing</p>
  
  <pre><code>   E=zeros(n);
     W=spalloc(n^2,n^2,4*n^2);
  
     for i=1:n^2
  
      E(i)=1;
       tmp=imwarp(E,tform);
      E(i)=0;
  
      W(:,i)=tmp(:);
  
     end
  </code></pre>
  
  <p>but this is brute force and slow. </p>
  
  <p>The routine <a href="http://www.mathworks.com/matlabcentral/fileexchange/44669-func2mat-convert-linear-function-to-matrix" rel="nofollow">FUNC2MAT</a> is somewhat more optimal in that it uses the loop to compute/gather the sparse entry data I,J,S of each column W(:,i). Then, after the loop, it uses this to construct the overall sparse matrix. It also offers the option of using a PARFOR loop. However, this is still slower than I would like.</p>
  
  <p>Can anyone suggest more speed-optimal alternatives? </p>
  
  <p><strong>EDIT:</strong></p>
  
  <p>For those uncomfortable with my claim that imwarp(A,tform) is linear w.r.t. A, I include the demo script below, which tests that the superposition property is satisfied for random input images and tform data. It can be run repeatedly to see that the <code>nonlinearityError</code> is always small, and easily attributable to floating point noise.</p>
  
  <pre><code>tform=affine2d(rand(3,2));
   %tform=projective2d(rand(3));
  
  fun=@(A) imwarp(A,tform,'cubic');
  
  I1=rand(100); I2=rand(100);
  c1=rand; c2=rand;
  
  
  LHS=fun(c1*I1+c2*I2);  %left hand side
  RHS=c1*fun(I1)+c2*fun(I2); %right hand side
  
  linearityError = norm(LHS(:)-RHS(:),'inf')
  </code></pre>