It is a bit difficult to follow your code.
In my additions, these lines:
idx = find(contourData(1,:) == contourLevel);
clen = contourData(2,idx);
find and return the indices of the starting colums of the contours, and the lengths of each contour (in index units), corresponding to the contour matrix conventions. The code then uses contourf too fill the middle patch, and a loop:
axis([min(x1) max(x2) min(y1) max(y2)])
cm = 'rb';
for k = 1:numel(idx)
idxrng = (1:clen(k)-1)+idx(k);
yval = min(y1)*(k==1) + max(y2)*(k==2);
patch([contourData(1,idxrng), flip(contourData(1,idxrng))], [contourData(2,idxrng), ones(size(contourData(2,idxrng)))*yval],cm(k))
end
to fill the others.
Try this —
% Data-Example
X = [8, 12, 16, 20, 24, 28, 32];
Y = [8, 20, 32, 44, 56, 68, 80];
Z1 = [0.2795, 0.2799, 0.2814, 0.2817, 0.2800, 0.2809, 0.2808;
0.2786, 0.3040, 0.3093, 0.3113, 0.3130, 0.3131, 0.3118;
0.2824, 0.2908, 0.2975, 0.2997, 0.2992, 0.2922, 0.2930;
0.2864, 0.2847, 0.2871, 0.2798, 0.2753, 0.2714, 0.2686;
0.2850, 0.2776, 0.2802, 0.2699, 0.2684, 0.2680, 0.2681;
0.2851, 0.2830, 0.2732, 0.2684, 0.2678, 0.2684, 0.2688;
0.2857, 0.2824, 0.2747, 0.2688, 0.2668, 0.2690, 0.2687];
figure
% Given contourlevel
contourLevel = 0.28;
% Get contourlines
contourData = contourf(X, Y, Z1, [1 1]*contourLevel);
colormap(eye(3));
idx = find(contourData(1,:) == contourLevel);
clen = contourData(2,idx);
% figure
hold on
for k = 1:numel(idx)
idxrng = (1:clen(k)-1)+idx(k);
% plot(contourData(1,idxrng), contourData(2,idxrng), '.-')
[x1(k),x2(k)] = bounds(contourData(1,idxrng));
[y1(k),y2(k)] = bounds(contourData(2,idxrng));
end
axis([min(x1) max(x2) min(y1) max(y2)])
cm = 'rb';
for k = 1:numel(idx)
idxrng = (1:clen(k)-1)+idx(k);
yval = min(y1)*(k==1) + max(y2)*(k==2);
patch([contourData(1,idxrng), flip(contourData(1,idxrng))], [contourData(2,idxrng), ones(size(contourData(2,idxrng)))*yval],cm(k))
end
hold off
return
% % Get X-Y-Data of contourlines
% k = 1;
% allContourPoints = [];
% while k < size(contourData, 2)
% level = contourData(1, k);
% numPoints = contourData(2, k);
% points = contourData(:, k+1:k+numPoints);
%
% if level == contourLevel
% allContourPoints = [allContourPoints; points(1, :)', points(2, :)'];
% end
%
% k = k + numPoints + 1;
% end
%
% % Remove Duplicates
% allContourPoints = unique(allContourPoints, 'rows');
%
% % Extend X and Y (in preparation to merge contourline-Data and Z1-Data)
% extendedX = sort(unique([X(:); allContourPoints(:, 1)]));
% extendedY = sort(unique([Y(:); allContourPoints(:, 2)]));
% [gridX, gridY] = meshgrid(extendedX,extendedY);
% % Extend Z1 by interpolation
% extendedZ1 = griddata(X, Y, Z1, gridX, gridY, 'linear');
% % Make sure that Z-Data of contourlines are = contourlevel
% for i = 1:size(allContourPoints, 1)
% [~, closestXIdx] = min(abs(gridX(1, :) - allContourPoints(i, 1)));
% [~, closestYIdx] = min(abs(gridY(:, 1) - allContourPoints(i, 2)));
% disp(['Original:',num2str(extendedZ1(closestYIdx, closestXIdx))])
% extendedZ1(closestYIdx, closestXIdx) = contourLevel;
% end
%
% % Creat Binary MAtrix of extendedZ1
% extendedBinaryZ1 = extendedZ1 <= contourLevel;
%
% % Find connected areas in extendedBinaryZ1
% [B, L] = bwboundaries(extendedBinaryZ1, 'noholes');
%
% % Show contourf and Boundary-Data
% figure
% % plot contourf at contourlevel
% contourf(X,Y,Z1,[contourLevel,contourLevel])
% hold on
% % plot Boundary-Data
% for k = 1:length(B)
% b = B{k};
%
% x_values = extendedX(b(:, 2));
% y_values = extendedY(b(:, 1));
%
% switch k
% case 1
% x_color = 'rX';
% case 2
% x_color = 'gX';
% case 3
% x_color = 'bX';
% otherwise
% x_color = 'mX';
% end
% plot(x_values, y_values, x_color, 'MarkerSize', 10, 'LineWidth', 2);
% end
.
