function [DeltaX, R, Sigma] = build_feature(part, n_neighbours)
    % Input:
    %   part            3 x n_part      coordinates of particles
    %   n_neighbours    1 x 1           number of neighbours to use in
    %                                   building topological feature set
    %
    % Output:
    %   DeltaX_hat      3 x n_neighbours x n_part
	% the relative position of n_neighbours neighbouring particles,
	% in a coordinate system which is aligned and scaled with 
	% the principle axes of variation of the neighbours
	%
	% i.e. positions of neighbours up to am arbitrary rotation and scale
    % 
    % for a set of neighbours X = [x1, x2, x3, ... xn]
    %   relative to a particle at x = 0
    % the local coordinate system R is defined such that
    %   [R, Sigma] = eig(X * X.')
	% we then define 
	%	X^ = inv(sqrt(Sigma)) * R.' * X
	% so that
	%	X^ * X^' = inv(sqrt(Sigma)) * R.' * X * X.' * R * inv(sqrt(Sigma))
    %            = inv(sqrt(Sigma)) * R.' * R * Sigma * R.' R * inv(...)
	%		     = eye(3)
	
    %% get indices of nearest neighbours
    n_part  = size(part, 2); 
    [i_nn,r]= knnsearch(part.', part.', 'K', n_neighbours+1);
    % exclude self
	i_nn = i_nn(:, 2:end);
	r    = r(:,2:end);
	
    %% get displacement to nearest neighbours
    DeltaX  = reshape(part(:,i_nn.'), [3 n_neighbours n_part]) ...
            - reshape(part, [3 1 n_part]);
       
%     %% define rotation invariant coordinate system
%     % based on singular value decomposition
%     % for n = 1 : n_part
%     %   [U,sig] = svd(DeltaX(:,:,n))
%     % end
%     
%     % compute singular value decomposition of displacement vectors
%     DeltaX_var	= mtimesx(DeltaX, 'n', DeltaX, 't');
%     % calculate eigenvectors of DeltaX * DeltaX.', ordering eigenvalues in
%     % descending order of magnitude
%     [Sigma, R]   = ndsyev(single(DeltaX_var), 's-');
% 	Sigma		= double(Sigma);
% 	R			= double(R);
% 	% apply transformation: rotation and scale
% 	DeltaX_rot = mtimesx(R, 't', DeltaX, 'n');
% 	%DeltaX_rot = DeltaX_rot ./ sqrt(reshape(Sigma, [3 1 n_part]));
% 	% check: DeltaX_rot * DeltaX_rot' = eye(3)
% 	%DeltaX_rot_var = mtimesx(DeltaX_rot, 'n', DeltaX_rot, 't');
end
