How can I smooth a surf with many spikes?

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.

[DGM]

How is this different from the prior question?

https://www.mathworks.com/matlabcentral/answers/1839283-how-to-smooth-surf?s_tid=prof_contriblnk

[Alexandra Roxana]

I left the entire code here. And there was more a general question.

[KSSV]

Did you try increasing M and N?

[Alexandra Roxana]

For M=60 and N=80, there are even more spikes.

[Star Strider]

It might be easier to help with this if we knew what you want to do. There might be better and more efficient ways to code it.

[Alexandra Roxana]

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]

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.

Answer 1

[Image Analyst]

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]

Thank you very much! The removal color program compiled. What values could I change so it can look the way it should? 0 and 255?

The grayscale one showed only two plots. I get the error:

Error using medfilt2

Expected input number 1, A, to be two-dimensional.

Error in medfilt2>parse_inputs (line 106)

validateattributes(a, ...

Error in medfilt2 (line 48)

[a, mn, padopt] = parse_inputs(varargin{:});

Error in salt_and_pepper_noise_removal_grayscale1 (line 43)

medianFilteredImage = medfilt2(noisyImage, [3 3]);

[Image Analyst]

@Alexandra Roxana I just tried them both and they work beautifully. You must have altered them somehow to adapt them to your images instead of the MATLAB demo images. But you forgot to attach your altered code and images. If you do that I might be able to fix it.

[Alexandra Roxana]

I only changed the directory. But I can attach the files.

[Image Analyst]

You need to supply your actual image array. The 2-D array of numbers. You can't give it a 3-D-ish rendering of a surface plot of the image. It needs the image itself, not a screenshot of a rendering of it. So I guess that's sol, rather than the saved screenshot of it plotted with surf().

[Alexandra Roxana]

I did something like this. What values could I change to make it look more appropriate? I'm not quite proficient with MATLAB.

    figure(2)
    set(gcf, 'Position', get(0,'Screensize')); 
    noisyRGB = imnoise(sol,'salt & pepper', 0.05);
    imshow(noisyRGB);

[Image Analyst]

Here is sol displayed as an image rather than a surface rendering:

Exactly what pixels do you want filtered or modified in intensity, and why?

[Alexandra Roxana]

@Image Analyst I would like the spikes mostly from the red part to be not so obvious, that means more close values like the ones from the close neighboring sections.

[Image Analyst]

Did you try a median filter, or even better, the salt and pepper demo? You don't need to create a demo image with noise like I did. You already have your noisy data in sol. Don't create more noise in it. Just use what I did with the noisyImage

noisyImage = sol;

[Alexandra Roxana]

I did this:

 figure(2)
 noisyImage=sol;
 medianFilteredImage = medfilt2(noisyImage, [3 3]);
 imshow(medianFilteredImage);

[Alexandra Roxana]

and it looks like this

Is it how it's supposed to look? I want it to be coloured, the blue parts to be intact.

[Image Analyst]

I don't know how it's supposed to look. I'm asking you.

From what I can see in your surface plot of the data it looks like the blue parts are the left and right two columns of the matrix. Do you want to replace the filtered matrix with the original matrix in those columns?

filteredMatrix(:, 1:2) = originalMatric(:, 1:2);
filteredMatrix(:, end-1:end) = originalMatric(:, end-1:end);

Your data is not inherently colored. It's not a true color RGB image. You can apply a colormap to pseudocolor the matrix, but you need to specify the colormap. Like what gray level gets mapped to what RGB color. I have no idea what that would be. What is this array anyway? What does it represent? Why do you not know what it "should" look like? If you don't know what it should look like then why do you think it's wrong as it is? Maybe it's fine. What are you going to do with it anyway? What would be the next step if you could filter it somehow?

[Alexandra Roxana]

@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]

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?

Answer 2

[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

https://uk.mathworks.com/help/images/ref/imfilter.html

and I would suggest to use same size and replicate padding to avoid boundary problems.

[Alexandra Roxana]

I did this:

img=surf(X,Y,sol',C);
 h=[1 1 1;1 1 1;1 1 1]/9;
imfilter(img,h)

But it doesn't work. So, in my case, it works only if I upload the image?

[Constantino Carlos Reyes-Aldasoro]

no no, when you do img=surf what you are doing is saving the handle of the actual figure (with axes and all) into img. You need to filter your data, if the data is "sol" try

sol_filtered = imfilter(sol,h);

then

surf(X,Y,sol_filtered,C)

if you want to smooth more, then you would need to change the kernel, but try this first.

[DGM]

Note that by default, imfilter() zero-pads the array boundaries, so you'll get edge effects there.

[Alexandra Roxana]

Thank you! I did this

    h=[1 1 1;1 1 1;1 1 1]/9;
    sol_filtered = imfilter(sol,h);
    surf(X,Y,sol_filtered',C)
    colormap([1 0 0; 0 0 1]);

[Alexandra Roxana]

And it looks something like this

What should I do to have the boundaries quite the same and the plot to have almost the same structure?

[Image Analyst]

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]

@Image Analyst Thank you very much for your answers and explanations! Very much appreciated.