How can I smooth a surf with many spikes?

10 visualizaciones (últimos 30 días)
Alexandra Roxana
Alexandra Roxana el 1 de Nov. de 2022
Comentada: Alexandra Roxana el 5 de Nov. de 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.
  7 comentarios
Alexandra Roxana
Alexandra Roxana el 1 de Nov. de 2022
Editada: Alexandra Roxana el 1 de Nov. de 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
Alexandra Roxana el 1 de Nov. de 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.

Iniciar sesión para comentar.

Respuesta aceptada

Image Analyst
Image Analyst el 1 de Nov. de 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 comentarios
Alexandra Roxana
Alexandra Roxana el 4 de Nov. de 2022
Editada: Alexandra Roxana el 4 de Nov. de 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
Image Analyst el 4 de Nov. de 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?

Iniciar sesión para comentar.

Más respuestas (1)

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

Iniciar sesión para comentar.

Productos


Versión

R2017a

Community Treasure Hunt

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

Start Hunting!

Translated by