lundi 14 septembre 2015

Nearest neighborS search by Binary Partition Tree

There are a lot of excellent references on Internet, one of them is:

The missing part: some code to play around. Let's try to make something generic. What to we need to make a class partitionable ?

    public interface IIsPartionableByBT {
        Int32 Id { get; }

        float BTX { get; }
        float BTY { get; }
        float BTZ { get; }

        float DistTo(IIsPartionableByBT c);
    }

    public delegate float BTCooProjector(IIsPartionableByBT coo);

We also need Cells/Nodes for the tree:

    public class MbtBSPNode<T> where T :  IIsPartionableByBT {
        public MbtBSPNode<T> Parent { get; set; }
        public T Cell { get; set; }
        public MbtBSPNode<T> Left { get; set; }
        public MbtBSPNode<T> Right { get; set; }

        public ICollection<MbtBSPNode<T>> Cells { get; set; }

        public UInt32 Depth { get { if (Parent == null) return 0; return Parent.Depth + 1; } }
    }

And the tree himself:

    public class MbtBSPTree<T> where T : class, IIsPartionableByBT {
        public MbtBSPNode<T> Root { get; private set; }

        public UInt16 SplitDepthBy { get; private set; }

        public UInt16 LimitDepthTo { get; private set; }

        public MbtBSPTree(IEnumerable<T> l, UInt16 splitDepthBy = 0, UInt16 limitDepthTo = 0) {
            if (SplitDepthBy == 0) {
                if (l.Min(x => x.BTZ) == l.Max(x => x.BTZ)) {
                    SplitDepthBy = 2;
                } else {
                    SplitDepthBy = 3;
                }
            } else {
                SplitDepthBy = splitDepthBy;
            }
            LimitDepthTo = limitDepthTo;

            Root = BuildBSPTree(l, 0);
        }

        public BTCooProjector DefProjector(UInt16 axis) {
            switch (axis) {
                case 0:
                    return (n) => n.BTX;
                case 1:
                    return (n) => n.BTY;
                case 2:
                    return (n) => n.BTZ;
                default:
                    throw new Exception("valeur de split de profondeur non gérée");
            }
        }

        public MbtBSPNode<T> BuildBSPTree(IEnumerable<T> l, UInt16 depth, MbtBSPNode<T> p = null) {
            if (l.Count() == 0)
                return null;
            UInt16 axis = (UInt16)(depth % SplitDepthBy);

            MbtBSPNode<T> n = new MbtBSPNode<T>();
            n.Parent = p;

            if (LimitDepthTo > 0 && depth == LimitDepthTo) {
                n.Cells = new List<MbtBSPNode<T>>();
                foreach (T cc in l) {
                    n.Cells.Add(new MbtBSPNode<T> { Cell = cc });
                }
            } else {
                Double med = Double.MinValue;
                List<T> ll = new List<T>();
                List<T> lr = new List<T>();

                BTCooProjector proj = DefProjector((UInt16)(depth % SplitDepthBy));
                med = l.Median(x => proj(x));
                foreach (T c in l) {
                    if (proj(c) < med) {
                        ll.Add(c);
                    } else {
                        if (proj(c) == med) {
                            if (n.Cell == null)
                                n.Cell = c;
                            else
                                ll.Add(c);
                        } else {
                            lr.Add(c);
                        }
                    }
                }

                if (depth <= 2) {
                    //this will lead to 8 threads
                    Task<MbtBSPNode<T>> tl = Task<MbtBSPNode<T>>.Factory.StartNew(() =>
                        BuildBSPTree(ll, (UInt16)(depth + 1), n));
                    n.Right = BuildBSPTree(lr, (UInt16)(depth + 1), n);
                    //synchro thread
                    n.Left = tl.Result;
                } else {
                    n.Left = BuildBSPTree(ll, (UInt16)(depth + 1), n);
                    n.Right = BuildBSPTree(lr, (UInt16)(depth + 1), n);
                }
            }

            return n;
        }

#if DEBUG
        public Int32 ScannedNodeNumber { get; set; }
#endif

        public IEnumerable<MbtBSPNode<T>> GetNearestNodes(T cel, float radius) {
#if DEBUG
            ScannedNodeNumber = 0;
#endif
            MbtBSPNode<T> current;
            Queue<MbtBSPNode<T>> q = new Queue<MbtBSPNode<T>>();
            q.Enqueue(Root);

            while (q.Count > 0) {
                current = q.Dequeue();
                BTCooProjector proj;

                if (current.Cell != null) {
#if DEBUG
                    ScannedNodeNumber++;
#endif
                    if (cel.DistTo(current.Cell) <= radius && cel.Id != current.Cell.Id) {
                        yield return current;
                    }

                    proj = DefProjector((UInt16)(current.Depth % SplitDepthBy));

                    if (proj(cel) - radius <= proj(current.Cell)) {
                        if (current.Left != null)
                            q.Enqueue(current.Left);
                    }
                    if (proj(cel) + radius > proj(current.Cell)) {
                        if (current.Right != null)
                            q.Enqueue(current.Right);
                    }
                }

                if (current.Cells != null) {
                    foreach (MbtBSPNode<T> cc in current.Cells) {
#if DEBUG
                        ScannedNodeNumber++;
#endif
                        if (cel.Id != cc.Cell.Id && cel.DistTo(cc.Cell) <= radius)
                            yield return cc;
                    }
                }

            }
        }

One can find a median algorithm in here: