離散（点群）データからCADデータ（数式）への変換

はじめに

弊社で１５年以上前にボクセル等の離散幾何を製造業向けにモデリング、加工シミュレーションで応用できないかを実証研究を行ったことがあります。その後３Dスキャンデータ（点群）データを設計データ（CADデータ）に変換できないかを具体的に検討する機会があり、経験が役に立った記憶があります。実際複数のメーカー（機器メーカー、3Dプリンタメーカー）からの問い合わせもありこの技術は役立つのでは無いかと思い業務後の時間で色々取り組みました。
 
結論から言うとうまく行かなかったという話になるのですがお問い合わせのあった企業様向けに何とか提案できないかと思い、同様の技術を謳う企業に問い合わせて「こういうこと、できない？」と聞いたところかなりの高額でしたができるという企業がありました。

無償ベンチマークを依頼したところ結果は無残でした。右が私の方法でCADデータ化したやったやつです。

無償ベンチマークを依頼したところ結果は無残でした。右が私の方法でCADデータ化したやったやつです。元形状への忠実性は私の方が良い感じですが設計データ、CADデータとしてみた場合はNGでしょう。左はフィレット部分を頑張った形跡が見られますがやはりNGですが形状特徴をうまく切り分けているようですがもう一歩という感じでしょうか。

やはり元形状への忠実性は私の方が良い感じです。左はフィレット部分を頑張った形跡が見られますが所々隙間が生じています。これだと実用的ではないので聞いてみたところ回答は「想定している前提条件を満たしていない」ということでした。理想的なデータならうまくできるということでしょう。
 
それで依頼元と話した結果どっちも使えるレベルかに至っていないという回答でした。実際、依頼元からは「やはり、、、無理なのかなぁ」と言われました。
 
求めているのは前処理とか専門的な処理を必要とせずに測定データを放り込んだら後は自動的に設計データ（CADデータ）に変換できるシステムだったので結論としては当然な結果でした。設計データ（CADデータ）に変換するということは測定データ（ランダムな点群）から曲面の数式を作るということになります。
 
以下はその当時（確か2017,18年頃）の苦戦の記録です。当然ですが未解決課題がそのまま残っています。

真ん中はメッシュデータで点群から表面になるべき三角形のパッチを生成したメッシュデータです。右側の表示はCADデータに変換した結果をCADで表示している図です。 どうやって実現するか当時私が考えた基本的な攻略案は以下の通りです。

基本的な攻略案

•       点群からMesh生成 陰関数フィッティングを利用
•       MeshをParameter化(UV座標付加)
•       $${Y=(x,y,z)}$$と$${X=(u,v)}$$から曲面Sを計算
$${min(YーS(X))}$$を解いてSを求める
 

RBF(radial basis function)によるMesh再構成

一旦、点群データに対して陰関数をフィットするということを思いつきました。CADで使われる曲面を直接フィットするには点座標(X,Y,Z)とパラメータUV座標が必要ですが陰関数なら点座標(X,Y,Z)だけで先ずはフィットできるというメリットがあります。

陰関数フィットにはRBF(radial basis function)が最適なことは全く別のプロジェクトでお客様と一緒に使い倒した経験もあり最も優秀かつ簡素なアルゴリズムということを経験的に知っていたのでそれを使うことにしました。

こうして$${f(ｘ,ｙ,ｚ)}$$にフィットができると陰関数なので$${f(x,y,z)=0}$$となる$${（ｘ,ｙ,ｚ）}$$が曲面の表面になっているはずです。例えば半径ｒの球面であれば

Jules Bloomenthal's "marching cubes" polygonizer
Bloomenthal（ An Implicit Surface Polygonizer ）によるボクセルからメッシュを生成するアルゴリズムでメッシュ化が可能です。つまり陰関数の表面ボクセルを作ってメッシュ化します。
http://www.unchainedgeometry.com/jbloom/papers.html
Public Domain Polygonizer

かなり忠実に表面を無駄なく覆うメッシュが作成されることが分かります。

ピン角の問題

世の中そんなにあまくないです。メッシュ化した結果として問題が明確になります。アルゴリズムの特性上表面ボクセルの大きさが表面状態をボカシてしまうのです。つまり、本来ピン角だった所がガタガタ、あるいは丸まってしまっています。これは工業製品のようなデータでは完全にNGです。

marching cubesの限界はボクセルの分解能に依存して絶対にピン角情報が欠落するということです。ただ、課題としてf(x,y,z)にフィットした時点でそもそもピン角情報が欠落していたとしたら？
これはmarching cubesの限界だといって放置できません。しかし、幸運なことにコーナー、ピン角、エッジを形状分類してフィッティングするという優れた方法があることが分かりました。MPU法というアルゴリズムです。

Multi-level Partition of Unity

コーナー、ピン角、エッジを形状分類してフィッティング関数を分け（Local Shape Functions）て陰関数に合成する

b2はB-spline関数、 A,B,C,D,E,Fを最小二乗法で計算して決定します。
簡単に説明すると局所的にフィットする陰関数を変えて最後にそれらを合成するアルゴリズムになっています。つまり、MPU法の特徴の一つはそれぞれのセルにフィッティングする曲面の特徴形状を推定してそれぞれに適合するローカル形状陰関数を指定するところにあります。

JU, T., LOSASSO, F., SCHAEFER, S., AND WARREN, J. 2002. Dual contouring of
hermite data. ACM Transactions on Graphics 21, 3 (July), 339–346. Proceedings of ACM SIGGRAPH 2002.
 
OHTAKE, Y., AND BELYAEV, A. G. 2002. Dual/primal mesh optimization for polygonized
implicit surfaces. In 7th ACM Symposium on Solid Modeling and Applications,
171–178.
BLOOMENTHAL, J. 1988. Polygonization of implicit surfaces. Computer Aided
Geometric Design 5, 4, 341–356.
 
KOBBELT, L. P., BOTSCH, M., SCHWANECKE, U., AND SEIDEL, H.-P. 2001.
Feature-sensitive surface extraction from volume data. In Proceedings of SIGGRAPH
2001, ACM Press / ACM SIGGRAPH, Computer Graphics Proceedings,
Annual Conference Series, 57–66.

ローカルフィッティングの様子は図のようになります。

フィットしてない部分（正確に言うと点が乗っていない部分）は０値ではないので陰関数曲面としては自動的に対象外、つまりカットされていくので問題ありません。
 
陰関数のメリットはこういった合成（集合演算）が異常とも思えるくらい楽になります。イメージしやすい例としてまたまた球面に登場してもらいます。

ともかく点の座標のx,y,zそれぞれ2乗して足して半径の二乗を引き算して０だったら球の表面ということです。2個の球面を合成するイメージを考えてみると分かりやすいです。

丁度球の表面が赤線の部分です。そして中（青い部分）は$${f(x,y,z)}$$を計算すると負になっていて外側（薄緑）は正になっています。つまり０に成る部分（赤線）がそれぞれの表面となります。そして2個の球面の足し算が旨く球面の合成（集合演算：和演算）になっていることが分かります。
もう少し正確に言えば球f1と球f2で
　${f(x,y,z) = min( f1(x,y,z), f2(x,y,z)) = 0}$ という陰関数です。

これは球に限らず任意の陰関数表現された曲面に適用できます。これは数値計算上は異常に平易でパラメトリックで表現された形状であればこれとは逆に平易な式では書けないほど難しいものになると思います。

論文ではピン角もちゃんと生成すると言っているが、、、ピン角が出ない？

当然と言えば当然です。$${f(x,y,z)}$$にフィットした時点でそもそもピン角情報が欠落していなかったとしてもmarching cubesの限界の問題は解決していません。つまり、本当にピン角が再生されているのか確認できないわけで結局は可視化においてもメッシュ化も変更する必要があります。そしてまたまた幸運なことにちゃんとそこを考えてくれている人が世の中には居ました。拡張marching cubesです。EMC（Extended Marching Cubes）Dual Contouringというアルゴリズムです。早速実装して試してみました。

左がmarching cubes、右がMPUと拡張marching cubes を使った結果です。一目瞭然で拡張marching cubes はちゃんとエッジ、ピン角が正確に推定できています。素晴らしい！！
 
 
しかし、やはりというかそんなに甘くなったです。

非多様体を生成しやすくなりむしろ問題は深刻になっていました。ただこれは単に私の実装が悪いからかも知れないです。

問題はフィットの時点で悪さしているのかメッシュ化の時点で悪さしてしまったのか、あるいは両方なのか？この時点では良く分かりません。じゃ、どうしたら良いのか？これは最終手段ではありますがレイトレするしかないです。レイトレならメッシュ化を回避して直接f(x,y,z)=0をレンダリングできるのでf(x,y,z)=0に問題があるのかは正確に確認できます。
Raytraceing

レイトレなら関数通りに可視化されるので勝手に角が丸められることは無いですが、レイトレは計算コストが高すぎます。今ならリアルタイムでできるかも（？）
結局レイトレしてみました。エッジ、ピン角、コーナーも正確に推定できていることは間違いありませんし、非多様体を生成していないことも確認できました。

ではメッシュ化の問題はどうするのか？根本的な問題は残っていますがユーザー負担は増えますが実験したところスキャン等の解像度を上げてサンプリングする測定点数を増やすことでエッジ、ピン角、コーナーの問題はある程度は緩和することができます。気持ちよく解決はしていませんが、ここで止まらずに次のステップを確認することにします。

メッシュ化できていても情報は三角形の頂点座標（ｘ、ｙ、ｚ）しかありません。これにＵＶ座標を追加する必要があります。いわゆるテクスチャマッピングです。

なぜこれが必要かと言えば最終的にはＣＡＤの曲面（数式）Ｓを得たいからです。この最適化を計算すればＳが得られるはずです。
実は後になってこんな面倒なことをしなくても良かったのではないか？と思う部分でもあります。
Meshは三角形の集合体ですがこれをうまく四角形の集合体に出来ればもっと楽だったかもしれません。

パラメタライゼーション（Parameterization）

テクスチャマッピングを自動でやるにはパラメタライゼーションという技術を使いますがどういう技術か？と言えば３Ｄ形状に綺麗にラッピングする技術です。３Ｄ形状に正方形の方眼紙を被せて切れ目や穴、折り目や隙間なく綺麗にラッピングします。そうすると３Ｄ形状にぴったりと張り付いた正方形の方眼紙の座標がＵＶ座標になります。これは数学的には平面に位相同型にすることです。トポロジ数学（位相数学）知っている方ならこう言うはずです。絶対に無理じゃん！球体とか。

それとこういうやつ（うさぎ（スタンフォードバニーと言います））、絶対に切れ目とか入れずに平面に展開できないんだから絶対に無理。
であればうまく切れ込みを入れてやれば平面に引き延ばせるんじゃないの？そうしたら綺麗にラッピングできるからそれでやり切れるよね？創意工夫です。
 
ところが話はそう簡単では無いです。球体とかスタンフォードバニーは切れ目を入れれば大丈夫ですがこれは球面と同相だからです。

例えばドーナッツのような形状

この形は切れ込みを２個入れると平面に展開可能です。世の中色んな３D形状があるけど与えられた３D形状に幾つ切れ込みを入れたら平面に展開可能なのか？という問題になります。これは浮き輪の数が幾つかを数え上げる問題と等価なので位相の種数が分かれば幾つ切れ込みを入れたら良いのかが分かります。

位相（トポロジー）

ここで少しMeshの位相に関して幾つか基本的な事実を上げると、
オイラー標数 χ に対して種数を g としたとき、閉曲面　χ = 2 − 2g 、b 個の境界を持つ曲面　
χ = 2 − 2g − b　が成り立っている。
オイラー標数は頂点数 V = q0, 辺の数 E = q1, 面の数 F = q2　を使って
 χ= q0- q1+  q2 という形にも言い変えられます。
 

種数 gを元にMeshを平面と位相同型に切開手術を施して、頂点とUV空間と対応付けを行います。対応付けにはmesh parameterizationを実施します。この時UV空間からの面積が均一になるようにストレッチ計算も行い、ひずみエネルギーを最小化する最適化問題として解くことで可能な限り綺麗にラッピングできるはずです。

Mesh cut

Meshを切開手術
基本的にはMeshに切れ目を入れる。

切れ目を入れることで平面に展開できるようなります。

メッシュを対象に計算でやるとここまで綺麗には切れ目を入れられないですが可能になります。切れ目を計算で求めた結果はこうなりました。

結構汚い感じの切れ目ですがちゃんと平面（４方形）に展開できます。

mesh parameterization

これはひたすら連立方程式を解く問題に帰着します。解法には特異値分解を使用しました。そして境界条件として４方形の境界を(0,0)->(1,0)->(1,1)->(0,1)というUV値で拘束することで正方形の目盛りはかならず０～１になります。この拘束はちょうど切れ目のラインと一致させるということです。

実験結果が以下です。計算過程で数値エラーを起こすことがしばしばあったのでやり方を変えました。連立方程式はBi-CGSTABで計算して微小な誤差がＵＶ空間では致命的になるので座標等全ての浮動小数点データは４倍精度（Quadruple precision floating point number）で計算できるように改良しました。

これでようやくゴールに近づきました。よく見ると怪しい箇所がありますが、、、

Surface fitting

Parameterizationで座標とUV座標がマッチングしたら

を計算します。CADデータにしたいのでSは自由度の高いNURBS曲面を想定しました。

これは数理最適化問題に帰着しますし制約条件が無いのでわりと簡単に解ける問題です。直接解くと計算が終わらないのでパッチ毎に計算することも出来ます。これはUVが確定しているので領域分けが簡単になります。UV空間は２Dの四角形なので

こんな感じに分割して計算できます。
注意が必要なのはそのまま計算するとつなぎ目で連続性が失われてしまいます。実際にはこの先が旨く行くのか分からなかったのでこの時点ではこの計算はやりませんでした。実際、やらなくて良かったと思って居ます。というのも別の問題が後に出てきてここで足止めしてしまい時間を浪費することは避けられました。

曲面のたわみ/しわ　問題

実装して実験した結果はコントロールポイントが近接し過ぎて（縮退）してしまうため「しわ」が発生してしまいました（左）。

数学的は正しいはずですが数値としてコントロールポイントが近接したため曲面の状態が連続性を失ってしまいました。無限精度で計算できるのであればこの問題は無いはずですがそんなのは無理です。そこで創意工夫です。コントロールポイントが近接し過ぎてしまわないようにモデルをスケーリングしたところこの問題は解消しました。（右は100倍にスケーリング）コントロールポイントが近接してしまうのはParameterizationで十分なストレッチができないことで境界付近にパラメータのたわみや引き延ばしが発生し曲面生成でそのまま再現されていました。

完全に解消はできていませんが最終的には

こんな感じです。右はCADデータ(IGES)になっているのでCADに取り込むことが可能になっています。実際、CADで取り込んで表示したときのキャプチャです。

大きな残課題

Meshの曲面化は概ね実証できましたが曲面品質を上げるにはMeshの段階でsegmentationを行い、各セグメントパッチに対して曲面化できれば品質は向上すると思いました。
 
実際、手作業でやってみると品質は向上していることが一目でわかります。

できればこれを自動的にやりたいです。手作業でやると

このような非多様体を作りこんでしまい、Parameterizationでエラーになってしまうことがしばしばあります。そこでCGAL（https://www.cgal.org/）を利用しました。

特徴的な部分をパーツにするようなsegmentationを実施できる点が素晴らしい。セグメンテーションをやってみると品質の差は歴然です。

しわやすじが無くなって綺麗になっていることが分かります。つまり十分なストレッチができていることになります。そして各セグメントをCAD曲面化すると

完全に問題が解消していることが分かります。つまりコントロールポイントが近接して縮退するような事態は回避できていることが確認できます。
 
しかし、同時に大きな問題にも突き当たります。もともとは滑らかな３Dモデルだったものが複数の曲面で構成され曲面間の連続性を失ってしまいました。つまり、複数の曲面の接続境界で連続になるようにコントロールポイントを修正する縫い合わせが必要になってしまいました。
 
 
その他の課題としては工業製品の場合は特に部品を認識しアセンブリ毎のセグメントパッチができれば曲面化と同時にアセンブリ構造もある程度は構築可能になると思います。
 
課題を整理すると
（１）メッシュ化で非多様体を生成しない頑健な方法は？
（２）コントロールポイントが近接して縮退させない方法は？
（３）部品を認識しアセンブリ毎のセグメンメーションする頑健な方法は？
（４）頑健なコントロールポイントを修正する縫い合わせ方法は？
 
おそらくどの課題も世の中で解決済だと思いますが問題なのは「頑健性」です。理想的なデータであれば理論が示す期待通りの結果にはなると思いますが問題は理想的なデータをユーザーが入力するとは限らない点です。そして、

これらは私には無理そうです。なのでこの課題の段階で諦めました。

（１）に関しては手作業で編集するという方法があります。しかし、それ以前にMeshのトポロジ崩壊が問題になります。手作業でのMeshのトポロジ修正は大変な肉体労働です。これはやって見ると分かります。いくつか代表的な例を挙げます。

参考までに実施した検証ケースをお見せします。点群データからCADデータへの変換例です。

これらは全てCADデータ(IGES)になっています。かなり複雑な形状もCADデータ化はできています。ただし、課題点はまだ残っているため実用上の価値は殆どないでしょう。