function rfIm = getRFIm(xRF, yRF, varargin)
    
    if ~isempty(varargin)
        normFactor = varargin{1};
    else
        normFactor = max(abs(([xRF(:);yRF(:)])));
    end
    % Normalize the tuning curves.
    xRFArray = repmat(xRF / normFactor, fliplr(size(yRF)));
    yRFArray = repmat(yRF' / normFactor, size(xRF));
    % In case of saturation set to one.
    xRFArray(xRFArray >= 1) = 1;
    yRFArray(yRFArray >= 1) = 1;  
    xRFArray(xRFArray <= -1) = -1;
    yRFArray(yRFArray <= -1) = -1;  
    
    % Split arrays into nonnegative and nonpositive ones.
    [negYRFArray, posYRFArray] = rectifyArray(yRFArray);
    [negXRFArray, posXRFArray] = rectifyArray(xRFArray);
    
%%
    % Define image of negative RF, product of negative gives positive.
    negRF = negXRFArray .* negYRFArray;
    negIm = getSameSignRFIm(negRF, [60 80 0]);    
    % Define image of positive RF.
    posRF = posXRFArray .* posYRFArray;
    posIm = getSameSignRFIm(posRF, [160 40 0]);
    % Define image of positive X and negative Y RF, and vice versa.
    posNegIm = getSameSignRFIm(posXRFArray .* abs(negYRFArray), [0 0 0]);
    negPosIm = getSameSignRFIm(abs(negXRFArray) .* posYRFArray, [0 0 0]);
    rfIm = negIm + posIm + posNegIm + negPosIm;
    
end
%%
    
function [nonPosArray, nonNegArray] = rectifyArray(array)
    nonPosArray = array;
    nonPosArray(nonPosArray > 0) = 0;
    nonNegArray = array;
    nonNegArray(nonNegArray < 0) = 0;
end

function rfIm = getSameSignRFIm(rfArray, hclColor)
    % For hcl many colrs have artifacts in low luminances, the
    % previous choice of colors provides near optimal solution to the
    % problem, and tries to be colorblind-friendly.
    % Assign hcl values, first luminance.
    rfIm(:, :, 3) = rfArray * 100;
    % Now make the rest with hue and chroma.
    rfIm(:, :, 1) = rfArray ./ rfArray * hclColor(1);
    rfIm(:, :, 2) = rfArray ./ rfArray * hclColor(2);
    % Remove the artifacts of division by zero.
    rfIm(isnan(rfIm)) = 0;
    for i = 1: size(rfIm, 1)
        for j = 1: size(rfIm, 2)
            rfIm(i, j, :) = hcl2rgb(rfIm(i, j, 1), rfIm(i, j, 2), ...
                                    rfIm(i, j, 3)) / 255;
        end
    end
    % Set zero values to black to avoid conversion problems.
    rfIm(rfArray == 0) = 0;
end