How can I smooth a surf with many spikes?

15 views (last 30 days)
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.
  7 Comments
Alexandra Roxana
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.

Sign in to comment.

Accepted Answer

Image Analyst
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.
  13 Comments
Image Analyst
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?

Sign in to comment.

More Answers (1)

Constantino Carlos Reyes-Aldasoro
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.
  7 Comments
Alexandra Roxana
Alexandra Roxana on 5 Nov 2022
@Image Analyst Thank you very much for your answers and explanations! Very much appreciated.

Sign in to comment.

Products


Release

R2017a

Community Treasure Hunt

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

Start Hunting!