function bsstp = bsstensorproduct(bssvals)
%BSSTENSORPRODUCT - Compute the tensor product at a point.
%
% INPUTS:
%
%    bssvals:
%       The B-spline basis function values for all of the overlappiong
%       basis functions in each dimension that are used to compute the
%       tensor product.
%
%       This is an M by NDIM matrix where M is the order of the B-spline
%       series and NDIM is the number of dimensions.
%
% OUTPUTS:
%
%    bsstp:
%
%       An M^NDIM length array containing the column oriented tensor
%       product derived from the input matrix.
%
%

%
% Copyright 2014-2018 Biometrics Northwest LLC
%

%
% Modified: 2014-11-25
% Author  : Kevin R. Gehringer
%
%    Replaced an outer product with a vector of ones and a multiplication
%    with a diagonal matrix with an outer product.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

%
% THIS FUNCTION DOES NO ERROR CHECKING. IT IS AN INTERNAL ROUTINE THAT
% SHOULD NOT BE CALLED DIRECTLY.
%

[m,ndim] = size(bssvals);

bsstp      = zeros(m^ndim,1);

bsstp(1:m) = bssvals(:,1);
mp         = m;
for i = 2:ndim
   tmp          = bsstp(1:mp)*bssvals(:,i)';
%   mpm           = mp*m;
   mp           = mp*m;
   bsstp(1:mp) = tmp(:);
%   bsstp(1:mpm) = reshape(repmat(bsstp(1:mp),1,m).*repmat(bssvals(:,i)',mp,1),mpm,1);
%   mp = mpm;
end

