function [U,XI,YI] = takeo_asap(V,F,b,bc)
  % TAKEO_ASAP Solve "As-Similar-As-Possible" according to "As-Rigid-As-Possible
  % Shape Manipulation" by Igarashi et al., this is their "First Step". This is
  % equivalent up to factor of 2.5 to lscm
  %
  % U = takeo_arap(V,F,b,bc)
  %
  % Inputs:
  %   V  #V by 2 list of rest domain positions
  %   F  #F by 3 list of triangle indices into V
  %   b  #b list of indices of constraint (boundary) vertices
  %   bc  #b by 2 list of constraint positions for b
  % Outputs:
  %   U  #V by 2 list of new positions
  %   X1  #F*3 by 1 list of coefficents, so that V(I,:) can be described in
  %     terms of V(J,:) and V(K,:), where I,J,K are the indices of each
  %     triangle
  %   Y1  #F*3 by 1 list of coefficents, see above
  %
  % Note about X1,Y1:
  % V(I,:) = 
  %   (V(J,:) + [XI XI].* (V(K,:) - V(J,:)) + [YI -YI].* fliplr(V(K,:)-V(J,:)))
  %
  % See also: arap, takeo_arap, lscm
  %


  % We will build the quadratic system matrix per triangle
  % First we need for each triangle to write each vertex in terms of the other
  % two vertices.

  [XI,YI] = relative_coordinates(V,F);
  % reshape coordinate into single tall column
  XI = reshape(XI,size(F,1)*3,1);
  YI = reshape(YI,size(F,1)*3,1);

  % number of vertices
  n = size(V,1);
  % only works in 2D
  assert(size(V,2) == 2)
  % number of triangles
  nt = size(F,1);
  

  % Indices of each triangle vertex, I, and its corresponding two neighbors, J
  % and K
  I = [F(:,1);F(:,2);F(:,3)];
  J = [F(:,2);F(:,3);F(:,1)];
  K = [F(:,3);F(:,1);F(:,2)];
  % rename indices to so that I1 and I2 and so on index into vertex positions
  % as single column. Namely V(I,1) = V(I1) and V(I,2) = V(I2) etc.
  I1 = I;
  I2 = I+n;
  J1 = J;
  J2 = J+n;
  K1 = K;
  K2 = K+n;

  % Construct quadratic system matrix each triplet from II,JJ,VV contains an
  % set of elements per vertex per triangle in the quadratic energy summation
  II = [ ...
    J1
    J2
    J1
    K1
    J1
    K2
    J2
    K1
    J2
    K2
    K1
    K2
    I1
    J1
    I1
    J2
    I2
    J1
    I2
    J2
    I1
    I2
    I1
    K1
    I1
    K2
    I2
    K1
    I2
    K2
    ];
  JJ = [
    J1
    J2
    K1
    J1
    K2
    J1
    K1
    J2
    K2
    J2
    K1
    K2
    J1
    I1
    J2
    I1
    J1
    I2
    J2
    I2
    I1
    I2
    K1
    I1
    K2
    I1
    K1
    I2
    K2
    I2
    ];
  VV = [ ...
    (1-2*XI+XI.^2+YI.^2)
    (1-2*XI+XI.^2+YI.^2)
    (XI-XI.^2-YI.^2)
    (XI-XI.^2-YI.^2)
    (YI)
    (YI)
    (-YI)
    (-YI)
    (XI-XI.^2-YI.^2)
    (XI-XI.^2-YI.^2)
    (XI.^2+YI.^2)
    (XI.^2+YI.^2)
    (-1+XI)
    (-1+XI)
    (YI)
    (YI)
    (-YI)
    (-YI)
    (-1+XI)
    (-1+XI)
    ones(nt*3,1)
    ones(nt*3,1)
    (-XI)
    (-XI)
    (-YI)
    (-YI)
    (YI)
    (YI)
    (-XI)
    (-XI)
    ];
  % Assembly matrix
  A = sparse(II,JJ,VV,2*n,2*n);

  % solve
  U = min_quad_with_fixed(A,zeros(2*n,1),[b b+n],bc(:));
  % reshape into columns
  U = reshape(U,n,2);

end 

%E = ...
%  sum(sum(((U(J,:) + [XI XI].* (U(K,:) - U(J,:)) + [YI -YI].* fliplr(U(K,:)-U(J,:))) - U(I,:)).^2,2));
%  E
%
%
%E = ...
%  sum( ...
%      (U(J,1) + XI.* U(K,1) - XI.*U(J,1) + YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      (XI.* U(K,1) - XI.*U(J,1) + YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%  
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      (- XI.*U(J,1) + YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%  
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      XI.^2.*U(J,1).^2 + ...
%      -2.*XI.*U(J,1).*YI.* U(K,2) + 2.*XI.*U(J,1).*YI.*U(J,2) + 2.*XI.*U(J,1).*U(I,1) + ...
%      (YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%  
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      XI.^2.*U(J,1).^2 + ...
%      -2.*XI.*U(J,1).*YI.* U(K,2) + 2.*XI.*U(J,1).*YI.*U(J,2) + 2.*XI.*U(J,1).*U(I,1) + ...
%      YI.^2.*U(K,2).^2 + ...
%      - 2.*YI.^2.*U(J,2).*U(K,2) - 2.*YI.*U(I,1).*U(K,2) + ...
%      (- YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%  
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      XI.^2.*U(J,1).^2 + ...
%      -2.*XI.*U(J,1).*YI.* U(K,2) + 2.*XI.*U(J,1).*YI.*U(J,2) + 2.*XI.*U(J,1).*U(I,1) + ...
%      YI.^2.*U(K,2).^2 + ...
%      - 2.*YI.^2.*U(J,2).*U(K,2) - 2.*YI.*U(I,1).*U(K,2) + ...
%      YI.^2.*U(J,2).^2 + ...
%      + 2.*YI.*U(I,1).*U(J,2) + ...
%      U(I,1).^2+...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%  
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      - 2.*XI.*U(J,1).^2 + ...
%      XI.^2.*U(J,1).^2 + ...
%      YI.^2.*U(J,1).^2 + ...
%      ...
%      U(J,2).^2 + ...
%      YI.^2.*U(J,2).^2 + ...
%      - 2.*XI.*U(J,2).^2 + ...
%      XI.^2.*U(J,2).^2 + ...
%      ...
%      - 2.*YI.*U(J,1).*U(J,2) +...
%      2.*XI.*YI.*U(J,1).*U(J,2) + ...
%      2.*YI.*U(J,1).*U(J,2) + ...
%      - 2.*XI.*YI.*U(J,1).*U(J,2) + ...
%      ...
%      2.*XI.*U(J,1).*U(K,1) + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + ...
%      - 2.*YI.^2.*U(J,1).*U(K,1) + ...
%      ...
%      2.*YI.*U(J,1).*U(K,2) + ...
%      -2.*XI.*YI.*U(J,1).* U(K,2) + ...
%      2.*XI.*YI.*U(J,1).*U(K,2) + ...
%      ...
%      - 2.*XI.*YI.*U(J,2).*U(K,1) + ...
%      - 2.*YI.*U(J,2).*U(K,1) + ...
%      2.*XI.*U(J,2).*YI.* U(K,1) + ...
%      ...
%      - 2.*YI.^2.*U(J,2).*U(K,2) + ...
%      2.*XI.*U(J,2).*U(K,2) + ...
%      -2.*XI.^2.*U(J,2).*U(K,2) + ...
%      ...
%      - 2.*U(I,1).*U(J,1) + ...
%      2.*XI.*U(I,1).*U(J,1) + ...
%      - 2.*YI.*U(I,2).*U(J,1) + ...
%      ...
%      2.*YI.*U(I,1).*U(J,2) + ...
%      ...
%      - 2.*U(I,2).*U(J,2) + ...
%      2.*XI.*U(I,2).*U(J,2) + ...
%      ...
%      XI.^2.*U(K,1).^2 + ...
%      YI.^2.*U(K,1).^2 + ...
%      ...
%      YI.^2.*U(K,2).^2 + ...
%      XI.^2.*U(K,2).^2 + ...
%      ...
%      U(I,1).^2+...
%      ...
%      U(I,2).^2 + ...
%      ...
%      2.*XI.*YI.*U(K,1).*U(K,2) + ...
%      - 2.*XI.*YI.*U(K,1).*U(K,2) + ...
%      ...
%      - 2.*XI.*U(I,1).*U(K,1) + ...
%      ...
%      - 2.*YI.*U(I,1).*U(K,2) + ...
%      ...
%      2.*YI.*U(I,2).*U(K,1) + ...
%      ...
%      - 2.*XI.*U(I,2).*U(K,2) + ...
%  0 ...
%  );
%
%
%E = ...
%  sum( ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,1).^2 + ...
%      ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,2).^2 + ...
%      ...
%      (2*XI-2*XI.^2-2*YI.^2).*U(J,1).*U(K,1) + ...
%      ...
%      (2.*YI).*U(J,1).*U(K,2) + ...
%      ...
%      (-2*YI).*U(J,2).*U(K,1) + ...
%      ...
%      (2*XI-2*XI.^2-2*YI.^2).*U(J,2).*U(K,2) + ...
%      ...
%      (-2+2*XI).*U(I,1).*U(J,1) + ...
%      ...
%      (2.*YI).*U(I,1).*U(J,2) + ...
%      ...
%      (-2*YI).*U(I,2).*U(J,1) + ...
%      ...
%      (-2+2*XI).*U(I,2).*U(J,2) + ...
%      ...
%      (XI.^2+YI.^2).*U(K,1).^2 + ...
%      ...
%      (XI.^2+YI.^2).*U(K,2).^2 + ...
%      ...
%      (1).*U(I,1).^2+...
%      ...
%      (1).*U(I,2).^2 + ...
%      ...
%      (-2*XI).*U(I,1).*U(K,1) + ...
%      ...
%      (-2*YI).*U(I,1).*U(K,2) + ...
%      ...
%      (2*YI).*U(I,2).*U(K,1) + ...
%      ...
%      (-2*XI).*U(I,2).*U(K,2) + ...
%  0 ...
%  );
%
%E = ...
%  sum( ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,1).^2 + ...
%      ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,2).^2 + ...
%      ...
%      2*(XI-XI.^2-YI.^2).*U(J,1).*U(K,1) + ...
%      ...
%      2*(YI).*U(J,1).*U(K,2) + ...
%      ...
%      2*(-YI).*U(J,2).*U(K,1) + ...
%      ...
%      (2*XI-2*XI.^2-2*YI.^2).*U(J,2).*U(K,2) + ...
%      ...
%      2*(-1+XI).*U(I,1).*U(J,1) + ...
%      ...
%      2*(YI).*U(I,1).*U(J,2) + ...
%      ...
%      2*(-YI).*U(I,2).*U(J,1) + ...
%      ...
%      2*(-1+XI).*U(I,2).*U(J,2) + ...
%      ...
%      (XI.^2+YI.^2).*U(K,1).^2 + ...
%      ...
%      (XI.^2+YI.^2).*U(K,2).^2 + ...
%      ...
%      (1).*U(I,1).^2+...
%      ...
%      (1).*U(I,2).^2 + ...
%      ...
%      2*(-XI).*U(I,1).*U(K,1) + ...
%      ...
%      2*(-YI).*U(I,1).*U(K,2) + ...
%      ...
%      2*(YI).*U(I,2).*U(K,1) + ...
%      ...
%      2*(-XI).*U(I,2).*U(K,2) + ...
%  0 ...
%  );
