A Gaussian Filter for Plate Flatness Evaluation System with 3-D Scanner

Synopsis:

LIDAR (light detection and ranging) system was applied to a plate flatness evaluation system. Plate flatness surfaces are reconstructed from many points generated by LIDAR with a smoothing spline method. We defined a smoothing spline functional with sampling measure weights. The equivalent number of parameters defined on this functional does not depend on the distributions of samples. The approximation of the equivalent number of parameters is derived when the number of samples becomes infinity. This approximation greatly reduced the calculation time needed to estimate the optimal smoothing. The smoothing spline calculation cost was so high that new algorithms (FMM: fast multi-pole method) were introduced and we developed the smoothing engine, which was applied to practical problems. The engine generated clear surfaces and was robust to various dirty points cloud.

1. 緒言

近年，3次元レーザースキャナー(3Dスキャナー)と呼ばれる測定装置が汎用的，低コストとなり，3Dスキャナーを活用した鋼板の曲り/平坦度を測定する方法が造船工場等で適用されつつある。鉄鋼業で3Dスキャナーを活用する場合，オンライン平坦度/形状計のように移動している鋼板を測定することはできないが，一方で，ポータブル性に優れ，オンサイトでの測定，および定点固定された鋼板の形状測定には，適切な測定装置と考えられる。

3Dスキャナーを鋼板の平坦度，形状測定に適用する問題点は，3Dスキャナーから出力される数百万点の大規模点群データの処理である。鋼板表面に射影された点群データはグリッド状とならず，粗密のある規則的でない分布となり，曲面化処理が複雑となる。個々の点データには，測定誤差が含まれ，精度を高めるためには不規則分布点群に対する平滑化処理が必要となる。このような問題に対して，平滑化スプライン法を適用することを考える。

平滑化スプライン法は，誤差を含むサンプリング値からノンパラメトリックな曲線あるいは曲面を推定するための回帰法として知られており，様々な理工学分野で用いられている。平滑化スプラインは理論背景が明確な位相補償ガウシアンフィルターであり，ノイズ除去のための信号処理^9,13)，画像処理での画像復元¹⁸⁾およびノイズ除去¹⁰⁾，磁場，重力場に対する逆問題⁶⁾，医療データの統計処理⁸⁾および3Dスキャナー等を用いて測定されたデータからの形状曲面生成^12,14)で応用されている。

平滑化スプライン法は，平滑化パラメータを制御することで，任意の平滑化された回帰関数を計算することが可能である。平滑化パラメータを小さくしすぎると，回帰関数が測定データとフィッティングしすぎ(過適合)となり，回帰関数がギザギザとなる。逆に平滑化パラメータを大きくしすぎると，平滑化されすぎてしまい，回帰関数から重要な情報が得られなくなる。このことから，測定対象形状，データに依存して，ほどよく適合した回帰関数，また適切な平滑化パラメータが存在することがわかる。適切な平滑化パラメータを求める手法として，一般化交差検証法^2,3)(GCV)が知られているが，フル正方逆行列演算が必要となり，点群数が多くなると計算時間が膨大となり，適切な平滑化パラメータを求めることが実質不可能となる。

曲面生成のための平滑化スプラインの計算には，非常に計算コストがかかり，大規模点群処理には不適であったが，近年，計算コストの壁を突破できる可能性を持つ新たなアルゴリズムが開発されている。Beatson and Newsa⁴⁾は，多重極展開法(Fast Multi-pole Method, FMM)を用いた高速回帰モデル評価手法を提案し，劇的に計算時間を低減させることが可能となった。Beatsonら⁵⁾は，近似的なCardinal Spline(局所台の性質を持つスプライン)を用いた係数求解のための前処理(Preconditioning)と，GMRES(Generalized Minimal Residual Method)を用いてスプライン係数を高速に計算する手法を提案した。

本論文で述べる曲面生成方法には，過去に開発された平滑化スプライン法に対して独自の改良も加えた。実測された点群データに対応するために，サンプリング測度重み付き平滑化スプライン汎関数を定義する^16,17)。この汎関数を用いることで，大規模点群に対する適切な平滑化パラメータを高速かつ実用的レベルで求めることができる。3Dスキャナーで得られた大規模な実データの平滑化処理を行い，本論文で提案した手法が実用的に有効であることを示す。

2. 平滑化スプライン回帰

平滑化3次スプライン汎関数は次式のように定義される。   

$$\Pi =\sum _{i=1}^m\frac{1}{2}{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2+\gamma {\int }_0^L\frac{1}{2}{\left(\frac{d^{2}f}{d{x}^2}\right)}^{2}dx$$(1)

ここで，mはサンプリング数，Lはサンプリング領域の長さ，x⁽ⁱ⁾_M(i=1,…,m)は0≤x^(l)_M<…<x^(m)_M≤Lのようになるi番目のサンプリング点の位置，z⁽ⁱ⁾_M(i=1,…,m)はi番目のサンプリング点の値，fは回帰曲線，そしてγは平滑化パラメータである。式(1)の汎関数に変分原理を適用すると，   

$$\begin{array}{c}\delta \Pi ={\int }_0^L\left\{\gamma \frac{d^{4}f}{d{x}^4}+\sum _{i=1}^m\left(f-z_M^{\left(i\right)}\right)\text{Dirac}\left(x-x_M^{\left(i\right)}\right)\right\}\delta fdx\\ +\gamma {\frac{d^{2}f}{d{x}^2}\delta \left(\frac{df}{dx}\right)|}_{x=0}^{x=L}-\gamma {\frac{d^{3}f}{d{x}^3}\delta f|}_{x=0}^{x=L}=0\end{array}$$(2)

ここで，Dirac(x)はディラックのデルタ関数である。式(2)から，平滑化3次スプラインの微分方程式   

$$\gamma \frac{d^{4}f}{d{x}^4}+\sum _{i=1}^m\left(f-z_M^{\left(i\right)}\right)\text{Dirac}\left(x-x_M^{\left(i\right)}\right)=0$$(3)

と境界条件   

$$\frac{d^{2}f\left(0\right)}{d{x}^2}=\frac{d^{2}f\left(L\right)}{d{x}^2}=\frac{d^{3}f\left(0\right)}{d{x}^3}=\frac{d^{3}f\left(L\right)}{d{x}^3}=0$$(4)

が得られる。式(3)の微分方程式を，式(4)の境界条件のもとで解くと，回帰曲線は，   

$$f\left(x\right)=c_1+c_{2}x+\sum _{i=1}^m\frac{d_M^{\left(i\right)}}{12}{\left|x-x_M^{\left(i\right)}\right|}^3$$(5)

となり，また回帰曲線の拘束条件は，   

$$\sum _{i=1}^m{d}_M^{\left(i\right)}=\sum _{i=1}^m{x}_M^{\left(i\right)}{d}_M^{\left(i\right)}=0$$(6)

となる。ここで，c₁，c₂およびd⁽ⁱ⁾_M(i=1,…,m)は未知のパラメータである。

式(5)を式(3)へ代入すると，   

$$f\left(x_M^{\left(i\right)}\right)+\gamma d_M^{\left(i\right)}=z_M^{\left(i\right)}\left(i=1,\cdots ,m\right)$$(7)

となる。式(7)と式(6)を用いると，次式の連立一次方程式が得られる。   

$$\left[\begin{array}{cc}\left[A_M\right]+\gamma \left[I_{m\times m}\right]& \left[P_M\right]\\ {}^T\left[P_M\right]& \left[O_{2\times 2}\right]\end{array}\right]\left\{\begin{array}{c}\left\{d_M\right\}\\ \left\{c\right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{z_M\right\}\\ \left\{O_2\right\}\end{array}\right\}$$(8)

ここで，[I_m×m]はm行m列の単位行列，[O_2×2]は2行2列の零行列，{O₂}は2行1列の零ベクトルであり，そして，   

$$\begin{array}{c}{}^T\left[A_M\right]=\left[\begin{array}{ccc}{}^T\left[a_M\left(x_M^{\left(1\right)}\right)\right]& \cdots & {}^T\left[a_M\left(x_M^{\left(m\right)}\right)\right]\end{array}\right]\\ \left[a_M\left(x\right)\right]=\left[\begin{array}{ccc}\frac{1}{12}{\left|x-x_M^{\left(1\right)}\right|}^3& \cdots & \frac{1}{12}{\left|x-x_M^{\left(m\right)}\right|}^3\end{array}\right]\end{array}$$(9)

  

$$\begin{array}{c}{}^T\left[P_M\right]=\left[\begin{array}{ccc}{}^T\left[p\left(x_M^{\left(1\right)}\right)\right]& \cdots & {}^T\left[p\left(x_M^{\left(m\right)}\right)\right]\end{array}\right]\\ \left[p\left(x\right)\right]=\left[\begin{array}{cc}1& x\end{array}\right]\end{array}$$(10)

  

$${}^T\left\{d_M\right\}=\left[\begin{array}{ccc}{d}_M^{\left(1\right)}& \cdots & d_M^{\left(m\right)}\end{array}\right]$$(11)

  

$$\begin{array}{c}{}^T\left\{c\right\}=\left[\begin{array}{cc}{c}_1& c_2\end{array}\right]\\ {}^T\left\{z_M\right\}=\left[\begin{array}{ccc}{z}_M^{\left(1\right)}& \cdots & z_M^{\left(m\right)}\end{array}\right]\end{array}$$(12)

である。式(8)を解くことによって，未知のパラメータ{c}および{d_M}を求めることができる。式(9)，(10)，(11)および(12)を式(5)へ代入すると，マトリクス表現による回帰曲線は次式となる。   

$$f\left(x\right)=\left[a_M\left(x\right)\right]\left\{d_M\right\}+\left[p\left(x\right)\right]\left\{c\right\}$$(13)

未知のパラメータを式(13)に代入することにより，回帰曲線は求められる。

3. 一般化交差検証法と情報量規準

一般化交差検証法(GCV)は，最適な平滑化パラメータを推定するための最も有名な方法のひとつである。GCVでは，評価関数   

$$V\left(\gamma \right)=\sum _{i=1}^m{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2/{\left(1-\frac{1}{m}\text{trace}\left(\left[H\right]\right)\right)}^2$$(14)

を極小化する平滑化パラメータが，最適な平滑化パラメータとする。[H]はハット行列と呼ばれ，次式で定義される。   

$$\left\{f\left(\left\{x_M\right\}\right)\right\}=\left[H\right]\left\{z_M\right\}$$(15)

ここで，   

$${}^T\left\{f\left(\left\{x_M\right\}\right)\right\}=\left[\begin{array}{ccc}f\left(x_M^{\left(1\right)}\right)& \cdots & f\left(x_M^{\left(m\right)}\right)\end{array}\right]$$

式(8)から求められる未知のパラメータ{c}と{d_M}を代入して得られる式(13)の回帰曲線を，式(15)の定義式に代入すると，平滑化3次スプラインのハット行列は次式となる。   

$$\left[H\right]=\left[\begin{array}{cc}\left[A_M\right]& \left[P_M\right]\end{array}\right]{\left[\begin{array}{cc}\left[A_M\right]+\gamma \left[I_{m\times m}\right]& \left[P_M\right]\\ {}^T\left[P_M\right]& \left[O_{2\times 2}\right]\end{array}\right]}^{-1}\left[\begin{array}{c}\left[I_{m\times m}\right]\\ \left[O_{2\times m}\right]\end{array}\right]$$

式(14)のGCVの評価関数は，情報量規準と関連があることが知られている。特に，情報量規準の中で最も有名で単純な赤池情報量規準¹⁾   

$$AIC=-2LL+2k$$(16)

は，サンプリング数が大きければ，GCVの評価関数と等価である。LLは最大対数尤度であり，次式で定義される。   

$$LL=-\frac{m}{2}\log \left(\sum _{i=1}^m{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2\right)$$

また，kは回帰モデルのパラメータ数あるいは自由度である。式(14)のGCVの評価関数は，単調増加関数を用いて情報量規準   

$$\begin{array}{c}I{C}_{GCV}=m\log \left(V\left(\gamma \right)\right)=m\log \left(\sum _{i=1}^m{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2\right)-m\log \left({\left(1-\frac{1}{m}\text{trace}\left(\left[H\right]\right)\right)}^2\right)\\ =-2LL-m\log \left({\left(1-\frac{1}{m}\text{trace}\left(\left[H\right]\right)\right)}^2\right)\end{array}$$(17)

に変換される。もしサンプリング数が無限大であれば，式(17)の情報量規準は次式となる。   

$$AIC=-2LL+2\text{trace}\left(\left[H\right]\right)$$(18)

式(16)と式(18)を比較すると，trace([H])は，平滑化スプラインに対するパラメータ数とみなすことができる。等価パラメータ数(ENOP)は，次式で定義される。   

$$k_{GCV}=\text{trace}\left(\left[H\right]\right)$$(19)

サンプリング数が大きくなると，ENOPの計算に時間がかかる。なぜなら，ENOPは 行 例の逆行列を含んでいるからである。平滑化パラメータが変わるごとに，逆行列を計算しなければいけない。De Nicolaoら⁷⁾は，ENOPの高速計算手法を開発した。この計算法を曲線モデルに適用することは可能であるが，曲面モデル，あるいは高次モデルへの適用は言及されていない。ENOPの近似式はサンプリング長さのパラメータを含んでおり，ゆえにサンプリング長さに影響される。ENOPは回帰モデルの複雑性のようなものを示唆しているので，ENOPはサンプリング点の分布に依存してはいけない。サンプリング点の分布に依存しない平滑化スプラインを定義しなければいけない。

4.サンプリング測度重み付き平滑化スプライン回帰

重み付き平滑化l次スプライン汎関数は，次式で定義される。   

$$\Pi =\sum _{i=1}^m\frac{1}{2}{\omega }_M^{\left(i\right)}{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2+\gamma {\int }_0^L\frac{1}{2}{\left(\frac{d^{\left(l+1\right)/2}f}{d{x}^{\left(l+1\right)/2}}\right)}^{2}dx$$(20)

ここで，ω⁽ⁱ⁾_M(i=1,…,m)はi番目のサンプリング点の重みであり，lはl=1，3，5，…となるようなスプライン次数である。重みが次式で定義されると考える。   

$$\begin{array}{c}{\int }_0^{L}g\left(x\right)dx\approx \sum _{i=1}^m{\omega }_M^{\left(i\right)}g\left(x_M^{\left(i\right)}\right)\\ \sum _{i=1}^m{\omega }_M^{\left(i\right)}=L,{\omega }_M^{\left(i\right)}\ge 0\left(i=1,\cdots ,m\right)\end{array}$$(21)

ここで，g(x)は任意の関数である。式(21)は数値積分公式である。もしg(x)が折れ線内挿関数ならば，重みは次式となる。   

$${\omega }_M^{\left(i\right)}=\frac{x_M^{\left(i+1\right)}-x_M^{\left(i-1\right)}}{2}\left(i=1,\cdots ,m\right)$$

ここで，x⁽⁰⁾_M=x⁽¹⁾_Mおよびx_M^(m+1)=x^(m)_Mとする。

式(20)の汎関数に変分原理を適用すると，   

$$\begin{array}{c}\delta \Pi ={\int }_0^L\left\{{\left(-1\right)}^{\left(l+1\right)/2}\gamma \frac{d^{l+1}f}{d{x}^{l+1}}+\sum _{i=1}^m{\omega }_M^{\left(i\right)}\left(f-z_M^{\left(i\right)}\right)\text{Dirac}\left(x-x_M^{\left(i\right)}\right)\right\}\delta fdx\\ +\sum _{j=1}^{\frac{l+1}{2}}{\left(-1\right)}^{j-1}\gamma {\frac{d^{\left(l+1\right)/2+j-1}f}{d{x}^{\left(l+1\right)/2+j-1}}\delta \left(\frac{d^{\left(l+1\right)/2-j}f}{d{x}^{\left(l+1\right)/2-j}}\right)|}_{x=0}^{x=L}=0\end{array}$$(22)

式(22)から，重み付き平滑化l次スプラインの微分方程式   

$${\left(-1\right)}^{\left(l+1\right)/2}\gamma \frac{d^{l+1}f}{d{x}^{l+1}}+\sum _{i=1}^m{\omega }_M^{\left(i\right)}\left(f-z_M^{\left(i\right)}\right)\text{Dirac}\left(x-x_M^{\left(i\right)}\right)=0$$(23)

と境界条件   

$$\frac{d^{\left(l+1\right)/2+j-1}f\left(0\right)}{d{x}^{\left(l+1\right)/2+j-1}}=\frac{d^{\left(l+1\right)/2+j-1}f\left(L\right)}{d{x}^{\left(l+1\right)/2+j-1}}=0\left(j=1,\cdots ,\frac{l+1}{2}\right)$$(24)

が得られる。式(23)の微分方程式を境界条件(24)のもとで解くと，回帰曲線は，   

$$f\left(x\right)=\sum _{j=1}^{\left(l+1\right)/2}{c}_j{x}^{j-1}+\sum _{i=1}^m\frac{d_M^{\left(i\right)}}{2\times \left(l\right)!}{\left|x-x_M^{\left(i\right)}\right|}^l$$(25)

となり，また拘束条件は，   

$$\sum _{i=1}^m{\left(x_M^{\left(i\right)}\right)}^{j-1}{d}_M^{\left(i\right)}=0\left(j=1,\cdots ,\frac{l+1}{2}\right)$$(26)

となる。

式(25)を式(23)へ代入すると，   

$$f\left(x_M^{\left(i\right)}\right)+{\left(-1\right)}^{\left(l+1\right)/2}\frac{\gamma }{{\omega }_M^{\left(i\right)}}{d}_M^{\left(i\right)}=z_M^{\left(i\right)}\left(i=1,\cdots ,m\right)$$(27)

式(27)と式(26)を用いると，次式の連立一次方程式が得られる。   

$$\left[\begin{array}{cc}\left[A_M\right]+{\left(-1\right)}^{\left(l+1\right)/2}\gamma {\left[{\Omega }_M\right]}^{-1}& \left[P_M\right]\\ {}^T\left[P_M\right]& \left[O_{\left(l+1\right)/2\times \left(l+1\right)/2}\right]\end{array}\right]\left\{\begin{array}{c}\left\{d_M\right\}\\ \left\{c\right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{z_M\right\}\\ \left\{O_{\left(l+1\right)/2}\right\}\end{array}\right\}$$(28)

ここで，   

$$\begin{array}{c}\left[{\Omega }_M\right]=\left[\begin{array}{ccc}{\omega }_M^{\left(1\right)}& & 0\\ & \ddots & \\ 0& & {\omega }_M^{\left(m\right)}\end{array}\right]\\ \left[a_M\left(x\right)\right]=\left[\begin{array}{ccc}\frac{1}{2\times \left(l\right)!}{\left|x-x_M^{\left(1\right)}\right|}^l& \cdots & \frac{1}{2\times \left(l\right)!}{\left|x-x_M^{\left(m\right)}\right|}^l\end{array}\right]\end{array}$$(29)

  

$$\left[p\left(x\right)\right]=\left[\begin{array}{cccc}1& x& \cdots & x^{\left(l+1\right)/2-1}\end{array}\right]$$(30)

  

$${}^T\left\{c\right\}=\left[\begin{array}{ccc}{c}_1& \cdots & c_{\left(l+1\right)/2}\end{array}\right]$$(31)

式(28)を解くことで，未知のパラメータ{c}と{d_M}を求めることができる。式(29)，(30)，(11) および(31)を式(5)へ代入すると，マトリクス表現の回帰曲線は式(13)となる。未知のパラメータを式(13)に代入することで，回帰曲線が求められる。

重み付き平滑化スプラインに対するGCVは，評価関数   

$$V\left(\gamma \right)=\sum _{i=1}^m{\omega }_M^{\left(i\right)}{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2/{\left(1-\frac{k_{GCV}}{m}\right)}^2$$(32)

を極小化する平滑化パラメータが最適な平滑化パラメータであると定義される。重み付き平滑化スプラインのハット行列は，次式となる。   

$$\left[H\right]=\left[\begin{array}{cc}\left[A_M\right]& \left[P_M\right]\end{array}\right]{\left[\begin{array}{cc}\left[A_M\right]+{\left(-1\right)}^{\left(l+1\right)/2}\gamma {\left[{\Omega }_M\right]}^{-1}& \left[P_M\right]\\ {}^T\left[P_M\right]& \left[O_{\left(l+1\right)/2\times \left(l+1\right)/2}\right]\end{array}\right]}^{-1}\left[\begin{array}{c}\left[I_{m\times m}\right]\\ \left[O_{\left(l+1\right)/2\times m}\right]\end{array}\right]$$(33)

5. 等価パラメータ数の近似式

もしサンプリング数が無限大ならば，式(23)のサンプリング測度重み付き平滑化スプラインの微分方程式は次式となる。   

$${\left(-1\right)}^{\left(l+1\right)/2}\gamma \frac{d^{l+1}f}{d{x}^{l+1}}+f-z_M=0$$(34)

ここで，z_Mはサンプリング値の関数である。ラプラス変換を式(34)に適用すると，   

$$f\left(s\right)=h\left(s\right)z_M\left(s\right)$$(35)

ここで，sはラプラス変数であり，h(s)は次式で与えられる伝達関数である。   

$$h\left(s\right)=\frac{1}{{\left(-1\right)}^{\left(l+1\right)/2}\gamma s^{l+1}+1}$$(36)

式(35)と式(15)のハット行列の定義式を比較すると，伝達関数 はハット行列と同じ機能をしていることがわかる。フーリエ変換を式(36)の伝達関数に適用すると，次式の周波数応答関数が得られる。   

$$h\left(\omega \right)=\frac{1}{{\left(-1\right)}^{\left(l+1\right)/2}\gamma {\left(j\omega \right)}^{l+1}+1}=\frac{1}{\gamma {\omega }^{l+1}+1}$$(37)

ここで，jは虚数単位であり，ωは角周波数である。周波数応答関数h(ω)は，平滑化スプラインが位相補償ガウシアンフィルターとして機能していることを示している。

三角関数に対する自由度kを   

$$k=L/L_{\omega }$$(38)

として定義する。ここで，L_ωは次式で与えられる半波長である。   

$$L_{\omega }=\pi /\omega$$(39)

式(39)を式(38)へ代入すると，   

$$k=L\omega /\pi$$(40)

式(40)を式(37)へ代入すると，自由度に対する応答関数h(k)が次式で与えられる。   

$$h\left(k\right)={\left(\gamma {\left(\frac{\pi }{L}\right)}^{l+1}{k}^{l+1}+1\right)}^{-1}$$

ENOPがハット行列の固有値の総和であることを意味している式(19)のENOPの定義式に従うと，無限サンプリング数に対するENOPは，   

$$k_A\left(\gamma \right)={\int }_0^{\infty }h\left(k\right)dk={\int }_0^{\infty }{\left(\gamma {\left(\frac{\pi }{L}\right)}^{l+1}{k}^{l+1}+1\right)}^{-1}dk$$(41)

として定義される。ここで，k_AはENOPの近似値となる。変数変換   

$$\kappa =\frac{\pi }{L}{\gamma }^{\frac{1}{l+1}}k$$

を式(41)の積分に適用すると，ENOPの近似式k_Aは次式となる。   

$$k_A\left(\gamma \right)=\frac{L}{\pi }{\gamma }^{-\frac{1}{l+1}}{\int }_0^{\infty }\frac{1}{{\kappa }^{l+1}+1}d\kappa$$(42)

さらに，変数変換   

$$1-\xi =\frac{1}{{\kappa }^{l+1}+1}$$

を式(42)の積分に適用すると，ENOPの近似式k_Aは次式となる。   

$$\begin{array}{c}{k}_A\left(\gamma \right)=\frac{L}{\pi }{\gamma }^{-\frac{1}{l+1}}\times \frac{1}{l+1}{\int }_0^1{\xi }^{\frac{1}{l+1}-1}{\left(1-\xi \right)}^{\frac{l}{l+1}-1}d\xi \\ =\frac{L}{\pi }{\gamma }^{-\frac{1}{l+1}}\times \frac{1}{l+1}{\rm B}\left(\frac{1}{l+1},\frac{l}{l+1}\right)=\frac{L}{\pi \left(l+1\right)}{\rm B}\left(\frac{1}{l+1},\frac{l}{l+1}\right){\gamma }^{-\frac{1}{l+1}}\end{array}$$

ここで，Bはオイラーのベータ関数である。ENOPの近似式k_Aは単純な式となり，サンプリング点の分布に依存しない。ベータ関数値は，解析的となる。例えば，スプライン次数が3のとき，ベータ関数値は次式となる。   

$${\rm B}\left(1/4,3/4\right)=\sqrt{2}\pi$$

そして，ENOPの近似式は次式となる。   

$$k_A\left(\gamma \right)=\frac{\sqrt{2}L}{4}{\gamma }^{-\frac{1}{4}}$$(43)

以上のENOPの近似式(43)の導出は形式的である。式(42)はCraven and Wahba²⁾およびGolubら³⁾中でも見ることができ，理論詳細はこれらの論文に記載されている。

Fig.1は，スプライン次数が3，サンプリング数が101，サンプリング領域の長さが1.0であるときの，式(43)から得られるENOPの近似値k_Aと式(19)から得られるENOPの厳密値k_GCVを示している。ENOPの厳密値と比較するために，ENOPの近似式にバイアスが加えられている。経験的には，バイアス値は(l+1)/4となる。しかし，バイアスは，GCVあるいは情報量規準においては重要ではない。ENOPの近似値k_Aは，サンプリング数の半分(約50)より小さなENOPの領域では，ENOPの厳密値k_GCVによく一致しているが，サンプリング数の半分より大きなENOPの領域では，ENOPの厳密値k_GCVに一致しなくなる。しかし，サンプリング数の半分より大きなENOPは使用されるべきでないと普通いわれている。もし最適なGCVあるいは情報量規準により決定されるENOPがサンプリング数の半分より大きくなるならば，サンプリング数を増やすべきである。

Fig. 1.

 Equivalent number of parameters.

6.サンプリング測度重み付き平滑化薄板スプライン回帰

重み付き平滑化薄板スプライン(TPS)の汎関数は，次式で定義される。   

$$\begin{array}{c}\Pi =\sum _{i=1}^m\frac{1}{2}{\omega }_M^{\left(i\right)}{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)},y_M^{\left(i\right)}\right)\right)}^2+\gamma {\int }_{\Omega }\frac{1}{2}{\left({\nabla }^{2}f\right)}^{2}d\Omega \\ \nabla =\left[\begin{array}{cc}\frac{\partial }{\partial x}& \frac{\partial }{\partial y}\end{array}\right]\end{array}$$(44)

ここで，Ωはサンプリング領域の面積，(x⁽ⁱ⁾_M，y⁽ⁱ⁾_M)(i=1,…,m)はi番目のサンプリング点の位置である。重みが次式で定義されると考える。   

$$\begin{array}{c}{\int }_0^{\Omega }g\left(x,y\right)d\Omega \approx \sum _{i=1}^m{\omega }_M^{\left(i\right)}g\left(x_M^{\left(i\right)},y_M^{\left(i\right)}\right)\\ \sum _{i=1}^m\left[{\omega }_M^{\left(i\right)}\right]=\Omega ,{\omega }_M^{\left(i\right)}\ge 0\left(i=1,\cdots ,m\right)\end{array}$$(45)

ここで，g(x, y)は任意の関数である。式(45)は数値積分公式を示している。例えば，重みは数値求積法^11,15)から求められる。ただし，この求積法では，サンプリング測度ω⁽ⁱ⁾_M(i=1,…,m)の正値を保証できない。もしサンプリング測度が負となる場合には，対応するサンプリング点を除去する等の工夫が必要である。

変分原理を式(44)の汎関数に適用すると，   

$$\begin{array}{c}\delta \Pi ={\int }_{\Omega }\left\{\gamma {\nabla }^{4}f+\sum _{i=1}^m{\omega }_M^{\left(i\right)}\left(f\left(x_M^{\left(i\right)},y_M^{\left(i\right)}\right)-z_M^{\left(i\right)}\right)\text{Dirac}\left(x-x_M^{\left(i\right)},y-y_M^{\left(i\right)}\right)\right\}\delta fd\Omega \\ -\gamma {\int }_{\Gamma }\left(\nabla \left({\nabla }^{2}f\right)\cdot \mu \right)\delta fd\Gamma +\gamma {\int }_{\Gamma }{\nabla }^{2}f\left(\delta \left(\nabla f\right)\cdot \mu \right)d\Gamma =0\end{array}$$(46)

ここで，Γは領域Ωの境界であり，μは境界Γの単位法線ベクトルである。式(46)から，平滑化TPSの偏微分方程式は，   

$$\gamma {\nabla }^{4}f+\sum _{i=1}^m{\omega }_M^{\left(i\right)}\left(f\left(x_M^{\left(i\right)},y_M^{\left(i\right)}\right)-z_M^{\left(i\right)}\right)\text{Dirac}\left(x-x_M^{\left(i\right)},y-y_M^{\left(i\right)}\right)=0$$(47)

となり，境界条件は，   

$${\nabla }^{2}f=\nabla \left({\nabla }^{2}f\right)\cdot \mu =0\text{on}\Gamma$$

となる。一般的に，偏微分方程式(47)の解は，解析的ではないが，もし領域Ωが無限大であるならば，偏微分方程式(47)の解は単純になり，回帰曲面は，   

$$f\left(x,y\right)=c_1+c_{2}x+c_{3}y+\sum _{i=1}^m\frac{d_M^{\left(i\right)}}{8\pi }{r}_M^{\left(i\right)}{\left(x,y\right)}^2\log \left(r_M^{\left(i\right)}\left(x,y\right)\right)$$(48)

となり，拘束条件は，   

$$\sum _{i=1}^m{d}_M^{\left(i\right)}=\sum _{i=1}^m{x}_M^{\left(i\right)}{d}_M^{\left(i\right)}=\sum _{i=1}^m{y}_M^{\left(i\right)}{d}_M^{\left(i\right)}=0$$

となる。ここで，   

$$r_M^{\left(i\right)}\left(x,y\right)=\sqrt{{\left(x-x_M^{\left(i\right)}\right)}^2+{\left(y-y_M^{\left(i\right)}\right)}^2}\left(i=1,\cdots ,m\right)$$

式(48)のマトリクス表現による回帰曲面は，次式となる。   

$$f\left(x,y\right)=\left[a_M\left(x,y\right)\right]\left\{d_M\right\}+\left[p\left(x,y\right)\right]\left\{c\right\}$$

ここで，   

$$\begin{array}{c}\left[a_M\left(x,y\right)\right]=\left[\begin{array}{ccc}\frac{1}{8\pi }{r}_M^{\left(1\right)}{\left(x,y\right)}^2\log \left(r_M^{\left(1\right)}\left(x,y\right)\right)& \cdots & \frac{1}{8\pi }{r}_M^{\left(m\right)}{\left(x,y\right)}^2\log \left(r_M^{\left(m\right)}\left(x,y\right)\right)\end{array}\right]\\ \left[p\left(x,y\right)\right]=\left[\begin{array}{ccc}1& x& y\end{array}\right]\\ {}^T\left\{c\right\}=\left[\begin{array}{ccc}{c}_1& c_2& c_3\end{array}\right]\end{array}$$

未知のパラメータ{c}および{d_M}は，スプライン次数l=3を代入すると式(28)から求められる。GCVの評価関数は式(32)であり，ハット行列は式(33)である。

もしサンプリング数が無限大ならば，式(47)のサンプリング測度重み付き平滑化TPSの偏微分方程式は次式となる。   

$$\gamma {\nabla }^{4}f+f-z_M=0$$(49)

ラプラス変換を式(49)に適用すると，伝達関数は，次式で与えられる。   

$$h\left(s_1,s_2\right)=\frac{1}{\gamma {\left(s_1{}^2+s_2{}^2\right)}^2+1}$$(50)

ここで，s₁とs₂はラプラス変数である。フーリエ変換を式(50)の伝達関数に適用すると，次式の周波数応答関数が得られる。   

$$h\left({\omega }_1,{\omega }_2\right)=\frac{1}{\gamma {\left({\left(j{\omega }_1\right)}^2+{\left(j{\omega }_2\right)}^2\right)}^2+1}=\frac{1}{\gamma {\left({\omega }_1{}^2+{\omega }_2{}^2\right)}^2+1}$$(51)

ここで，ω₁とω₂は角周波数である。式(51)から，サンプリング測度重み付き平滑化スプライン回帰は位相補償ガウシアンフィルターとなる。領域Ωが長方形(L₁×L₂)とし，式(40)の角周波数と自由度の関係式を式(51)の周波数応答関数に代入すると，   

$$h\left(k_1,k_2\right)={\left(\gamma {\left({\left(\frac{\pi }{L_1}{k}_1\right)}^2+{\left(\frac{\pi }{L_2}{k}_2\right)}^2\right)}^2+1\right)}^{-1}$$

ここで，k₁とk₂は自由度である。一次元自由度空間に対する式(41)のENOPの定義を二次元自由度空間に拡張すると，平滑化TPSのENOPの近似式は，次式となる。   

$$k_A\left(\gamma \right)={\int }_0^{\infty }{\int }_0^{\infty }h\left(k_1,k_2\right)d{k}_{1}d{k}_2={\int }_0^{\infty }{\int }_0^{\infty }{\left(\gamma {\left({\left(\frac{\pi }{L_1}{k}_1\right)}^2+{\left(\frac{\pi }{L_2}{k}_2\right)}^2\right)}^2+1\right)}^{-1}d{k}_{1}d{k}_2$$(52)

変数変換   

$$\begin{array}{c}r={\left(\frac{\pi }{L_1}{k}_1\right)}^2+{\left(\frac{\pi }{L_2}{k}_2\right)}^2\\ \tan \theta =\frac{L_1{k}_2}{L_2{k}_1}\end{array}$$

を式(52)の積分に適用すると，ENOPの近似式は，次式となる。   

$$k_A\left(\gamma \right)=\frac{L_1{L}_2}{2{\pi }^2}{\int }_0^{\frac{\pi }{2}}d\theta {\int }_0^{\infty }\frac{1}{\gamma r^2+1}dr=\frac{L_1{L}_2}{4\pi }\times \frac{1}{2}{\rm B}\left(\frac{1}{2},\frac{1}{2}\right){\gamma }^{-\frac{1}{2}}=\frac{L_1{L}_2}{8}{\gamma }^{-\frac{1}{2}}$$(53)

式(53)のENOPの近似式は，長方形領域にのみ適用される。任意の領域に適用可能なように一般化されたENOPの近似式は，次式で与えられる。   

$$k_A\left(\gamma \right)=\frac{\Omega }{8}{\gamma }^{-\frac{1}{2}}$$(54)

汎関数が，   

$$\Pi =\sum _{i=1}^m\frac{1}{2}{\omega }_M^{\left(i\right)}{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)}\right)\right)}^2+\gamma {\int }_{\Omega }\frac{1}{2}{\left({\nabla }^{n+1}f\right)}^{2}d\Omega$$

として定義される平滑化n次多重調和スプラインのENOPの近似式は，次式となる。   

$$k_A\left(\gamma \right)=\frac{\Omega }{4\pi \left(n+1\right)}{\rm B}\left(\frac{1}{n+1},\frac{n}{n+1}\right){\gamma }^{-\frac{1}{n+1}}$$

Fig.2は，平滑化TPSが用いられ，サンプリング数が21*21=441であり，サンプリング領域が[0,1]×[0,1]の正方形であるときの，式(54)から得られるENOPの近似値と式(19)から得られるENOPの厳密値をしめしている。ENOPの厳密値と比較できるように，バイアスがENOPの近似式に加えられている。ENOPの近似値k_Aは，サンプリング数の半分(約220)より小さなENOPの領域で，ENOPの厳密値k_GCVとよく一致している。

Fig. 2.

 Equivalent number of parameters.

7.3Dスキャナーとサンプリング測度重み付き平滑化TPSによる鋼板表面形状推定

3Dスキャナーを用いて測定した鋼板表面のサンプリング値をFig.3に示す。鋼板の長さは5.475 m，幅は2.143 mである。x軸は鋼板の圧延方向であり，目盛り値は実際の値を10倍している。y軸は鋼板の幅方向である。10倍とした理由は，幅方向と比較して長手方向の形状の波長は短く，長手方向と幅方向の形状の波長をほぼ同じとするためである。サンプリング点は不均等に配置されており，サンプリング数は25691である。3DスキャナーはFARO社製Photon 120であり，測定値は約±2 mmの測定誤差(カタログ値)を含んでいる。Fig.3のサンプリング値の三角形メッシュ補間によるDEM (Digital Elevation Model)はFig.4である。鋼板表面のDEMはギザギザしている。測定対象は圧延された鋼板であるので，このようなギザギザは実際に存在していない。このギザギザは3Dスキャナーの測定誤差であるので，平滑化スプライン法を用いてギザギザを除去することを考える。

Fig. 3.

 Plate surface points cloud.

Fig. 4.

 Plate surface DEM.

サンプリング点数が多くなると(3000以上)，計算コストにおいて式(8)のシステムを解くことが困難となる。本ケースでは，Beatson and Newsa⁴⁾が提案したFMM(Fast Multipole Method)による高速計算法で式(48)を評価し，GMRES(Generalized Minimal Residual Method)により式(8)のシステムを解いた。GMRESで使用する行列の前処理(Preconditioning)では，Beatsonら⁵⁾が提案した近似的なCardinal Functionを用いる手法を参考とした。

本ケースでは，ベイズ情報量規準(Bayesian Information Criterion)   

$$BI{C}_A=m\log \left(\sum _{i=1}^m{\omega }_M^{\left(i\right)}{\left(z_M^{\left(i\right)}-f\left(x_M^{\left(i\right)},y_M^{\left(i\right)}\right)\right)}^2\right)-\log m\times m\log \left({\left(1-\frac{k_A}{m}\right)}^2\right)$$

を極小とする平滑化パラメータを最適な平滑化パラメータとした。本ケースでは，サンプリング数が多いためGCV法を適用することが実質的に不可能である。Fig.5に平滑化パラメータの探索結果を示す。横軸は平滑化パラメータ，縦軸はベイズ情報量規準である。点は探索中のベイズ情報量規準であり，白丸点は探索された最適な平滑化パラメータに対応するベイズ情報量規準である。最適な平滑化パラメータ値は5.09×10⁻³であり，対応するENOPは199.9であった。Fig.6に最適な平滑化パラメータで推定された鋼板の曲面を示す。Fig.5の計算結果は，従来の人による計測結果とほぼ一致することを確認した。

Fig. 5.

 Results of Bayesian information criterion.

Fig. 6.

 Result of optimal regression for plate surface.

本計算で使用した計算機の仕様は，2.66 GHz Intel(R)Core(TM)2 Quad CPU，3.25 GB RAM，Windows XP x86であり，アプリケーションの開発環境はMicrosoft Visual C++ 2005である。Fig.6の曲面を計算するための計算時間は9.8秒であり，最適な平滑化パラメータを求めるための計算時間は226秒であった。

製造条件が同じであれば，板のサイズが異なっていても製造された鋼板上に現れる波の波長分布はほぼ同じとなる。鋼板上に現れる波の代表的な半波長をL_ωとおき，式(38)と式(43)を用いると，平滑化スプライン曲線の平滑化パラメータは，   

$$\gamma =\frac{1}{64}{L}_{\omega }{}^4$$(55)

となり，鋼板の半波長のみに依存することがわかる。式(55)から，鋼板の半波長が変わらなければ(鋼板の製造条件が変わらなければ)，鋼板毎で適切な平滑化パラメータを計算する必要はなく，同じ平滑化パラメータを用いてもよい。Fig.7のように，一度求められた適切な平滑化パラメータを用いて複数の鋼板の曲面を推定することができる。

Fig. 7.

 Results of regressions for plate surfaces.

8. 結論

本論文では，3Dスキャナーで得られた鋼板表面の大規模点群データから平坦度曲面を再構成するための手法を提案し，サンプリング測度重み付き平滑化スプライン法における等価パラメータ数を高速に近似計算する方法，および実データを用いた平坦度曲面の計算を行い，以下の成果を得た。

(1)サンプリング測度重み付き平滑化スプライン法，およびその情報量規準を提案した。

(2)サンプリング測度重み付き平滑化スプラインにおいて，サンプリング数を無限とした場合のシステム伝達関数，周波数応答関数を求めた。

(3)周波数応答関数から等価パラメータ数の理論解を提案し，本理論解が実システムの等価パラメータ数の近似値であることを示した。

(4)等価パラメータ数の高速近似計算法と情報量規準を用いることで，高速に最適な平滑化パラメータを計算できることを示した。

(5)実際の大規模なサンプリング値からも，最適な平滑化曲面を効率的，現実的に求められることがわかり，工業的応用が可能であることを明らかにした。

文献

• 1)   H.  Akaike: IEEE Trans. Autom. Control, 19(1974), 716.
• 2)   P.  Craven and  G.  Wahba: Numer. Math., 31(1979), 377.
• 3)   G.H.  Golub,  M.  Heath and  G.  Wahba: Technometrics, 21(1979), 215.
• 4)   R.K.  Beatson and  G.N.  Newsa: Comput. Math. Appl., 24(1992), 7.
• 5)   R.K.  Beatson,  J.B.  Cherrie and  C.T.  Mouat: Adv. Comput. Math., 11(1999), 253.
• 6)   E.  Haber and  D.  Oldenburg: Comput. Geosci., 4(2000), 41.
• 7)   G.  De Nicolao,  G.  Ferrari-Trecate and  G.  Sparacino: Automatica, 36(2000), 1733.
• 8)   X.  Lin,  G.  Wahba,  D.  Xiang,  F.  Gao,  R.  Klein and  B.  Klein: Ann. Stat., 28(2000), 1570.
• 9)   G.  De Nicolao and  G.  Ferrari-Trecate: IEEE Trans. Neural Netw., 12(2001), 228.
• 10)   A.  Gtchev and  K.  Egiazarian: ISPA 2001 Proc. 2nd Int. Symp., IEEE, Pula, (2001), 580.
• 11)  青江信一郎:観測量解析方法，特開2002-365047, (2002).
• 12)   J.C.  Carr,  R.K.  Beatson,  J.B.  Cherrie,  T.J.  Mitchell,  W.R.  Fright,  B.C.  McCallum and  T.R.  Evans: Proc. 28th Annu. Conf. Computer Graphics and Interactive Technique, ACM, New York, (2001), 67.
• 13)   M.  Unser and  T.  Blu: IEEE Trans. Signal Proces., 53(2005), 2146.
• 14)   N.A.  Gumerov and  R.  Duraiswami: J. Comput. Phys., 215(2006), 363.
• 15)   A.  Sommariva and  M.  Vianello: Computing, 76(2006), 295.
• 16)   S.  Aoe and  G.M.  Hulbert: 9th National Cong. on Computational Mechanics, USACM, San Francisco, (2007).
• 17)  青江信一郎，植竹規人，林宏優:平面形状の測定方法，特開2012-37313, (2012).
• 18)   N.A.  Gumerov and  R.  Duraiswami: SIAM J. Sci. Comput., 29(2007), 1876.