How to find point(s) of intersection between two scatter surfaces and the zero plane?

Please help, I have searched a lot for this inquiry, but couldn't find anything to solve my issue;

I have two surfaces (green and red in image), and the zero plane (black). I need an algorithm that will find the point(s) where all three surfaces meet. I was able to create a line (black and blue in image) at surface interfaces, but was unable to find distinct points of intersection. The surfaces are created with following code

x1=eCLdata(:,1);
  y1=eCLdata(:,2);
  Z1=eCLdata(:,3);
  Z2=eCLdata(:,4);
  Z3=zeros(400,1);

    xlin=linspace(min(x1),max(x1),50);
    ylin=linspace(min(y1),max(y1),50);

    [X,Y]=meshgrid(xlin,ylin);
    f1=scatteredInterpolant(x1,y1,Z1);
    Ze=f1(X,Y);
    f2=scatteredInterpolant(x1,y1,Z2);
    Zp=f2(X,Y);
    f3=scatteredInterpolant(x1,y1,Z3);
    Z0=f3(X,Y);

So the surfaces are stored in three arrays with x,y,z values. The data was generated from very complicated procedures, so I can't work directly with the data source, I just have the scatter data to work with.

Please help?

Erin

Answer 1

[John D'Errico]

Hint: contour

1. That is, use contour to identify the level surface at 0 for each surface. The result will be a set of piecewise linear line segments. It will be a curve in the x-y plane, so the set of points where that surface is zero. That is what a level surface is, as well, what contour produces. As I recall, contourc just returns the level set, without producing a plot.

You use contourc to produce only the zero level surface by a call like this:

C = contourc(X, Y, Z,[0,0]);

This will give you two sets of curves.

2. Once you have the two sets of curves, then use a tool that will find the intersections of those curves. I like Doug Schwarz's utility, intersections, found on the file exchange.

fast-and-robust-curve-intersections

So, how would this work? Easy enough.

For example, find intersections of the surfaces, and the plane at z==0.

z1fun = @(x,y) x.^2 + 2*y.^2 - 3;
z2fun = @(x,y) sin(2*x - x.*y);

xvec = -5:.01:5;
yvec = xvec;
[x,y] = meshgrid(xvec,yvec);

z1 = z1fun(x,y);
z2 = z2fun(x,y);

contour(xvec,yvec,z1,[0 0],'r')
hold on
contour(xvec,yvec,z2,[0 0],'b')

As you can see, there should be exactly six intersections. We can root them out of the countours themselves, from:

C1 = contourc(xvec,yvec,z1,[0 0]);
C2 = contourc(xvec,yvec,z2,[0 0]);

[Erin Pratt]

John Thank you very much for the response.

I have done just as you had mentioned, but I have another additional issue.

If you can see the image, the Doug Schwarz's utility has come back with many false intersections. The heavy black and blue lines are the zero contours from the image above, but the intersection line has picked many points off of the black line that are nowhere near the blue line. Any ideas?

Thanks!

[John D'Errico]

Your image link was broken as you created it, so I fixed it.

My conjecture is: When contour sees an array with complicated discontinuous contours, it still returns them as one long array, packed together. In between each group of segments, there is an extra column that tells you the contour level, as well as the number of points in that segment. I should have gone into this in my answer. For example, in the case of the contours for C2, that I generated, we see this:

size(C2)
ans =
           2       17296

C2(:,1)
ans =
     0
    16

C2(:,18)
ans =
     0
   124

So the first contour block as 16 points in it. They are stored in columns 2:17. Then column 18 tells me there are 124 columns coming up, containing the second contour block, etc. This continues for 17296 columns, so there are MANY contours that were identified.

The solution is probably to replace those columns with NaNs, since intersections can deal with NaNs as a curve separator. Something like this:

k = 1;
C2copy = C2;
while k < size(C2,2);
  n = C2(2,k);
  C2copy(:,k) = NaN;
  k = k + n + 1;
end

When I do that on the contours created from the two surfaces I created in my answer, it does find 6 solutions, that match those I see in the plot.

[xint,yint] = intersections(C2copy(1,:),C2copy(2,:),C1(1,2:end),C1(2,2:end));

[xint,yint]
ans =
         -1.71552291085084         0.168721842093352
         -1.05785622319556        -0.969772893538178
                         0         -1.22473469387755
                         0          1.22473469387755
          1.05785622319556        -0.969772893538178
          1.71552291085084         0.168721842093352

Note that I only had to repair the C2 contours, since the C1 contour was known to be a simple ellipse.

[Erin Pratt]

John

Wow!!! That really did it. Now I only see the real intersections.

Your awesome, thanks!