You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: vl_tps.m
    path of file: ./repos/2016_person_re-ID/vlfeat/toolbox/imop
    
    the code of the file until where you have to start completion: function [phi,S] = vl_tps(x1,x2,Y)
% VL_TPS  Compute the thin-plate spline basis
%   PHI=VL_TPS(X1,X2,Y) returns the basis PHI of a thin-plate spline
%   (TPS) defined on the domain X1,X2 with control points Y.
%
%   X1 and X2 are MxN matrices specifying the grid vertices.  When
%   warping images, these usually correspond to image pixels.
%
%   Y is a 2xK matrix specifying the control points, one per
%   column. Ofthen Y is a subset of the domain X1,X2, but this is not
%   required.
%
%   PHI is a (K+3)xNxM matrix, with one layer per basis element. Each
%   basis element is a function of the domain X1,X2.
%
%   [PHI,S] = VL_TPS(X1,X2,Y) additionally returns the stiffness matrix S
%   of the TPS.
%
%   See also: VL_WTPS(), VL_HELP().

    
function [phi,S] = vl_tps(x1,x2,Y)
% VL_TPS  Compute the thin-plate spline basis
%   PHI=VL_TPS(X1,X2,Y) returns the basis PHI of a thin-plate spline
%   (TPS) defined on the domain X1,X2 with control points Y.
%
%   X1 and X2 are MxN matrices specifying the grid vertices.  When
%   warping images, these usually correspond to image pixels.
%
%   Y is a 2xK matrix specifying the control points, one per
%   column. Ofthen Y is a subset of the domain X1,X2, but this is not
%   required.
%
%   PHI is a (K+3)xNxM matrix, with one layer per basis element. Each
%   basis element is a function of the domain X1,X2.
%
%   [PHI,S] = VL_TPS(X1,X2,Y) additionally returns the stiffness matrix S
%   of the TPS.
%
%   See also: VL_WTPS(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

X = [x1(:)';x2(:)'] ;

K = size(Y,2) ;
Q = size(X,2) ;
U = vl_tpsu(Y,Y) ;
L = [[ones(1,K); Y], zeros(3) ; U, ones(K,1), Y'] ;
invL = inv(L) ;

tmp = vl_tpsu(Y,X) ;
phi = invL * [ ones(1,Q) ; X(1,:) ; X(2,:) ; tmp ] ;

[M,N] = size(x1) ;
phi = reshape(phi,K+3,M,N) ;

if nargout > 1
  % See Bookstein; note that here the terms are re-arranged a bit
  invLn = invL(1:K, end