      1function [pre] = MLS_precompute(p_i,q_i,V,alpha,size_of_img,type)
      2% PRECOMPUTATION
      3
      4%% common for all  MLS tranformation methods
      5% for each eval point v all control points need to be looped
      6
      7% new variables
      8pre.w_i = zeros(size(p_i,1),size(V,1));
      9pre.sum_w_i = zeros(size(V,1),1);
     10pre.p_ast = zeros(size(V,1),2);
     11pre.Similarity = zeros(2*size(p_i,1),2*size(V,1));
     12%pre.Rigid = zeros(2*size(p_i,1),2*size(V,1));
     13
     14%pre = cell(size_of_img);
     15for eval = 1:size(V,1)
     16    v = V(eval,:);
     17
     18    % weights w_i = 1/length_of_vect(p_i-v)^(2alpha)
     19    w_i = 1 ./ (sqrt(sum((p_i-v).^2,2)).^(2*alpha))
     20
     21    % p? and q? are weighted centroids.
     22    ss_w_i = size(w_i)
     23    ss_s_w_i = size(sum(w_i))
     24    ss_pi=size(p_i)
     25    sum_w_i = sum(w_i);
     26
     27    p_ast = sum([w_i,w_i].*p_i, 1) ./ sum_w_i;
     28    ss_p_ast = size(p_ast)
     29
     30    % p^_i for substitution
     31    p_hat_i = p_i - repmat(p_ast,size(p_i,1),1);
     32
     33
     34
     35    if strcmp(type,'Affine') || strcmp(type,'All')
     36        pre{v(1),v(2)}.Affine = precompute_Affine(p_hat_i,p_ast,w_i,v);
     37    end
     38    if strcmp(type,'Similarity') || strcmp(type,'Rigid') || strcmp(type,'All')
     39        %[pre{v(1),v(2)}.Similarity, pre{v(1),v(2)}.Rigid] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v);
     40    end
     41    % needed for q
     42    %pre{v(1),v(2)}.w_i = w_i;
     43    %pre{v(1),v(2)}.sum_w_i = sum_w_i;
     44
     45    % needed for Rigid deform
     46    %pre{v(1),v(2)}.p_ast = p_ast;
     47
     48
     49    % new pre
     50    pre.w_i(:,eval) = w_i;
     51    pre.sum_w_i(eval,1) = sum_w_i;
     52    pre.p_ast(eval,:) = p_ast;
     53
     54    [pre.Similarity(:,2*eval-1:2*eval), a] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v);
     55end
     56end

     57
     58%% Aj - scalar for each control point p at given evaluation point v
     59% -> for each v vector (Aj_1;Aj_2;Aj_3;...)
     60function Affine = precompute_Affine(p_hat_i,p_ast,w_i,v)
     61% Affine transform precompute
     62% Single scalar for each control point p -> vector
     63% ([2,2]-[1,1])*([1;1]*1*[1,1])*[1;1]
     64s = zeros(2);
     65for i = 1:length(w_i)
     66    s = s + (p_hat_i(i,:))' * w_i(i) * p_hat_i(i,:);
     67end
     68inv_sum = 1/(s(1)*s(4)-s(3)*s(2)) * [s(4) -s(3); -s(2) s(1)];
     69inv_sum = -inv(s); %added -
     70Affine = zeros(length(w_i),1);
     71for j = 1:length(w_i)
     72    Affine(j) = (v-p_ast)*inv_sum*(p_hat_i(j,:)');
     73end
     74end
     75
     76%% Ai - matrix for each control point p at given evaluation point v
     77% mu_s scalar for each evaluation point v

     78function [Similarity,Rigid] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v)
     79% Similarity transform precompute
     80
     81mu_s = 0;
     82for i = 1:length(w_i)
     83    mu_s = mu_s + w_i(i) * p_hat_i(i,:) * p_hat_i(i,:)';
     84end
     85
     86v_p_ast = (v-p_ast);
     87perpend = (-1)*(perpend_vec(v-p_ast));
     88Matrix_T = [(v-p_ast); (-1)*(perpend_vec(v-p_ast))]';
     89
     90%Similarity = zeros(2,2,length(w_i));
     91%Rigid = zeros(2,2,length(w_i));
     92
     93Similarity = zeros(2*length(w_i),2);
     94Rigid = zeros(2*length(w_i),2);
     95for i = 1:length(w_i)
     96    A_i_temp = w_i(i) ...
     97        * [p_hat_i(i,:) ; (-1)*perpend_vec(p_hat_i(i,:))]...
     98        * Matrix_T;
     99
    100    Similarity(2*i-1:2*i,:) = A_i_temp / mu_s;
    101    Rigid(2*i-1:2*i,:) = A_i_temp;
    102end
    103end
    104%% support function for precompute_SimilarityAndRigid
    105% (x,y) -> (-y,x)
    106function vec = perpend_vec(vec)
    107temp = vec(1);
    108vec(1) = (-1)*vec(2);
    109vec(2) = temp;
    110end
    111
    112
    113