Hi, I want to how can we convert the contour geometry to stl file?

Hello,

I would like to know how can we convert 2D-geometry I obtained from isocontour of the signed distance function to stl.file.

This is because I want to generate triangular mesh for the inside of the isocontour using "generateMesh" function.

For example, I need to mesh the inside of the circular domain shown in contour.png. The isocontour is plottted using the code below.

Plus, I would be happy if you could inform me how to conduct this meshing without stl.file.

clc; clear; close all;

% Generating simple SDF
[x, y] = meshgrid(linspace(-2, 2, 100), linspace(-2, 2, 100));
sdf = sqrt(x.^2 + y.^2) - 1;  

% Extracting 0-level isocontour of the SDF.
isoValue = 0;  
contour(x(1,:), y(:,1), sdf, [isoValue isoValue]); 

Answer 1

[Star Strider]

Try something like this —

% clc; clear; close all;

% Generating simple SDF
[x, y] = meshgrid(linspace(-2, 2, 100), linspace(-2, 2, 100));
sdf = sqrt(x.^2 + y.^2) - 1;  

% Extracting 0-level isocontour of the SDF.
isoValue = 0;  
c = contour(x(1,:), y(:,1), sdf, [isoValue isoValue]); 
axis('equal')

x = c(1,2:end).';                                                                       % First Row Of ‘c’  Are The ‘X’ Coordinates, Delete The First Element From Each
y = c(2,2:end).';                                                                       % Second Row Of ‘c’ Are The ‘Y’ Coordinates, Delete The First Element From Each
DT = delaunay(x,y);                                                                     % Delaunay Triangulation

Warning: Duplicate data points have been detected and removed.
Some point indices will not be referenced by the triangulation.

DT = 198×3

74 71 65 26 24 37 15 46 24 65 62 61 196 96 15 115 108 107 78 76 83 96 196 115 165 196 174 15 21 18

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
triplot(DT, x, y)                                                                       % Plot Result
axis('equal','padded')

TR = triangulation(DT,x,y);                                                         % Create ‘triangulation’ Object

Warning: Some input points are not referenced by the triangulation.

TR = 
  triangulation with properties:

              Points: [201x2 double]
    ConnectivityList: [198x3 double]

stlwrite(TR, 'STL_Test.stl', 'text')                                                % Write Triangulation Object To STL File

TRr = stlread('STL_Test.stl')                                                       % Check: Read Triangulation Object From STL File

TRr = 
  triangulation with properties:

              Points: [200x3 double]
    ConnectivityList: [198x3 double]

boundaryEdges = freeBoundary(TRr).';                                                % Get Boundary Edges

figure
triplot(TRr)                                                                        % Plot ‘triangulation’ Objeect
hold on
plot(TRr.Points(boundaryEdges,1), TRr.Points(boundaryEdges,2), '-r', 'LineWidth',2) % Plot Boundary
hold off
axis('equal','padded')

.

[ArxIv]

Hi, thank you for your detailed answer!

If possible, could you help me for the case with changed topology ?

For example, what can I do if I only want the specific domain with -0.5<=SDF<=0 (domain between red and purple boundary) to be meshed finally ?

clc; clear; close all;
% Generating simple SDF
[x, y] = meshgrid(linspace(-2, 2, 100), linspace(-2, 2, 100));
sdf = sqrt(x.^2 + y.^2) - 1;  
% Extracting 0-level isocontour of the SDF.
isoValue = 0;  

c1 = contour(x(1,:), y(:,1), sdf, [isoValue isoValue],"r"); 
hold on;
c2 = contour(x(1,:), y(:,1), sdf, [isoValue-0.5 isoValue-0.5]);

[Star Strider]

My pleasure!

That is an extension of my original code, hoever it is not possible to detect and specifically plot the inner boundary. The isInterior function should work here, accorrding to the documentation on it, however it doesn’t. II can’t get it to work on either the original triiangulation or the one read from the file.

I consider this a bug. Anyway, extracting the inner contour from the imported file is not possible.

Thee featureEdges function works (in the sense that it at least returns a result, as opposed to isInterior that completely fails), however it is not appropriate for these data.

I get the impression that this problem may not be appropriate for a STL file, since it usually extracts and stores surface data, and this is 2D rather than 3D.

This is the best I can do with this —

% clc; clear; close all;
% Generating simple SDF
[x, y] = meshgrid(linspace(-2, 2, 100), linspace(-2, 2, 100));
sdf = sqrt(x.^2 + y.^2) - 1;  
% Extracting 0-level isocontour of the SDF.
isoValue = 0;  

[c{1},h{1}] = contour(x(1,:), y(:,1), sdf, [isoValue isoValue],"r"); 
hold on;
[c{2},h{2}] = contour(x(1,:), y(:,1), sdf, [isoValue-0.5 isoValue-0.5]);
hold off
axis('equal')

% c{1}
% c{2}

Levels{1} = h{1}.LevelList;
Levels{2} = h{2}.LevelList;

% return

for k1 = 1:numel(Levels)
    idxv = find(c{k1} == Levels{k1});
    for k2 = 1:numel(idxv)
        lenv = c{k1}(2,idxv(k2));
        xc{k1,k2} = c{k1}(1,2:lenv);
        yc{k1,k2} = c{k1}(2,2:lenv);
    end
end

xc

xc = 2x1 cell array
    {[      0.1417 0.1414 0.1010 0.0606 0.0202 -0.0202 -0.0606 -0.1010 -0.1414 -0.1417 -0.1818 -0.2222 -0.2626 -0.3030 -0.3133 -0.3434 -0.3838 -0.4163 -0.4242 ... ] (1x200 double)}
    {[0.1844 0.1818 0.1414 0.1010 0.0606 0.0202 -0.0202 -0.0606 -0.1010 -0.1414 -0.1818 -0.1844 -0.2222 -0.2626 -0.2645 -0.3030 -0.3201 -0.3434 -0.3631 -0.3838 ... ] (1x96 double)}

yc

yc = 2x1 cell array
    {[-0.9899 -0.9900 -0.9949 -0.9982 -0.9998 -0.9998 -0.9982 -0.9949 -0.9900 -0.9899 -0.9833 -0.9750 -0.9649 -0.9530 -0.9495 -0.9392 -0.9234 -0.9091 -0.9055 ... ] (1x200 double)}
    {[ -0.4646 -0.4658 -0.4796 -0.4897 -0.4963 -0.4996 -0.4996 -0.4963 -0.4897 -0.4796 -0.4658 -0.4646 -0.4478 -0.4255 -0.4242 -0.3975 -0.3838 -0.3631 -0.3434 ... ] (1x96 double)}

xd = cat(2,xc{:}).';                                                                % First Row Of ‘c’  Are The ‘X’ Coordinates, Delete The First Element From Each
yd = cat(2,yc{:}).';                                                                % Second Row Of ‘c’ Are The ‘Y’ Coordinates, Delete The First Element From Each
DT = delaunayTriangulation(xd,yd);                                                               % Delaunay Triangulation

figure
triplot(DT, xd, yd)                                                                 % Plot Result
axis('equal','padded')

% return

% TR = triangulation(DT,xd,yd);                                                       % Create ‘triangulation’ Object

stlwrite(DT, 'STL_Test.stl', 'text')                                                % Write Triangulation Object To STL File

TRr = stlread('STL_Test.stl');                                                      % Check: Read Triangulation Object From STL File
TRr

TRr = 
  triangulation with properties:

              Points: [296x3 double]
    ConnectivityList: [390x3 double]

boundaryEdges = freeBoundary(TRr).';                                                % Get Boundary Edges
% Lvii = isInterior(TRr.Points)
% Lvii = isInterior(TRr.ConnectivityList)
% Lvii = isInterior(DT);
% FE = featureEdges(TRr, pi/10);

figure
triplot(TRr)                                                                        % Plot ‘triangulation’ Objeect
hold on
plot(TRr.Points(boundaryEdges,1), TRr.Points(boundaryEdges,2), '-r', 'LineWidth',2) % Plot Boundary
% plot(TRr.ConnectivityList(Lvii,1), TRr.Points(LVii,2), '-r', 'LineWidth',2) % Plot Boundary
% triplot(TRr.ConnectivityList(Lvii,:),TRr.Points(:,1),TRr.Points(:,2)) 
hold off
axis('equal','padded')

I experimented with every function I can find to extract the inner contour. None of them work on the imported data.

.

[ArxIv]

Thank you for your help.

I really appreciate it!

[Star Strider]

As always, my pleasure!

[DGM]

I'm going to use a more complicated example scenario, but the original example with the conical surface is included and does work.

%% choose your data and levels

% Generating simple SDF
%[x, y] = meshgrid(linspace(-2,2,100)); % original annulus example
%sdf = sqrt(x.^2 + y.^2) - 1;  
[x y sdf] = peaks(100); % multiple regions, dissimilar lengths

% pick two specific levels
%contourlevels = [-0.5 0]; % for first example
contourlevels = [1.5 2.5]; % for second example

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get contours for two specific levels
figure(1)
[C,hc] = contour(x(1,:),y(:,1),sdf,contourlevels); 
axis equal; grid on

% get vertex sequences from contours
% Vc is a cell array of [x y] vertex lists describing individual curves,
% each column in Vc contains the curves for a particular level
for kl = 1:numel(contourlevels) % index over each level
    startidx = find(C(1,:) == contourlevels(kl)); % start of each curve
    for kc = 1:numel(startidx) % index over each curve in this level
        blockstart = startidx(kc); % the starting index for this block
        lenv = C(2,blockstart); % the length of this curve + 1
        Vc{kc,kl} = C(:,blockstart+(2:lenv)).'; % C is [level, x; length, y]
    end
end

% reshape the cell array to a single column
% ideally, we should first flip the vertex lists in each column
% such that the interior curves are CW and exterior curves are ccw
% but we don't necessarily know which is which.
% delaunayTriangulation() seems to figure it out without ordering
% but it might not always work as expected.
Vc = Vc(:);

% construct the edge lists (constraints) in the same direction
% the vertex lists in Vc don't need to be the same length
% i'm only presuming these curves are supposed to be closed
% if they're not closed, they will be.
sz = cellfun(@(x) size(x,1),Vc); % length of each vertex list
szc = [0; cumsum(sz)];
E = zeros(0,2);
for k = 1:numel(Vc)
    v = 1+szc(k):szc(k+1);
    thisE = [v; circshift(v,-1)].';
    E = [E; thisE];  %#ok<*AGROW>
end

% TRIANGULATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do a constrained triangulation to generate the triangles
T = delaunayTriangulation(cell2mat(Vc),E); % use the constraints
F = T.ConnectivityList(isInterior(T),:); % eliminate extraneous faces
V = T.Points;

% write it
T = triangulation(F,V);
stlwrite(T,'test.stl')

% display it using patch()
figure
patch('faces',F,'vertices',V,'facecolor',[1 1 1]*0.8,'edgecolor','k');
axis equal; grid on

My interpretation of the task might be taking some liberties. I'm not sure whether there was a good reason to be using multiple contour objects, so I just used one. I'm assuming that we're only dealing with two contour levels at a time. I also don't know whether we're supposed to just assume that the contours are always closed, or how open contours should be handled. I simply didn't handle them in any special way. I also don't see how we can assume which curve is interior without prior knowledge of convexity. I just let delaunayTriangulation() try to figure it out, but that might not always work.

I bet there's an even more concise way to do this specifically for contourf plots, but that's a blind guess based on an extremely confused understanding of how some of the undocumented descendant graphics objects work. Might be totally wrong about how they work, and in the best case, it would be a fragile, undocumented solution.

[ArxIv]

Very insightful. Thank you very much!