Unexpected result at interpolating values from scattered data

The goal is to create gridded data from scattered data. My scattered data (sample: XS1 and XS2) have [x,y,z] values and appear as multiple lines. To generate gridded data, I tried to interpolate the scattered data on a pre-generated grid using scatteredInterpolant(x,y,z). Unfortunately, the result is not as expected, as some of the sample points do not seem to be used in the interpolation.

Fig. 1: Plot of sample data XS1 and XS2

Fig 2: Unexpected result, most of XS1-data seems to be unused during interpolation.

Example of my interpolation problem:

%% Data
XS1 = [191.569000000000	9974.91400000000	22.8200000000000
191.556100000000	9974.93660000000	22.8200000000000
191.561900000000	9974.97310000000	22.8200000000000
191.628200000000	9975.32030000000	22.7400000000000
191.654200000000	9975.45940000000	22.6700000000000
191.667900000000	9975.54300000000	22.6000000000000
191.671700000000	9975.72320000000	22.5600000000000
191.650300000000	9975.84250000000	22.5100000000000
191.627500000000	9975.88500000000	22.5100000000000
191.556000000000	9975.94840000000	22.4900000000000
191.487700000000	9976.18690000000	22.4700000000000
191.338800000000	9976.31250000000	22.4600000000000
191.308400000000	9976.34640000000	22.4700000000000
191.170600000000	9976.50180000000	22.4500000000000
191.164400000000	9976.51000000000	22.4600000000000
191.154000000000	9976.52600000000	22.4500000000000
191.151000000000	9976.64480000000	22.4500000000000
191.047000000000	9976.78200000000	22.4400000000000
191.023800000000	9976.81970000000	22.4500000000000
190.958600000000	9976.87300000000	22.4600000000000
190.946300000000	9976.90080000000	22.4500000000000
190.941500000000	9976.96480000000	22.4500000000000
190.906000000000	9977.02150000000	22.4600000000000
190.896000000000	9977.08220000000	22.4400000000000
190.889300000000	9977.19400000000	22.4400000000000
190.874100000000	9977.29160000000	22.4600000000000
190.826200000000	9977.47420000000	22.4400000000000
190.657500000000	9977.96680000000	22.4600000000000
190.628300000000	9978.06250000000	22.4800000000000
190.571800000000	9978.27970000000	22.4800000000000
190.467700000000	9978.70870000000	22.5100000000000
190.446800000000	9978.82710000000	22.5400000000000
190.423700000000	9978.99860000000	22.5500000000000
190.426500000000	9979.22390000000	22.5800000000000
190.439800000000	9979.33960000000	22.6100000000000
190.501600000000	9979.60210000000	22.6400000000000
190.576800000000	9979.81980000000	22.6900000000000
190.619600000000	9979.91880000000	22.7400000000000
190.691400000000	9980.06460000000	22.7700000000000
190.770000000000	9980.19600000000	22.8400000000000
190.821100000000	9980.23690000000	22.9000000000000
190.949700000000	9980.25050000000	23.0100000000000
191.039300000000	9980.30000000000	23.0200000000000
191.153900000000	9980.30530000000	23.0300000000000
191.158500000000	9980.32710000000	23.1200000000000
191.159100000000	9980.34640000000	23.1400000000000
191.183000000000	9980.34770000000	23.2300000000000];

XS2 = [210.629000000000	9975.80200000000	22.8000000000000
210.602400000000	9975.77940000000	22.8100000000000
210.473700000000	9975.69270000000	22.7800000000000
210.454900000000	9975.67190000000	22.8200000000000
210.437500000000	9975.64880000000	22.7300000000000
210.378800000000	9975.52460000000	22.7100000000000
210.343700000000	9975.42220000000	22.7100000000000
210.325000000000	9975.36220000000	22.6800000000000
210.292600000000	9975.26140000000	22.6700000000000
210.275200000000	9975.20210000000	22.6300000000000
210.116400000000	9974.65730000000	22.4900000000000
209.950700000000	9974.07280000000	22.2900000000000
209.898400000000	9973.90170000000	22.2700000000000
209.870100000000	9973.80770000000	22.2400000000000
209.781500000000	9973.54460000000	22.2100000000000
209.600500000000	9972.98720000000	22.1700000000000
209.513200000000	9972.71100000000	22.2400000000000
209.363600000000	9972.22120000000	22.3100000000000
209.318800000000	9972.06310000000	22.3700000000000
209.266600000000	9971.78050000000	22.4400000000000
209.257900000000	9971.68210000000	22.5500000000000
209.256000000000	9971.57920000000	22.7100000000000
209.258800000000	9971.50120000000	22.7200000000000
209.265500000000	9971.40990000000	22.8100000000000
209.270700000000	9971.35170000000	22.7500000000000
209.276500000000	9971.29730000000	22.9000000000000
209.286200000000	9971.15370000000	23.0400000000000
209.179600000000	9971.12910000000	23.1700000000000
209.088800000000	9971.09510000000	23.1700000000000
209.049700000000	9971.09770000000	23.2400000000000];

%% Generate gridded data
n = 50;
[XSix,XSiy] = meshgrid(linspace(max(XS1(:,1)),min(XS2(:,1)),n),linspace(max(XS1(:,2)),min(XS2(:,2)),n));

%% Intepolate values
F = scatteredInterpolant([XS1(:,1);XS2(:,1)],[XS1(:,2);XS2(:,2)],[XS1(:,3);XS2(:,3)]);
F.ExtrapolationMethod = 'none';
F.Method = 'linear';
XSivq = F(XSix,XSiy);

%% Plot data

close all
figure()
hold on
plot3(XS1(:,1),XS1(:,2),XS1(:,3),'-','DisplayName','XS1','linewidth',2);
plot3(XS1(:,1),XS1(:,2),XS1(:,3),'.','HandleVisibility','off','markersize',5);

plot3(XS2(:,1),XS2(:,2),XS2(:,3),'-','DisplayName','XS2','linewidth',2);
plot3(XS2(:,1),XS2(:,2),XS2(:,3),'.','HandleVisibility','off','markersize',5);

plot3(XSix(:,35),XSiy(:,35),XSivq(:,35),'-','DisplayName','interpolated');
plot3(XSix(:,35),XSiy(:,35),XSivq(:,35),'.','HandleVisibility','off','LineWidth',5);

plot3(XSix,XSiy,XSivq,'.','HandleVisibility','off','LineWidth',5);

view(-105,32)
legend
xlabel('X'), ylabel('Y'), zlabel('Z')

As can be seen, most of the values at XS1 are not used during the interpolation. I hope someone can help me with this problem.

[dpb]

Just guessing without really digging into it, I'd guess

F.ExtrapolationMethod = 'none';

is prevents having any of those other values and the root cause is that one of the two interpolating variables is outside the dually-bounded area.

[Martin]

Thank you very much for your reply. I have tested your hint, unfortunately it does not solve the problem.

Answer 1

[Karim]

You could try another approach, by first fitting a curve trough the points of the two lines. Then evaluating those line alone equal distances. And using those points you can interpolate manually in a simple loop. See below for a demonstration

% load the data for the curves (only to make the answer a bit shorter, this
% holds the same data as the XS1 and XS2 grid in the OP's question)
load('X_data.mat') 

% grid for the interpolation
n = 50;

% flip the direction of the curve
% this is done so that we can perform the interpolation in a consistent direction
XS2up = flipud(XS2); 

% fit a curve trough the points, and evaluate alone "n" equal distances 
pt1 = interparc(n,XS1(:,1),XS1(:,2),XS1(:,3),'linear');
pt2 = interparc(n,XS2up(:,1),XS2up(:,2),XS2up(:,3),'linear');

% intialize some placeholders for the interpolated grid
GridX = zeros(n,n);
GridY = zeros(n,n);
GridZ = zeros(n,n);

% fill in the first and last set of points (i.e. the curves themself)
GridX(:,1) = pt1(:,1); GridX(:,end) = pt2(:,1);
GridY(:,1) = pt1(:,2); GridY(:,end) = pt2(:,2);
GridZ(:,1) = pt1(:,3); GridZ(:,end) = pt2(:,3);

% perform the interpolation
for i = 1:(n-1)
    GridX(:,i) = pt1(:,1) + (pt2(:,1)-pt1(:,1)) * i/n;
    GridY(:,i) = pt1(:,2) + (pt2(:,2)-pt1(:,2)) * i/n;
    GridZ(:,i) = pt1(:,3) + (pt2(:,3)-pt1(:,3)) * i/n;
end

% plot the data
figure
hold on
plot3(XS1(:,1),XS1(:,2),XS1(:,3),'r','linewidth',2);
plot3(XS2(:,1),XS2(:,2),XS2(:,3),'b','linewidth',2);
surf(GridX,GridY,GridZ,'EdgeAlpha',0.1)
view([50,12])
xlabel('X'), ylabel('Y'), zlabel('Z')

% make an extra view, since we can't see data in 3D using the online editor
figure
hold on
plot3(XS1(:,1),XS1(:,2),XS1(:,3),'r','linewidth',2);
plot3(XS2(:,1),XS2(:,2),XS2(:,3),'b','linewidth',2);
surf(GridX,GridY,GridZ,'EdgeAlpha',0.1)
view([71,33])
xlabel('X'), ylabel('Y'), zlabel('Z')

Note: this answer makes use of the following file exchange submission: interparc - File Exchange - MATLAB Central (mathworks.com)

[Martin]

Thank you for your new approach. It does what I was looking for, great job!