function hRFEnsembleAx = plotRFEnsemble(xFilters, yFilters, splitSign, doSmoothing, rfInds, plotImage, hue, useSubplots)
%% hRFEnsembleAx = plotRFEnsemble(xFilters, yFilters, splitSign, doSmoothing, rfInds, plotImage, hue, useSubplots)
%% Get an array with all xFilters and all yFilters, get RF image from the
% cartesian product and plot the ensemble image of all RFs.
% fi

% createFullScreenFig;

normFactor = 0.2626;%[];
if ~exist('splitSign', 'var')
    splitSign = 0;
end

if ~exist('doSmoothing', 'var')
    doSmoothing = 0;
end

if exist('rfInds', 'var')
    if ~isempty(rfInds)
        xFilters = xFilters(rfInds, :);
        yFilters = yFilters(rfInds, :);
    end
end

if ~exist('plotImage', 'var')
    plotImage = 1;
    useSubplots = 0;
elseif plotImage == 0 && ~exist('useSubplots', 'var')
    useSubplots = 1;
end

if ~exist('hue', 'var')
    hue = [];
end

if doSmoothing
    xFilters = smoothdata(xFilters, 2, 'gaussian', floor(size(xFilters, 2) * 0.15));
    yFilters = smoothdata(yFilters, 2, 'gaussian', floor(size(xFilters, 2) * 0.15));
end

% This adds a nan colum for each cell to introduce a break in case of
% plotting as a single row trace several RFs.
if ~plotImage
    xFilters = [xFilters nan * ones(size(xFilters, 1), 1)];
    yFilters = [yFilters nan * ones(size(yFilters, 1), 1)];
end


nRFs = size(yFilters, 1);
% Assuming equal length tuning curves for x and y.
nXPixelsPerRF = size(xFilters, 2);
nYPixelsPerRF = size(yFilters, 2);

nRFsPerRow = ceil( sqrt( nRFs) );
nRFsPerCol = ceil(nRFs / nRFsPerRow);

xFilters = [xFilters; ones(nRFsPerCol * nRFsPerRow - nRFs, size(xFilters, 2)) * nan];
yFilters = [yFilters; ones(nRFsPerCol * nRFsPerRow - nRFs, size(yFilters, 2)) * nan];

nRFs = size(yFilters, 1);
alphaMarker = 0.7;
% nRFsPerRow = floor( sqrt( nRFs) );
% while mod( nRFs, nRFsPerRow )
%     nRFsPerRow = nRFsPerRow - 1;
% end
% nRFsPerCol = ceil(nRFs / nRFsPerRow);

if plotImage
    if splitSign
        if ~isempty(normFactor)
            allRFsCell = arrayfun(@(x) getRFIm(xFilters(x, :), yFilters(x, :), normFactor), 1:nRFs, 'UniformOutput', false);
        else
            allRFsCell = arrayfun(@(x) getRFIm(xFilters(x, :), yFilters(x, :)), 1:nRFs, 'UniformOutput', false);
        end
    else
        allRFsCell = arrayfun(@(x) yFilters(x, :)' * xFilters(x, :), 1:nRFs, 'UniformOutput', false);
    end
    allRFsArray = cell2mat(reshape(allRFsCell, [], nRFsPerRow));
    hAllRFsIm = imagesc(allRFsArray, ZeroCenteredBounds(allRFsArray, [0 100]));
    hRFEnsembleAx = gca;
    hRFEnsembleAx.XTick = linspace(0.5, nXPixelsPerRF * nRFsPerRow + 0.5, nRFsPerRow + 1);
    hRFEnsembleAx.YTick = linspace(0.5, nYPixelsPerRF * nRFsPerCol + 0.5, nRFsPerCol + 1);
    
    hRFEnsembleAx.XAxis.Visible = 'off';
    hRFEnsembleAx.YAxis.Visible = 'off';
    hRFEnsembleAx.Layer = 'top';
    
    if splitSign
        set(hRFEnsembleAx, 'xgrid', 'on', 'ygrid', 'on', 'gridlinestyle', '-', 'xcolor', 'w', 'ycolor', 'w');
    else
        set(hRFEnsembleAx, 'xgrid', 'on', 'ygrid', 'on', 'gridlinestyle', '-', 'xcolor', 'k', 'ycolor', 'k');
    end
    pbaspect([nXPixelsPerRF * nRFsPerRow, nYPixelsPerRF * nRFsPerCol, 1])
    axis tight
    % axis square
    colormap(brewermap([], 'RdBu'))
    cptFile =  'greenPurple'; %'PRGn_11'; %
    colormapG2P = cptcmap(cptFile, 'ncol', 256);
    colormap(colormapG2P);
else
    if useSubplots
        for iRow = 1: nRFsPerCol
            for jCol = 1: nRFsPerRow
                rfInd = nRFsPerRow * (iRow - 1) + jCol;
                hSubAx(rfInd) = subplot_tight(nRFsPerCol, nRFsPerRow, rfInd, [0 0]);
                hXLines = plot(xFilters(rfInd, :), '-', 'LineWidth', 2, 'MarkerSize', 6);
                hold on;
                hYLines = plot(yFilters(rfInd, :), '-', 'LineWidth', 2, 'MarkerSize', 6);
                hRFEnsembleAx = gca;
                if ~isempty(hue)
                    % Fill in the values of hue to match the extra NaN
                    % filters used to make a square plot.
                    hue = [hue ones(1, size(xFilters, 1) - numel(hue))];
                    colorMap = flipud(colormap(inferno));
                    colorMap = colorMap(ceil(size(colorMap, 1) * 0.2): end, :);
                    colors = colorMap(round(((hue - min(hue)) / range(hue) * (size(colorMap, 1) - 1)) + 1), :);
                    hXLines.Color = colors(rfInd, :);
                    hYLines.Color = alphaMarker * colors(rfInd, :) + (1 - alphaMarker) * [1 1 1];
                    hXLines.MarkerFaceColor = colors(rfInd, :);
                    hYLines.MarkerFaceColor = alphaMarker * colors(rfInd, :) + (1 - alphaMarker) * [1 1 1];                    
                    hXLines.MarkerEdgeColor = 'none';
                    hYLines.MarkerEdgeColor = 'none';
                else
                    [hXLines.Color] = deal('k');
                    [hYLines.Color] = deal('r');
                end
                hRFEnsembleAx.XAxis.Visible = 'off';
                hRFEnsembleAx.YAxis.Visible = 'off';
                axis tight
                minY = min([xFilters(:); yFilters(:)]);
                maxY = max([xFilters(:); yFilters(:)]);
                ylim([minY, maxY])
                pbaspect([nXPixelsPerRF, nYPixelsPerRF, 1])
            end
        end
        hRFEnsembleAx = hSubAx;
    else
        vertPlotOffset = max(range(xFilters(:)), range(yFilters(:)));
        hXLines = plot(reshape(xFilters', nXPixelsPerRF * nRFsPerRow, []) + ...
                    repmat((1 : nRFsPerCol), nXPixelsPerRF * nRFsPerRow, 1) * ...
                    vertPlotOffset);
        hold on;
        hYLines = plot(reshape(yFilters', nYPixelsPerRF * nRFsPerRow, [])  + ...
                 repmat((1 : nRFsPerCol), nYPixelsPerRF * nRFsPerRow, 1) * ...
                 vertPlotOffset);
        hRFEnsembleAx = gca;
        hRFEnsembleAx.XTick = linspace(0.5, nXPixelsPerRF * nRFsPerRow + 0.5, nRFsPerRow + 1);
        hRFEnsembleAx.YTick = linspace(vertPlotOffset, vertPlotOffset * nRFsPerCol, nRFsPerCol);
        [hXLines.Color] = deal('k');
        [hYLines.Color] = deal('r');
        set(hRFEnsembleAx, 'xgrid', 'on', 'ygrid', 'on', 'gridlinestyle', '-', 'xcolor', 'k', 'ycolor', 'k');
        hRFEnsembleAx.XAxis.Visible = 'off';
        hRFEnsembleAx.YAxis.Visible = 'off';
        hRFEnsembleAx.Layer = 'top';
        pbaspect([nXPixelsPerRF * nRFsPerRow, nYPixelsPerRF * nRFsPerCol, 1])
        axis tight
    end
end

% title('Cartesian product (X \times Y) of white noise spatial filters')



% axis square

end