Matlab code for generating some shapes using signed distance

Question

SAMUEL AYINDE 2018-10-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/423827-matlab-code-for-generating-some-shapes-using-signed-distance

回答： SAMUEL AYINDE 2019-4-23

Picture_1.jpg

In order to get a good understanding of my question, please, you may first run the matlab code below. This code can generate the shape for rectangle, ellipse and circle using signed distance if you uncomment the portion corresponding to each shape.

Please, I need a code that can give the shapes in the attached picture (Picture_1.jpg) using signed distance. Thank you so much.

   %%%Signed distance code for generating shapes
      clear all;
      close all;
      clc;
      DomainWidth=2;
      DomainHight=1;
      ENPC=40;   
      ENPR=80;
      EW = DomainWidth / ENPR;       %  The width of each finite element.
      EH = DomainHight / ENPC;        %   The hight of each finite element.
      M = [ ENPC + 1 , ENPR + 1 ];
      [ x,y ] = meshgrid( EW * [ -0.5 : ENPR + 0.5 ] , EH * [ -0.5 : ENPC + 0.5 ]);
      [ FENd.x, FENd.y, FirstNdPCol] = MakeNodes(ENPR,ENPC,EW,EH);
      LSgrid.x = x(:); LSgrid.y = y(:);       % The coordinates of Level Set grid 1
      cx = DomainWidth/2;
      cy= DomainHight/2;
      a = cx;
      b = 0.8*cy;
      %%Generate circle
        tmpPhi=  sqrt ( ( LSgrid . x - cx  ) .^2 + ( LSgrid . y  - cy  ) .^2 ) - DomainHight/2; 
       LSgrid.Phi = -((tmpPhi.')).';
       %%Generate ellipse
      %   tmpPhi=   ( (LSgrid . x - cx)/a  ) .^2 + (( LSgrid . y  - cy)/(0.8*cy)  ) .^2  - 1;
      %  LSgrid.Phi = -((tmpPhi.')).';
      %%Generate rectangle
      % lower = [cx - 0.5 * DomainWidth, cy - 0.25 * DomainHight];
      % upper = [cx + 0.5 * DomainWidth, cy + 0.25 * DomainHight];
      % Phi11 =max(LSgrid.x  - upper(1), lower(1) - LSgrid . x );
      % Phi11 =max(Phi11,LSgrid.y   - upper(2));
      % Phi11 =max(Phi11,lower(2) - LSgrid . y );
      % LSgrid.Phi = -Phi11;
        FENd.Phi = griddata( LSgrid.x, LSgrid.y, LSgrid.Phi, FENd.x, FENd.y, 'cubic');
    figure(10)
    % Figure of the full contour plot of the signed distance of the shape
         contourf( reshape( FENd.x, M), reshape(FENd.y , M), reshape(FENd.Phi, M));
          axis equal;  grid on;drawnow; colorbar;
     figure(11)
    % Figure of the scaled contour plot showing only the shape
          contourf( reshape( FENd.x, M), reshape(FENd.y , M), reshape(FENd.Phi, M), [0 0] );
          axis equal;  grid on;drawnow; colorbar;
   figure(12)
    %Figure of the signed distance function 
        h3=surface(x, y, reshape(-LSgrid.Phi , M + 1));  view([37.5  30]);  axis equal;  grid on; 
        set(h3,'FaceLighting','phong','FaceColor','interp', 'AmbientStrength',0.6); light('Position',[0 0 1],'Style','infinite'); colorbar;  
       %%Code for creating nodes    
      function [NodesX, NodesY, FirstNdPCol] = MakeNodes(EleNumPerRow,EleNumPerCol,EleWidth,EleHight)
      [ x , y ]= meshgrid( EleWidth * [ 0 : EleNumPerRow ], EleHight * [0 : EleNumPerCol]);
      FirstNdPCol = find( y(:) == max(y(:)));
      NodesX = x(:); NodesY = y(:);
      end

6 个评论
显示 4更早的评论隐藏 4更早的评论

jonas 2018-10-13

Not sure if I can, but I will give it a try. Anyway, if the question is clear then usually someone will answer eventually :)

SAMUEL AYINDE 2018-10-13

These are other pictures

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

jonas 2018-10-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/423827-matlab-code-for-generating-some-shapes-using-signed-distance#answer_341349

编辑：jonas 2018-10-14

在 MATLAB Online 中打开

I had to read a bit about "signed distance fields" and came across this blogpost

https://blogs.mathworks.com/graphics/2016/03/07/signed-distance-fields/

Have your read it? Should not be too hard to implement for your shapes.

I would attempt the following more 'general' algorithm.

describe the shape with a set of x-y-coordinates
create a mesh
find the closest distance to each point in the mesh to the set of x-y-coordinates

If I understand correctly, that is what the "signed distance field" describe, i.e. the closest distance to a shape from any point in the domain.

This is the code for a single horizontal line from [0,0.5] to [1,0.5].

% The "shape"
x=0:0.01:1;
y=ones(size(x)).*.5;
% The "domain"
[X,Y] = meshgrid(0:0.01:1,0:0.01:1);
% Closest point between domain and shape
[X,Y] = meshgrid(0:0.01:1,0:0.01:1);
[id,Z] = dsearchn([x',y'],[X(:),Y(:)]);
scatter3(X(:),Y(:),Z,[],Z)

The only thing you have to change is the sign of those points that are inside of the "shape". This is however fairly simple, using the inpolygon function. Let's try one for the sinusoidal shape:

% Build polygons
x=0:0.01:1;
Y{1} = 0.2.*sin(x/0.05)+0.3
Y{2} = 0.2.*sin(x/0.05)+0.7
py = [Y{1},fliplr(Y{2})] 
px = [x,fliplr(x)];
% Define domain
[X,Y] = meshgrid(0:0.01:1,0:0.01:1);
% Calculate distance
[~,Z] = dsearchn([px',py'],[X(:),Y(:)]);
[in]=inpolygon(X(:),Y(:),px,py)
Z(in) = -Z(in)
Z = -Z;
% Plot
surf(X,Y,reshape(Z,size(X)),'edgecolor','none');hold on
p = polyshape(px,py)
plot(p)
axis([0 1 0 1])

.

There seems to be some edge-effect, in particular to the left. Other than that, it looks quite OK.

14 个评论
显示 12更早的评论隐藏 12更早的评论

jonas 2018-10-14

编辑：jonas 2018-10-14

在 MATLAB Online 中打开

Happy that it works for you! The nice thing about this approach is that its general. All you have to do is to be creative in the way you build your polygon at the start of the code.

% Build linear piecewise polygon
% Corners of top top saw shape
xt=0:1/10:1;
yt=[repmat([0 0.1],1,5) 0]
%Interpolate between corners
xi = 0:0.005:1;
yi = interp1(xt,yt,xi)
% build polygon
px = [xi,fliplr(xi)]
py = [yi+0.5,fliplr(-yi+0.5)];
p = polyshape(px,py)
plot(p)
%%Same code as before............
% Define domain
[X,Y] = meshgrid(0:0.01:1,0:0.01:1);
% Calculate distance
[~,Z] = dsearchn([px',py'],[X(:),Y(:)]);
[in]=inpolygon(X(:),Y(:),px,py)
Z(in) = -Z(in)
Z = -Z;
% Plot
surf(X,Y,reshape(Z,size(X)),'edgecolor','none');hold on
p = polyshape(px,py)
plot(p)
axis([0 1 0 1])

All you have to do is learn to build polygons, which is very easy. The only tricky part is that you have to build a closed shape, which is why I have used fliplr() at some places (flip the vector horizontally).

It is also important to interpolate, because dsearchn looks for the closest point. Higher resolution on the mesh as well as higher resolution on the boundary of the polygon will give better results.

SAMUEL AYINDE 2018-10-15

编辑：SAMUEL AYINDE 2018-10-15

在 MATLAB Online 中打开

@jonas, I appreciate you. Please, how can I make those two cases into three-dimensional. I need a three dimensional shape in signed distance function of the shapes in the picture. For instance, the matlab code below will compute the three-dimensional signed distance for sphere. Thank you so much.

close all;
clear all;
clc;
DomainWidth=2;
DomainHight=1;
DomainLenght=1;
ENPC=40;
ENPR=80;
ENPL=40;
nelx = ENPR;
nely = ENPC;
nelz=ENPL;
EW = DomainWidth / ENPR;       %  The width of each finite element.
EH = DomainHight / ENPC;        %   The hight of each finite element.
EL = DomainLenght/ENPL;           %   The length of each finite element.
M = [ ENPC + 1 , ENPR + 1 , ENPL+1];
[ x ,y,z ] = meshgrid( EW * [ -0.5 : ENPR + 0.5 ] , EH * [ -0.5 : ENPC + 0.5 ], EL * [ -0.5 : ENPL + 0.5 ]);
[ FENd.x, FENd.y,FENd.z, FirstNdPCol ] = MakeNodes(ENPR,ENPC,ENPL,EW,EH, EL);
LSgrid.x = x(:); LSgrid.y = y(:); LSgrid.z = z(:);
cx = DomainWidth / 200 * [ 33.33  100  166.67  0   66.67  133.33  200  33.33  100  166.67  0   66.67  133.33  200  33.33  100  166.67];
cy = DomainHight / 100 * [   0     0     0     25   25      25     25    50    50    50    75    75     75     75    100   100   100]; 
cz = DomainLenght/ 100 * [   0     0     0     25   25      25     25    50    50    50    75    75     75     75    100   100   100]; 
for i = 1 : length( cx )
       tmpPhi( :, i ) = - sqrt ( ( LSgrid . x - cx ( i ) ) .^2 + ( LSgrid . y  - cy ( i ) ) .^2 + ( LSgrid . z  - cz ( i ) ) .^2 ) + DomainHight/8;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
end;
LSgrid.Phi =  -(max(tmpPhi.')).';                          % The level set function value on each Level Set grid
FENd.Phi = griddata( LSgrid.x, LSgrid.y,LSgrid.z, LSgrid.Phi, FENd.x, FENd.y,FENd.z,'natural');
 figure(110)
    isosurface( reshape( FENd.x, M), reshape(FENd.y , M),  reshape(FENd.z , M),reshape(FENd.Phi , M), 0);
%%Code to make nodes  
function [NodesX, NodesY, NodesZ,FirstNdPCol] = MakeNodes(EleNumPerRow,EleNumPerCol,EleNumPerLength,EleWidth,EleHight,EleLength)
[ x , y,z ]= meshgrid( EleWidth * [ 0 : EleNumPerRow ], EleHight * [0 : EleNumPerCol], EleLength * [0 : EleNumPerLength]);
FirstNdPCol = find( y(:) == max(y(:)));
NodesX = x(:); NodesY = y(:); NodesZ = z(:);
end

jonas 2018-10-15

编辑：jonas 2018-10-15

在 MATLAB Online 中打开

OKay. I'm just going to give you an idea of how this works. Then you will have to solve it for your shapes.

First, define a set of scattered points describing a sphere.

[Xt,Yt,Zt] = sphere(200);
Xt=Xt./4;
Yt=Yt./4;
Zt=Zt./4;

Define the domain, in 3D

[Xc,Yc,Zc] = meshgrid(-1:0.01:1,-1:0.01:1,-1:0.01:1);

Calculate the shortest distance between every point in the domain and a point on your sphere.

[~,Vc] = dsearchn([Xt(:),Yt(:),Zt(:)],[Xc(:),Yc(:),Zc(:)]);

Now, since this is a sphere, the radius of the isosurface at V = 0 should be equal to the radius of our sphere. All points inside of the sphere are defined by V < 0. Try drawing a negative isosurface.

isosurface(Xc,Yc,Zc,reshape(Vc,size(Xc)),-0.1);hold on

As you can see, nothing is drawn because the output of dsearchn is always positive. Now we have to find the points inside of the sphere and change their signs. To do this, we first perform a triangulation to define a shell.

tri = delaunayn([Xt(:) Yt(:) Zt(:)]);

Next, we find points inside of the shell.

tn = tsearchn([Xt(:) Yt(:) Zt(:)], tri, [Xc(:) Yc(:) Zc(:)]); 
IsInside = ~isnan(tn);

Finally, change sign on the points inside of the shell

Vc(IsInside) = -Vc(IsInside);

Try plotting again

isosurface(Xc,Yc,Zc,reshape(Vc,size(Xc)),-0.1);hold on

It takes some time, but the result is a nice little sphere.

For further questions you should post another thread. I'm always writing my answers so that others can learn from them in the future, would they come across this thread. No one reads this far down in the comments though.

jonas 2018-10-15

编辑：jonas 2018-10-15

My pleasure! I will update my previous answer with an image of the results if the run finishes :)

SAMUEL AYINDE 2018-10-16

Thank you so much, @jonas!

请先登录，再进行评论。

Answer 2

SAMUEL AYINDE 2019-4-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/423827-matlab-code-for-generating-some-shapes-using-signed-distance#answer_371775

New Doc 2019-04-23 00.41.24_1.jpg

Please, can someone help me with the code to generate the signed distance for this geometry? Thank you so much.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Matlab code for generating some shapes using signed distance

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

14 个评论
显示 12更早的评论隐藏 12更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Matlab code for generating some shapes using signed distance

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

14 个评论 显示 12更早的评论隐藏 12更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

14 个评论
显示 12更早的评论隐藏 12更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论