function lyngby_batch_bmr(filenames, varargin)

% lyngby_batch_bmr    - Batch for Binomial mixture ROC curve
%
%       lyngby_batch_bmr(filenames, 'PropertyName', 'PropertyValue')
%     
%       Input:    filenames   ANALYZE filenames as cell array strings
% 
%       Property: LambdaCommon [ Replications | {ReplicationsMethods} ]
%                 NumberOfThresholds [ {8} | an integer ]
%                 Path         { {''} ] filepath for result files
%                 Prefix       [ {lyngby_bmr} | a string ] Filename
%                              prefix for the result files.
%                 VolMask      'Vol' structure with volume mask
%
%       This functions is a batch for the binomial mixture estimate of
%       the ROC-curve (receiver operating characteristics curve).
%
%       Input is filenames for summary images of replicated
%       experiments. Filenames should be organized in a cell array of
%       string with replications in rows. It is possible to have
%       several methods and this should be organized in columns of the
%       'filenames' input argument. 
%
%       'LambdaCommon' determines whether the mixing coefficient in
%       the binomial mixture model (lmabda) should be common only across
%       replications or common both across replications and methods.
%
%       'NumberOfThresholds' determine the number of thresholds
%       applied to the summary images. Each threshold will correspond
%       to a point on the ROC-curve. This number should be larger than
%       two. 
%
%       The output of this function is a text file with results and a
%       EPS-file with a plot of the ROC-curves. 
%
%       Ref: Genovese, Noll, Eddy, MRM 38:497+, 1997.
%
%       Example:
%         % A specific data set
%         Runs = [ 1 2 3 8 9 10 14 15 ]; 
%         Methods = { 'fir_strength' 'fg_zscore' 'cc_energy' 'ica_image1'};
%         for m = 1:length(Methods)
%           method = Methods{m};
%           k = 1;
%           for r = 1:length(Runs) 
%             filenames{r,m} = sprintf('psn_jerom25_run%02d_%s', Runs(r), Methods{m});
%             k = k + 1;
%           end
%         end
%         vmask = lyngby_read_analyze(filenames{1,1})
%         lyngby_batch_bmr(filenames, 'VolMask', vmask);
%         !gv lyngby_bmr.eps
%                 
%       See also LYNGBY, LYNGBY_PDE_BINMIX, LYNGBY_READ_ANALYZE,
%                LYNGBY_READ_VOLUMES. 
%                      
% $Id: lyngby_batch_bmr.m,v 1.5 2004/03/19 09:11:45 fnielsen Exp $


    % Check input arguments
    if nargin < 1
      error('Too few input arguments');
    end
    if ~iscellstr(filenames)
      error('Input argument for ''filenames'' should be cellstrings');
    end
    if size(filenames,1) < 4
      error(['The number of rows in filenames should be < 4, ' ...
	    'ie, there should be more the 3 replications']);
    end
    if mod(length(varargin), 2)
      error([ 'Properties should contain an equal number arguments. ' ...
	    'Both property name and property value'])
    end
    
    
    % Default values
    filename_path     = '';
    filename_prefix   = 'lyngby';
    lambdaCommon      = 'replicationsmethods';
    numberOfThresholds = 8;
    sColor = 'rgbkcmy';
    sMark  = '.ox+*sdv^';
    Vmask  = [];         
    
    
    
    % Read properties
    n = 1;
    while n < length(varargin)
      arg = lower(varargin{n}); 

      if strcmp(arg, 'lambdacommon')
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg) & ismember(arg, { 'replications' ...
	      'replicationsmethods' })
	  lambdaCommon = arg;
	else
	  error('The argument with ''LambdaCommon'' should be a string');
	end

      elseif strcmp(arg, 'path')
	n = n + 1;
	arg = varargin{n};
	if isstr(arg) 
	  filename_path = arg;
	else
	  error('The argument with ''Path'' should be a string');
	end

      elseif strcmp(arg, 'prefix')
	n = n + 1;
	arg = varargin{n};
	if isstr(arg) 
	  filename_prefix = arg;
	else
	  error('The argument with ''Prefix'' should be a string');
	end

      elseif strcmp(arg, 'volmask') 
	n = n + 1;
	arg = varargin{n};
	if isstruct(arg) 
	  Vmask = arg;
	else
	  error('Argument with ''VolMask'' should be a structure');
	end
	
      elseif strcmp(arg, 'numberofthresholds')
	n = n + 1;
	arg = varargin{n};
	if isnumeric(arg) & length(arg)==1 & arg>2
	  numberOfThresholds = arg;
	else
	  error('Argument with ''NumberOfThresholds'' should be an integer');
	end

      
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end
    
    

    % Sizes 
    replications = size(filenames, 1);
    methods      = size(filenames, 2); 
    

    switch lambdaCommon 
      
    case 'replications'
      
      % Iterate over methods
      ThetaAppend = [];
      for m = 1:methods
	
	% Load data into a 'mat' structure
	M = lyngby_read_volumes(filenames(:,m), 'Volmask', Vmask);
	
	iThresholds = round(linspace(1, size(M.matrix,2)-1, ...
	    numberOfThresholds)) * replications; 
	[Xsort, i] = sort(M.matrix(:));
	
	clear N 
	for n = 1:length(iThresholds)
	  
	  iThreshold = iThresholds(n);
	  threshold = Xsort(end-iThreshold);

	  Xthreshold = zeros(size(M.matrix));
	  Xthreshold(M.matrix>threshold) = 1;
    
	  [N(n,:),x] = hist(sum(Xthreshold), 1:replications);
	end
	
	% Estimate binomial mixture model
	Theta(:,:,m) = lyngby_pde_binmix(N);
	ThetaAppend = [ ThetaAppend Theta(:,:,m) ];
	
      end

    case 'replicationsmethods'
      
      % Iterate over methods
      j = 1;
      for m = 1:methods
	
	% Load data into a 'mat' structure
	M = lyngby_read_volumes(filenames(:,m), 'Volmask', Vmask);
	
	iThresholds = round(linspace(1, size(M.matrix,2)-1, ...
	    numberOfThresholds)) * replications; 
	[Xsort, i] = sort(M.matrix(:));
	
	for n = 1:length(iThresholds)
	  
	  iThreshold = iThresholds(n);
	  threshold = Xsort(end-iThreshold);

	  Xthreshold = zeros(size(M.matrix));
	  Xthreshold(M.matrix>threshold) = 1;
    
	  [N(j,:),x] = hist(sum(Xthreshold), 1:replications);
	  j = j + 1;
	end
      end

      % Estimate binomial mixture model
      ThetaAppend = lyngby_pde_binmix(N);
      Theta = permute(reshape(ThetaAppend, numberOfThresholds, ...
	  methods, 3), [1 3 2]); 

      
    otherwise
      error('Internal error (lambdaCommon)')
    end

    % Save parameters.
    save(fullfile(filename_path, sprintf('%s_bmr.txt', ...
	filename_prefix)), '-ascii', 'ThetaAppend') 

   
    % Construct a plot
    hFig = figure('visible', 'off');
    plot([ 0 1 ], [ 0 1], 'linestyle', '--', 'linewidth', 2)
    axis([0 1 0 1])
    xlabel('P_I, Estimated false positive rate')
    ylabel('P_A, Estimated true positive rate')
    for m=1:methods
      hold on
      c = sColor(mod(m-1,length(sColor))+1);
      h(m) = plot(Theta(:,1,m), Theta(:,2,m), [ '-' c ], 'linewidth', 2);
      plot(Theta(:,1,m), Theta(:,2,m), [ '.' c ], 'markersize', 40);
      hold off
    end

    print(hFig, '-depsc', fullfile(filename_path, ...
	sprintf('%s_bmr.eps', filename_prefix)));



    








