function bsstpidx = bsstensorindex(bssidx)
%BSSTENSORINDEX - Compute the tensor product indices at a point.
%
% INPUTS:
%
%    bssidx:
%       The B-spline basis function indices 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:
%
%    bsstpidx:
%
%       An M^NDIM by NDIM length array containing the column oriented tensor
%       product indices derived from the input matrix.
%
%

%
% Copyright 2014-2018 Biometrics Northwest LLC
%

%
% Modified: 2014-11-25
% Author  : Kevin R. Gehringer
%
%    Updated the comments to corectly specify the dimensions of the output
%    as M^NDIM by NDIM rather than M^NDIM by 1.
%
%    Replaced outer products with vectors of ones with calls to REPMAT.
%    Since REPMAT is a builtin function now it is faster than the outer
%    product for replicating data.
%
% 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(bssidx);

bsstpidx        = zeros(m^ndim,ndim);
bsstpidx(1:m,1) = bssidx(:,1);
mp               = m;
for i = 2:ndim
   mpm                   = mp*m;
   tmpi                  = repmat((1:mp)',1,m);
   bsstpidx(1:mpm,1:i-1) = bsstpidx(tmpi,1:i-1);
   tmpm                  = repmat(bssidx(:,i)',mp,1);
   bsstpidx(1:mpm,i)     = tmpm(:);
   mp                    = mpm;
end
