Interpolation of two-dimensional within a given space window

Dear all

I am attaching a dataset, called "Points.txt", where the first and second columns represent the spatial coordinates along the x- and y-th directions, respectively, and the third column the value of a variable at each listed pair of points (x,y). I was wondering which could be a good way to increase the number of points mainly inside the hexagonal shape within the current listed points are restricted to:

I have tried the following:

pgon=nsidedpoly(6,'Center',[0 0],'SideLength',1);
%% Let's rearrange the data associated to the interlayer isotropic exchange, Jinter
% Let's load the preliminar data
Jinter=cell2mat(struct2cell(load(strcat(pwd,'\Preliminar-Data\Jinter.mat'))));
% Let's construct the unit vectors of the in-plane unit cell
x=unique(Jinter(:,1));
y=unique(Jinter(:,2));
a1_points=[max(x) 0];
a2_points=[-0.5 max(y)];
a1_modulus=sqrt(dot(a1_points,a1_points));
a2_modulus=sqrt(dot(a2_points,a2_points));
a1=a1_points./a1_modulus;
a2=a2_points./a2_modulus;
% Let's construct the supercell for Jinter
diff_y=max(y)-min(y);
counter=0;
for dy=[-1 0 1]
    for i=1:length(Jinter(:,1))
        if (Jinter(i,1:2)++dy*[0 1]*diff_y)~=any(Jinter(:,1:2))
            counter=counter+1;
            Jinter_extended(counter,1:2)=Jinter(i,1:2)+dy*[0 1]*diff_y;
            Jinter_extended(counter,3)=Jinter(i,3);
        end
    end
end
counter=length(Jinter_extended(:,1));
for i=1:length(Jinter(:,1))
    if (Jinter(i,1:2)+[a1(1)+abs(a2(1)) abs(a2(2))])~=any(Jinter_extended(:,1:2))
        counter=counter+1;
        Jinter_extended(counter,1:2)=Jinter(i,1:2)+[a1(1)+abs(a2(1)) abs(a2(2))];
        Jinter_extended(counter,3)=Jinter(i,3);
    end
end
counter=length(Jinter_extended(:,1));
for i=1:length(Jinter(:,1))
    if (Jinter(i,1:2)+[a1(1)+abs(a2(1)) -abs(a2(2))])~=any(Jinter_extended(:,1:2))
        counter=counter+1;
        Jinter_extended(counter,1:2)=Jinter(i,1:2)+[a1(1)+abs(a2(1)) -abs(a2(2))];
        Jinter_extended(counter,3)=Jinter(i,3);
    end
end
counter=length(Jinter_extended(:,1));
for i=1:length(Jinter(:,1))
    if (Jinter(i,1:2)+[-(a1(1)+abs(a2(1))) abs(a2(2))])~=any(Jinter_extended(:,1:2))
        counter=counter+1;
        Jinter_extended(counter,1:2)=Jinter(i,1:2)+[-(a1(1)+abs(a2(1))) abs(a2(2))];
        Jinter_extended(counter,3)=Jinter(i,3);
    end
end
counter=length(Jinter_extended(:,1));
for i=1:length(Jinter(:,1))
    if (Jinter(i,1:2)+[-(a1(1)+abs(a2(1))) -abs(a2(2))])~=any(Jinter_extended(:,1:2))
        counter=counter+1;
        Jinter_extended(counter,1:2)=Jinter(i,1:2)+[-(a1(1)+abs(a2(1))) -abs(a2(2))];
        Jinter_extended(counter,3)=Jinter(i,3);
    end
end
clear a1 a2 a1_modulus a2_modulus a1_points a2_points Jinter diff_y x y dy
% Let's plot it
u17=figure(17)
scatter(Jinter_extended(:,1),Jinter_extended(:,2),[],Jinter_extended(:,3),'o','filled');
colormap(gca,bluewhitered(256));
clim([min(Jinter_extended(:,3)) max(Jinter_extended(:,3))]);
clr17=colorbar;
box on;
pbaspect([1 6/4 1]);
clear u17 clr17
% Now, let's interpolate data
x_extended=linspace(min(Jinter_extended(:,1)),max(Jinter_extended(:,1)));
y_extended=linspace(min(Jinter_extended(:,2)),max(Jinter_extended(:,2)));
Jinter_interpolated=RegularizeData3D(Jinter_extended(:,1),Jinter_extended(:,2),Jinter_extended(:,3),x_extended,y_extended,'interp','bicubic','smoothness',1e-4);
x_interpolation=linspace(-1,1,201);
y_interpolation=linspace(-sqrt(3)/2,sqrt(3)/2,201);
[X,Y]=meshgrid(x_interpolation,y_interpolation);
interpolated_Jinter=interp2(x_extended,y_extended,Jinter_interpolated,X,Y,'cubic');
clear X Y Jinter_interpolated x_extended y_extended Jinter_extended
% Let's plot it
u18=figure(18)
uimagesc(x_interpolation,y_interpolation,interpolated_Jinter);
colormap(gca,bluewhitered(256));
% hold on
axis xy;
% plot(pgon);
clim([min(min(interpolated_Jinter)) max(max(interpolated_Jinter))]);
clr18=colorbar;
box on;
pbaspect([1 sqrt(3)/2 1]);
clear u18 clr18
% Now, let's rearrange and save the data
x_intermediate=linspace(-100,100,201);
y_intermediate=linspace(-100,100,201);
counter=0;
for i=1:length(x_interpolation)
    for j=1:length(y_interpolation)
        counter=counter+1;
        Jinter_save(counter,1)=x_intermediate(i);
        Jinter_save(counter,2)=y_intermediate(j);
        Jinter_save(counter,3)=interpolated_Jinter(i,j);
    end
end
writematrix(Jinter_save,strcat(pwd,'\Interpolated-Data\Interpolated_J_Inter.txt'),'Delimiter','space');
clear interpolated_Jinter Jinter_save x_interpolation y_interpolation counter i j x_intermediate y_intermediate

where the "Jinter.mat" is exactly the same dataset as in "Points.txt". The problem with this is that, for example, if I take a look at the generated "Interpolated_J_Inter.txt", the value for (0,0) differs notably.

Any ideas?

[Umar]

Hi @Richard Wood,

Try using a different interpolation method, such as 'linear' or 'nearest', to see if it produces more accurate results.

interpolated_Jinter = interp2(x_extended, y_extended, Jinter_interpolated, X, Y, 'linear');

For guidance on linear or nearest, please refer to https://www.mathworks.com/help/matlab/ref/interp2.html

Also, try experimenting with different parameters, such as the number of points in the x_interpolation and y_interpolation vectors, to achieve a finer resolution in the interpolated data. Remember to save the modified code and run it to generate the updated "Interpolated_J_Inter.txt" file. Check the values at (0,0) to see if they are closer to the expected values.

I hope this helps! Let me know if you have any further questions.

Answer 1

[KSSV]

T = readtable('https://in.mathworks.com/matlabcentral/answers/uploaded_files/1744881/Points.txt');
x = T.(1) ;
y = T.(2) ;
z = T.(3) ; 
dx = 0.001 ;  dy = 0.001;    % change it to desired 
[X,Y] = meshgrid(min(x):dx:max(x),min(y):dy:max(y)) ; 
Z = griddata(x,y,z,X,Y) ; 

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

%
idx = boundary(x,y) ; 
[in,on] = inpolygon(X,Y,x(idx),y(idx)) ; 
Z(~(in | on)) = NaN ; 

h = pcolor(X,Y,Z);
h.EdgeColor = 'none' ;

[Umar]

@Richard Wood,

@KSSV has shared his input, it seems that he is familiar with unique function but did want to use it to ensure the integrity of data but if you choose to get rid of warning message as it appears in @KSSV code,then this is what modified code would look like.

x = T.(1); y = T.(2); z = T.(3);
dx = 0.001; dy = 0.001;
% Remove duplicate points
[uniquePoints, ia, ~] = unique([x, y], 'rows');
x = uniquePoints(:, 1);
y = uniquePoints(:, 2);
z = z(ia);
[X, Y] = meshgrid(min(x):dx:max(x), min(y):dy:max(y));
Z = griddata(x, y, z, X, Y);
idx = boundary(x,y) ;
[in,on] = inpolygon(X,Y,x(idx),y(idx)) ;
Z(~(in | on)) = NaN ;
h = pcolor(X,Y,Z);
h.EdgeColor = 'none' ;

For guidance on this function, please refer to,

https://www.mathworks.com/help/matlab/ref/double.unique.html

@KSSV, thanks for your contribution, really appreciated.