Matlab code for generating some shapes using signed distance

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

[jonas]

What is the purpose of this? Do you just want to draw the shapes? If so, better use polyshape.

[SAMUEL AYINDE]

No, I am not just drawing shapes. I want the signed-distance field. Take the ellipse or any of the shape for example, if you plot the figure with the code below, you can have a better understanding of what I am talking about.

figure(10)
    contourf( reshape( FENd.x, M), reshape(FENd.y , M), reshape(FENd.Phi, M),);
    axis equal;  grid on;drawnow;colorbar;

But if you plot it with the code below, you have removed the negative part, that is the reason you will have the shape on white background. What I need is not just a shape. Rather, I need a signed distance field that can describe or generate those shapes. Thank you so much.

 figure(10)
      contourf( reshape( FENd.x, M), reshape(FENd.y , M), reshape(FENd.Phi, M), [0 0] );
      axis equal;  grid on;drawnow;colorbar

[jonas]

Well, it's not clear from your figure. You should rather show this figure:

[SAMUEL AYINDE]

Alright. Thank you so much. You can generate that figure with the code below. I will appreciate it if you can help.

figure(11)
    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;

[jonas]

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]

These are other pictures

Answer 1

[jonas]

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.

[jonas]

@Walter: I see that you edited the post at some point while I was probably still editing myself. There is a risk your edit therefore got reverted.

[SAMUEL AYINDE]

@jonas, thank you so much. This is quite insightful. You gave an insight into case 3. How can I go about the case 1 and case 2 in the attached Picture_1.jpg? Thank you so much.

[jonas]

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]

@jonas, thank you so much.

[jonas]

My pleasure! You got some code and I learned what a "signed distance field" is.

win-win

[SAMUEL AYINDE]

@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]

Hello again! My old code should work for 3D as well. You just have to define the "polygon" in 3d space. The function inpolygon does not work for 3d, but check this thread for an alternative.

https://se.mathworks.com/matlabcentral/answers/421672-how-to-built-a-3d-binary-volume-array-from-set-of-x-y-z-coordinates

See this very simplistic example for a single line in a 3d domain

% Define shape
xt=0.2:0.01:0.8;
yt=repmat(0.5,1,length(xt));
zt=repmat(0.5,1,length(xt));

% Define domain
[X,Y,Z] = meshgrid(0:0.01:1,0:0.01:1,0:0.01:1);

% Calculate distance
[~,V] = dsearchn([xt',yt',zt'],[X(:),Y(:),Z(:)]);

% Plot
isosurface(X,Y,Z,reshape(V,size(X)),0.2);

[SAMUEL AYINDE]

@jonas,Thank you so much. Can you please explain what that thread is talk about and how it relates to our case? It's like the person want to make a 3D binary image.

[jonas]

You can use the above code to get a unsigned (?) distance, ie distances are positive by default. This was solved before by changing the sign of all pts inside the polygon. The inpolygon function does not work in 3d, but you can use the method in that link to Calculate what pts lie inside a 3d shell, e.g. A sphere. I can show you tonight, in mobile atm.

Maybe I misunderstood your question. Is the shape in 2d or 3d?

[SAMUEL AYINDE]

The shape is in 3d. Thank you so much.

[jonas]

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.

[SAMUEL AYINDE]

jonas, thank you so much.

[jonas]

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

[SAMUEL AYINDE]

Thank you so much, @jonas!

Answer 2

[SAMUEL AYINDE]

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