function regstr = landreg(y,time,ctau)

%function regstr = landreg(x,t,lmk)
%
%Landmark registration (using linear interpolation)
%
%INPUT:
%   x:   data matrix (n x m; NaN's acceptable)
%   t:   time grid (1 x m)
%   lmk: landmarks (n x p; NaN's acceptable)
%
%OUTPUT:
%   Struct with the following fields:
%       t:  same as input
%       xw: registered data (n x m)
%       mu: registered mean (1 x m)
%       w:  warping functions (n x m)
%       wi: inverse warping functions (n x m)
%
%NOTE: Make sure that the (non-missing) landmarks are in a strictly increasing
%       order or you'll get an error message from INTERP1Q.
%      The same goes for the input time grid T, of course. 
%
%No external routines are used.


% Checking consistency of input variables
[n,m] = size(y);
if length(time)~=m
    disp('Dimensions of X and T are incompatible')
    regstr = [];
    return
elseif size(time,1)==m,
    time = time';
end
if size(ctau,1)~=n
    disp('Dimensions of X and LMK are incompatible')
    regstr = [];
    return
end
if max(ctau(:))>time(m) | min(ctau(:))<time(1)
    disp('Some landmarks are outside time range')
    regstr = [];
    return
end

% Mean of non-missing landmarks
p = size(ctau,2);
ctaum = zeros(1,p);
for j=1:p
    ctaum(j) = mean(ctau(isfinite(ctau(:,j)),j));
end

% Registration by linear interpolation
warp = zeros(n,m);
warpinv = zeros(n,m);
freg = zeros(n,m);
for i = 1:n
    jj = isfinite(ctau(i,:));
    warp(i,:) = interp1q([time(1) ctaum(jj) time(m)]',[time(1) ctau(i,jj) time(m)]',time')';
    warpinv(i,:) = interp1q([time(1) ctau(i,jj) time(m)]',[time(1) ctaum(jj) time(m)]',time')';
    jjj = isfinite(y(i,:));
    freg(i,:) = interp1q(time(jjj)',y(i,jjj)',warp(i,:)')';
end
f0 = mean(freg);

regstr = struct('t',time,'xw',freg,'mu',f0,'w',warp,'wi',warpinv);