27 views (last 30 days)

I wrote a Matlab script that computes about a hundred of contour lines which all have the same general shape but are very slightly different in dimensions (and they can sometimes intersect). You can imagine concentric circles, except that the contours don't have a regular geometric shape and can often overlap each other. I would like to extract the contour of the 'blank space' that remains inside of this chaos:

The contour lines come from Matlab's "ContourMatrix" with this code:

[~,hcontour] = contour(x,y,z,'LevelList',k);

c = hcontour.ContourMatrix;

I then wrote a loop that selects the relevant contour lines and stores them in a matrix where the first three lines contain the X,Y and Z coordinates for the first contour line, then the next three lines contain the coordinates for the second line etc...

When I plot the outcome, the problem is that no single line is entirely contained within all the others, so I can't just compute the area of all the lines and see which is the smallest. I need to compute the boundary that only contains blank space and doesn't cross any of the contour lines.

Thanks in advance!

Wick
on 3 May 2018

Just saw that you had posted a response that the noise in the small x isn't a concern. Well, I've got an overkill of overkill solutions for you then. But it's quite robust.

I've interpolated the contours to add points between the dots. I then made a 2D histogram on a fine resolution grid to make a binary decision as to where there are points. I smoothed that binary function because contour won't work on constant values in Z. I extract the contours at the edges of the binary values and determine which contour is the inside one.

For the data you sent, it takes about 3 seconds to run when the 'interp_scale' variable is set to 1000 (not including plotting time). That's a resolution of 0.11 m or so. When I cranked 'interp_scale' to 10,000 it took 90 seconds to run but that's 0.011 m resolution so probably overkill. (about 60% worse resolution in 'y').

With that said, I'm comfortable that this is pretty bullet-proof. If you find that the contour on the inside isn't smooth you can decrease the multiplier on x_bins and y_bins to generate more interpolated points per 2D grid.

I'm sure there are ways to speed this up by only considering a small section of the overall region at a time but this brute-force solution works.

Good luck!

clearvars

tic

figure(1);

clf

W = load('testset.txt');

XX = W(1:3:end,:);

YY = W(2:3:end,:);

ZZ = W(3:3:end,:);

interp_scale = 300; % additional points for every point in the original contours

x_bins = 4 * interp_scale; % about 4x interp_scale is as high as you can go without having breaks in the histogram.

y_bins = 1 * interp_scale;

XX = [XX NaN*ones(size(XX,1),1)]; % breaks the segments apart to interpolate in a minute

YY = [YY NaN*ones(size(YY,1),1)];

ZZ = [ZZ NaN*ones(size(ZZ,1),1)];

X = reshape(XX',numel(XX),1);

Y = reshape(YY',numel(YY),1);

% this would be much more elegant if you only performed it for X > 50 on a contour-by-contour basis

X = interp1(X,linspace(1,numel(X),interp_scale*numel(X)),'pchip')';

Y = interp1(Y,linspace(1,numel(Y),interp_scale*numel(Y)),'pchip')';

X = X(~isnan(X));

Y = Y(~isnan(Y));

% axis([105 126 -27 -16])

x_edges = linspace(0,450,x_bins + 1);

y_edges = linspace(-80,80,y_bins + 1);

x_centers = 0.5*(x_edges(1:end-1) + x_edges(2:end));

y_centers = 0.5*(y_edges(1:end-1) + y_edges(2:end));

N = histcounts2(X,Y,x_edges,y_edges);

N(N>0) = 1;

% there's probably a built-in smoothing algortihm that's prettier than

% this, but this works

NB = zeros(size(N,1)+2,size(N,2)+2);

NB(2:end-1,2:end-1) = N;

NB(2:end-1,2:end-1) = N + 0.25*NB(1:end-2,2:end-1) + 0.25*NB(3:end,2:end-1) + ...

0.25*NB(2:end-1,1:end-2) + 0.25*NB(2:end-1,3:end);

N = NB(2:end-1,2:end-1);

if interp_scale <= 1000 % if there are few enough elements we'll plot the intermediate stages

subplot(221)

h=plot(XX,YY,'+k',X,Y,'.');

set(h(1),'MarkerSize',7);

set(h(2),'MarkerSize',.1);

subplot(222)

h = pcolor(x_centers,y_centers',N');

set(h,'EdgeColor','none');

axis([105 126 -27 -16])

subplot(223)

[C, Ch] = contour(x_centers, y_centers, N',[0.5,0.5]);

C3 = Ch.ContourMatrix;

subplot(224)

else

C3 = contourc(x_centers, y_centers, N',[0.75,0.75]); % use contourc if you're not plotting

end

% break out the contours into x,y curves

jj = 1;

kk = 1;

while jj < size(C3,2)

Z = C3(1,jj);

L = C3(2,jj);

P(kk).X = C3(1,jj+1:jj+L); %#ok<*SAGROW>

P(kk).Y = C3(2,jj+1:jj+L); % one P for each contour curve

jj = jj + 1 + L;

kk = kk + 1;

end

% we have many contours so we have to find which one is the closest.

distance_to_profile = zeros(length(P));

for jj = 1:length(P)

% (200,0) is in the middle of the opening. Shortest distance of any

% profile to it will be the one we want.

distance_to_profile(jj) = min((P(jj).X-200).^2 + P(jj).Y.^2);

end

[~,index] = min(distance_to_profile);

profile_x = P(index).X;

profile_y = P(index).Y;

toc

plot(XX,YY,'.k',profile_x,profile_y)

Wick
on 1 May 2018

Edited: Wick
on 2 May 2018

The solution is to find a convex hull that encompasses all the points from all the contours. An example using the 'peaks' function is shown.

figure(1)

clf

%#ok<*AGROW>

clearvars

W = peaks(25);

[xx,yy] = meshgrid(linspace(-1,1,25));

subplot(121);

% this draws the contours we're going to encompass

[~,hcontour] = contour(xx,yy,W,-3:2:5);

C3 = hcontour.ContourMatrix;

hold on

% this while loop allows me to build a matrix of all the X,Y pairs form all

% the contours. There is not (to my knowledge) a built-in way to do this.

jj = 1;

X = [];

Y = [];

while jj < size(C3,2)

Z = C3(1,jj);

L = C3(2,jj);

X = [X C3(1,jj+1:jj+L)];

Y = [Y C3(2,jj+1:jj+L)];

jj = jj + 1 + L;

end

plot(X,Y,'.k','MarkerSize',8)

hold off

dt = delaunayTriangulation(X',Y');

[Outside_Index, Av] = convexHull(dt);

Outside_X = dt.Points(Outside_Index,1);

Outside_Y = dt.Points(Outside_Index,2);

subplot(122)

plot(X,Y,'.k',Outside_X,Outside_Y,'-r')

Edited to remove the 'column' function calls. They're not necessary in this demo.

Wick
on 2 May 2018

sorry, my column function is my own utility that I didn't realize was in there. It's just:

function out = column(x)

out = reshape(x,numel(x),1);

return

I misread the original request. I thought you were looking for the outer envelope, not the inner one.

I don't know of a built-in way to do that. However, I have a quick hack that would work better for your results than my demo.

Pick a point that is "inside" the final result and convert the (x,y) pairs into polar coordinates about that point as the origin. invert the radius for each point (r = 1./r) and calculate new Cartesian coordinates. The exterior convex hull of that set of points should be the same as the interior of the original set once you convert them back to the original coordinate system.

Edit: The idea of using inverse radius doesn't work. The interior envelope curvature relative to the origin may change such that you'd need a concave section of the voronoi diagram which convexHull can't find.

Adam Danz
on 1 May 2020

FYI, there are several function on the file exchange the extract the controur line coordinates.

Opportunities for recent engineering grads.

Apply Today
## 6 Comments

## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563480

⋮## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563480

## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563496

⋮## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563496

## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563502

⋮## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563502

## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563823

⋮## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563823

## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563906

⋮## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563906

## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563929

⋮## Direct link to this comment

https://la.mathworks.com/matlabcentral/answers/394374-how-can-i-compute-a-contour-bounded-by-the-combination-of-multiple-contour-lines#comment_563929

Sign in to comment.