What is the best and quick way to find the intersections points set of two 3d surfaces?

Hi,

I am looking for the easiest best and quick way to find the intersections points of two surfaces in which the points for each curve be reported separately. As an example here there are two surfaces defined in specific domains which have two intersections. I am looking for the points on each one of these two intersection curves in seperated groups. e.g.

Curve 1:

[ x0 y0; x1 y1; x2 y2;...; xn yn]

Curve 2:

[ x0 y0; x1 y1; x2 y2;...; xn yn]

syms x y 

f = (x-1)*exp(-y^3-x^2);
g = .5-sqrt(x*exp(-x^2-y^2));
ezsurf(f,[0,2],[-2,2]);
hold on
ezsurf(g,[0,2],[-2,2]);

[Torsten]

f = @(x,y)(x-1).*exp(-y.^3-x.^2);
g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
h = @(x,y) f(x,y)-g(x,y);
fimplicit(h)

[Mehdi]

As I insisted in my question I am looking for a closed method that only reports the points of each curve separately. Your suggested method just plotted the f-g which is useless for me. In Maple there is a plots:-implicitplot command which groups the data of each curve seperately without plotting something. I am looking for similar command in Matlab, if any.

Answer 1

[Matt J]

You can to download this FEX submission,

https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates?s_tid=srchtitle

f = @(x,y)(x-1).*exp(-y.^3-x.^2);
g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
h = @(x,y) real(f(x,y)-g(x,y));

[X,Y]=deal(linspace(-6,6,1000));

[~,result]=getContourLineCoordinates(  contourc(X,Y,h(X,Y'),[0,0])   )

result = 1×2 cell array
    {996×2 double}    {327×2 double}

[Mehdi]

Interesting Method. Thanks

[Mehdi]

I mostly do not use @(x,y) in my codes, use syms instead. But here faced error. Do you know how to cope with it?

clear
syms x y 
f = (x-1)*exp(-y^3-x^2);
g = .5-sqrt(x*exp(-x^2-y^2));
h = real(f-g);
[X,Y]=deal(linspace(-6,6,1000));
[~,result]=getContourLineCoordinates(  contourc(X,Y,h(X,Y'),[0,0])   )

Index in position 1 is invalid. Array indices must be positive integers or logical values.

Error in indexing (line 1075)
R_tilde = builtin('subsref',L_tilde,Idx);

[Torsten]

Do you know how to cope with it?

By using Matt's code with function handles instead of symbolic expressions for f, g and h.

[Walter Roberson]

Your x and y are scalar symbolic variables,so your f, g, and h are scalar symbolic formula.

You try to invoke h(X,Y') which is an attempt to index the scalar symbolic formula. You have not defined a symbolic function.

You could define a symbolic function, such as

h(x,y) = real(f-g);

and then your h(X,Y') would be valid in itself. But it would give a symblic result, and contourc() only permits numeric inputs. So you would have to double() the result of calling the symbolic function h

[Matt J]

You could also convert to a non-symbolic version of the function, e.g.,

syms x y 
f = (x-1)*exp(-y^3-x^2);
g = .5-sqrt(x*exp(-x^2-y^2));
h = real(f-g)

h = 

 
hnum=str2func(vectorize(matlabFunction(h)))

hnum = function_handle with value:
    @(x,y)real(exp(-x.^2-y.^3).*(x-1.0))+real(sqrt(x.*exp(-x.^2-y.^2)))-1.0./2.0

[Torsten]

Working with real(f-g) might produce (x,y) pairs (with x < 0) that satisfy real(f(x,y)-g(x,y)) = 0, but that are not intersection points of the surfaces.

Choosing a positive search interval for x (e.g. linspace(0,6,1000) ) should be a better solution.

[Matt J]

Or maybe rewrite as below.

f = @(x,y)(x-1).*exp(-y.^3-x.^2);
g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
h = @(x,y) f(x,y)-g(x,y);

[X,Y]=deal(linspace(-6,6,1000));
Z=h(X,Y');
Z(imag(Z)~=0)=nan;

[~,result]=getContourLineCoordinates(  contourc(X,Y,Z,[0,0])   )

[Mehdi]

Thanks for the answers.

I want to put these codes in a loop which changes f sequentially. There are situations where there is not intersection between two surfaces and the program faced error. Is there any way to do if there is no intersection the program returns void ({}) or writes 'No intersection exists' instead of stopping the program?

clc
clear
for i=1:10:50
f = @(x,y)(x-i).*exp(-y.^3-x.^2);
g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
h = @(x,y) real(f(x,y)-g(x,y));
[X,Y]=deal(linspace(-6,6,1000));
[~,result]=getContourLineCoordinates(  contourc(X,Y,h(X,Y'),[0,0])   );
A{i}=result;
end

[Matt J]

for i=1:10:50
    
    f = @(x,y)(x-i).*exp(-y.^3-x.^2);
    g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
    h = @(x,y) f(x,y)-g(x,y);
    
    [X,Y]=deal(linspace(-6,6,1000));
    Z=h(X,Y');
    Z(imag(Z)~=0)=nan;
    
    M=contourc(X,Y,Z,[0,0]);
    if isempty(M), continue; end
    
    [~,result]=getContourLineCoordinates(  M  );
    A{i}=result;
    
end

[Mehdi]

As I pointed out in my question I want to find a way when there is no intersection between surfaces the program returns void ({}) or writes something such as 'No intersection exists'. Your suggested way returns anything ( {} or other thing when there is no iintersection) it just keeps running. I want to get A{7}={} or {NAN}.

[Walter Roberson]

if isempty(M)
    A{i} = Nan;
else
    [~,result]=getContourLineCoordinates(  M  );
    A{i}=result;
end

You should have been able to handle that yourself.

[Matt J]

@Mehdi Naturally, I assumed you had preallocated A. Why would I assume otherwise?

A=cell(1,numel(1:10:50)); %Preallocate

for i=1:10:50
    
    f = @(x,y)(x-i).*exp(-y.^3-x.^2);
    g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
    h = @(x,y) f(x,y)-g(x,y);
    
    [X,Y]=deal(linspace(-6,6,1000));
    Z=h(X,Y');
    Z(imag(Z)~=0)=nan;
    
    M=contourc(X,Y,Z,[0,0]);
    if isempty(M), continue; end
    
    [~,result]=getContourLineCoordinates(  M  );
    A{i}=result;
    
end

A

A = 1×5 cell array
    {1×2 cell}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

[Torsten]

Just to see what's going on.

f = @(x,y)(x-1).*exp(-y.^3-x.^2);
g = @(x,y) 0.5-sqrt(x.*exp(-x.^2-y.^2));
h = @(x,y) f(x,y)-g(x,y);
[X,Y]=deal(linspace(-6,6,1000));
Z=h(X,Y');
Z(imag(Z)~=0)=nan;
[~,result]=getContourLineCoordinates(  contourc(X,Y,Z,[0,0])   )

result = 1×2 cell array
    {996×2 double}    {327×2 double}

hold on
plot(result{1}(:,1),result{1}(:,2))
plot(result{2}(:,1),result{2}(:,2))
hold off

function [contourTable, contourArray] = getContourLineCoordinates(cm)
if ishandle(cm)
    cm = cm.ContourMatrix;
end

% Set up while loop
cmSize = size(cm,2);   	% number of columns in ContourMatrix
cmWindow = [0,0];      	% [start,end] index of moving window
contourArray = {};     	% Store the (x,y) coordinates of each contour line
% Extract coordinates of each contour line
while cmWindow(2) < cmSize
    cmWindow(1) = cmWindow(2) + 1;
    cmWindow(2) = cmWindow(2) + cm(2,cmWindow(1)) + 1;
    contourArray(end+1) = {cm(:,cmWindow(1):cmWindow(2)).'};  %#ok<AGROW>
end

% Separate the level, count, and coordinates.
level = cellfun(@(c)c(1,1),contourArray).';
numCoord = cellfun(@(c)c(1,2),contourArray).';
contourArray = cellfun(@(c)c(2:end,:),contourArray,'UniformOutput',false);

% Sort by level (just in case Matlab doesn't)
[~,sortIdx] = sort(level);

% Create a table with combined coordinates from all levels and grouping variable
levelsRep = cell2mat(arrayfun(@(v,n)repmat(v,n,1),level(sortIdx),numCoord(sortIdx),...
    'UniformOutput',false));
group = cell2mat(arrayfun(@(v,n)repmat(v,n,1),(1:numel(level)).',numCoord,...
    'UniformOutput',false));
contourTable = array2table([levelsRep, group, vertcat(contourArray{sortIdx})],...
    'VariableNames',{'Level','Group','X','Y'});

end

[Mehdi]

where is the problem?

clc
clear
syms eta__2 zeta__2
wxy2 =(6*((3*eta__2)/2 - (5*eta__2^3)/2)*((6*zeta__2)/6 - (6*zeta__2^3)/6))/6 - (7*((3*eta__2)/2 - (5*eta__2^3)/2)*((7*zeta__2^2)/7 - 7/7))/6 - (6*((6*zeta__2)/9 ));
wxy3 =wxy2- ((2*eta__2^5)/9);
wxy22= @(zeta__2, eta__2) wxy2;
wxy33= @(zeta__2, eta__2) wxy3;  
    wxy23 = @(zeta__2, eta__2) real(wxy33(zeta__2, eta__2)-wxy22(zeta__2, eta__2));
    Hvs2 = 0;
    [XX,YY]=deal(linspace(-6,6,1000));
    ZZ=wxy23(XX,YY');
    ZZ(imag(ZZ)~=0)=nan;
    Cotr=contourc(XX,YY,ZZ,[0,0]);

Error using contourc
Input arguments for contourc must be of type 'double'.

[Torsten]

clc
clear
wxy22= @(zeta__2, eta__2) (6*((3*eta__2)/2 - (5*eta__2.^3)/2)*((6*zeta__2)/6 - (6*zeta__2.^3)/6))/6 - (7*((3*eta__2)/2 - (5*eta__2.^3)/2)*((7*zeta__2.^2)/7 - 7/7))/6 - (6*((6*zeta__2)/9 ));;
wxy33= @(zeta__2, eta__2) wxy22(zeta__2, eta__2)- ((2*eta__2.^5)/9);  
    wxy23 = @(zeta__2, eta__2) real(wxy33(zeta__2, eta__2)-wxy22(zeta__2, eta__2));
    Hvs2 = 0;
    [XX,YY]=deal(linspace(-6,6,1000));
    ZZ=wxy23(XX,YY');
    ZZ(imag(ZZ)~=0)=nan;
    Cotr=contourc(XX,YY,ZZ,[0,0]);

[Mehdi]

I still have a problem when the functions are defined through summations( loops ) such as below:

clc
clear
%syms eta__2 zeta__2
wxy2=0;
wxy3=0;
for i=1:3
    for j=1:3
        wxy2=wxy2+ @(zeta__2, eta__2) legendreP(i, zeta__2)*legendreP(j, eta__2);
        wxy3=wxy2+  @(zeta__2, eta__2) ((zeta__2)^3)*legendreP(i, zeta__2)*legendreP(j, eta__2);
    end
end

Operator '+' is not supported for operands of type 'function_handle'.

wxy23 = @(zeta__2, eta__2) real(wxy3(zeta__2, eta__2)-wxy2(zeta__2, eta__2));
    [XX,YY]=deal(linspace(-0.999,0.999,1000));
    ZZ=wxy23(XX,YY');
    ZZ(imag(ZZ)~=0)=nan;
    Cotr=contourc(XX,YY,ZZ,[0,0]);

[Torsten]

You defined XX and YY as linspace(-0.999,0.999,1000).

So calculate ZZ for all values of XX x YY in a loop without using function handles for that.

And remove the syms line - it is not needed if you make numerical calculations.

[Mehdi]

"Operator '+' is not supported for operands of type 'function_handle'" error occurs before the calculation of ZZ. How to get rid of this error first?

[Torsten]

By not using function handles, as I wrote.