# How can I smooth a surf with many spikes?

15 views (last 30 days)
Alexandra Roxana on 1 Nov 2022
I'm having problems to smooth a surf for a FEM program. Any help would be much appreciated. I'm leaving the code below:
function [] = MEF3()
clc
close all
%date
a=0;
b=2;
c=0;
d=4;
M=40; %no of points on X
N=60; %no of points on Y
hx=(b-a)/(M+1);
hy=(d-c)/(N+1);
Nod=zeros((M+2)*(N+2),3); %node matrix
Tr=zeros(2*(M+1)*(N+1),3); %triangles matrix
sigma=2;
rhos=1024;
rho0=1000;
Lc=0.01;
vc=0.1;
mu=0.004;
Re=rho0*vc*Lc/mu;
T0=253;
alphaf=207; alphafast=T0*alphaf;
g=9.81*Lc/vc^2;
E=882000;
k=4900000; kast=k/(rhos^2);
alphas=6.1; alphasast=T0*alphas;
tau=1;
dt=tau/10;
cf=3617;
cs=3306;
Ec=vc^2/(cs*T0);
epsilon=1/10;
fs=230; fsast=230*(Lc/rhos);
Qf=0;
Qs=0;
kf=0.52;
ks=0.46;
Prf=(mu*cf)/kf;
Prs=(mu*cs)/ks*(rhos/rho0);
chi3=1/(Re*Prf)+(cs/cf)*(rhos/rho0)*1/(Re*Prs);
iter=3;
kt=0; %total nodes
ki=0; %internal nodes
x=zeros(1,M+2);
y=zeros(1,N+2);
%creating nodes matrix
for j=1:N+2
y(j)=c+hy*(j-1);
for i=1:M+2
x(i)=a+hx*(i-1);
kt=kt+1;
Nod(kt,1)=x(i);
Nod(kt,2)=y(j);
if(i==1)||(i==M+2)||(j==1)||(j==N+2)
Nod(kt,3)=0;
else
ki=ki+1;
Nod(kt,3)=ki;
end
end
end
%display(Nod)
%creating triangles matrix
kTr=0;
for j=1:N+1
for i=1:M+1
nod=(j-1)*(M+2)+i;
if (((i==M+1)&&(j==1)) || ((i==1)&&(j==N+1)))
kTr=kTr+1;
Tr(kTr,1)=nod;
Tr(kTr,2)=nod+M+2;
Tr(kTr,3)=nod+1;
kTr=kTr+1;
Tr(kTr,1)=nod+1;
Tr(kTr,2)=nod+M+2;
Tr(kTr,3)=nod+M+3;
else
kTr=kTr+1;
Tr(kTr,1)=nod;
Tr(kTr,2)=nod+M+2;
Tr(kTr,3)=nod+M+3;
kTr=kTr+1;
Tr(kTr,1)=nod;
Tr(kTr,2)=nod+M+3;
Tr(kTr,3)=nod+1;
end
end
end
%stiffness matrix and free term
solold=zeros(2*ki,1);
A=zeros(2*ki);
L=zeros(2*ki,1);
for t=1:(2*(M+1)*(N+1))
n1=Tr(t,1);
n2=Tr(t,2);
n3=Tr(t,3);
x1=Nod(n1,1);
y1=Nod(n1,2);
x2=Nod(n2,1);
y2=Nod(n2,2);
x3=Nod(n3,1);
y3=Nod(n3,2);
M1=[x1 y1 1; x2 y2 1; x3 y3 1];
de=det(M1);
ariat=abs(de)/2;
C=inv(M1); %barycentric matrix
a1=C(1,1);
b1=C(2,1);
c1=C(3,1);
a2=C(1,2);
b2=C(2,2);
c2=C(2,3);
a3=C(1,3);
b3=C(2,3);
c3=C(3,3);
%elementary stiffness matrix and elementary free term
Ae=zeros(6);
Le=zeros(6,1);
coef_a=[a1 a2 a3]; coef_b=[b1 b2 b3]; coef_c=[c1 c2 c3];
for k=1:3
for i=1:3
Ae(k,i)=...
Ae(k,i+3)=...
Ae(k+3,i)= ...
Ae(k+3,i+3)=...;
end
end
%Le is written as an approximation
Le(1)=...;
Le(2)=...;
Le(3)=...;
Le(4)=...;
Le(5)=...;
Le(6)=...;
n=[n1 n2 n3];
for k=1:3
if(Nod(n(k),3)~=0) %internal node
for i=1:3
if (Nod(n(i),3)~=0) %internal node
A(Nod(n(k),3),Nod(n(i),3))=...
A(Nod(n(k),3),Nod(n(i),3))+Ae(k,i);
A(Nod(n(k),3),Nod(n(i),3)+ki)=...
A(Nod(n(k),3),Nod(n(i),3)+ki)+ Ae(k,i+3);
A(Nod(n(k),3)+ki,Nod(n(i),3))=...
A(Nod(n(k),3)+ki,Nod(n(i),3))+Ae(k+3,i);
A(Nod(n(k),3)+ki,Nod(n(i),3)+ki)=...
A(Nod(n(k),3)+ki,Nod(n(i),3)+ki)+Ae(k+3,i+3);
end
end
L(Nod(n(k),3))=L(Nod(n(k),3))+Le(k);
L(Nod(n(k),3)+ki)=L(Nod(n(k),3)+ki)+Le(k+3);
end
end
end
alpha=inv(A)*L; %u approximated in the internal nodes
sol=zeros((M+2),(N+2));
for k=1:N*M
i=rem(k,M);
if i==0
i=M;
end
j=(k-i)/M+1;
solnew(i,j)=alpha(k);
end
figure(1)
sol(2:M+1,2:N+1)=solnew;
[X,Y]=meshgrid(a:hx:b,c:hy:d);
C = 1 + (X <= 0.25 | X >= 1.75);
surf(X,Y,sol',C);
colormap([1 0 0; 0 0 1]);
hold on
plot(0,0,'ko','MarkerSize',10,'MarkerFaceColor','black')
text(0,0,'(0,0)','Fontsize',20,'Color','k', 'Horiz','left','Vert','top')
plot(2,0,'k>','MarkerSize',10,'MarkerFaceColor','black')
text(2,0,'x_1','Fontsize',20,'Color','k', 'Horiz','left','Vert','top')
plot(0,4,'k<','MarkerSize',10,'MarkerFaceColor','black')
text(0,4,'x_2','Fontsize',20,'Color','k', 'Horiz','left','Vert','top')
xlabel('cm')
ylabel('cm')
zlabel('longitudinal fluid velocity (cm/s)')
end
Here is the result.
Alexandra Roxana on 1 Nov 2022
The solution is something like () and I have to plot only the first component of , but the surface plot looks bad so I chose to plot the entire solution.

Image Analyst on 1 Nov 2022
Try my attached modified median filter. Adapt as needed. Basically you identify spikes somehow, like thresholding (globally or locally) or whatever method you like. For example you could smooth the image with imfilter or conv2 and then look where the difference between smoothed and original is more than some value, then threshold it.
Then you replace the data in the spike locations with the values of the median filtered image from those locations.
Image Analyst on 4 Nov 2022
I've asked this before I believe but how big is too big? What threshold do you want to use?
Also, a point has 8 neighbors. Do you want to compare a point's value with the point value MOST different or compare it with the average of the 8 neighbors?

If you just want to smooth the values, you could use a filter, from a small uniform filter like [1 1 1;1 1 1;1 1 1]/9 to a Gaussian would help. The best function for this is imfilter
and I would suggest to use same size and replicate padding to avoid boundary problems.
Alexandra Roxana on 5 Nov 2022
@Image Analyst Thank you very much for your answers and explanations! Very much appreciated.

### Categories

Find more on Red in Help Center and File Exchange

R2017a

### Community Treasure Hunt

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

Start Hunting!