function G = unrollgm(W,N,H,V,D)
% UNROLLGM Unroll a grid CRF for data creating a graphical model
%
% G = UNROLLGM(W,N,H,V) where W is the weight struct for a gcrf, N
% is the node data matrix, H and V are the horizontal and vertical
% edge data matrices, respectively. G is a graphical model struct.

        
% Jerod Weinman
% jerod@acm.org
  
%	$Id: unrollgm.m 1 2009-05-07 13:41:00Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2004, 2008, 2009
 
%    This file is part of Matlab GridCRF.
%
%    Matlab GridCRF 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.
%
%    Matlab GridCRF 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 Matlab GridCRF.  If not, see <http://www.gnu.org/licenses/>.


  
  if nargin<5 
	D = 1;
  end;
  
  L = size(W.node,2);
  
  sz = gridsize(N,H);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate adjacency matrix

  K = nodenum(sz(1), sz(2));
  E = gadjmat(K);

  S = cell(size(E));
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate node energies

  % Note: notes are ordered in the grid with their standard matlab indices
  
  % Perform standard multiplication (weight dot prod)
  P = N*W.node(1:end-1,:) + ...
	  ones(sz(1)*sz(2),1)*W.node(end,:);
  
  % Transpose for inferencer
  P = P';
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate horizontal energies

  % H:  NxD
  % We: DxL^2
  % Uh: NxL^2
  We = reshape(W.hEdge,[L*L size(W.hEdge,3)])';

  % Weight dot product
  Uh = H*We(1:end-1,:) + ones(size(H,1),1)*We(end,:);
  
  [R,M] = ind2sub([sz(1) sz(2)-1],1:size(H,1));
  
  for e=1:length(R)

	k1 = K(R(e),M(e));
	k2 = K(R(e),M(e)+1);
	
	S{k1,k2} = reshape(Uh(e,:)',[L L]);
	
  end;
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate vertical energies

  % V: NxD
  % We: DxL^2
  % Uv: NxL^2

  We = reshape(W.vEdge,[L*L size(W.vEdge,3)])';
	  
  % Weight dot product
  Uv = V*We(1:end-1,:) + ones(size(V,1),1)*We(end,:);

  [R,M] = ind2sub([sz(1)-1 sz(2)],1:size(V,1));

  for e=1:length(R)
	
	k1 = K(R(e),M(e));
	k2 = K(R(e)+1,M(e));
	  
	S{k1,k2} = reshape(Uv(e,:)',[L L]);
	
  end;
  
  
  if D
	
	% Note: We could reduce the computation here by checking to see
    % if the CRF is symmetric and/or isotropic
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Calculate default horizontal energies
	
	% Node Data: NxD
	% Wex: DxL^2
	% Un: NxL^2
	We1 = reshape(W.dhEdge(:,:,:,1), ...
				  [L*L size(W.dhEdge,3)])';
	
	We2 = reshape(W.dhEdge(:,:,:,2), ...
				  [L*L size(W.dhEdge,3)])';

	% Weight dot product
	Un1 = N*We1;
	Un2 = N*We2;
	

	[R,M] = ind2sub([sz(1) sz(2)-1],1:size(H,1));
	
	for e=1:length(R)
	  
	  k1 = K(R(e),M(e));
	  k2 = K(R(e),M(e)+1);
	
				
	  S{k1,k2} = S{k1,k2} + reshape((Un1(k1,:)+Un2(k2,:))',[L L]);
	  
	end;
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Calculate default vertical energies

	% Node Data: NxD
	% Wex: DxL^2
	% Un: NxL^2
	We1 = reshape(W.dvEdge(:,:,:,1), ...
				  [L*L size(W.dvEdge,3)])';
	
	We2 = reshape(W.dvEdge(:,:,:,2), ...
				  [L*L size(W.dvEdge,3)])';

	% Weight dot product
	Un1 = N*We1;
	Un2 = N*We2;
	

	[R,M] = ind2sub([sz(1)-1 sz(2)],1:size(V,1));
  
	for e=1:length(R)
	  
	  k1 = K(R(e),M(e));
	  k2 = K(R(e)+1,M(e));
		
	  S{k1,k2} = S{k1,k2} + reshape((Un1(k1,:)+Un2(k2,:))',[L L]);

	end;
	
  end; % if D
  
  G = struct('numNodes', {K(end)}, ...
			 'adjmat', {E}, ...
			 'numLabels', {L}, ...
			 'nodeEnergies',{P}, ...
			 'edgeEnergies',{S} );
