function A = cellndgrid(D)
% CELLNDGRID Generation of arrays for N-D functions and interpolation.
%   A = CELLNDGRID(D) transforms the domain specified by a cell of vectors
%   x1, x2, x3, etc.  into a cell of arrays X1,X2,X3, etc. that can be
%   used for the evaluation of functions of N variables and N-D
%   interpolation.  The i-th dimension of the output array Xi are
%   copies of elements of the vector xi.
%
%   For example, to evaluate the function  x2*exp(-x1^2-x2^2-x3^3) over the
%   range  -2 < x1 < 2,  -2 < x2 < 2, -2 < x3 < 2,
%
%       [x1,x2,x3] = ndgrid(-2:.2:2, -2:.25:2, -2:.16:2);
%       z = x2 .* exp(-x1.^2 - x2.^2 - x3.^2);
%       slice(x2,x1,x3,z,[-1.2 .8 2],2,[-2 -.2])
%

% Jerod Weinman 
% jerod@acm.org


% $Id: cellndgrid.m 1 2009-05-07 13:41:00Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2008, 2009, 2011
 
%    This file is part of MatLibFactorGraph.
%
%    MatLibFactorGraph is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    MatLibFactorGraph is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with MatLibFactorGraph.  If not, see <http://www.gnu.org/licenses/>.

if nargin==0 
  error('cellndgrid:NotEnoughInputs', 'Not enough input arguments.'); 
end

nin = length(D);
nout = nin;

for i=length(D):-1:1,
  D{i} = full(D{i}); % Make sure everything is full
  siz(i) = numel(D{i});
end

if length(siz)<nout, siz = [siz ones(1,nout-length(siz))]; end

A = cell(1,nout);
for i=1:nout,
  x = D{i}(:); % Extract and reshape as a vector.
  s = siz; s(i) = []; % Remove i-th dimension
  x = reshape(x(:,ones(1,prod(s))),[length(x) s]); % Expand x
  A{i} = permute(x,[2:i 1 i+1:nout]); % Permute to i'th dimension
end
