東京都の点群データをUnityに読み込もう！

まず、これを見ていただきたい。
https://www.metro.tokyo.lg.jp/tosei/hodohappyo/press/2023/09/01/16.html

一部地域の点群データが公開されました。

こういった感じの結構精度の高いデータになります。

点の数が多すぎて、画像のように見えるかもしれませんが、手前の方をよく見ると点々になっていることが判ります。

これは、Unityに読み込んでみたくなりますね。

点群データについて、良く調べてみると.lasというファイルフォーマットでした。
なので、ネットでlasファイルの仕様について検索します。
そうすると、下記のようなpdfが見つかりました。
https://www.asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf

これを元に、インポーターを書きました。

LasImporter.cs

using System.Collections;
using System.Collections.Generic;
using System.IO;
using System.Text;
using UnityEngine;

public class LASImporter
{
    public void Import(string filename)
    {
        // Seekを多用しているため、FileStreamよりも、全bytes読み取ってからMemoryStreamのほうが速い(10倍ぐらい)
#if false
        using (FileStream fs = new FileStream(filename, FileMode.Open, FileAccess.Read))
        {
            Import(fs);
        }
#else
        var bytes = File.ReadAllBytes(filename);
        using (MemoryStream ms = new MemoryStream(bytes))
        {
            Import(ms);
        }
#endif
    }

    void Import(Stream stream)
    {
        BinaryReader br = new BinaryReader(stream);
        var fileSignature = br.ReadBytes(4);
        var fileSourceID = br.ReadUInt16();
        var globalEncoding = br.ReadUInt16();
        var guidData1 = br.ReadUInt32();
        var guidData2 = br.ReadUInt16();
        var guidData3 = br.ReadUInt16();
        var guidData4 = br.ReadBytes(8);
        var versionMajor = br.ReadByte();
        var versionMinor = br.ReadByte();
        var systemIdentifier = Encoding.UTF8.GetString(br.ReadBytes(32));
        Debug.Log(systemIdentifier);
        var generatingSoftware = Encoding.UTF8.GetString(br.ReadBytes(32));
        Debug.Log(generatingSoftware);
        var fileCreationDayOfYear = br.ReadUInt16();
        var fileCreationYear = br.ReadUInt16();
        var headerSize = br.ReadUInt16();
        Debug.Log("headerSize:" + headerSize);
        var offsetToPointData = br.ReadUInt32();
        Debug.Log("offsetToPointData:" + offsetToPointData);
        var numberOfVariableLengthRecords = br.ReadUInt32();
        Debug.Log("numberOfVariableLengthRecords:" + numberOfVariableLengthRecords);
        var pointDataFormatID = br.ReadByte();
        Debug.Log("pointDataFormatID:" + pointDataFormatID);
        var pointDataRecordLength = br.ReadUInt16();
        Debug.Log("pointDataRecordLength:" + pointDataRecordLength);
        var numberOfPointRecords = br.ReadUInt32();
        Debug.Log("numberOfPointRecords:" + numberOfPointRecords);
        var numberOfPointByReturn = new uint[5];
        for(int i=0;i<5;++i)
        {
            numberOfPointByReturn[i] = br.ReadUInt32();
            Debug.Log("numberOfPointByReturn[" + i + "]:" + numberOfPointByReturn[i]);
        }
        var xScaleFactor = br.ReadDouble();
        var yScaleFactor = br.ReadDouble();
        var zScaleFactor = br.ReadDouble();
        Debug.Log(xScaleFactor + "," + yScaleFactor + "," + zScaleFactor);
        var xOffset = br.ReadDouble();
        var yOffset = br.ReadDouble();
        var zOffset = br.ReadDouble();
        Debug.Log(xOffset + "," + yOffset + "," + zOffset);
        var maxX = br.ReadDouble();
        var minX = br.ReadDouble();
        var maxY = br.ReadDouble();
        var minY = br.ReadDouble();
        var maxZ = br.ReadDouble();
        var minZ = br.ReadDouble();
        Debug.Log(maxX + "," + maxY + "," + maxZ);
        Debug.Log(minX + "," + minY + "," + minZ);

        Debug.Log("Position:" + br.BaseStream.Position);
        Debug.Log("Length:" + br.BaseStream.Length);

        if (pointDataFormatID == 3)
        {
            br.BaseStream.Seek(offsetToPointData, SeekOrigin.Begin);
            var mesh = ReadPointData_Format3(br, (int)numberOfPointRecords,
                xScaleFactor,
                yScaleFactor,
                zScaleFactor,
                xOffset - minX,
                yOffset - minY,
                zOffset - minZ,
                pointDataRecordLength
                );
            GameObject o = new GameObject("las");
            var mf = o.AddComponent<MeshFilter>();
            mf.sharedMesh = mesh;
            var mr = o.AddComponent<MeshRenderer>();
        }
        else
        {
            Debug.LogError("Not support pointDataFormatID:" + pointDataFormatID);
        }

        Debug.Log("Last Position:" + br.BaseStream.Position);
    }

    Mesh ReadPointData_Format3(BinaryReader br,
        int numPoints,
        double xScaleFactor,
        double yScaleFactor,
        double zScaleFactor,
        double xOffset,
        double yOffset,
        double zOffset,
        int pointDataRecordLength
        )
    {
        var from = br.BaseStream.Position;
        Mesh mesh = new Mesh();
        Vector3[] points = new Vector3[numPoints];
        Color32[] colors = new Color32[numPoints];
        int[] indices = new int[numPoints];
        for (int i = 0; i < numPoints; ++i)
        {
            br.BaseStream.Seek(from + i * pointDataRecordLength, SeekOrigin.Begin);
            var x = br.ReadInt32(); // 4
            var y = br.ReadInt32(); // 8
            var z = br.ReadInt32(); // 12
            var intensity = br.ReadUInt16();    // 14
            var pack1 = br.ReadByte();  // 15
            var classification = br.ReadByte(); // 16
            var scanAngleRank = br.ReadByte();  // 17
            var userData = br.ReadByte(); // 18
            var pointSourceID = br.ReadUInt16();    // 20
            var gpsTime = br.ReadDouble();  // 28
            var red = br.ReadUInt16();  // 30
            var green = br.ReadUInt16();    //32
            var blue = br.ReadUInt16(); // 34

            points[i] = new Vector3(
                (float)(x * xScaleFactor + xOffset),
                (float)(z * zScaleFactor + zOffset),
                (float)(y * yScaleFactor + yOffset)
                );
            colors[i] = new Color(
                red / (float)ushort.MaxValue,
                green / (float)ushort.MaxValue,
                blue / (float)ushort.MaxValue,
                intensity / (float)ushort.MaxValue
                );
            indices[i] = i;
        }
        mesh.indexFormat = UnityEngine.Rendering.IndexFormat.UInt32;
        mesh.SetVertices(points);
        mesh.SetColors(colors);
        mesh.SetIndices(indices, MeshTopology.Points,0);

        return mesh;
    }
}

Debug.Logが残ってたりしますが、勝手に消してください。
これは、すごく適当にインポータを書いたので、ポイントデータFormat3にしか対応していません。
今回ダウンロードしたデータが、Format3だったためです。
lasには、ポイントデータFormat0～3があり、しっかりlas読み込みをしたいなら、先ほどのpdfをちゃんと読んで対応したら良いと思いますが、今回は適当にやりました。

上記のコードだけだと、インポートする処理だけが書かれたものなので、Unityのトップメニューから読み込めるようにしましょう。
下記のようなコードを、Editorフォルダに入れましょう。

Editor/LASImportMenu.cs

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEditor;

public class LASImportMenu 
{
    [MenuItem("LAS/LAS Import", false, 0)]
    public static void Import()
    {
        string filename = EditorUtility.OpenFilePanel("Import LAS", "", "las");
        if (string.IsNullOrEmpty(filename) == false)
        {
            LASImporter importer = new LASImporter();
            importer.Import(filename);
        }
    }
}

そうすると、Unityのトップメニュからーlasファイル選択して、読み込めるようになります。
マテリアル設定はしていないので、最初はシェーダーエラー紫の状態で出ると思います。
頂点カラーシェーダーの付いたマテリアルを設定してください。(Unlit/VertexColorなど)

感想

最初にリンクを張ったLAS仕様のpdfによると、pointDataFormatIDが３のときは、１点の情報量は34バイトで固定だとおもうのだけど、今回いくつかダウンロードしたlasファイルの内のいくつかは36バイトになっているものがあった。
そのため、１点ごとにSeekをするはめになり重くなっている。
FileStreamでSeekすると重い気がする。
高速化を考えて、インポート最初の処理を変えてみた。この部分。

        // Seekを多用しているため、FileStreamよりも、全bytes読み取ってからMemoryStreamのほうが速い(10倍ぐらい)
#if false
        using (FileStream fs = new FileStream(filename, FileMode.Open, FileAccess.Read))
        {
            Import(fs);
        }
#else
        var bytes = File.ReadAllBytes(filename);
        using (MemoryStream ms = new MemoryStream(bytes))
        {
            Import(ms);
        }
#endif

FileStreamで読んだときは、40.6821334秒だったのが、byte[]を読んでからのMemoryStreamだと4.4661267秒になった。

頂点座標は、YとZを入れ替えて読み取っているが、これは右手座標系Zアップから、左手座標系Yアップに変換のため。
LAS仕様書には、右手座標系とかZアップは書かれてないけど、たぶんこれであっていると思う。