# How can I smooth a surf with many spikes?

8 views (last 30 days)

Show older comments

Alexandra Roxana
on 1 Nov 2022

Commented: Alexandra Roxana
on 5 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.

##### 7 Comments

### Accepted Answer

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
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?

### More Answers (1)

Constantino Carlos Reyes-Aldasoro
on 1 Nov 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 Comments

Image Analyst
on 5 Nov 2022

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!