matthiasm@8: function [falseAlarmRate, detectionRate, area, th] = plotROC(confidence, testClass, col, varargin)
matthiasm@8: % You pass the scores and the classes, and the function returns the false
matthiasm@8: % alarm rate and the detection rate for different points across the ROC.
matthiasm@8: %
matthiasm@8: % [faR, dR] = plotROC(score, class)
matthiasm@8: %
matthiasm@8: %  faR (false alarm rate) is uniformly sampled from 0 to 1
matthiasm@8: %  dR (detection rate) is computed using the scores.
matthiasm@8: %
matthiasm@8: % class = 0 => target absent
matthiasm@8: % class = 1 => target present
matthiasm@8: %
matthiasm@8: % score is the output of the detector, or any other measure of detection.
matthiasm@8: % There is no plot unless you add a third parameter that is the color of
matthiasm@8: % the graph. For instance:
matthiasm@8: % [faR, dR] = plotROC(score, class, 'r')
matthiasm@8: %
matthiasm@8: % faR, dR are size 1x1250
matthiasm@8: 
matthiasm@8: if nargin < 3, col = []; end
matthiasm@8: [scale01] = process_options(varargin, 'scale01', 1);
matthiasm@8: 
matthiasm@8: S = rand('state');
matthiasm@8: rand('state',0);
matthiasm@8: confidence = confidence + rand(size(confidence))*10^(-10);
matthiasm@8: rand('state',S)
matthiasm@8: 
matthiasm@8: ndxAbs = find(testClass==0); % absent
matthiasm@8: ndxPres = find(testClass==1); % present
matthiasm@8: 
matthiasm@8: [th, j] = sort(confidence(ndxAbs));
matthiasm@8: th = th(fix(linspace(1, length(th), 1250))); 
matthiasm@8: 
matthiasm@8: cAbs = confidence(ndxAbs);
matthiasm@8: cPres = confidence(ndxPres);
matthiasm@8: for t=1:length(th)
matthiasm@8:   if length(ndxPres) == 0
matthiasm@8:     detectionRate(t) = 0;
matthiasm@8:   else
matthiasm@8:     detectionRate(t)  = sum(cPres>=th(t)) / length(ndxPres);
matthiasm@8:   end
matthiasm@8:   if length(ndxAbs) == 0
matthiasm@8:     falseAlarmRate(t) = 0;
matthiasm@8:   else
matthiasm@8:     falseAlarmRate(t) = sum(cAbs>=th(t)) / length(ndxAbs);
matthiasm@8:   end
matthiasm@8:   
matthiasm@8:   %detectionRate(t)  = sum(confidence(ndxPres)>=th(t)) / length(ndxPres);
matthiasm@8:   %falseAlarmRate(t) = sum(confidence(ndxAbs)>=th(t)) / length(ndxAbs);
matthiasm@8:   %detections(t)     = sum(confidence(ndxPres)>=th(t));
matthiasm@8:   %falseAlarms(t)    = sum(confidence(ndxAbs)>=th(t));
matthiasm@8: end
matthiasm@8: 
matthiasm@8: area = sum(abs(falseAlarmRate(2:end) - falseAlarmRate(1:end-1)) .* detectionRate(2:end));
matthiasm@8: 
matthiasm@8: if ~isempty(col)
matthiasm@8:     h=plot(falseAlarmRate, detectionRate, [col '-']);
matthiasm@8:     %set(h, 'linewidth', 2);
matthiasm@8:     e = 0.05;
matthiasm@8:     if scale01
matthiasm@8:       axis([0-e 1+e 0-e 1+e])
matthiasm@8:     else
matthiasm@8:       % zoom in on the top left corner
matthiasm@8:       axis([0-e 0.5+e 0.5-e 1+e])
matthiasm@8:     end
matthiasm@8:     grid on
matthiasm@8:     ylabel('detection rate')
matthiasm@8:     xlabel('false alarm rate')
matthiasm@8: end