Plotting level sets of a function on a triangulated surface.

I want to plot the level sets of one function g(x,y,z) as curves on an isosurface f(x,y,z)=c. The following code comes close, but I need help getting over the finish line.

[x,y,z] = meshgrid(-1:.025:1); 
f = x.^2 +y.^2 + z.^2;
g = x.^2 - y.^2;
c = 0.8
[faces,verts,colors] = isosurface(x,y,z,f,c,g);
patch('Vertices',verts,'Faces',faces,'FaceVertexCData',colors,...
    'FaceColor','interp','EdgeColor','none')
view(3)
axis equal
colormap(turbo(10))

In this example code, we can see the level sets as the boundaries between the colored regions, but what I really want is

• Curves on a solid-colored sphere
• The ability to choose which values of g(x,y,z) get plotted
• The ability to choose line styles, colors, and thicknesses

[Roy Goodman]

A comment on my original question. The issue boils down to computing the intersection of implicitly defined surfaces. The blog "Mike on MATLAB Graphics" discussed this issue in a post in 2015.

User TimeCoder turned the ideas discussed in this post into a File Exchange submission Intersection of 2 surfaces.

Jaroslaw Tuszynski submitted another solution to the same problem to File Exchange a few years earlier.

Between these two packages and Mathieu's accepted answer, I should be able to solve my problem.

Still, If I could extract the borders between the different colors in the colored image above, I would be done. I wonder how to do this..

Answer 1

[Mathieu NOE]

hello

this is maybe not yet the perfect solution, but ....you can get your curves as 3D points , isolated in separate clusters (thank you dbscan : DBSCAN Clustering Algorithm - File Exchange - MATLAB Central )

NB that I had to modify a bit the PlotClusterinResult.m file (attached)

now remains to make it a bit prettier...

here we have slected level = 0.2 and we get two clusters of scattered data (that need now to be transformed into a smooth closed curve)

[x,y,z] = meshgrid(-1:.025:1); 
f = x.^2 +y.^2 + z.^2;
g = x.^2 - y.^2;
c = 0.8;
[faces,verts,colors] = isosurface(x,y,z,f,c,g);
figure,
patch('Vertices',verts,'Faces',faces,'FaceVertexCData',colors,...
    'FaceColor','interp','EdgeColor','none')
view(3)
axis equal
hold on 
x = verts(:,1);
y = verts(:,2);
z = verts(:,3);

colormap(turbo(10));

level = 0.2; % select level 
tol = 0.01;
% extract points that are within tolerance
ind = abs(colors - level)<tol;
xs = x(ind);
ys = y(ind);
zs = z(ind);

%% Run DBSCAN Clustering Algorithm
epsilon=0.2;
MinPts=10;
X = [xs ys zs];
IDX=DBSCAN(X,epsilon,MinPts);


%% Plot Results
PlotClusterinResult(X, IDX);
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);
colorbar('vert')

[Mathieu NOE]

Improvement is here !

from the raw scattered data we need a few more things to do before we get a nice closed curve

• put back the samples in an ordered manner using a Greedy Nearest‐Neighbor Loop : for this job you need the attached function : reorder_data.m
• smooth the result - I opted for the attached function : smooth_3d_closed_curve.m
• these functions are now called from inside the PlotClusterinResult.m function (updated version in attachment) and also , in the mean time, you can pass the line properties using params as show in the new code below

result obtained so far for level = 0.2.

the sphere is shown in transparency mode for better "contour" line vizualisation

[x,y,z] = meshgrid(-1:.025:1); 
f = x.^2 +y.^2 + z.^2;
g = x.^2 - y.^2;
c = 0.8;
[faces,verts,colors] = isosurface(x,y,z,f,c,g);
figure,
patch('Vertices',verts,'Faces',faces,'FaceVertexCData',colors,...
    'FaceColor','interp','EdgeColor','none', 'FaceAlpha', 0.5)
view(3)
axis equal
hold on 
x = verts(:,1);
y = verts(:,2);
z = verts(:,3);

colormap(turbo(10));

level = 0.2; % select level 
tol = 0.01;
% extract points that are within tolerance
ind = abs(colors - level)<tol;
xs = x(ind);
ys = y(ind);
zs = z(ind);

%% Run DBSCAN Clustering Algorithm
epsilon=0.2;
MinPts=10;
X = [xs ys zs];
IDX=DBSCAN(X,epsilon,MinPts);


%% Plot Results
params.color = 'k';
params.linewidth = 3;

PlotClusterinResult(X, IDX, params);
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);
colorbar('vert')

[Roy Goodman]

Thanks, Mathieu,

I downloaded tried to run your code, but got the following error message

Unrecognized function or variable 'euclidean'.

Error in DBSCAN (line 22)

D=euclidean(X,X);

There doesn't seem to be a function euclidean in any of the m-files you uploaded.

[Mathieu NOE]

hello Roy

sorry for that , here you are !

FYI, the original dbscan code available on the Fex uses pdist2.m but that requires you to have the corresponding toolbox (stats or machine learning)

I found this euclidean function as a good alternative (there may be others)

have a nice day

[Roy Goodman]

That looks great. Thanks. Of course I uploaded a minimal example. If I have any trouble with the real problem, I will let you know.

[Mathieu NOE]

hello again

tx for accepting my answer

BTW, I found a faster alternative to the euclidean function I sent you before

with euclidean_distance.m it goes at least 10 times faster, can be valuable for large data sets !

DBSCAN code must be updated : (really minor update)

    % compute distance
%     D = pdist2(X,X); % method 1 (original code)
%     D = euclidean(X,X);  % alternative code (ok but slow)
    D = euclidean_distance(X,X);  % fastest alternative code 

Answer 2

[Roy Goodman]

Here is the solution I eventually came up with. It uses the intersection of two surfaces code that I reference in my previous comment.

function contoursOfGonF(f,g,xyzbox,gLevels)
% Plot the surface f(x,y,z)=0 and plot level contours of g(x,y,z) on that
% surface
% 
% Input parameters
% f - the function to be plotted as a surface
% g - the function whose contours will be laid on top
% xyzbox - the plotting region [xmin xmax ymin ymax zmin zmax]
% gLevels - if a vector, a set of g-values to plot
%         - if a scalar, the number of evenly spaced g-values to plot

set(gcf,'Visible','off')
[h1, hel1] = imsurf(f, xyzbox);
xyz=h1.Vertices;
x=xyz(:,1);y=xyz(:,2);z=xyz(:,3);

gg=g(x,y,z);
gMin = min(gg); gMax=max(gg);
if isscalar(gLevels)
    nG=gLevels;
    gLevels=linspace(gMin,gMax,nG+2); gLevels=gLevels(2:end-1);
else
    nG = length(gLevels);
    if all(gLevels>gMax) || all(gLevels<gMin)
        error('no g levels on the plotted surface')
    end
end
x=cell(nG,1);y=cell(nG,1);z=cell(nG,1);

hold on;axis equal;

for k = 1:nG
    gf = @(x,y,z)g(x,y,z)-gLevels(k);
    [~, hel2] = imsurf(gf, xyzbox);
    h = intercurve(hel1, hel2);
    x{k}=h.XData;y{k}=h.YData;z{k}=h.ZData;
end

clf
[h1,~]=imsurf(f,xyzbox);
hold on
for k=1:nG
    plot3(x{k},y{k},z{k},'k','LineWidth',2)
end
set(gcf,'Visible','on')
axis equal
set(h1,'facecolor',.8*[1 1 1])
set(h1,'FaceAlpha',.9)
camlight;

An example call

contoursOfGonF(@(x,y,z)x.^2+y.^2+z.^2-1,@(x,y,z)x.^2-y.^2,1.1*[-1 1 -1 1 -1 1],.4*(-2:2));

returns the following image

[Mathieu NOE]

hello @Roy Goodman

that is very good indeed !

thumbs up !!