function muboot = bootmlreg(x,regstr,q,B)
% function muboot = bootmlreg(x,regstr,q,B)
%
% Bootstrapped Maximum Likelihood Registration (as of 8/22/05)
%
% INPUT:
%   x:          (n x m)     Data matrix
%   regstr:     (struct)    Output of MLREG
%   q:          (scalar)    Number of variance components
%   B:          (scalar)    Number of bootstrap samples
%
% OUTPUT:
%   muboot:     (B x m)     Bootstrapped estimators of the mean
%

t = regstr.t;
tt = regstr.tout;
[n,m] = size(x);
mout = length(tt);
a = t(1);
b = t(m);

% Variance component estimation
z = regstr.xregsm-repmat(regstr.mu,n,1);
[phi,d] = eigs(z'*z/n,q);
l = sqrt(diag(d))*(b-a)/m;
phi = phi/sqrt((b-a)/m);
z = z*phi*((b-a)/m)*diag(1./l);

% Warped means and components
mug = zeros(n,m);
phig = zeros(n,m,q);
for i = 1:n
    mug(i,:) = interp1(tt,regstr.mu,regstr.g(i,:),'cubic');
    for k = 1:q
        phig(i,:,k) = interp1(tt,phi(:,k),regstr.g(i,:),'cubic');
    end
end
ep = zeros(n,m);
for k = 1:q
    ep = ep + l(k)*repmat(z(:,k),1,m).*phig(:,:,k);
end
ep = x-mug-ep;

% Bootstrap
muboot = zeros(B,mout);
xx = zeros(n,m);
for r = 1:B
    disp(['Bootstrap sample ' num2str(r) ' of ' num2str(B)])
    % Resampling
    i1 = randsample(1:n,n,1);
    i2 = randsample(1:n,n,1);
    i3 = randsample(1:n,n,1);
    i4 = randsample(1:m,m,1);
    xx = zeros(n,m);
    for k = 1:q
        xx = xx + l(k)*repmat(z(i2,k),1,m).*phig(i1,:,k);
    end
    xx = mug(i1,:) + xx + ep(i3,i4);
    % Bootstrapped estimator
    regstr1 = mlreg(xx,regstr.t,regstr.tout,regstr.theta0,regstr.tau, ...
        regstr.K,regstr.h,'mean',20,-1);
    muboot(r,:) = regstr1.mu;
end