# How can I smooth a surf with many spikes?

4 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
Edited: Alexandra Roxana on 1 Nov 2022
The variations in the graph are big, thus the spikes and I was wondering if it could be smoothed. I wonder if the data are chosen wrong. This is a TFSI problem and I'm solving it with FEM. The blue part represents the elastic structure and the red one the fluid structure.
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.
Alexandra Roxana on 4 Nov 2022
Edited: Alexandra Roxana on 4 Nov 2022
@Image Analyst I'm afraid this doesn't help either. I have done this and is starting to look somehow similar with what I would like to achieve.
figure(2)
h = [-2 0 2]/3;
sol_filtered = imfilter(sol,h,'conv');
surf(X,Y,sol_filtered',C)
colormap([1 0 0; 0 0 1]);
I don't know how to code it but I would like that when the difference between 2 points is too big and a spike appears, to compute the average, not to process the image. I don't know if that makes any sense but this is what I would like. There are some spikes left though in this image and the boundary doesn't quite look the same.
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.
Image Analyst on 5 Nov 2022
That filter simply blurs the image. Both "good", non-spiky data wil get smoothed, as well as isolated spikes, and edge pixels. If you want original "good" data to be retained/untouched, then you need to segment it to locate ONLY the spikes and then replace ONLY the spikes. I showed you how to do that in my Answer.
Alexandra Roxana on 5 Nov 2022
@Image Analyst Thank you very much for your answers and explanations! Very much appreciated.

### Categories

Find more on Image Filtering 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!