Scatteredinterpolant (linear) from symmetric data does not produce symmetric isolines

I utilised the scatteredinterpolant function to generate an interpolation from which I derived isolines. My data are situated on a non-equispaced grid of points that are symmetric with respect to the midline. In fact, I only have data on one side, so I placed data on the other side by copying the values from the original data set, exploiting the symmetry.

When I generate iso-lines using the NATURAL method, the isolines have a symmetric trend. If I use the simplest method, i.e. LINEAR, the iso-lines are not symmetric, which does not make sense.

Do you know what may cause this issue and how to sort it out?

Answer 1

[Bruno Luong]

scatteredinterpolant all methods is based on triangulation of the data points. Even the data points are symmetric, the triangulation migh not.

[ES_Thorny]

Thanks for you answer. Is it possible to force a symmetric triangulation in some way?

[John D'Errico]

@ES_Thorny

No. You cannot enforce a symmetric triangulation in scatteredInterpolant. That would require you to be able to supply the triangulation itself to scatteredInterpolant, and this is not an option.

Could you create a triangulation that has the property of symmetry that you desire, and using other interpolation tools perform the interpolation? The answer is, well, yes, in theory. It would require some serious effort on your part, and it would force you to learn a lot about triangulation methods. Is it worth the effort? That is, just to make a subtly different picture that you prefer, when as you found, other methods were better able to handle your problem?

[Bruno Luong]

The easiest workaround is to symmetrize the interpolation result if you query points are symmetric

[ES_Thorny]

Yes, I thought about it, but how can it be done? I need to understand what is the output of scatteredinterpolant ... and how I can ducplicate it on the other side ...

[Walter Roberson]

Map the query location...

QX = X;
mask = QX > MidPoint;
QX(mask) = 2 * MidPoint - QX(mask);

[Bruno Luong]

@ES_Thorny post your code (that can run) if you want to get any specific help. The output is scattered is the z coordinates of the query points (xq,yq) assumed you want to interpolated 2D function. It's described in the doc page. If you have specific question just ask.

[ES_Thorny]

    x_PP=[...]  % x coordinates
    y_PP=[...]  % x coordinate
    PP_i=[...] % P data
    
    x_PP_sym=x_PP; 
    y_PP_sym=-y_PP; 
    PP_i_sym=PP_i;
    
    % Adding other points symmetrically with respect to y=0 
    
    x_PP_sym=x_PP; 
    y_PP_sym=-y_PP; 
    PP_i_sym=PP_i;
    
    % remove duplicate points
    
    MM=[x_PP; y_PP; PP_i]';
    MM_u=unique(MM,'rows');
    x_PP=MM_u(:,1)';
    y_PP=MM_u(:,2)';
    P_v=MM_u(:,3)';
    
    % ---------------
    % Scattered Interpolant
    % ---------------
    
    % Grid 
    x0=min(x_PP);
    x1=max(x_PP);
    y0=min(y_PP);
    y1=max(y_PP);
    
    N=50; 
    xll=linspace(x0,x1,N) ;  
    yll=linspace(y0,y1,N) ; 
    
    [X,Y] = meshgrid(xll,yll) ;
    P=[x_PP,y_PP];
    V=PP_i;
    
    F = scatteredInterpolant(x_PP',y_PP',P_v'); 
    F.Method ='linear';
    F.ExtrapolationMethod ='none'; 
    Z=F(X,Y);
    
    % Contour plot
    
    clevels=[-95:10:-5 0 5:10:95];
    
    [ch0,ch]=contour(X,Y,Z,clevels);
    ch.LineWidth = 2;
    axis equal
    xlabel('x (mm)')
    ylabel('y (mm)');
    colormap(jet)
    Cbar=colorbar;
    Cbar.Label.String = 'Pressure (kPa)';
    caxis([clevels(1) clevels(end)])
    
    

[ES_Thorny]

Probably I have to build the conour levels on half grid and then duplicate the query point later.

[Bruno Luong]

    x_PP_sym=x_PP; 
    y_PP_sym=-y_PP; 
    PP_i_sym=PP_i;
    
    % Adding other points symmetrically with respect to y=0 
    
    x_PP_sym=x_PP; 
    y_PP_sym=-y_PP; 
    PP_i_sym=PP_i;

You create twice the variable and do not use it anywhere later. Clearly Your code is NOT correct.

[Bruno Luong]

This demo shows how to enforce the interpolation output to be symmetric, as it is an a priori expectation characteristic, whatever the input might be and or breaking symmetric due to scattering interpolation algorithm.

    [x_PP,y_PP,PP_i] = peaks(20);
    
    x_PP=x_PP(:);
    y_PP=y_PP(:);
    P_v=PP_i(:);
    
    % ---------------
    % Scattered Interpolant
    % ---------------
    
    % Grid 
    x0=min(x_PP);
    x1=max(x_PP);
    y0=min(y_PP);
    y1=-y0; %max(y_PP); % Change by Bruno
    
    N=50; 
    xll=linspace(x0,x1,N) ;  
    yll=linspace(y0,y1,N) ; 
    
    [X,Y] = meshgrid(xll,yll) ;
    P=[x_PP,y_PP];
    V=PP_i;
    
    F = scatteredInterpolant(x_PP,y_PP,P_v); 
    F.Method ='linear';
    F.ExtrapolationMethod ='none'; 
    Z=F(X,Y);
    Z=(Z+flip(Z,1))/2; % Add by Bruno, symmetrize the Z with to Y
    
    % Contour plot
    
    clevels=linspace(-4,4,10)

clevels = 1x10

-4.0000 -3.1111 -2.2222 -1.3333 -0.4444 0.4444 1.3333 2.2222 3.1111 4.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    
    [ch0,ch]=contour(X,Y,Z,10);
    ch.LineWidth = 2;
    axis equal
    xlabel('x (mm)')
    ylabel('y (mm)');
    colormap(jet)
    Cbar=colorbar;
    Cbar.Label.String = 'Pressure (kPa)';
    caxis([clevels(1) clevels(end)])
    yline(0)