function [SCM] = calc_strain_cross_matrix(cpts, cpts_d, REGS, TRI)
% function [SCM] = calc_strain_cross_matrix(cpts, cpts_d, REGS)
%
% A function that takes two sets of vertices, cpts and cpts_d together with
% a cell array that has indices to which vertices make up a particular
% region of the organ. For each region we calculate the best linear
% transform that takes one set of vertices to the other. We then return a
% matrix that contains the rotation matrices for each of the regions.
%
% Inputs:
% cpts - before the ungrowth Nx2 matrix
% cpts_d - after the ungrowth Nx2 matrix
% REGS - cell array of indices to the vertices for each region
%
% Outputs:
% SCM - a 3x3xR matrix where each 3x3 element is a rotation matrix
%
% Dr. A. I. Hanna (2007)

SCM = zeros(3,3, length(REGS),1);
for i=1:length(REGS)
    % Get the points for each region
    % the points before ungrowth
    pts_ind = unique(TRI(REGS{i},:));
    
    X = cpts(pts_ind,:);
    % the points after growth
    Y = cpts_d(pts_ind,:);
    
    % Normalize for translation
    X = X - ones(size(X,1),1)*mean(X);
    Y = Y - ones(size(Y,1),1)*mean(Y); 
    
    % Order the points in a clockwise fashion so we can interpolate along
    % edges
  %  [X, ind] = order_xy_pts(X);
  %  Y = Y(ind,:);
    
    % interpolate along the edges
    X = interp_edges(X);
    Y = interp_edges(Y);
    
    % we are looking for the T that sends the points before the ungrowth to the points after
    % the ungrowth, effectively this is XT = Y
    
    % calculate the transformation matrix
    T = inv(X'*X)*(X'*Y);

    [U, S, V] = svd(T);
    T = V*S*U';
    
    de = Y - X*T; 
    e = mean(sqrt(sum(de.^2,2)));
    
    % Take the SVD of this transformation matrix in order to extract the
    % rotation angles
    [U, S, V] = svd(T);
    
    % take the difference of thse angles to recalculate the roation of the
    % regions
    SCM(1:2,1:2,i) = V*U';
end
return;
%%%
%
%%%
function P = interp_edges(pts)
P = [];
pts = [pts; pts(1,:)];
for i=1:size(pts,1)-1
   p1 = pts(i,:);
   p2 = pts(i+1,:);
   midpt = (p1+p2)/2;
   P = cat(1, P, p1);
   P = cat(1, P, midpt);
end
return;