% contourf_overlay_demo.m
% =========================================================================
% NOTE: The functionality described by this demo script is obsolete for
%       MatLab versions prior to 2014b and only relevant for newer versions
%       because of reduced functionality in these versions.
%
% This script describes how to extract selected patches created by contourf
% to overlay these patches in a second contourf plot (e.g., in the case of
% two related quantities, Q1 and Q2, but you only want to plot Q2 in regions
% where Q1 > Q1_min, WITHOUT changing the shape of the contourf plot of Q2
% =========================================================================

% To see all required steps at once go to end of script (line 186)

close all;

% define your grid
i_grd = -50:5:50;
ni = length(i_grd);
x_grd = repmat(i_grd, ni,  1);
y_grd = repmat(i_grd', 1, ni);

% define your data/quantities (Q1, Q2)
% here: all data normalized and positive => range = [0,1]
quantity = cell(1,2);
quantity{1} = abs(x_grd.*y_grd./max(max(abs(x_grd.*y_grd))));
quantity{2} = abs((x_grd+y_grd)./max(max((x_grd+y_grd))));

% define your colorbar ticks based on your data
cTicks = 0:0.05:1;

% select Q1 range to be overlaid in Q2 plot
Q1_range = [0.15, 0.2];

% pause between plotting of individual rows of subplots?
doPause = false;

% create figure filling the left half of the screen
screenSize = get(groot, 'ScreenSize');
hFig = figure('Position', screenSize.*[1,1,0.5,1]);

%% 1st row of subplots 
% plot the original data in panels 1 and 2
for iq = 1:2
    subplot(4,2,iq,'align');
    % create contourf plot    
    contourf(x_grd, y_grd, quantity{iq}, cTicks);
    colorbar;
    colormap(parula(length(cTicks)-1));
    title(sprintf('Contourf Q%i\noriginal data', iq));
end

if doPause, pause, end

%% 2nd row of subplots
% plot the same data but set Q2 to NaN in regions where Q1 is within a
% defined data range, and plot contour lines of Q1 to illustrate region
% which is meant to be removed
% set data mask
mask = quantity{1} >= Q1_range(1) & quantity{1} < Q1_range(2);
% define epsilon value for finding contour lines
eps = 1.e-8;
for iq = 1:2
    subplot(4,2,iq+2,'align');
    plotData = quantity{iq};
    if  iq==1 %quantity 1
        % plot filled contours of Q1 AND store information in 'hc' structure
        [~, hc] = contourf(x_grd, y_grd, plotData, cTicks);
        % get contour lines of regions to be removed (using 'hc')
        nCLines = 0; % line counter
        CLines = {}; % cell to be filled with lines' (x,y)-pairs
        for ir = 1:2
            % hc.ContourMatrix is a 2-by-N Matrix with the contour level and the
            % number of entries (M) in the 1st and 2nd row, respectively, followed
            % by the M (x,y)-pairs for each contour line
            ic = find(abs(hc.ContourMatrix(1,:) - Q1_range(ir))<eps);
            if  isempty(ic)
                if  (Q1_range(ir)~=cTicks(1)) && (Q1_range(ir)~=cTicks(end))
                    close(hFig);
                    error(['Selected contour line does not exist in Q1 plot: %4.2f.\n', ...
                           'Choose another contour range.'], Q1_range(ir));
                else
                    continue
                end
            end
            for i = 1:length(ic)
                nCLines = nCLines + 1;
                nc = hc.ContourMatrix(2,ic(i)); % number of (x,y)-pairs to be read
                CLines{nCLines} = hc.ContourMatrix(:,ic(i)+(1:nc));
            end
        end
    else % quantity 2
        plotData(mask) = NaN; % set masked Q2 data to NaN
        contourf(x_grd, y_grd, plotData, cTicks);
    end
    colorbar;
    colormap(parula(length(cTicks)-1));
    % plot contour lines of Q1 to indicate regions that are supposed to be removed
    hold on
    for i = 1:nCLines
        hp = plot(CLines{i}(1,:), CLines{i}(2,:), 'Color', 'r', 'LineWidth', 2);
    end
    if  iq==1
        title(sprintf('Contourf Q%i\nwith {\\color{red}contour lines: %4.2f \\leq Q1 < %4.2f}', ...
                      iq, Q1_range));
    else
        title(sprintf(['Contourf Q%i\nwith Q2 = NaN where %4.2f \\leq Q1 < %4.2f, and \n', ...
                       '{\\color{red}corresponding contour lines from Q1}'], iq, Q1_range));
        xPos = 0.5*(min(x_grd(:))+max(x_grd(:)));
        yPos = min(y_grd(:))+0.05*(max(y_grd(:))-min(y_grd(:)));
        hText = text(xPos, yPos, 'Obviously not so nice ...');
        set(hText, 'FontSize', 12, 'FontWeight', 'bold', 'FontName', 'Arial', ...
            'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', ...
            'BackgroundColor', 'w', 'EdgeColor', ' k', 'LineWidth', 1.5);
    end
end

if doPause, pause, end

%% 3rd row of subplots
% plot the original data for Q1 and Q2 and overlay the raw information of 
% the excluded Q1 patch(es) in the Q2 plot 
% NOTE: In newer MatLab (post-2014a), individual contourf patches cannot 
%       be accessed as patches anymore. Instead they are stored as 3-by-M
%       matrices describing a series of triangular shapes that define the
%       contour patches. However, the series must be subdivided into the
%       individual patches, which is a bit tedious and not documented
%       in MatLab ...
%       The extraction of (A) the basic information and (B) the individal patches
%       is done in the following.
% =========================================================================

% (A) EXTRACTION of basic information

% STEP 1: do contourf for Q1 and store the plot information in a
%         variable. Here, 'hc' is used.
subplot(4,2,5,'align');
% create contourf plot    
[~,hc] = contourf(x_grd, y_grd, quantity{1}, cTicks);
colorbar;
caxis(cTicks([1, end]));
colormap(parula(length(cTicks)-1));
hold on
for i = 1:nCLines
    plot(CLines{i}(1,:), CLines{i}(2,:), 'Color', 'r', 'LineWidth', 2);
end
title(sprintf('Contourf Q%i\nwith {\\color{red}contour lines: %4.2f \\leq Q1 < %4.2f}', iq, Q1_range));

% STEP 2: extract the information on the triangular shapes
%   'drawnow' this ensures that the information needed for the next commands
%   are available at the time of their execution.
drawnow;
%   find the indices of the levels that were excluded
%   'hc.LevelList' contains the contour levels
iExcl = hc.LevelList >= Q1_range(1) & hc.LevelList < Q1_range(2);
%   Now we first extract the raw data needed for our overlay patch
vData = {hc.FacePrims(iExcl).VertexData};
%   This is undocumented MatLab - so, what does it do?
%   'hc.FacePrims(i).VertexData':
%     - contains a 3-by-M matrix of (x,y,z)-pairs describing a series of
%       triangular shapes (i.e., a line)
%     - in 2D, row 3 contains only zeros

%   this is the color we're going to use
cols = [.95, .95, .95];

% now do the contour plot for Q2 and overlay the triangular shapes
subplot(4,2,6,'align');
% create contourf plot    
contourf(x_grd, y_grd, quantity{2}, cTicks);
colorbar;
caxis(cTicks([1, end]));
colormap(parula(length(cTicks)-1));
hold on
for i = 1:nCLines
    hp = plot(CLines{i}(1,:), CLines{i}(2,:), 'Color','r', 'LineWidth', 2);
end
for i = 1:length(vData)
    plot(vData{i}(1,:), vData{i}(2,:), 'Color', cols, 'LineWidth', 1.5);
end
title(sprintf(['Contourf Q2\nwith {\\color{red}contour lines: %4.2f \\leq Q1 < %4.2f}', ...
               '\n and overlay of series of triangular shapes'], Q1_range));
hText = text(xPos, yPos, 'Still some work to do ...');
set(hText, 'FontSize', 12, 'FontWeight', 'bold', 'FontName', 'Arial', ...
    'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', ...
    'BackgroundColor', 'w', 'EdgeColor', ' k', 'LineWidth', 1.5);

%% 4th row of subplots
% plot the original data for Q1 and Q2 and overlay the selected Q1 patch(es)
% in the Q2 plot
% NOTE: The following code describes all steps needed to overlay the selected
%       Q1 contour patch(es) in the Q2 contour plot
% =========================================================================

% (B) EXTRACTION of patches from raw data and overlay

% STEP 1: do contourf for Q1 and store the plot information in a
%         variable. Here, 'hc' is used.
subplot(4,2,7,'align');
% create contourf plot    
[~,hc] = contourf(x_grd, y_grd, quantity{1}, cTicks);
colorbar;
caxis(cTicks([1, end]));
colormap(parula(length(cTicks)-1));
hold on
for i = 1:nCLines
    plot(CLines{i}(1,:), CLines{i}(2,:), 'Color', 'r', 'LineWidth', 2);
end
title(sprintf('Contourf Q%i\nwith {\\color{red}contour lines: %4.2f \\leq Q1 < %4.2f}', iq, Q1_range));

% STEP 2: extract the information on the triangular shapes
%   'drawnow' this ensures that the information needed for the next commands
%   are available at the time of their execution.
drawnow;
%   find the indices of the levels that were excluded
%   'hc.LevelList' contains the contour levels
iExcl = hc.LevelList >= Q1_range(1) & hc.LevelList < Q1_range(2);
%   Now we first extract the raw data needed for our overlay patch
vData = {hc.FacePrims(iExcl).VertexData};
sData = {hc.FacePrims(iExcl).StripData};
%   This is undocumented MatLab - so, what does it do?
%   'hc.FacePrims(i).VertexData':
%     - contains a 3-by-M matrix of (x,y,z)-pairs describing a series of
%       triangular shapes (i.e., a line)
%     - in 2D, row 3 contains only zeros
%   'hc.FacePrims(i).StripData':
%     - contains a vector of indices pointing to starting points of continuous
%       series of triangles to fill the i-th patch
%     - the last value of 'StripData' is 1 larger than the number of columns
%       in the corresponding 'VertexData'

% STEP 3: separate the basic information in VertexData into the information
%         on the individual triangles needed to fill each selected patch
nPatches = sum(iExcl);                        % number of patches
nTriangles = NaN(1,nPatches);                 % number of triangles per patch
triangles = struct('n', NaN(1,nPatches), ...
                   'x', [], 'y', []);         % set of (x,y) points for each triangle of all selected patches

for ip = 1:nPatches
    sData{ip} = single(sData{ip});
    % Conversion of sData to 'single' integer values is required to ensure 
    % usability in combination with other integer values
    nSeries = (length(sData{ip}) - 1) * 2;
    % '-1' because last entry of sData is "ending point of last triangle + 1"
    % '*2' because you need to fill a series of 'lower' and 'upper' triangles
    triangles(ip).x = [];
    triangles(ip).y = [];
    for it = 2:2:nSeries
        ind = sData{ip}(0.5*it):sData{ip}(0.5*it+1)-1; % indices of continuous series of (x,y)-points
        for i = 1:length(ind)
            % get 'lower' triangle's points
            if  ind(i)+2<=ind(end)
                triangles(ip).x = [triangles(ip).x, vData{ip}(1,ind(i)+[0:2,0])'];
                triangles(ip).y = [triangles(ip).y, vData{ip}(2,ind(i)+[0:2,0])'];
            end
            % get 'upper' triangles
            if  ind(i)+3<=ind(end)
                triangles(ip).x = [triangles(ip).x, vData{ip}(1,ind(i)+[1:3,1])'];
                triangles(ip).y = [triangles(ip).y, vData{ip}(2,ind(i)+[1:3,1])'];
            end
        end
    end
    triangles(ip).n = size(triangles(ip).x,2);
end

%   this is the color we're going to use
cols = [.95, .95, .95];

% now do the contour plot for Q2 and overlay the triangular shapes
subplot(4,2,8,'align');
% create contourf plot    
contourf(x_grd, y_grd, quantity{2}, cTicks);
colorbar;
caxis(cTicks([1, end]));
colormap(parula(length(cTicks)-1));
hold on
% draw the Q1 patches by drawing the individual triangles
for ip = 1:nPatches
    for it = 1:triangles(ip).n
        if sum(isnan(triangles(ip).x(:,it)))>0, continue, end
        hp = patch(triangles(ip).x(:,it),triangles(ip).y(:,it), 'k');
        set(hp, 'FaceColor', cols, 'EdgeColor', 'none');
    end
end
% draw the Q1 contour lines
for i = 1:nCLines
    hp = plot(CLines{i}(1,:), CLines{i}(2,:), 'Color','r', 'LineWidth', 2);
end

title(sprintf(['Contourf Q2\nwith {\\color{red}contour lines: %4.2f \\leq Q1 < %4.2f}', ...
               '\n selected Q1 patch(es) as overlay'], Q1_range));
hText = text(xPos, yPos, 'Just want we wanted!');
set(hText, 'FontSize', 12, 'FontWeight', 'bold', 'FontName', 'Arial', ...
    'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', ...
    'BackgroundColor', 'w', 'EdgeColor', ' k', 'LineWidth', 1.5);