function [output binarizedContour] = CORF(I, sigma, t)
% CORF Contour detection based on a computational model of a simple cell.
%
% VERSION 22/04/2012
% CREATED BY: George Azzopardi and Nicolai Petkov, University of Groningen,
%             Johann Bernoulli Institute for Mathematics and Computer Science, Intelligent Systems
%
%If you use this script please cite the following paper:
%   George Azzopardi and Nicolai Petkov, "A CORF computational model of a
%     simple cell that relies on LGN input outperforms the Gabor function
%     model", 2012, DOI: 10.1007/s00422-012-0486-6
% 
%   CORF achieves orientation selectivity by combining the output - at certain 
%   positions with respect to the center of the CORF operator - of center-on 
%   and center-off difference of Gaussians (DoG) functions by a weighted geometric mean. 
%
%   CORF takes as input:
%      I -> intensity image
%      sigma -> the standard deviation of the outer Gaussian function of the 
%               DoG operator
%      t -> high threshold used for hysteresis thresholding
%
%   CORF returns:
%      output -> maximum superposition of CORF responses that correspond to 12 orientations
%      binarizedContour -> A contour map that is obtaiend by first thinning
%                          output and then performs hysteresis thresholding 
%                          with a high threshold t and a low threshold that 
%                          is a fraction 0.5 of the given t.
%
%   Example: [o bc] = CORF(imread('rino.pgm'),2.5,0.3);
%
%   The image rino.pgm is taken from the RuG data set of 40 images of
%   natural scenes, which can be downloaded from:
%   http://www.cs.rug.nl/~imaging/databases/contour_database/contour_database.html

% configure CORF operator
operator = configureCORF(sigma,0.5);

% Set number of orientations
noriens = 12;
orienslist = 0:180/noriens:359;

% Preprocessing
inputImage = double(I);
if ndims(inputImage) == 3
    inputImage = double(rgb2gray(inputImage));
end
if max(inputImage(:)) > 1
    inputImage = inputImage ./ 255;
end

% Apply difference of Gaussians function.
padwidth = ceil(max(operator.params.rho));
DoGresponse(:,:,1) = applyDoG(inputImage, 0, padwidth, operator.params.sigma, operator.params.sigmaRatio, 0);        
DoGresponse(:,:,2) = -DoGresponse(:,:,1);            
DoGresponse(:,:,1) = DoGresponse(:,:,1) .* (DoGresponse(:,:,1) > 0);
DoGresponse(:,:,2) = DoGresponse(:,:,2) .* (DoGresponse(:,:,2) > 0);

% Apply the CORF oeprator at different orientations
data = [];
output = zeros(size(DoGresponse,1),size(DoGresponse,2),noriens);
for rot = 1:length(orienslist)
    rotatedOperator = operator;
    rotatedOperator.tuple(4,:) = rotatedOperator.tuple(4,:) + orienslist(rot)*pi/180;    
    [output(:,:,rot) data] = getCORFresponse(DoGresponse,rotatedOperator,data);
end
output = output(padwidth+1:end-padwidth,padwidth+1:end-padwidth,:);

% Merge output of all orientations by using maximum superposition
[output, oriensMatrix] = calc_viewimage(output,1:size(output,3), orienslist*pi/180);   

% Perform thinning
thinningOutput = thinning(output, oriensMatrix, 2);     

%Perform hysteresis thresholding
t = t * max(thinningOutput(:));
binarizedContour = hysthresh(thinningOutput, t, 0.5*t);

function [output data] = getCORFresponse(DoGresponse,operator,data)

sz = [size(DoGresponse,1),size(DoGresponse,2)];

if isempty(data)
    data.params = [];   
    data.tupleOutput = cell(0);
    data.location = [];
    data.paramsindex = 1;
    data.current.tupleOutput = zeros(sz(1),sz(2),size(operator.tuple,2)); 
end

% The weight vector is requried for the computation of the weighted
% geometric mean at the bottom of the function.
weightVector = zeros(1,size(operator.tuple,2));
sgm = max(operator.tuple(3,:)) / 3;

for i = 1:size(operator.tuple,2)                                   
    polarity = operator.tuple(1,i);
    rho = operator.tuple(3,i);
    phi = operator.tuple(4,i);
    
    weightVector(i) = (exp(-rho^2/(2*sgm*sgm)));    
    
    mem = find(ismember(data.params,[polarity rho],'rows'),1); 
    if ~isempty(mem)
        % This tuple output is obtained by appropriate shifting of another
        % tuple output that was computed for the same values of polarity,
        % sigma and rho.
        [col row] = pol2cart(phi,rho);
        shiftrow = -(data.location{mem}(1)-fix(row));
        shiftcol = data.location{mem}(2)-fix(col);
        data.current.tupleOutput(:,:,i) = circshift(data.tupleOutput{mem},[shiftrow,shiftcol]);                                                        
    else
        DoG = DoGresponse(:,:,polarity+1);                    
        data.params(data.paramsindex,:) = [polarity rho];        
        [col row] = pol2cart(phi,rho);

        r = (operator.params.d0 + operator.params.alpha*rho)/2;
        if r > 0
            smoothfilter = fspecial('gaussian',round([2*r+1,2*r+1]),r/3);
            data.current.tupleOutput(:,:,i) = conv2(DoG,smoothfilter,'same');
            data.current.tupleOutput(:,:,i) = circshift(data.current.tupleOutput(:,:,i),[fix(row),fix(-col)]);
        else
            data.current.tupleOutput(:,:,i) = DoG;
        end
        data.current.tupleOutput(:,:,i) = data.current.tupleOutput(:,:,i) .^ weightVector(i);
        
        % We use the following variables to reuse the computations obtained
        % for the same values of parameters: polarity, sigma and rho.
        data.tupleOutput{data.paramsindex} = data.current.tupleOutput(:,:,i);
        data.location{data.paramsindex} = [round(row) round(col)];
        data.paramsindex = data.paramsindex + 1;
    end      
end

% compute the weighted geometric mean
output = prod(data.current.tupleOutput,3).^(1/sum(weightVector));

function operator = configureCORF(sigma, sigmaRatio)
% configureCORF: configure a CORF operator for the given parameters.

% The parameters alpha and d0 are set as reported in the above mentioned paper
operator.params.alpha = 0.9;
operator.params.d0 = 2;

operator.params.sigma = sigma;
operator.params.sigmaRatio = sigmaRatio;

% The following are the rho values that were used in the experiments
% reported in the above mentioned paper.
if sigma >= 1 && sigma < 2.5
    operator.params.rho = [14.38 6.9796 3.0310 1.4135];
elseif sigma >= 2.5 && sigma < 4
    operator.params.rho = [3.0515 6.1992 12.6488 24.62];
elseif sigma >= 4 && sigma <= 5
    operator.params.rho = [3.3021 4.7877 9.2467 18.08 34.43];
else
    error('The value of the parameter sigma is out of bounds');
end

maxRadius = ceil(max(operator.params.rho))+1;
stimulus = zeros((2*maxRadius));
stimulus(:,1:maxRadius) = 1;
center = [maxRadius maxRadius];

%Obtain the output of DoG function for the synthetic edge stimulus
DoGresponse(:,:,1) = applyDoG(stimulus, 0, maxRadius, sigma, sigmaRatio, 1);
DoGresponse(:,:,2) = -DoGresponse(:,:,1);
DoGresponse(:,:,1) = DoGresponse(:,:,1) .* (DoGresponse(:,:,1) > 0);
DoGresponse(:,:,2) = DoGresponse(:,:,2) .* (DoGresponse(:,:,2) > 0);

operator.tuple = [];
for r = 1:length(operator.params.rho)
    [polarity rho phi] = getTuple(DoGresponse,operator.params.rho(r),center);    
    operator.tuple = [operator.tuple [polarity; repmat(sigma,1,length(polarity)); rho; phi]];    
end

function [polarity rho phi] = getTuple(DoGresponse,radius,fp)              
% Determine the values of the polar coordinates (rho,phi)

if radius <= 0
    error('Parameter radius must be greater than 0');
end

phi = []; polarity = [];
x = 1:360;

for pol = 1:size(DoGresponse,3)       
    DoG = DoGresponse(:,:,pol);
    y = DoG(sub2ind(size(DoG),round(fp(1) + radius*cos(pi/2+x*pi/180)),round(fp(2) + radius*sin(pi/2+x*pi/180))));
       
    % Threshold low values
    y(y < 0.1*max(DoGresponse(:))) = 0;   
    y = round(y*1000)/1000;
    
    BW     = bwlabel(imregionalmax(y));
    npeaks = max(BW(:));    

    for i = 1:npeaks
        phi(end+1) = mean(x(BW == i)) * pi/180;
        polarity(end+1) = pol - 1;            
    end     
end
rho = repmat(radius,1,length(phi));

function output = applyDoG(inputImage,polarity,padwidth,sigma,sigmaRatio,crop)

% pad input image
paddedInputImage = padarray(inputImage,[padwidth padwidth],'both','symmetric');

% create DoG operator
sz = size(inputImage) + padwidth + padwidth;    
g1 = fspecial('gaussian',sz,sigma);
g2 = fspecial('gaussian',sz,sigma*sigmaRatio);
if polarity == 1
    DoG = g2 - g1;  
elseif polarity == 0
    DoG = g1 - g2;
else
    error('Polarity must be either 0 (on) or 1 (off)');
end

% compute DoG
output = fftshift(ifft2(fft2(DoG,sz(1),sz(2)) .* fft2(paddedInputImage)));

if crop == 1
    output = output(padwidth+1:end-padwidth,padwidth+1:end-padwidth);
end

function [result, oriensMatrix] = calc_viewimage(matrices, dispcomb, theta)
% VERSION 14/05/04
% CREATED BY: M.B. Wieling and N. Petkov, Groningen University,
%             Department of Computer Science, Intelligent Systems
%
% CALC_VIEWIMAGE: calculates the maximum-superposition of all the matrices stored
% in MATRICES (according to the L-infinity norm). It uses only the matrices for
% which the index is entered in DISPCOMB, e.g. if DISPCOMB contains the values
% 1,2,4: only the first, second and fourth matrix contained in MATRICES are used
% for the superposition. This method also calculates the orientationmatrix (ORIENSMATRIX) 
% which stores the maximum orientation response of each point in the resulting matrix
% (RESULT) - for use in CALC_THINNING. A progressbar of the calculations is also shown. 
%   CALC_VIEWIMAGE(MATRICES, DISPCOMB, THETA) 
%   calculates the single viewing image according to the following parameters
%     MATRICES - the matrices which hold all the convolutions for each orientation
%     DISPCOMB  - the indexes of each matrix (in MATRICES) which should be used for the
%                 superposition  
%     THETA - a list of all the orientations - to create ORIENSMATRIX

% initialize values
oriensMatrix = 0;
tmpMaxConv = -Inf;
result = -Inf;
cnt1 = 1;

if (size(dispcomb,2) == 1)
  result = matrices(:,:,dispcomb(1));
else

  % calculate the superposition (L-infinity norm)
  while (cnt1 <= size(dispcomb,2))
    % calculate the maximum orientation-response in each point (based on the absolute values)
    oriensMatrixtmp1 = (abs(matrices(:,:,dispcomb(cnt1))) > tmpMaxConv) .* theta(dispcomb(cnt1));
    oriensMatrixtmp2 = (abs(matrices(:,:,dispcomb(cnt1))) <= tmpMaxConv) .* oriensMatrix;
    oriensMatrix = oriensMatrixtmp1 + oriensMatrixtmp2;
    tmpMaxConv = max(abs(matrices(:,:,dispcomb(cnt1))), tmpMaxConv);
   
    % calculate the superposition
    result = max(result,abs(matrices(:,:,dispcomb(cnt1))));
    cnt1 = cnt1 + 1;
  end
end
 
function result = thinning(matrix, oriensMatrix, method)
% VERSION 27/02/05
% CREATED BY: M.B. Wieling and N. Petkov, Groningen University,
%             Department of Computer Science, Intelligent Systems
%
% THINNING: reduces the edges to a width of 1 pixel. This is done 
%           by looking at the two pixels next to the current pixel.
%           the neighbourpixels are determined by the orientation of 
%           the current point (this is stored in the exact same location
%           as in the matrix ORIENSMATRIX). A progressbar of the
%           calculations is also shown. 
%   THINNING(MATRIX, ORIENSMATRIX, METHOD) thins the edge
%      MATRIX - the matrix which should be thinned
%      ORIENSMATRIX - the matrix which holds for every pixel of MATRIX the
%                     orientation (in the same position)
%      METHOD - the method of thinning: METHOD == 1: simple method, just
%               take the value of the nearest pixels to compare with.
%               e.g. if the orientation = 20 degrees, the points S & N are chosen
%               to compare to, if the orientation = 25 degrees the points SW & NE
%               are chosen (see below, P = current point)
%                             NW  N  NE
%                             W   P   E
%                             SW  S  SE
%               METHOD == 2: the values to compare with are calculated using
%               interpolation based on the surrounding pixels (NW, N, NE, E, SE, S, SW, W) 
%
% corrected an error: <pi := <=pi

% create the coordinate-system
h = size(matrix,1); % height of image
w = size(matrix,2); % width of image
mx = max(h,w);
[xcoords, ycoords] = meshgrid(1:mx);
xcoords = xcoords(1:h, 1:w);
ycoords = ycoords(1:h, 1:w);

% set every value between 0 and pi
oriensMatrix = mod(oriensMatrix, pi);

% add a border of zeros to 'matrix' and 'oriensMatrix' to ease calculations.
matrixB(h+2,w+2) = 0;
matrixB(2:h+1,2:w+1) = matrix(:,:);
oriensMatrixB(h+2,w+2) = 0;
oriensMatrixB(2:h+1,2:w+1) = oriensMatrix(:,:);

if (method == 1) % simple thinning
  result = 0;
  for I=2:h+1 % the rows
    for J=2:w+1 % the columns
      orien = oriensMatrixB(I,J);
      
      % calculate dx and dy - this is correct as can be seen
      % in a picture
      dx = (orien < (3/8)*pi) - (orien > (5/8)*pi); 
      dy = ((orien > (1/8)*pi) & (orien <= (1/2)*pi)) - ((orien > (1/2)*pi) & (orien < (7/8)*pi));

      % normally the same pixels would be checked if (dy and dx > 0) and (dy and dx < 0). because
      % different pixels should be checked with a gabor-orientation of 45 and 135 this difference should
      % be checked
      if (dy < 0) & (dx < 0)
       result(I,J) = ((matrixB(I,J) >= matrixB(I+dy, J+dx)) & (matrixB(I,J) >= matrixB(I-dy, J-dx))) * matrixB(I,J);  
      else
       result(I,J) = ((matrixB(I,J) >= matrixB(I-dy, J+dx)) & (matrixB(I,J) >= matrixB(I+dy, J-dx))) * matrixB(I,J);  
      end  
    end
  end
else % linear thinning
  %hb = waitbar(0,'Applying linear thinning, please wait ... (Step 6/7)'); % display a progressbar
  result = 0;
  for I=2:h+1 % the rows
    for J=2:w+1 % the columns
      orien = oriensMatrixB(I,J);
      % get the values of the surrounding pixels
      north = matrixB(I-1, J); 
      northeast = matrixB(I-1, J+1); 
      east = matrixB(I, J+1);
      southeast = matrixB(I+1, J+1);
      south = matrixB(I+1, J);
      southwest = matrixB(I+1, J-1);
      west = matrixB(I, J-1);
      northwest = matrixB(I-1, J-1);
    
      % calculate the value of the points in one line (using interpolation)
      if (orien <= (1/4)*pi)
          fraction = orien/((1/4)*pi);
          pnt1 = (1-fraction) * east + (fraction) * northeast;
          pnt2 = (1-fraction) * west + (fraction) * southwest;     
      elseif (orien <= (1/2)*pi)
          fraction = (orien-(1/4)*pi)/((1/4)*pi);
          pnt1 = (1-fraction) * northeast + (fraction) * north;  
          pnt2 = (1-fraction) * southwest + (fraction) * south;
      elseif (orien <= (3/4)*pi)
          fraction = (orien-(1/2)*pi)/((1/4)*pi);
          pnt1 = (1-fraction) * north + (fraction) * northwest;
          pnt2 = (1-fraction) * south + (fraction) * southeast;
      elseif (orien <= pi)
          fraction = (orien-(3/4)*pi)/((1/4)*pi);
          pnt1 = (1-fraction) * northwest + (fraction) * west;
          pnt2 = (1-fraction) * southeast + (fraction) * east;
      else
        orien
      end
      result(I,J) = ( (matrixB(I,J) >= pnt1) & (matrixB(I,J) >= pnt2) ) * matrixB(I,J);
    end
    %waitbar(I/h); % update the progressbar
  end
  %close(hb);
end

% removing the borders
result = result(2:h+1, 2:w+1); 

function bw = hysthresh(im, T1, T2)
% HYSTHRESH - Hysteresis thresholding
%
% Usage: bw = hysthresh(im, T1, T2)
%
% Arguments:
%             im  - image to be thresholded (assumed to be non-negative)
%             T1  - upper threshold value
%             T2  - lower threshold value
%
% Returns:
%             bw  - the thresholded image (containing values 0 or 1)
%
% Function performs hysteresis thresholding of an image.
% All pixels with values above threshold T1 are marked as edges
% All pixels that are adjacent to points that have been marked as edges
% and with values above threshold T2 are also marked as edges. Eight
% connectivity is used.
%
% It is assumed that the input image is non-negative
%
% Author: Peter Kovesi   
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk @ csse uwa edu au   http://www.csse.uwa.edu.au/~pk   
%
% December 1996  - Original version
% March    2001  - Speed improvements made (~4x)

%
% A stack (implemented as an array) is used to keep track of all the
% indices of pixels that need to be checked.
% Note: For speed the number of conditional tests have been minimised
% This results in the top and bottom edges of the image being considered to
% be connected.  This may cause some stray edges to be propagated further than 
% they should be from the top or bottom.
%

if (T2 > T1 | T2 < 0 | T1 < 0)  % Check thesholds are sensible
  error('T1 must be >= T2 and both must be >= 0 ');
end

[rows, cols] = size(im);    % Precompute some values for speed and convenience.
rc = rows*cols;
rcmr = rc - rows;
rp1 = rows+1;

bw = im(:);                 % Make image into a column vector
pix = find(bw > T1);        % Find indices of all pixels with value > T1
npix = size(pix,1);         % Find the number of pixels with value > T1

stack = zeros(rows*cols,1); % Create a stack array (that should never
                            % overflow!)

stack(1:npix) = pix;        % Put all the edge points on the stack
stp = npix;                 % set stack pointer
for k = 1:npix
    bw(pix(k)) = -1;        % mark points as edges
end


% Precompute an array, O, of index offset values that correspond to the eight 
% surrounding pixels of any point. Note that the image was transformed into
% a column vector, so if we reshape the image back to a square the indices 
% surrounding a pixel with index, n, will be:
%              n-rows-1   n-1   n+rows-1
%
%               n-rows     n     n+rows
%                     
%              n-rows+1   n+1   n+rows+1

O = [-1, 1, -rows-1, -rows, -rows+1, rows-1, rows, rows+1];

while stp ~= 0            % While the stack is not empty
    v = stack(stp);         % Pop next index off the stack
    stp = stp - 1;
    
    if v > rp1 & v < rcmr   % Prevent us from generating illegal indices
			    % Now look at surrounding pixels to see if they
                            % should be pushed onto the stack to be
                            % processed as well.
       index = O+v;	    % Calculate indices of points around this pixel.	    
       for l = 1:8
	   ind = index(l);
	   if bw(ind) > T2   % if value > T2,
	       stp = stp+1;  % push index onto the stack.
	       stack(stp) = ind;
	       bw(ind) = -1; % mark this as an edge point
	   end
       end
    end
end

bw = (bw == -1);            % Finally zero out anything that was not an edge 
bw = reshape(bw,rows,cols); % and reshape the image