function [B,D] = marginals(C,X,G,V)
% MARGINALS Calculate probabilities of data 
%
% B = MARGINALS(C,X) where X is a feature image. B contains the
% corresponding probabilities with labels along the third dimension.
%
% B = MARGINALS(C,X,G) where G indicates whether to exponentiate
% the probabilities (G==1) [default] or return them in logspace (G==0).
%
% B = MARGINALS(C,X,G,V) where if V==1 the energy is computed only
% over valid areas, if V==0 [default] the returned result is the same size
% (rows and cols) as the feature image. If V==2 the result is the same size,
% but invalid areas have NaN values.
  
% SDMAXENT.MARGINALS  
    
% Jerod Weinman
% jerod@acm.org
% (c) 2008
  
% $Id: marginals.m 1 2009-05-07 13:41:00Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2007, 2008, 2011
 
%    This file is part of Matlab SDMaxEnt.
%
%    Matlab SDMaxEnt 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 SDMaxEnt 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 SDMaxEnt.  If not, see <http://www.gnu.org/licenses/>.
      
  if nargin<3
	G = 1;
  end;
  
  M = 0; % Switch to zero invalid areas
  
  if nargin<4 || V~=1
    M = exist('V','var') && V==2;
    V = 'same';
  else
    V = 'valid';
  end;
  

  
  szW = size(C.params.weights);
  sz = size(X);

  L = szW(4); % Number of labels
  H = szW(3); % Number of image channels
  
  %  B = zeros([sz(1:2)-szW(1:2)+1 L]);
  %B = zeros([sz(1:2) L]);

  if all(szW(1:2)>=12)
    dofft = 1;
    Xf = funstack( fftextend2stack(X,szW(1:2)), 'fft2' );
  else
    dofft = 0;
  end;
  
  for c=1:L
	
	uu = 0;
	
	if ~(c==L) || anyv(C.params.weights(:,:,:,L))

	  for h=1:H

        if dofft
          Wf = fft2( fftextend2(C.params.weights(:,:,h,c),sz(1:2)) );
          uu = uu + real(fftcrop2(ifft2(Xf(:,:,h).*Wf),sz(1:2),szW(1:2)));
        else
          uu = uu + conv2(X(:,:,h),C.params.weights(:,:,h,c),V);
        end;
        
      end;
      
    end;
	
    if c==0
      B = zeros([size(uu) L]);
    end;
    
	B(:,:,c) = uu+C.params.bias(c);
	
  end;
  
  mx = max(B,[],3);
  lgZ = log(sum(exp(B-mx(:,:,ones(1,L))),3)) + mx;
  
  B = B - lgZ(:,:,ones(1,L));
  
  if G
	B = exp(B);
  end;
  
  D = ceil(get(C,'windowSize')/2);

  if M
    B(1:D(1)-1,:,:) = NaN;
    B(end-D(1)+1:end,:,:) = NaN;
    B(:,1:D(2)-1,:) = NaN;
    B(:,end-D(2)+1:end,:) = NaN;
  end;
