You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
What is the best and quick way to find the intersections points set of two 3d surfaces?
9 views (last 30 days)
Show older comments
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]);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1118660/image.png)
2 Comments
Torsten
on 6 Sep 2022
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
on 7 Sep 2022
Edited: Mehdi
on 7 Sep 2022
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.
Accepted Answer
Matt J
on 7 Sep 2022
Edited: Matt J
on 8 Sep 2022
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}
19 Comments
Mehdi
on 7 Sep 2022
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);
Walter Roberson
on 7 Sep 2022
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
on 7 Sep 2022
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 =
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1119690/image.png)
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
on 7 Sep 2022
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.
Mehdi
on 8 Sep 2022
Moved: Matt J
on 8 Sep 2022
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
on 8 Sep 2022
Edited: Matt J
on 8 Sep 2022
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
on 8 Sep 2022
Edited: Mehdi
on 8 Sep 2022
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
on 8 Sep 2022
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
on 8 Sep 2022
Edited: Matt J
on 8 Sep 2022
@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
on 8 Sep 2022
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1120260/image.png)
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
on 16 Sep 2022
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'.
Input arguments for contourc must be of type 'double'.
Torsten
on 17 Sep 2022
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
on 17 Sep 2022
Edited: Mehdi
on 17 Sep 2022
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]);
More Answers (0)
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)