function [output, mse, psnr] = LiRandhawa_DirectionSelectedOutput(interpH, interpV, omV, omH, X, Xcfa_NEG)
% [output, mse, psnr] = LiRandhawa_DirectionSelectedOutput
%                               (interpH,interpV, omV, omH, X, Xcfa_NEG)
%
% Construct final demosaicked output image using Li & Randhawa orientation maps
% to select the vertical or horizontal interpolation values of the 
% green channel.
% Written by Matt Golub for Psych 221, Winter Quarter, 2008.
%
%   INPUTS:    interpH -  An interpolated RGB image formed by LiHandhawa_interp.m
%                         where the green channel was interpolated across
%                         horizontal rows.
%              interpV -  An interpolated RGB image formed by LiHandhawa_interp.m
%                         where the green channel was interpolated down
%                         vertical columns.
%              omV -      An orientation map containing ones at pixels where vertically 
%                         interpolated values should be used.
%              omH -      An orientation map containing ones at pixels where horizontally 
%                         interpolated values should be used.  This may simply be double(~omV).
%              X -        the original RGB image before the CFA array was
%                         applied.  Scaled from 0-1.
%              Xcfa_NEG - The RGB image X after the application of a Bayer
%                         color filter array (CFA) with -1 values replacing
%                         filtered values. This input can be produced by bayerCFA_NEG.m.
%
%   OUTPUTS:  output -    The final output image.
%             mse -       A 1x3 vector containing the mean-square error of
%                         the output compared to X.
%             psnr -      A 1x3 vector containing the peak signal-to-noise ratio of
%                         the output compared to X.

    [height width depth] = size(interpH);

    % The red and blue channels are the same in interpV and interpH
    output(:,:,1) = interpH(:,:,1);
    output(:,:,3) = interpH(:,:,3);

    % Here we fill in the green channel
    output(:,:,2) = Xcfa_NEG(:,:,2);
    for n = 1:height
        for m = 1:width
            if (output(n,m,2)==-1)  %Found a pixel who needs an interpolated value
                if (omV(n,m)==1)  %Vertical edge, use only vertical interpolation
                    output(n,m,2) = interpV(n,m,2);
                elseif (omH(n,m)==1) %Horizontal edge, use only vertical interpolation
                    output(n,m,2) = interpH(n,m,2);
                else %Use average between vertical and horizontal interpolations
                    output(n,m,2) = .5*(interpV(n,m,2) + interpH(n,m,2));
                end
            end
        end
    end


    mse = mean( mean( (255*(X-output)) .^2) );
    psnr = 10*log10(((max(max(output*255))).^2)./mse);

end

