      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    aaaa = (sqrt(sum(( p_i-repmat(v,size(p_i,1),1) ).^2,2))
     20    w_i = ones(size(p_i,1),1) ./ (sqrt(sum(( p_i-repmat(v,size(p_i,1),1) ).^2,2)).^(2*alpha))
     21
     22    % p? and q? are weighted centroids.
     23    ss_w_i = size(w_i)
     24    ss_s_w_i = size(sum(w_i))
     25    ss_pi=size(p_i)
     26    sum_w_i = sum(w_i);
     27
     28    p_ast = sum([w_i,w_i].*p_i, 1) ./ sum_w_i;
     29    ss_p_ast = size(p_ast)
     30
     31    % p^_i for substitution
     32    p_hat_i = p_i - repmat(p_ast,size(p_i,1),1);
     33
     34
     35
     36    if strcmp(type,'Affine') || strcmp(type,'All')
     37        %pre{v(1),v(2)}.Affine = precompute_Affine(p_hat_i,p_ast,w_i,v);
     38    end
     39    if strcmp(type,'Similarity') || strcmp(type,'Rigid') || strcmp(type,'All')
     40        %[pre{v(1),v(2)}.Similarity, pre{v(1),v(2)}.Rigid] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v);
     41    end
     42    % needed for q
     43    %pre{v(1),v(2)}.w_i = w_i;
     44    %pre{v(1),v(2)}.sum_w_i = sum_w_i;
     45
     46    % needed for Rigid deform
     47    %pre{v(1),v(2)}.p_ast = p_ast;
     48
     49
     50    % new pre
     51    pre.w_i(:,eval) = w_i;
     52    pre.sum_w_i(eval,1) = sum_w_i;
     53    pre.p_ast(eval,:) = p_ast;
     54
     55    [pre.Similarity(:,2*eval-1:2*eval), a] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v);
     56end
     57end
     58
     59%% Aj - scalar for each control point p at given evaluation point v
     60% -> for each v vector (Aj_1;Aj_2;Aj_3;...)
     61%function Affine = precompute_Affine(p_hat_i,p_ast,w_i,v)
     62% Affine transform precompute
     63% Single scalar for each control point p -> vector
     64% ([2,2]-[1,1])*([1;1]*1*[1,1])*[1;1]
     65
     66%s = zeros(2);
     67%for i = 1:length(w_i)
     68 %   s = s + (p_hat_i(i,:))' * w_i(i) * p_hat_i(i,:);
     69% end
     70% inv_sum = 1/(s(1)*s(4)-s(3)*s(2)) * [s(4) -s(3); -s(2) s(1)];
     71% inv_sum = -inv(s); %added -
     72% Affine = zeros(length(w_i),1);
     73% for j = 1:length(w_i)
     74%     Affine(j) = (v-p_ast)*inv_sum*(p_hat_i(j,:)');
     75% end
     76% end
     77
     78%% Ai - matrix for each control point p at given evaluation point v
     79% mu_s scalar for each evaluation point v
     80function [Similarity,Rigid] = precompute_SimilarityAndRigid(p_hat_i,p_ast,w_i,v)
     81% Similarity transform precompute
     82
     83mu_s = 0;
     84for i = 1:length(w_i)
     85    mu_s = mu_s + w_i(i) * p_hat_i(i,:) * p_hat_i(i,:)';
     86end
     87
     88v_p_ast = (v-p_ast);
     89perpend = (-1)*(perpend_vec(v-p_ast));
     90Matrix_T = [(v-p_ast); (-1)*(perpend_vec(v-p_ast))]';
     91
     92%Similarity = zeros(2,2,length(w_i));
     93%Rigid = zeros(2,2,length(w_i));
     94
     95Similarity = zeros(2*length(w_i),2);
     96Rigid = zeros(2*length(w_i),2);
     97for i = 1:length(w_i)
     98    A_i_temp = w_i(i) ...
     99        * [p_hat_i(i,:) ; (-1)*perpend_vec(p_hat_i(i,:))]...
    100        * Matrix_T;
    101
    102    Similarity(2*i-1:2*i,:) = A_i_temp / mu_s;
    103    Rigid(2*i-1:2*i,:) = A_i_temp;
    104end
    105end
    106%% support function for precompute_SimilarityAndRigid
    107% (x,y) -> (-y,x)
    108function vec = perpend_vec(vec)
    109temp = vec(1);
    110vec(1) = (-1)*vec(2);
    111vec(2) = temp;
    112end
    113
    114
    115