function neighbours = regionQuery(P, n, E)
    
    E2 = E^2;   
    [dim, Npts] = size(P);
    neighbours = [];
    
    for i = 1:Npts
        if i ~= n
            % Test if distance^2 < E^2 
            v = P(:,i)-P(:,n);
            dist2 = v'*v;
            if dist2 < E2 
               neighbours = [neighbours i];     
            end
        end
    end
    
end % of regionQuery