Быстрого вычисления матриц варпа

Question

Для фиксированного и заданного TForm, команда imwarp в Image Processing Toolbox

B = imwarp(A,tform) 

линейно относительно А, значит, существует некоторый разреженная матрица W, в зависимости от TForm но независимый от, таким образом, что выше, могут быть реализованы

что то же самое

B(:)=W*A(:) 

для всех А фиксированных известных размеров [N, N]. Мой вопрос в том, есть ли быстрые/эффективные варианты для вычисления W. Матричная форма необходима, когда мне нужна операция транспонирования W. '* B (:), или если мне нужно сделать W \ B (:) или аналогичную линейную алгебраическую вещи, которые я не могу сделать непосредственно через imwarp в одиночку.

Я знаю, что можно вычислить W столбец за столбцом, делая

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 

, но это грубая сила и медленно.

Процедура FUNC2MAT несколько более оптимальна тем, что она использует цикл для вычисления/сбора разреженных данных ввода I, J, S каждого столбца W (:, i). Затем, после цикла, он использует это для построения общей разреженной матрицы. Он также предлагает возможность использования петли PARFOR. Однако это все еще медленнее, чем хотелось бы.

Может ли кто-нибудь предложить более оптимальные по скорости альтернативы?

EDIT:

Для тех дискомфортных с моим утверждением, что imwarp (A, ТГогт) линейна w.r.t. A, я включаю демо-скрипт ниже, который проверяет, что свойство суперпозиции выполнено для случайных входных изображений и данных tform. Его можно многократно запускать, чтобы увидеть, что nonlinearityError всегда мал и легко связан с шумом с плавающей запятой.

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') 

Answer 1

Это на самом деле довольно просто:

W = sparse(B(:)/A(:)); 

Обратите внимание, что W не является уникальной, но эта операция, вероятно, производит самый редкий результат. Другим способом его расчета будет

W = sparse(B(:) * pinv(A(:))); 

, но это приводит к значительно менее разреженному (но все еще актуальному) результату.

Answer 2

Я построил коробление матрицы, используя поля оптического потока [и, v], и она работает хорошо для моего приложения

% this function computes the warping matrix 

% M x N is the size of the image 

function [ Fw ] = generateFwi(u,v,M,N) 

    Fw = zeros(M*N, M*N); 

    k =1; 

    for i=1:M 

    for j= 1:N 

     newcoord(1) = i+u(i,j); 
     newcoord(2) = j+v(i,j); 



     newi = newcoord(1); 


     newj = newcoord(2); 


     if newi >0 && newj >0 

     newi1x = floor(newi); 
     newi1y = floor(newj); 
     newi2x = floor(newi); 
     newi2y = ceil(newj); 
     newi3x = ceil(newi); % four nearest points to the given point 
     newi3y = floor(newj); 
     newi4x = ceil(newi); 
     newi4y = ceil(newj); 
     x1 = [newi,newj;newi1x,newi1y]; 
     x2 = [newi,newj;newi2x,newi2y]; 
     x3 = [newi,newj;newi3x,newi3y]; 
     x4 = [newi,newj;newi4x,newi4y]; 
     w1 = pdist(x1,'euclidean'); 
     w2 = pdist(x2,'euclidean'); 
     w3 = pdist(x3,'euclidean'); 
     w4 = pdist(x4,'euclidean'); 
     if ceil(newi) == floor(newi) && ceil(newj)==floor(newj) % both the new coordinates are integers 
      Fw(k,(newi1x-1)*N+newi1y) = 1; 
     else if ceil(newi) == floor(newi) % one of the new coordinates is an integer 
      w = w1+w2; 
      w1new = w1/w; 
      w2new = w2/w; 
      W = w1new*w2new; 
      y1coord = (newi1x-1)*N+newi1y; 
      y2coord = (newi2x-1)*N+newi2y; 
      if y1coord <= M*N && y2coord <=M*N 
      Fw(k,y1coord) = W/w2new; 
      Fw(k,y2coord) = W/w1new; 
      end 
      else if ceil(newj) == floor(newj) % one of the new coordinates is an integer 
      w = w1+w3; 
      w1 = w1/w; 
      w3 = w3/w; 
      W = w1*w3; 
      y1coord = (newi1x-1)*N+newi1y; 
      y2coord = (newi3x-1)*N+newi3y; 
      if y1coord <= M*N && y2coord <=M*N 
       Fw(k,y1coord) = W/w3; 
       Fw(k,y2coord) = W/w1; 
      end 
      else  % both the new coordinates are not integers 
       w = w1+w2+w3+w4; 
       w1 = w1/w; 
       w2 = w2/w; 
       w3 = w3/w; 
       w4 = w4/w; 
       W = w1*w2*w3 + w2*w3*w4 + w3*w4*w1 + w4*w1*w2; 
       y1coord = (newi1x-1)*N+newi1y; 
       y2coord = (newi2x-1)*N+newi2y; 
       y3coord = (newi3x-1)*N+newi3y; 
       y4coord = (newi4x-1)*N+newi4y; 
       if y1coord <= M*N && y2coord <= M*N && y3coord <= M*N && y4coord <= M*N 
       Fw(k,y1coord) = w2*w3*w4/W; 
       Fw(k,y2coord) = w3*w4*w1/W; 
       Fw(k,y3coord) = w4*w1*w2/W; 
       Fw(k,y4coord) = w1*w2*w3/W; 
       end 
      end 
      end 
     end 
     else 
     Fw(k,k) = 1; 
     end 
     k=k+1; 
    end 
    end 
end