KD-Treeを用いたRevitファイル間の要素比較をC#で実装する

こんにちは！
今週のAKARI Tech BlogはDX Solution事業部の池之上が担当します。

燈ではRevit周辺の開発を数多く行っています。
建築の実務で使われるモデルは要素数が数万、数十万単位になるため高速なアドインの実装には一工夫が必要です。
今回は異なるファイル間で要素同士を紐付ける際に用いたKD-Treeという手法をご紹介します。

背景
ナイーブな線形探索
KD-Treeとは？
実装
速度比較
まとめ
We're Hiring!

背景

BIM運用において修正前のモデルAと修正後のモデルB、あるいは意匠モデルと構造モデルの間で同じ位置にある特定の要素（例えば1Fにある同じ座標の柱など）を特定したいケースが多々あります。
しかしインスタンスのIdはファイル固有であるためモデルを跨いで共通のキーとして使うことができない場合があります。
結局のところ配置点のような空間情報を頼りに探すことになりますが、単純な全件ループ処理では数万の要素を抱える実プロジェクトでは処理時間が膨大になりフリーズしてしまうことも珍しくありません。
本記事ではこのボトルネックを解決するために計算幾何学のアルゴリズムである KD-Tree をC#で実装し、Revitアドインで高速に要素を特定する手法を紹介します。

ナイーブな線形探索

最も単純な方法は、モデルAの各要素に対して、モデルBの全要素を一つずつ走査する総当たり方式です。
例えば柱を検索するなら次のような処理が考えられます。

var targetPoints = new List<XYZ>
{
    new XYZ(0,0,0),
    new XYZ(10,0,0),
    new XYZ(20,0,0),
}; // 実際はモデルAの柱の配置点など
double threshold = 1; // 検索範囲(フィート)

using (var collector = new FilteredElementCollector(documentB))
{
    var detectedElements = collector
    .OfCategory(BuiltInCategory.OST_StructuralColumns)
    .WhereElementIsNotElementType()
    .Where(e => e.Location is LocationPoint);

    foreach (var point in targetPoints)
    {
        var nearbyElement = detectedElements
        .FirstOrDefault(e =>
        {
            var locPoint = e.Location as LocationPoint;
            return locPoint.Point.DistanceTo(point) < threshold;
        });
　　// nearByElementの処理
    }
}

今回はKD-Treeによりこれを高速化します。

KD-Treeとは？

KD-Tree（K-Dimensional Tree）は、多次元空間に存在する点データを効率よく管理・検索するためのデータ構造です。
2次元や3次元の座標データを扱う場面で特に有効で、空間内の近い点やある範囲に収まる点を高速に見つけることができます。
よく利用されるのは下記のような用途です。

最近傍点検索：

　ある座標に最も近い点を素早く探したい場合。

範囲検索：

　指定した半径や範囲内に存在する点をまとめて取得したい場合。

空間インデックスとして：

　3D CAD、ゲーム、地図アプリ、ロボットの経路探索など、空間上のオブジェクト管理や検索が必要な場合。

3次元のKD-Treeは各ノードごとにX軸・Y軸・Z軸などの座標軸で空間を分割し再帰的に木構造を作っていきます。
例えば
1段目はX軸で分割
2段目はY軸で分割
3段目はZ軸で分割
4段目はまたX軸…
というように分割軸を順番に切り替えながら構築します。
探索時に探索不要な部分を効率よく枝刈りできるのが特徴です。

実装

KD-Treeのライブラリは多数存在するのですがDLL競合に悩むこともあるため今回は自前実装を行ってみます。
Revitアドインのコマンドで使うことが目的であるため、要素の追加/削除やそれに伴う再構築を考慮しない簡易的なものとします。
またElementはカテゴリ次第で特定の配置点を持たないため、IElementWrapperを定義し要素ごとに独自の位置取得メソッドを実装することにします。
例えば点ベースファミリならそのままLocationPointを使えますし、線ベースファミリなら両端の中点を返すようにしておいて範囲検索後に厳密判定を行うことが考えられるでしょう。(ブレースのように中点で交差する要素がある場合にはさらに一工夫必要になります。)

internal interface IElementWrapper
{
    XYZ GetPosition();
}

internal class KDTree<T> where T : IElementWrapper
{
    private readonly Node _root;
    private readonly List<T> _failedItems = new List<T>();

    public KDTree(IEnumerable<T> items)
    {
        var pts = new List<NodeItem>();

        foreach (var it in items ?? Enumerable.Empty<T>())
        {
            try
            {
                var p = it.GetPosition();
                if (p == null) continue;
                pts.Add(new NodeItem { Item = it, Point = p });
            }
            catch
            {
                _failedItems.Add(it);
            }
        }

        _root = Build(pts, 0);
    }

    // 最近傍を1つ返す
    public T Nearest(XYZ query)
    {
        if (query == null || _root == null) return default;
        T bestItem = default;
        double bestDist2 = double.MaxValue;
        NearestSearch(_root, query, ref bestItem, ref bestDist2);
        return bestItem;
    }

    // 最大距離制限付きで最近傍を1つ返す
    public T Nearest(XYZ query, double maxDistance)
    {
        if (query == null || _root == null) return default;
        T bestItem = default;
        double bestDist2 = maxDistance * maxDistance;
        NearestSearch(_root, query, ref bestItem, ref bestDist2);
        return bestItem;
    }

    // 半径内のすべてを返す
    public IList<T> RadiusSearch(XYZ query, double radius)
    {
        var results = new List<T>();
        if (query == null || _root == null || radius < 0) return results;
        double r2 = radius * radius;
        RadiusSearch(_root, query, r2, results);
        return results;
    }

    private class Node
    {
        public NodeItem Item;
        public Node Left;
        public Node Right;
        public int Axis; // 0=x,1=y,2=z
    }

    private class NodeItem
    {
        public T Item;
        public XYZ Point;
    }

    private Node Build(List<NodeItem> pts, int depth)
    {
        if (pts == null || pts.Count == 0) return null;
        int axis = depth % 3;
        int cmp(NodeItem a, NodeItem b)
        {
            switch (axis)
            {
                case 0: return a.Point.X.CompareTo(b.Point.X);
                case 1: return a.Point.Y.CompareTo(b.Point.Y);
                default: return a.Point.Z.CompareTo(b.Point.Z);
            }
        }

        pts.Sort(cmp);
        int mid = pts.Count / 2;

        var node = new Node
        {
            Axis = axis,
            Item = pts[mid]
        };

        var leftList = pts.Take(mid).ToList();
        var rightList = pts.Skip(mid + 1).ToList();

        node.Left = Build(leftList, depth + 1);
        node.Right = Build(rightList, depth + 1);
        return node;
    }

    private void NearestSearch(Node node, XYZ query, ref T bestItem, ref double bestDist2)
    {
        if (node == null) return;
        double d2 = DistanceSquared(node.Item.Point, query);
        if (d2 < bestDist2)
        {
            bestDist2 = d2;
            bestItem = node.Item.Item;
        }

        int axis = node.Axis;
        double diff = GetCoordinate(query, axis) - GetCoordinate(node.Item.Point, axis);
        Node first = diff <= 0 ? node.Left : node.Right;
        Node second = diff <= 0 ? node.Right : node.Left;

        if (first != null) NearestSearch(first, query, ref bestItem, ref bestDist2);

        if (second != null && diff * diff < bestDist2)
        {
            NearestSearch(second, query, ref bestItem, ref bestDist2);
        }
    }

    private void RadiusSearch(Node node, XYZ query, double r2, List<T> results)
    {
        if (node == null) return;
        double d2 = DistanceSquared(node.Item.Point, query);
        if (d2 <= r2) results.Add(node.Item.Item);

        int axis = node.Axis;
        double diff = GetCoordinate(query, axis) - GetCoordinate(node.Item.Point, axis);

        if (diff <= 0)
        {
            if (node.Left != null) RadiusSearch(node.Left, query, r2, results);
            if (node.Right != null && diff * diff <= r2) RadiusSearch(node.Right, query, r2, results);
        }
        else
        {
            if (node.Right != null) RadiusSearch(node.Right, query, r2, results);
            if (node.Left != null && diff * diff <= r2) RadiusSearch(node.Left, query, r2, results);
        }
    }

    private static double GetCoordinate(XYZ p, int axis)
    {
        switch (axis)
        {
            case 0: return p.X;
            case 1: return p.Y;
            default: return p.Z;
        }
    }

    private static double DistanceSquared(XYZ a, XYZ b)
    {
        double dx = a.X - b.X;
        double dy = a.Y - b.Y;
        double dz = a.Z - b.Z;
        return dx * dx + dy * dy + dz * dz;
    }
}

実装のポイントは検索時の枝刈りです。
二乗距離のまま保持しているので少しわかりづらいかもしれません。diffはクエリ点から分割平面までの距離を符号付きで持っています。
分割平面までの距離が現在の暫定最小距離よりも大きい場合、反対側の空間には現在の暫定解より近い点が存在しないため検索をする必要がなくなります(枝刈り)。
これにより木構造の半分以上のノードを無視することができ計算量を削減できます。
RadiusSearchの場合も枝刈りの考え方は同様ですが検索半径が常に固定値になります。

private void NearestSearch(Node node, XYZ query, ref T bestItem, ref double bestDist2)
{
    if (node == null) return;
    double d2 = DistanceSquared(node.Item.Point, query);
    if (d2 < bestDist2)
    {
        bestDist2 = d2;
        bestItem = node.Item.Item;
    }

    int axis = node.Axis;
    double diff = GetCoordinate(query, axis) - GetCoordinate(node.Item.Point, axis);
    Node first = diff <= 0 ? node.Left : node.Right;
    Node second = diff <= 0 ? node.Right : node.Left;

    if (first != null) NearestSearch(first, query, ref bestItem, ref bestDist2);

    if (second != null && diff * diff < bestDist2)
    {
        NearestSearch(second, query, ref bestItem, ref bestDist2);
    }
}

3次元のままだと検索の流れが理解しづらいのですが、2次元で表現するとこのような感じでしょうか。
最短半径の円が分割平面をまたがる場合は左右ノードとも検索、またがらない場合はかからない側のノードの検索を省略できます。

速度比較

最後に速度比較をしてみましょう!
対象はRevitのサンプル構造モデルの梁としました。下記が比較結果です。

要素数: 1061
1. 線形探索: 2527ms
2. KDTree:
　構築時間: 22ms
　検索時間: 3ms
　合計時間: 25ms

結果として100倍程度の高速化になりました!
要素数が増えればさらに差は大きくなるでしょう。

まとめ

この記事ではRevitアドインにおけるKD-Treeを用いたファイル間の要素比較について紹介しました。
高機能な外部ライブラリを利用する方が良いケースもあると思いますが、今回のようなシンプルな実装であればアドインのDLL競合に悩まされることなくプロジェクトに組み込めるのが大きなメリットです。
BIMデータの活用が進むにつれ、数万単位の要素を扱うシーンは今後さらに増えていくでしょう。そうした際のパフォーマンス改善の一手としてこうしたアルゴリズムを取り入れてみてはいかがでしょうか。

We're Hiring!

燈ではRevitに明るいエンジニアを募集しております!

興味がある方はぜひカジュアル面談でお話ししましょう!
akariinc.co.jp