      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    aaa = repmat(v,size(p_i,1),1)
     20    aaaa = p_i-aaa
     21    w_i = ones(size(p_i,1),1) ./ (sqrt(sum(( p_i-repmat(v,size(p_i,1),1) ).^2,2)).^(2*alpha))
     22
     23    % p? and q? are weighted centroids.
     24    ss_w_i = size(w_i)
     25    ss_s_w_i = size(sum(w_i))
     26    ss_pi=size(p_i)
     27    sum_w_i = sum(w_i);
     28
     29    p_ast = sum([w_i,w_i].*p_i, 1) ./ sum_w_i;
     30    ss_p_ast = size(p_ast)
     31
     32    % p^_i for substitution
     33    p_hat_i = p_i - repmat(p_ast,size(p_i,1),1);
     34
     35
     36
     37    if strcmp(type,'Affine') || strcmp(type,'All')
     38        %pre{v(1),v(2)}.Affine = precompute_Affine(p_hat_i,p_ast,w_i,v);
     39    end
     40    if strcmp(type,'Similarity') || strcmp(type,'Rigid') || strcmp(type,'All')
     41        %[pre{v(1),v(2)}.Similarity, pre{v(1),v(2)}.Rigid] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v);
     42    end
     43    % needed for q
     44    %pre{v(1),v(2)}.w_i = w_i;
     45    %pre{v(1),v(2)}.sum_w_i = sum_w_i;
     46
     47    % needed for Rigid deform
     48    %pre{v(1),v(2)}.p_ast = p_ast;
     49
     50
     51    % new pre
     52    pre.w_i(:,eval) = w_i;
     53    pre.sum_w_i(eval,1) = sum_w_i;
     54    pre.p_ast(eval,:) = p_ast;
     55
     56    [pre.Similarity(:,2*eval-1:2*eval), a] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v);
     57end
     58end
     59
     60%% Aj - scalar for each control point p at given evaluation point v
     61% -> for each v vector (Aj_1;Aj_2;Aj_3;...)
     62%function Affine = precompute_Affine(p_hat_i,p_ast,w_i,v)
     63% Affine transform precompute
     64% Single scalar for each control point p -> vector
     65% ([2,2]-[1,1])*([1;1]*1*[1,1])*[1;1]
     66
     67%s = zeros(2);
     68%for i = 1:length(w_i)

     69 %   s = s + (p_hat_i(i,:))' * w_i(i) * p_hat_i(i,:);
     70% end
     71% inv_sum = 1/(s(1)*s(4)-s(3)*s(2)) * [s(4) -s(3); -s(2) s(1)];
     72% inv_sum = -inv(s); %added -
     73% Affine = zeros(length(w_i),1);
     74% for j = 1:length(w_i)
     75%     Affine(j) = (v-p_ast)*inv_sum*(p_hat_i(j,:)');
     76% end
     77% end
     78
     79%% Ai - matrix for each control point p at given evaluation point v
     80% mu_s scalar for each evaluation point v
     81function [Similarity,Rigid] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v)
     82% Similarity transform precompute
     83
     84mu_s = 0;
     85for i = 1:length(w_i)
     86    mu_s = mu_s + w_i(i) * p_hat_i(i,:) * p_hat_i(i,:)';
     87end
     88
     89v_p_ast = (v-p_ast);
     90perpend = (-1)*(perpend_vec(v-p_ast));
     91Matrix_T = [(v-p_ast); (-1)*(perpend_vec(v-p_ast))]';
     92
     93%Similarity = zeros(2,2,length(w_i));
     94%Rigid = zeros(2,2,length(w_i));
     95
     96Similarity = zeros(2*length(w_i),2);
     97Rigid = zeros(2*length(w_i),2);
     98for i = 1:length(w_i)
     99    A_i_temp = w_i(i) ...
    100        * [p_hat_i(i,:) ; (-1)*perpend_vec(p_hat_i(i,:))]...
    101        * Matrix_T;
    102
    103    Similarity(2*i-1:2*i,:) = A_i_temp / mu_s;
    104    Rigid(2*i-1:2*i,:) = A_i_temp;
    105end
    106end
    107%% support function for precompute_SimilarityAndRigid
    108% (x,y) -> (-y,x)
    109function vec = perpend_vec(vec)
    110temp = vec(1);
    111vec(1) = (-1)*vec(2);
    112vec(2) = temp;
    113end
    114
    115
    116