Matlab code for generating some shapes using signed distance

2018 10 13
2 个回答
回答已采纳
更新时间：2019 4 23
6 次查看（30 天）

0 个投票

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
jonas 2018-10-13
What is the purpose of this? Do you just want to draw the shapes? If so, better use polyshape.
SAMUEL AYINDE
SAMUEL AYINDE 2018-10-13
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
jonas 2018-10-13
Well, it's not clear from your figure. You should rather show this figure:
SAMUEL AYINDE
SAMUEL AYINDE 2018-10-13
编辑：SAMUEL AYINDE 2018-10-13
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
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 :)

请先登录，再进行评论。

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

 采纳的回答

jonas
jonas 2018-10-14
编辑：jonas 2018-10-14

0 个投票

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
jonas 2018-10-14
@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
SAMUEL AYINDE 2018-10-14
编辑：SAMUEL AYINDE 2018-10-14
@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
jonas 2018-10-14
编辑：jonas 2018-10-14
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.
jonas
jonas 2018-10-14
My pleasure! You got some code and I learned what a "signed distance field" is.
win-win
SAMUEL AYINDE
SAMUEL AYINDE 2018-10-15
编辑：SAMUEL AYINDE 2018-10-15
@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
jonas 2018-10-15
编辑：jonas 2018-10-15
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
SAMUEL AYINDE 2018-10-15
@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
jonas 2018-10-15
编辑：jonas 2018-10-15
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
SAMUEL AYINDE 2018-10-15
The shape is in 3d. Thank you so much.
jonas
jonas 2018-10-15
编辑：jonas 2018-10-15
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
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 :)

请先登录，再进行评论。

更多回答（1 个）

SAMUEL AYINDE
SAMUEL AYINDE 2019-4-23

0 个投票

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

类别

在 帮助中心 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by