3Dデータの入出力(ascii)についてまとめてみた

# はじめに
基本的に可視化されたメッシュデータを見たいときは、MeshLabとかCloudCompareを使っていますが、フリーソフトにやりたい機能がない、というときは自分でプログラムを組むことになります。今回はその基本のファイルの入出力について（意外と記事になかったので）まとめました。

## 対象
obj, off, ply, stl
plyとstlはバイナリ形式もありますが、今回はasciiのみ扱います。

## 前提
- 頂点の法線情報はなし
- 面はすべて三角形
- 色・テクスチャの情報はなし

## クラス
前提に沿った頂点と面のクラスを用意します。メッシュはこれらの配列で、頂点は座標情報、面は頂点配列のインデックスを保持する、という一般的なメッシュデータファイルの形式をイメージしています。シンプルすぎてあまり実用的なデータ構造ではないですが、とりあえずプログラム上に存在させることを考えます。
```cpp
class Vertex{
public:
    Vertex(double x, double y, double z): x(x), y(y), z(z) {}
    double x, y, z;
};

class Face{
public:
    Face(unsigned int v1_idx, unsigned int v2_idx, unsigned int v3_idx):
        v1_idx(v1_idx), v2_idx(v2_idx), v3_idx(v3_idx) {}
    unsigned int v1_idx, v2_idx, v3_idx;
};

std::vector<Vertex> vertices;
std::vector<Face> faces;
```

# 入力
## 共通
ファイルの入力には `ifstream` を用います。ファイルが読み込めなかったときの例外処理は用意しておきます。こんな感じ。
```cpp
std::ifstream fin(filename, std::ios::in);
if(fin.fail()){
    std::cerr << "couldn' t read the file." << std::endl;
    exit(EXIT_FAILURE);
}
```
 `std::getline` でファイルの情報を1行ずつ文字列として読み込むことができます。
```cpp
std::string reader;
std::getline(fin, reader);
```
基本的に頂点・面の情報はファイルの1行に1つ、スペース区切りで入ります。そのため文字列をスペースで分割できるものを用意しておきましょう。例えば `sstream` を使えば
```cpp
// reader = "0.1 0.2 0.3";
double x, y, z; 
std::stringstream ss;
ss << reader;
ss >> x >> y >> z;
```
みたいな方法が使えます。もちろん、スペースで区切られた`std::vector<std::string>` を返す関数を自作してもいいと思います（こっちの方が汎用性は高い）。

## obj
[前提](#前提)に沿ったobjファイルのフォーマットはこんな感じ（[参考](https://www.hiramine.com/programming/3dmodelfileformat/objfileformat.html)）。
```example.obj
v 0.1 0.2 0.3
v 0.4 0.5 0.6
...
f 1 2 3
f 4 5 6
...
```
頂点や面の数（≒ファイルの行数）は与えられないので、最後までファイルを読み込む必要があります。こんな感じで書けばEOFで止まってくれます。
```cpp
while(std::getline(fin, reader)){
    // 処理
}
```
あとは一番最初の文字がvかfかで判断すれば良いでしょう。ただし、.objのfaceのindexは1から始まるので注意が必要です。
```cpp
    std::stringstream reader_data;
    reader_data << reader;
    std::string head; double i1, i2, i3;
    reader_data >> head >> i1 >> i2 >> i3;
    
    // reader = "v 0.1 0.2 0.3"
    if(head == "v"){
        vertices.push_back(Vertex(i1, i2, i3));
    }
    
    // reader = "f 1 2 3";
    if(head == "f"){
        faces.push_back(Face(i1 - 1, i2 - 1, i3 - 1));
    }
```

## .off
[前提](#前提)に沿ったoffファイルのフォーマットはこんな感じ（[参考](https://en.wikipedia.org/wiki/OFF_(file_format))）。
```example.off
OFF
4 4 0
1.0 0.0 0.0
0.0 1.0 0.0
...
3 0 1 2
3 1 2 3
...
```
1行目はOFFファイルであることを表しています。NOFFと書かれていることもあり、この場合は頂点の法線ベクトルが横に入って `1.0 0.0 0.0 0.1 0.2 0.3` みたいになります。2行目は左から頂点数・面数・辺数を表していますが、辺の情報が入っているOFFファイルを私は見たことないです。  
先ほどのように `while` を使うと頂点と面の境目の判定がちょっとめんどくさいです。せっかく2行目で頂点数と面数が与えられているので、 `for` 文の方がやりやすそうです。
```cpp
std::getline(fin, reader);
assert(reader == "OFF");
std::getline(fin, reader);
std::stringstream header_data;
header_data << reader;
unsigned int num_verts, num_faces, num_edges;
header_data >> num_verts >> num_faces >> num_edges;

for(int i = 0; i < num_verts; i++){
    // 以下略
}
for(int i = 0; i < num_faces; i++){
    // 以下略
}
```
ちなみに、offファイルでは空行とコメントに気を付ける必要があります。というのも、objファイルは各行の1文字目がvかfでない場合は何もしないので自然にスルーできていましたが、offファイルはそのような識別ができません。コメント行は#で始まるので下記のようなものをそれっぽいところに散りばめて回避しましょう。
```cpp
while(reader.length() == 0 || reader[0] == '#') std::getline(fin, reader);
```

## ply
[前提](#前提)に沿ったplyファイルのフォーマットはこんな感じ（[参考](http://paulbourke.net/dataformats/ply/)）。
```example.ply
ply
format ascii 1.0
element vertex 8
property float x
property float y
property float z
element face 12
property list uchar int vertex_indices
end_header
0.0 0.0 10.0
10.0 0.0 10.0
...
3 0 1 2
3 3 2 1
...
```
offファイルに比べるとヘッダ行が長いですね。2行目はフォーマットがasciiであることを表します。 `element` で始まる行は頂点数や面数を取得できます。 `property` で始まる行は、データの型と名前が載っています。長い行のところ（property list unchar...）は、int型のvertex_indicesがlistとして入っていてそのlistの要素数はuchar型、みたいなことです。 `comment` で始まる行がコメントとして扱われます。そして、 `end_header` と書かれた行でヘッダーが終わりであることを示します（これ以降コメントが書かれることは想定していないと思います）。その後はOFFとほとんど同じですね。  
ヘッダの本質情報は `element` と `end_header` ですので、こんな感じで読み込めばいいでしょう。
```cpp
unsigned int num_verts, num_faces;
while(reader != "end_header"){
    std::getline(fin, reader);
    reader_data << reader;
    std::string i1, i2, i3;
    reader_data >> i1 >> i2 >> i3;
    if(i1 == "element"){
        if(i2 == "vertex") num_verts = std::stoul(i3);
        if(i2 == "face") num_faces = std::stoul(i3);
    }
}
for(int i = 0; i < num_verts; i++){
    // 以下略
}
for(int i = 0; i < num_faces; i++){
    // 以下略
}
```
<details><summary>propertyの型読み取り</summary><div>

基本的に頂点座標を `double` ・インデックスを `uint32` にしておけば、大は小を兼ねると思いますが、ヘッダに記述されたpropertyの型情報を厳密に守りたい場合は、テンプレートを用いるといいでしょう。
```cpp
template <typename T>
class Vertex{
public:
    Vertex(T x, T y, T z): x(x), y(y), z(z) {}
    T x, y, z;
};

template <typename T>
class Face{
public:
    Face(T v1_idx, T v2_idx, T v3_idx): 
        v1_idx(v1_idx), v2_idx(v2_idx), v3_idx(v3_idx) {}
    T v1_idx, v2_idx, v3_idx;
};

// ヘッダを読み込む。例えば頂点座標がfloat、頂点インデックスがushortだった
std::string vert_type, vert_idx_type;
while(reader != "end_header"){
    std::getline(fin, reader);
    reader_data << reader;
    std::string i1, i2, i3, i4, i5;
    reader_data >> i1 >> i2 >> i3 >> i4 >> i5;
    if(i1 == "property"){
        // ハードコーディングがひどいので、うまいこと修正してください
        if(i3 == "x") vert_type = i2;
        if(i5 == "vertex_indices") vert_idx_type = i4;
    }
}
if(vert_type == "float"){
    std::vector<Vertex<float>> vertices;
    if(vert_idx_type == "ushort"){
        std::vector<Face<unsigned short>> faces;
        // 以下略
    }
}
```
</div></details>

## stl
[前提](#前提)に沿ったstlファイルのフォーマットはこんな感じ（[参考](https://www.hiramine.com/programming/3dmodelfileformat/stlfileformat.html)）。
```example.stl
solid example
    facet normal  0.0  0.0  1.0
        outer loop
            vertex  0.0  0.0  10.0
            vertex  10.0  0.0  10.0
            vertex  0.0  10.0  10.0
        endloop
    endfacet
    facet normal  0.0  0.0  1.0
        outer loop
            vertex  10.0  10.0  10.0
            vertex  0.0  10.0  10.0
            vertex  10.0  0.0 10.0
        endloop
    endfacet
...
endsolid
```
STLは3Dプリンタ用としてもよく使われるフォーマットなのですが、今までの3つと比べるとかなり毛色が違うことがわかります。面の数だけ、その面の法線と構成する各頂点の座標が格納されています。このフォーマットは以下の点で冗長です。
- 面の法線が格納されている
    - 3点の座標 $\boldsymbol{v_1}, \boldsymbol{v_2}, \boldsymbol{v_3}$ があれば、$(\boldsymbol{v_2} - \boldsymbol{v_1}) \times (\boldsymbol{v_3} - \boldsymbol{v_1})$ を正規化して法線ベクトルを出せるので不要です。 方向は$\boldsymbol{v_1} \rightarrow \boldsymbol{v_2} \rightarrow \boldsymbol{v_3}$ まわりに右ネジを回すと進む方向が主流です。
- 頂点の座標情報が何回も出てくる
    - 1つの頂点が1つの面にしか属しないということはまずないでしょう。

とはいえ、ここではあくまで[共通](#共通)で作ったデータ構造に落とし込むことを考えましょう。 `vertices` は `Vertex` のvectorなので、単純な方法としては、線形探索して同じものが見つかったら（等価演算子をオーバーロードしておきましょう）そのindexを取得、なければverticesに挿入してそのサイズを取得してindexとする、という方法があります。
```cpp
// Vertexクラス内に追記
bool operator==(const Vertex &other) const {
    return x == other.x && y == other.y && z == other.z;
}

while(std::getline(fin, reader)){
    // 略
    ss >> head >> x >> y >> z;
    if(head == "vertex"){
        Vertex v1 = Vertex(x, y, z);
        std::getline(fin, reader);  // 2つ目の頂点
        // 略
        Vertex v2 = Vertex(v2x, v2y, v2z);
        std::getline(fin, reader);  // 3つ目の頂点
        // 略
        Vertex v3 = Vertex(v3x, v3y, v3z);

        bool v1_new = true; bool v2_new = true; bool v3_new = true;
        unsigned int v1_idx, v2_idx, v3_idx;
        for(int i = 0; i < vertices.size(); i++){
            if(vertices[i] == v1){
                v1_idx = i;
                v1_new = false;
            }
            // v2, v3も同様
        }
        if(v1_new){
            v1_idx = vertices.size();
            vertices.push_back(v1);
        }
        // v2, v3も同様
        faces.push_back(Face(v1_idx, v2_idx, v3_idx));
    }
}
```
この方法は、1回の挿入を $\mathcal{O}(1)$ でできますが、その前に重複がないかを頂点数 $N$ に対し線形検索$\mathcal{O}(N)$するので、面数$M$ でこの処理をすることを考えると全体の計算量は $\mathcal{O}(MN)$ です。重複を防ぐ挿入に $\mathcal{O}(\log N)$ ででき平衡二分木（set）を使い、vectorに直して二分探索（$\mathcal{O}(\log N)$）する方が高速です。

自作クラス型のsetを使うためには比較演算子をオーバーロードする必要があります。ちなみに $v_1<v_2$ と $v_2<v_1$ がどっちもfalseを返してきたときに、$v_1$ と $v_2$ が同じ値だという判定をするので、等価演算子のオーバーロードはなくていいらしいです。
```cpp
// Vertexクラス内に追記
bool operator<(const Vertex &other) const {
    return x < other.x || (x == other.x && y < other.y) || (x == other.x && y == other.y && z < other.z);
}
```
setは挿入をするたびにindexが変わるので、ファイルを2回走査して「頂点を重複がないように `vertices` に挿入」→「各面の各頂点のインデックスを探し `faces` に挿入」 します。1回で済む上手い方法があったら教えてください。
```cpp
std::set<Vertex> set_vertices;
// 挿入
while(std::getline(fin, reader)){
    std::stringstream ss;
    ss << reader;
    std::string head; double v1x, v1y, v1z;
    ss >> head >> v1x >> v1y >> v1z;
    if(head == "vertex"){
        Vertex v1 = Vertex(v1x, v1y, v1z);
        
        std::getline(fin, reader);  // 2つ目の頂点
        // 略
        Vertex v2 = Vertex(v2x, v2y, v2z);
        
        std::getline(fin, reader);  // 3つ目の頂点
        // 略
        Vertex v3 = Vertex(v3x, v3y, v3z);

        set_vertices.insert(v1);
        set_vertices.insert(v2);
        set_vertices.insert(v3);
    }
}

std::vector<Vertex> vertices_(set_vertices.begin(), set_vertices.end());
vertices = vertices_;
std::sort(vertices.begin(), vertices.end());

// 探索
while(std::getline(fin2, reader)){
        // 略
        ss >> head >> v1x >> v1y >> v1z;
        if(head == "vertex"){
            Vertex v1 = Vertex(v1x, v1y, v1z);
            std::vector<Vertex>::iterator v1_itr = std::lower_bound(vertices.begin(), vertices.end(), v1);
            unsigned int v1_id = std::distance(vertices.begin(), v1_itr);

            std::getline(fin2, reader);  // 2つ目の頂点
            // 略
            Vertex v2 = Vertex(v2x, v2y, v2z);
            std::vector<Vertex>::iterator v2_itr = std::lower_bound(vertices.begin(), vertices.end(), v2);
            unsigned int v2_id = std::distance(vertices.begin(), v2_itr);

            std::getline(fin2, reader);  // 3つ目の頂点
            // 略
            Vertex v3 = Vertex(v3x, v3y, v3z);
            std::vector<Vertex>::iterator v3_itr = std::lower_bound(vertices.begin(), vertices.end(), v3);
            unsigned int v3_id = std::distance(vertices.begin(), v3_itr);

            faces.push_back(Face(v1_id, v2_id, v3_id));

            printf("\r%07d", i);
            fflush(stdout);
            i += 1;
        }
    }
```
この方法（$ 2 \times \mathcal{O}(M \log N)$）で $N \fallingdotseq 200,000$ 、$M \fallingdotseq 400,000$ のメッシュを読み込むのに10秒弱かかりました（MeshLabでのインポートは3秒ほどでした）。objなどは$\mathcal{O}(M+N)$ で1秒弱なので結果としては妥当でしょうか。もう少し高速化できそうな気がしますね。

<details><summary>ちなみに</summary><div>
はじめは以下のようにsetからインデックスを取得することを考えたのですが、すごく時間がかかっていました。

```cpp
// 探索
while(std::getline(fin2, reader)){
    // 略
    ss >> head >> v1x >> v1y >> v1z;
    if(head == "vertex"){
        Vertex v1 = Vertex(v1x, v1y, v1z);
        std::getline(fin, reader);  // 2つ目の頂点
        // 略
        Vertex v2 = Vertex(v2x, v2y, v2z);
        std::getline(fin, reader);  // 3つ目の頂点
        // 略
        Vertex v3 = Vertex(v3x, v3y, v3z);

        std::set<Vertex>::iterator v1_itr = set_vertices.find(v1);
        unsigned int v1_idx = std::distance(set_vertices.begin(), v1_itr);
        // v2, v3も同様
        faces.push_back(Face(v1_idx, v2_idx, v3_idx));
    }
}
```
[こちら](https://cpprefjp.github.io/reference/iterator/distance.html)を見たところ、setのイテレータ間の距離を求めるには $\mathcal{O}(N)$ かかるらしく、それが原因でした。
</div></details>

# 出力
結構長くなってしまったのと、フォーマットが分かっていればやることはほとんど変わらないので、offだけにします（笑）。

## 共通
ファイルの出力には `ofstream` を用います。ファイルが読み込めなかったときの例外処理は用意しておきます。その後は標準出力とそっくりな書き方で書き込むことができます。
```cpp
std::ofstream fout(filename, std::ios::out);
if(fout.fail()){
    std::cerr << "couldn' t open the file." << std::endl;
    exit(EXIT_FAILURE);
}
```
## off
ヘッダ行を作ります。
```cpp
fout << "OFF" << std::endl;

int num_verts = vertices.size();
int num_triangles = triangles.size();
int num_edges = 0;
fout << num_verts << " " << num_triangles << " " << num_edges << std::endl;
```
頂点数、面数だけファイルに書き込みます。単純ですね。
```cpp
for(int i = 0; i < num_verts; i++){
    fout << vertices[i].x << " " << vertices[i].y << " " << vertices[i].z << std::endl;
}
for(int i = 0; i < num_triangles; i++){
    fout << 3 << " " << triangles[i].v1_idx << " " << triangles[i].v2_idx << " " << triangles[i].v3_idx << " " << std::endl;
}
```

# おわりに
## 使ったSTL
```cpp
#include <algorithm>
#include <assert.h>
#include <iostream>
#include <fstream>
#include <set>
#include <sstream>
#include <string>
#include <time.h>
#include <utility>
#include <vector>
```

## 参考にしたサイト
- [リファレンス - cpprefjp C++日本語リファレンス](https://cpprefjp.github.io/reference)
- [hiramine.com - 3Dモデルファイルフォーマット](https://www.hiramine.com/programming/3dmodelfileformat)
- [OFF (file format)) - Wikipedia](https://en.wikipedia.org/wiki/OFF_(file_format))
- [PLY - Polygon File Format](http://paulbourke.net/dataformats/ply/)

以上です。ここまで読んでいただきありがとうございました。