Outlier detection and removal of 3D vector fields

A filter for 3D vector fields. A single parameter is used in the code to detect the spurious vectors from 3D randomly distributed vector fil

function [ output_data ] = Outlier_detection( Point_coordinate, Displacement, Threshold )
% This code is used to detect spurious vectors from a set of 3D randomly distributed vectors.
%
% By Zhuang Cheng, City University of Hong Kong. Please cite the following paper when you use this code.
%
% Cheng, Z. and Wang, J., 2018. Quantification of the strain field of sands based on X-ray micro-tomography: A comparison between a grid-based method and a mesh-based method. Powder Technology.
%
%
%
% Point_coordinate (N*3) is the scatter points' coordinates, i.e., (x1, y1, z1), (x2, y2, z2),...

% Displacement (N*3) is the cooresponding displacement vector (Ui, Vi, Wi) of the scatter point (xi, yi, zi)

% Threshold is a parameter used to judge whether a vector is an outlier or not. Threshold is suggested to be within the range of 3~4. Generally, a lower Threshold leads to a larger number of outliers detected.

% output_data is a matrix with the dimension of M*6 (the first 3 columns are the coordinates of the scatter points that are not outliers. The following 3 columns are the corresponding displacements.)

%
% For detailed explanation of the code, please refer to:
%
% Cheng, Z. and Wang, J., 2018. Quantification of the strain field of sands based on X-ray micro-tomography: A comparison between a grid-based method and a mesh-based method. Powder Technology.
% https://www.sciencedirect.com/science/article/pii/S0032591018310854
% https://www.researchgate.net/publication/329548259_Quantification_of_the_strain_field_of_sands_based_on_X-ray_micro-tomography_A_comparison_between_a_grid-based_method_and_a_mesh-based_method
%
dt=delaunayTriangulation(Point_coordinate);
element=dt.ConnectivityList;% P*4, each row denotes the index (i.e., the row number of Point_coordinate) of the four vertex of a tetrahedron.
P1=element(:,1);
P2=element(:,2);
P3=element(:,3);
P4=element(:,4);
for i=1:1:size(Point_coordinate,1)
row1=find(P1==i);
neighbour1=element(row1,2:4);
n1=neighbour1(:);

row2=find(P2==i);
neighbour2=element(row2,3:4);
neighbour22=element(row2,1);
n2=neighbour2(:);
n22=neighbour22(:);
n222=union(n2,n22);

row3=find(P3==i);
neighbour3=element(row3,1:2);
neighbour33=element(row3,4);
n3=neighbour3(:);
n33=neighbour33(:);
n333=union(n3,n33);

row4=find(P4==i);
neighbour4=element(row4,1:3);
n4=neighbour4(:);

n=union(union(union(n1,n222),n333),n4);
vertex(i).neighbour=n; %Data recording the neighbour vertex of the centered vertex i.
disp=Displacement(n,:);
vertex(i).nei_disp=disp;%Data recording the displacement of the neighbour vertex
distance=repmat(Point_coordinate(i,:),length(n),1);
distance=distance-Point_coordinate(n,:);
vertex(i).nei_dist=sqrt((distance(:,1)).^2+(distance(:,2)).^2+(distance(:,3)).^2);%Data recording the distance between the center vertex and the neighbour vertex.

n1=[]; n2=[]; n22=[]; n222=[]; n3=[]; n33=[]; n333=[]; n4=[]; n=[];
neighbour1=[]; neighbour2=[]; neighbour22=[]; neighbour3=[]; neighbour33=[]; neighbour4=[];
row1=[]; row2=[]; row3=[]; row4=[];
disp=[]; distance=[];
end

for j=1:1:size(Point_coordinate,1)
dist=cat(1,vertex(j).nei_dist);
med_d=median(dist);
elapse=(sqrt(med_d^2+0.4)-med_d)/2; %Refer to the work of Duncan, J., Dabiri, D., Hove, J. and Gharib, M., 2010. Universal outlier detection for particle image velocimetry (PIV) and particle tracking velocimetry (PTV) data. Measurement Science and Technology, 21(5), p.057002.
for k=1:1:3
U0=Displacement(j,k);
A=U0/(elapse+med_d);
Displ=cat(1,vertex(j).nei_disp);
Ui=Displ(:,k);
B_array=Ui./(dist+elapse);
B=median(B_array);
r(k)=abs(A-B)/(median(abs(B_array-B))+elapse);
end
Displ=[]; Ui=[]; B_array=[];
if norm(r)>Threshold
outlier(j)=1;
else
outlier(j)=0;
end
dist=[];
end

output_data(:,1:3)=Point_coordinate(find(outlier==0),:);
output_data(:,4:6)=Displacement(find(outlier==0),:);

end