Hi Gerrit,
Happy new year and something wrong with my gmail, so I just notice your reply.
The code for the neighbor calculation you provide seems to simple and I decide to post a new topic to ask, but thanks anyway.
For Sionelli's code, I also debug it and get it to run. The code I used is as follow to plot it out
cs=crystalSymmetry('cubic'); % I change it to crystalSymmetry
ss=specimenSymmetry('triclinic'); % I change it to specimenSymmetry
ebsd=loadEBSD_generic('ebsd_file.txt', 'CS', cs, 'SS', ss, 'ColumnNames', {'Euler1', 'Euler2', 'Euler3', 'x','y'}, 'Bunge');
plot(ebsd('indexed'),ebsd('indexed').orientations);% To plot out the orientation map
the result is

However, even the noisy result I can not got, it is exactly as the same as the original.

I suspect maybe it is due to the different RD/TD definition between different software, as the stupid part is that I use INCA to do the scan, and the OIM could only receive HKL channel 5 data.
One question about the Sionelli's code, do you know why he choose the threshold as 30 at the beginning? It seems set as a big range to cover as much grains as possible, not like the definition like the misorientation between two grains, right?
By the way, the good news is that I have find out Sionelli's email and he is now working in Nottingham University, he replied in November to review his code but not receive any reply from him yet. If there any news from him, I will let you know right now.
The following is the boundary I plot by using my own code and the different axis definition between different software.

If you could figure anything wrong from here or new idea, please let me know.
Best regards,
Shuang