How can I compute a contour bounded by the combination of multiple contour lines?

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]

Antoine,

I've spent some time thinking about this problem because it irked me. The trouble is, you're going to have to define some threshold for what constitutes "inside" the profile.

When defining the exterior hull of a field of points it's easy to define "outside" because you can declare that no segment between any two points in the cloud crosses your outside hull.

On the inside, you can always cut across any interior shape and end up eliminating part of the area. In your figure above, if you put a vertical line across the region where it necks down you could cut off the tail to the left. That seems obviously wrong, but when you zoom in at the spacing-between-points level, how do you decide when you should make a high aspect ratio triangle that goes further out or skip across to the next point? I don't think there's a solution that doesn't involve some sort of threshold or smoothing to accomplish what you want.

If you're willing to smooth the set of points, I'd recommend using a 2D kernel estimator to approximate a distribution function that generates your points. You can then take a contour of that function (or of its gradient or del2) at whatever level you decide is a good threshold.

I have had good luck using Botev's kde2d code for such purposes. If you upload a set of example contour values I can help you run it through this code if you'd like.

[Antoine Bridet]

Thank you very much for your help Chad!

I think that you're right. I've been trying to implement different methods for the past couple of days but it isnt't as easy as defining the 'outside' for the reason that you mentionned.

Although not impossible, we should be careful with any kind smoothing of the data because the purpose of my program is to qualify prototypes.

I'll add more information about my most promissing attempt and a set of contour values tomorrow!

[Wick]

Actually, I think a kernel estimation may work well to define a loci of points to use as origins in the inverse radius idea below. The kernel should bleed into the shape you want to keep. Then you can cycle through the process below and find local boundaries for each starting point not worrying about whether part of the shape is cut off. That will be reached from other starting points.

If you send some data I'm sure I can whip that together in no time flat.

[Antoine Bridet]

Dear Chad,

I've attached a sample data set stored in a .txt file. Its dimensions are 75x1025, with the following structure:

• Row 1 -> X Data 1
• Row 2 -> Y Data 1
• Row 3 -> Z Data 1
• Row 4 -> X Data 2
• Row 5 -> Y Data 2
• etc...

Even though I've stored the Z values, they are irrelevant to our problem so I think that you can ignore them. My best attempt so far is to use Matlab's 'polyshape' function to transform the contour lines into polygons, and then the 'intersect' function which returns the common area between two polygons. It works, but it doesn't look like a very efficient solution.

[Wick]

Yes, if you do that one contour at a time it should work. Is computational time a concern? I doubt my solution is going to be quick.

Edit: Just took a look at some of those contours. The data for X < 50 is quite noisy and doesn't lend itself to the analysis you suggest. And the data from 200 < x < 300 is too sparse to use kde2d effectively.

Working on something else. It won't be fast either.

[Antoine Bridet]

It turns out that the computation time is quite reasonnable even in a "real life" situation. The data set I sent you has 25 lines, but even ith 120 lines the program needs about 2.5 seconds to run this section which is quite acceptable (considering that I'm using a relatively modest laptop).

Edit: To add a little bit of context, the data I provided represents lenghts in meters. The X-axis data goes from 0m to about 400m and the Y-axis data goes from about -60m to +60m. The noise present in the X-Axis between 0m and say 10m isn't much of a concern and is mostly a result of the way that the measurement is done.

Answer 1

[Wick]

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)

[Antoine Bridet]

Dear Chad,

Thanks a lot for your dedication! This really solves my problem and I do agree that it seems to be more robust. I am going to test this with different data sets and sizes and see how it behaves.

Cheers!

Answer 2

[Wick]

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.

[Antoine Bridet]

Thank you for your answer. A few remarks:

• I use a loop similar to yours in order to build a matrix of the X,Y pairs. I haven't found a built-in solution either.
• What is this ' column ' function? When I try to run your code, i get the following error:

Undefined function or variable 'column'.

• I've found a ' columns ' function, is this the one you were refering to? If so, I'll install the required Database Toolbox.

EDIT: I tried to run this code simply by removing the 'column' function and I think it still works. So the result is the smallest surface that contains all the points, which is interesting but not exactly what I'm after. I'm trying to find the largest surface that doesn't contain any point. I'll try to adapt your code and I'll get back to you but if you have any suggestion I'm also interested!

[Wick]

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]

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

Answer 3

[Antoine Bridet]

Even though it might not be the definitive answer to the problem, here is my attempt at a working solution:

j_storage = dlmread('testset.txt');

j_param1 = 1;
j_param2 = 1;

testpgon = figure('Name','Test Polygons',... % Parameters for the  figure
    'rend','opengl',... % WARNING : do NOT use painters as renderer!
    'OuterPosition',[0 0 1920 1200]);

xlim([-50;500])
ylim([-100;100])
hold off

while j_param1<=size(j_storage,1)
    j_Pgon{j_param2} = polyshape(j_storage(j_param1,:),j_storage(j_param1+1,:),'Simplify',true);

      if j_param2 == 2
          j_commonArea = intersect(j_Pgon{j_param2-1},j_Pgon{j_param2});
          plot(j_commonArea)
          text(380,55,num2str(j_param2),'Color','red','FontSize',24)
      elseif j_param2 > 2
          j_commonArea = intersect(j_Pgon{j_param2},j_commonArea);
          plot(j_commonArea)
          text(380,55,num2str(j_param2),'Color','red','FontSize',24)
      elseif j_param2 == 1
          plot(j_Pgon{j_param2})
          text(380,55,num2str(j_param2),'Color','red','FontSize',24)
      end

      xlim([0;400])
      ylim([-60;60])
      drawnow

      j_param1 = j_param1 + 3;
      j_param2 = j_param2 + 1;
  end

plot(j_commonArea)
text(380,55,num2str(j_param2),'Color','red','FontSize',24)
xlim([0;400])
ylim([-60;60])

Note that my program uses a dataset in an array called 'j_storage' so to run it on your machine, you'll want to run this code with the attached file in the same folder (it is the same as the one already provided in the discussion above).