function orMat = orientationMatrix(cfaIn)
% orMat = orientationMatrix(cfaIn)
%
% Computes a raw orientation matrix based on the 1995 Li & Randhawa paper.
% Written by Matt Golub for Psych 221, Winter Quarter, 2008.
%
% INPUTS:  cfaIn - The original image file after CFA filtering with
%                  BayerCFA.m.
% OUTPUTS: orMat - The output orientation matrix.  1's at pixel locations
%                  represent potential vertical edges.  0's represent 
%                  potential horizontal edges.

    [rows cols chans]=size(cfaIn);

    if chans~=3
        error('Input must be the RGB output of a CFA filter (NxMx3 matrix)');
    end

    flatMat = cfaIn(:,:,1) + cfaIn(:,:,2) + cfaIn(:,:,3);

    padBottom = flatMat;
    padBottom(end+1:end+2,:) = 0;

    padTop(3:2+rows,1:cols) = flatMat;
    padTop(1:2,:) = 0;

    padLeft(1:rows,3:2+cols) = flatMat;
    padLeft(1:rows,1:2) = 0;

    padRight = flatMat;
    padRight(1:rows,end+1:end+2) = 0;

    V = abs(padTop - padBottom);
    V = V(2:end-1,:);
    V(1,:) = 0;
    V(end,:) = 0;

    H = abs(padLeft - padRight);
    H = H(:,2:end-1);
    H(:,1) = 0;
    H(:,end) = 0;

    orMat = double(V<H); %Each pixel contains a 1 if there is a 
    %potential vertical edge (ie we should use vertical spline
    %interpolation only, horizontal interpolation will cause
    %blurring of the edge).

end