(57)【特許請求の範囲】
【請求項１】
３次元スキャナーからなるデータ値計測装置と、
サンプリング測度演算器、平滑化演算器および画像処理器をそなえる計算機と、
処理した画像を表示する出力装置
を有する平面形状測定装置を用いて、被検体表面の平面形状を測定するに際し、
平滑化演算器において、座標系（ｘ，ｙ）上に任意に分布するｍ個のサンプリング点
【数１】
0005540977.tif 000070
と、各サンプリング点に対応した誤差を含むデータ値
【数２】
0005540977.tif 000071
から曲面を求めるために平滑化スプライン法を用いるものとし、その際、下記式（１）で示される汎関数Π(ｆ)を極小化する曲面ｆ（ｘ，ｙ）を回帰曲面とすることを特徴とする平面形状の測定方法。
記
【数３】
0005540977.tif 000072
【数４】
0005540977.tif 000073
ｌ・・・スプライン次数
γ・・・平滑化パラメータ
式（１）の右辺積分項の被積分関数は、
ｌが奇数のとき、次式（２）
【数５】
0005540977.tif 000074
ｌが偶数のとき、次式（３）
【数６】
0005540977.tif 000075
で示される。
ただし、上記サンプリング測度〔数４〕は、次式（４）を満たすものとする。
【数７】
0005540977.tif 000076
ここで、Ａは全てのサンプリング点を含む座標系（ｘ，ｙ）上の領域（サンプリング領域）の面積である。
【請求項２】
平滑化パラメータγを決定するために用いる情報量規準の残差平方和ＲＳＳを、次式（５）とすることを特徴する請求項１に記載の平面形状の測定方法。
【数８】
0005540977.tif 000077
【請求項３】
平滑化パラメータγを決定するために用いる情報量規準のパラメータ数ｋを、次式（６）を用いて計算することを特徴する請求項２に記載の平面形状の測定方法。
【数９】
0005540977.tif 000078
ここで、Ｂ（ａ，ｂ）は、次式（７）で定義されるオイラーのベータ関数
【数１０】
0005540977.tif 000079
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、平面形状の測定方法に関するものである。
【背景技術】
【０００２】
従来から、例えば原料ヤードの原料山や厚板の平坦度、コークス炉の内壁面等の被検体の表面状態を把握するために、平滑化スプライン法を用いた平面形状の測定方法が使用されている。
【０００３】
例えば非特許文献１に開示されている平滑化スプライン法は、座標系（ｘ，ｙ）上に任意に分布するｍ個のサンプリング点
【数１】
0005540977.tif 000002
と、各サンプリング点に対応した誤差を含むデータ値
【数２】
0005540977.tif 000003
から回帰曲面ｆ（ｘ，ｙ）を推定するために用いられる方法である。
【０００４】
この平滑化スプライン法において、平滑化スプライン曲面ｆ（ｘ，ｙ）を定義する汎関数Πは、次式（８）で与えられる。
【数３】
0005540977.tif 000004
ここで、ｌはスプライン次数と呼ばれる所定の自然数、γは平滑化パラメータと呼ばれる所定の正の実数である。
上記式（８）の右辺積分項の被積分関数は、
ｌが奇数のときは、次式（９）
【数４】
0005540977.tif 000005
で、一方ｌが偶数のときは、次式（１０）
【数５】
0005540977.tif 000006
で示される。
【０００５】
ここに、平滑化スプライン曲面は、式（８）の汎関数を極小とする曲面として定義される。式（８）から、平滑化スプライン曲面の偏微分方程式は、次式（１１）となる。
【数６】
0005540977.tif 000007
ここで、Dirac（ｘ，ｙ）はディラックのデルタ関数である。式（１１）から、解析的に次式（１２）で示される平滑化スプライン曲面と次式（１３）で示される拘束条件式が求められる。
【数７】
0005540977.tif 000008
【数８】
0005540977.tif 000009
ここで、
【数９】
0005540977.tif 000010
【数１０】
0005540977.tif 000011
は、未知パラメータであり、また
【数１１】
0005540977.tif 000012
である。
平滑化スプライン曲面を示す式（１２）の第２項は、グリーン関数の和であり、また式（１２）の平滑化スプライン曲面は、Ｃ2l-1連続性をもつことが知られている。
【０００６】
以下、未知パラメータ〔数９〕、〔数１０〕を求める方法について述べる。式（１２）をマトリクス形式で表現すると次式（１５）となり、また式（１３）をマトリクス形式で表現すると次式（１６）となる。
【数１２】
0005540977.tif 000013
【数１３】
0005540977.tif 000014
ここで、｛Ｏ(l+1)(l+2)/2｝ は、(l+1)(l+2)/2行１列の零ベクトル、
【数１４】
0005540977.tif 000015
【数１５】
0005540977.tif 000016
【数１６】
0005540977.tif 000017
【数１７】
0005540977.tif 000018
である。式（１２）を式（１１）へ代入すると、次式（２１）となり、
【数１８】
0005540977.tif 000019
また式（１５）を式（２１）へ代入すると、次式（２２）となる。
【数１９】
0005540977.tif 000020
ここで、
【数２０】
0005540977.tif 000021
【数２１】
0005540977.tif 000022
【０００７】
式（２４）と式（１６）を用いると、次式（２５）に示す連立一次方程式が得られる。
【数２２】
0005540977.tif 000023
ここで、［ＩmXm］はｍ行ｍ列の単位行列、［Ｏ(l+1)(l+2)/2X(l+1)(l+2)/2］は(l+1)(l+2)/2行(l+1)(l+2)/2列の零行列である。
式（２５）を数値的に解くことで、未知パラメータ｛ｄM｝、｛ｃ｝を求めることができ、次式（２６）となる。
【数２３】
0005540977.tif 000024
そして式（２６）を式（１５）に代入すると、次式（２７）に示す平滑化スプライン曲面が求められる。
【数２４】
0005540977.tif 000025
【０００８】
次に、平滑化パラメータγを求める方法について述べる。
平滑化スプライン法において、平滑化パラメータを求めるためには、一般化交差検証法（Generalized Cross Validation: GCV法）が用いられる。ＧＣＶ法では、次式（２８）の評価関数を極小化する平滑化パラメータを最適な平滑化パラメータとする。
【数２５】
0005540977.tif 000026
ここで、ＲＳＳは次式（２９）で表される残差平方和であり、［Ｈ(γ)］はハット行列と呼ばれるもので、次式（３０）で定義される。
【数２６】
0005540977.tif 000027
【数２７】
0005540977.tif 000028
ここで、
【数２８】
0005540977.tif 000029
【数２９】
0005540977.tif 000030
【０００９】
式（２７）の平滑化スプライン曲面にサンプリング点群（｛ｘM｝,｛ｙM｝） を代入した式と、式（３０）のハット行列の定義式を比較することで、平滑化スプライン曲面のハット行列は次式（３３）となる。
【数３０】
0005540977.tif 000031
【００１０】
また、式（２８）のＧＣＶ法の評価関数は、情報量規準と関連があることが知られている。特に、赤池情報量規準は、サンプリング数が大きければ、ＧＣＶ法の評価関数と等価となる。なお、赤池情報量規準は、次式（３４）で定義される。
【数３１】
0005540977.tif 000032
ここで、ＬＬは次式（３５）で定義される最大対数尤度、ｋはモデルのパラメータ数あるいは自由度である。
【数３２】
0005540977.tif 000033
【００１１】
式（２８）のＧＣＶ法の評価関数は、単調増加関数を用いて次式（３６）で示される情報規準の評価関数ICGCVに変換される。
【数３３】
0005540977.tif 000034
サンプリング数ｍを無限大とすると、式（３６）の情報量規準の評価関数は次式（３７）となる。
【数３４】
0005540977.tif 000035
式（３４）と式（３７）を比較すると、式（３７）の右辺第二項のtrace([Ｈ(γ)])は、平滑化スプライン曲面に対するパラメータ数とみなすことができる。等価パラメータ数は、次式（３８）で定義される。
【数３５】
0005540977.tif 000036
【先行技術文献】
【特許文献】
【００１２】
【特許文献１】特開2002−365047号公報
【非特許文献】
【００１３】
【非特許文献１】「Smoothing Noisy Data with Spline Functions (Estimating the Correct Degree of smoothing by the Method of Generalized Cross-Validation), Numer. Math., 31(1979),pp. 377-403」
【発明の概要】
【発明が解決しようとする課題】
【００１４】
最適な平滑化パラメータを求めるためには、式（２８）の評価関数Ｖ(γ)を計算する必要がある。しかしながら、式（２８）で用いられる式（３３）のハット行列は、ｍ＋(l+1)(l+2)/2行ｍ＋(l+1)(l+2)/2列の逆行列を含んでおり、サンプリング数ｍが大きくなるほど、この逆行列の計算時間が長くなってしまい、実用的でなくなる。
【００１５】
本発明は、上記の問題を有利に解決するもので、従来よりも計算時間の短縮化を図り、従来と遜色のない誤差および滑らかさの下で曲面を把握ことが可能な、新しい平面形状の測定方法を提案することを目的とする。
【課題を解決するための手段】
【００１６】
以下、本発明の解明経緯について説明する。
図１に、本発明の実施に用いて好適な平面形状測定装置の概略を示す。
この測定装置は、３次元スキャナーからなるデータ値計測装置１と、計算機２と、処理した画像を表示する出力装置３を有しており、また計算機２は、サンプリング測度演算器４、平滑化演算器５および画像処理器６をそなえている。なお、７は、平滑化データを保存するための記憶装置である。
【００１７】
平滑化スプライン曲面を推定するには、平滑化パラメータを求める必要があり、そのためには、図２に示すようなフローチャートに従って、最適な平滑化パラメータや未知パラメータを求めるのであるが、従来の方式では、かかるパラメータを算出するのに長時間を要していたことは、前述したとおりである。
【００１８】
そこで、発明者らは、等価パラメータ数を短時間で計算する方法について考察した。
すなわち、本発明では、平滑化スプライン曲面を定義する汎関数として、式（８）の代わりに式（３９）に示す汎関数を用いるのである。
【数３６】
0005540977.tif 000037
ここで、
【数３７】
0005540977.tif 000038
は、サンプリング測度（サンプリング点の相当面積）である。
また、このサンプリング測度は、次式（４０）を満たすものとする。
【数３８】
0005540977.tif 000039
ここで、Ａは全てのサンプリング点を含む座標系（ｘ，ｙ）上の領域（サンプリング領
域）の面積である。
【００１９】
もし、サンプリング点がサンプリング領域内に均一に存在する場合には、サンプリング測度は次式（４１）で与えられる。
【数３９】
0005540977.tif 000040
一方、サンプリング点がサンプリング領域内に不均一に存在する場合には、例えば、特許文献１の方法でサンプリング測度を求めることができる。平滑化スプライン曲面は、式（３９）の汎関数を極小とする曲面として定義される。
式（３９）から、平滑化スプライン曲面の偏微分方程式は、次式（４２）となる。
【数４０】
0005540977.tif 000041
式（４２）から、解析的に式（１２）の平滑化スプライン曲面と式（１３）の拘束条件式が求められる。
【００２０】
次に、未知パラメータ〔数９〕、〔数１０〕を求める方法について述べる。式（１２）を式（４２）へ代入すると、次式（４３）となる。
【数４１】
0005540977.tif 000042
式（１５）を式（４３）へ代入すると、次式（４４）が得られる。
【数４２】
0005540977.tif 000043
ここで、［ΩM］は、次式（４５）で示される。
【数４３】
0005540977.tif 000044
【００２１】
式（４４）と式（１６）を用いると、次式（４６）の連立一次方程式が得られる。
【数４４】
0005540977.tif 000045
式（４６）を数値的に解くことで、未知パラメータ｛ｄM｝、｛ｃ｝を求めることができ、次式（４７）となる。
【数４５】
0005540977.tif 000046
式（４７）を式（１５）に代入すると、次式（４８）の平滑化スプライン曲面が求められる。
【数４６】
0005540977.tif 000047
【００２２】
次に、平滑化パラメータγを求める方法について述べる。残差平方和を次式（４９）、ハット行列を次式（５０）とすることで、ＧＣＶ法あるいは情報量規準を用いて最適な平滑化パラメータを求めることができる。
【数４７】
0005540977.tif 000048
【数４８】
0005540977.tif 000049
式（３８）で定義される等価パラメータ数の近似式を求める。サンプリング数ｍが無限大ならば、式（４２）の偏微分方程式は次式（５１）となる。
【数４９】
0005540977.tif 000050
ここで、ｚMはサンプリングデータ曲面である。式（５１）にラプラス変換を適用すると、次式（５２）となる。
【数５０】
0005540977.tif 000051
ここで、ｓ1は変数ｘに対応するラプラス変数、ｓ2は変数ｙに対応するラプラス変数、ｈ（ｓ1,ｓ2）は伝達関数であり、次式（５３）となる。
【数５１】
0005540977.tif 000052
【００２３】
式（５２）と式（３０）のハット行列の定義式を比較すると、伝達関数ｈ（ｓ1,ｓ2） は、ハット行列と同じ機能をしていることがわかる。フーリエ変換を式（５３）の伝達関数に適用すると、次式（５４）の周波数応答関数ｈ（ω1,ω2）が得られる。
【数５２】
0005540977.tif 000053
ここで、ｊは虚数単位、ω1は変数ｘに対応する角周波数、ω2は変数ｙに対応する角周波数である。式（５４）から、平滑化スプライン曲面がローパスフィルター、ガウシアンフィルターとして機能していることがわかる。
【００２４】
次に、三角関数に対する自由度ｋを次式（５５）として定義する。
【数５３】
0005540977.tif 000054
ここで、Ｌは三角関数が存在する領域の長さ、Ｌωは次式（５６）で与えられる三角関数の半波長である。
【数５４】
0005540977.tif 000055
ここで、ωは三角関数の各周波数である。
式（５６）を式（５５）へ代入すると、次式（５７）が求まる。
ｋ＝（Ｌ／π）・ω            ・・・  （５７）
サンプリング領域を長方形とすると、変数ｘに対応する角周波数ω1に対する自由度ｋ1は式（５８）となり、変数ｙに対応する角周波数ω2に対する自由度ｋ2は式（５９）となる。
ｋ1＝（Ｌ1／π）・ω1         ・・・  （５８）
ｋ2＝（Ｌ2／π）・ω2         ・・・  （５９）
ここで、Ｌ1はサンプリング領域のｘ軸方向の長さ、Ｌ2はサンプリング領域のｙ軸方向の長さである。
式（５８）と式（５９）を式（５４）の周波数応答関数に代入すると、次式（６０）となる。
【数５５】
0005540977.tif 000056
【００２５】
式（３８）の等価パラメータ数の定義から、等価パラメータ数はハット行列の固有値の総和である。したがって、無限サンプリング数に対する等価パラメータ数は、次式（６１）で定義される。
【数５６】
0005540977.tif 000057
式（６１）で定義される等価パラメータ数は、式（３８）で定義される等価パラメータ数の近似値となる。式（６２）と式（６３）の変数変換を、式（６１）右辺の積分に適用すると、等価パラメータ数の近似式は次式（６４）となる。
【数５７】
0005540977.tif 000058
【数５８】
0005540977.tif 000059
【数５９】
0005540977.tif 000060
ここで、ｒ，θは変数、Ｂ（ａ,ｂ）は次式（６５）で定義されるオイラーのベータ関数である。
【数６０】
0005540977.tif 000061
式（６４）の等価パラメータの近似式は、サンプリング領域が長方形である場合にのみに適用される。任意のサンプリング領域に適用可能なように一般化された等価パラメータ数の近似式は、次式（６６）となる。
【数６１】
0005540977.tif 000062
【００２６】
等価パラメータ数の近似値は、式（６６）を計算することにより求めることができ、式（６６）の計算時間は、式（３８）の計算時間よりも非常に短くなる。
また、式（６６）の等価パラメータ数の近似値は、サンプリング点の分布と無関係に求めることができる。
このことは、式（３９）の平滑化スプライン曲面を定義する汎関数に、サンプリング測度を導入したことに由来している。平滑化スプライン曲面もサンプリング点の分布に依存しなくなる。
本発明は、上記の知見に立脚して得られたものである。
【００２７】
すなわち、本発明の要旨構成は次のとおりである。
１．３次元スキャナーからなるデータ値計測装置と、
サンプリング測度演算器、平滑化演算器および画像処理器をそなえる計算機と、
処理した画像を表示する出力装置
を有する平面形状測定装置を用いて、被検体表面の平面形状を測定するに際し、
平滑化演算器において、座標系（ｘ，ｙ）上に任意に分布するｍ個のサンプリング点〔数１〕と、各サンプリング点に対応した誤差を含むデータ値〔数２〕から曲面を求めるために平滑化スプライン法を用いるものとし、その際、〔数３６〕で示される汎関数Π(ｆ)を極小化する曲面ｆ（ｘ，ｙ）を回帰曲面とすることを特徴とする平面形状の測定方法。
上記式（３９）において、〔数３７〕はサンプリング測度、ｌはスプライン次数、γは平滑化パラメータである。
式（３９）の右辺積分項の被積分関数は、ｌが奇数のとき〔数４〕で、一方ｌが偶数のときは〔数５〕で示される。
ただし、サンプリング測度〔数３７〕は、式（４０）を満たすものとする。
式（４０）において、Ａは全てのサンプリング点を含む座標系（ｘ，ｙ）上の領域（サンプリング領域）の面積である。
【００２９】
２．平滑化パラメータγを決定するために用いる情報量規準の残差平方和ＲＳＳを、〔数４７〕とすることを特徴する上記１に記載の平面形状の測定方法。
【００３０】
３．平滑化パラメータγを決定するために用いる情報量規準のパラメータ数ｋを、〔数６１〕を用いて計算することを特徴する上記２に記載の平面形状の測定方法。
上記式（６６）において、Ｂ（ａ，ｂ）は、〔数６０〕で定義されるオイラーのベータ関数である。
【発明の効果】
【００３１】
本発明によれば、適切な平滑化スプライン曲面を高速、かつ短時間で算出することができる。
【図面の簡単な説明】
【００３２】
【図１】本発明の実施に用いて好適な平面形状測定装置の概略図である。
【図２】平滑化パラメータの演算要領を示すフローチャートである。
【図３】サンプリング点とサンプリング測度との関係を示した図である。
【図４】平滑化パラメータγと等価パラメータ数ｋとの関係を示した図である。
【図５】式（７０）で求めたサンプリング値を示した図である。
【図６】推定された式（６９）により得られた平滑化スプライン曲面を示した図である。
【図７】式（６８）から求めた情報量規準の評価関数を示した図である。
【図８】計測データをｘｙｚ座標系でプロットした図である。
【図９】（θ，φ）投影面上にデータ点を投影したプロットを示した図である。
【図１０】計測データを間引きした、曲面推定のためのデータを示した図である。
【図１１】平滑化処理して推定された、平滑化スプライン曲面を示した図である。
【発明を実施するための形態】
【００３３】
以下、本発明を具体的に説明する。
まず、本発明の第一の実施形態について説明する。式（６６）の等価パラメータ数近似式と式（３８）の等価パラメータ数を比較する。スプライン次数ｌを１とすると、式（６６）の等価パラメータ数近似式は次式（６７）となる。
【数６２】
0005540977.tif 000063
【００３４】
図３に示すように、サンプリング領域を［０,１］ラ［０,１］の正方形とし、サンプリング点がサンプリング領域に均等に４４１個（ｘ軸方向に２１個ラｙ軸方向に２１個）存在する場合について考える。図３の高さ方向の軸は、特許文献１の計算法で求められたサンプリング点に対応するサンプリング測度の値である。サンプリング点が均等に配置されているので、サンプリング測度は、境界を除いてほぼ同じ値となる。図４に、式（６６）の等価パラメータ数近似式と、式（３８）の等価パラメータ数を比較したグラフを示す。図４の横軸は平滑化パラメータ、縦軸は等価パラメータ数である。  等価パラメータ数近似値と等価パラメータ数の間には、スプライン次数ｌに依存するバイアス差が存在する。本実施形態では、等価パラメータ数近似値に約２を加えると等価パラメータ数となる。しかしながら、式（３４）の情報量規準の評価式から、バイアス差は、最適な平滑化パラメータに影響を与えない。
【００３５】
等価パラメータ数近似値は、サンプリング数の半分より小さな等価パラメータ数の領域では、等価パラメータ数の厳密値とよく一致している。これに対し、サンプリング数の半分より大きな等価パラメータ数の領域では、等価パラメータ数の厳密値と一致しなくなる。等価パラメータ数近似式は、サンプリング数が無限大である場合の近似式であるため、等価パラメータ数がサンプリング数に近づいていくと一致しなくなる。サンプリング数の半分より大きな等価パラメータ数は使用されるべきでないと一般的にいわれており、もし推定された最適な等価パラメータ数がサンプリング数の半分よりも大きくなるようであれば、サンプリング数を増やすべきである。
【００３６】
次に、本発明の第二の実施の形態について説明する。
平滑化スプライン曲面を求めには、情報量規準の評価関数として、式（３６）の等価パラメータ数を等価パラメータ数近似値に置き換えた次式（６８）の評価関数を用いる。
【数６３】
0005540977.tif 000064
スプライン次数ｌは１とすると、平滑化スプライン曲面は次式（６９）となる。
【数６４】
0005540977.tif 000065
サンプリング領域、サンプリング点の配置は、第一の実施形態と同じものであり、図３に示したとおりである。サンプリング測度は、特許文献１に開示されている方法を用いた。サンプリング値を次式（７０）とする。
【数６５】
0005540977.tif 000066
ここで、ｚExact（ｘ,ｙ）は厳密曲面であり、式（７１）とする。
【数６６】
0005540977.tif 000067
但し、
【数６７】
0005540977.tif 000068
は、正規分布Ｎ（０，0.12）（平均値：０、分散：0.12）にしたがう計測誤差である。
推定される平滑化スプライン曲面が、厳密曲面ｚExact（ｘ,ｙ）に近づくほど、よい推定である。
【００３７】
図５に、式（７０）のサンプリング値
【数６８】
0005540977.tif 000069
を示す。
図６に、推定された式（６９）の平滑化スプライン曲面を示す。
図７に、式（６８）の情報量規準の評価関数を示す。
図７から、評価関数を極小化する最適な平滑化パラメータ値は1.0×10-5であり、最適な平滑化パラメータ値を式（６７）へ代入すると、最適な等価パラメータ数は約４５となる。
【００３８】
サンプリング数４４１に比べて最適な等価パラメータ数は小さいので、等価パラメータ数近似値を用いても問題はない。図６から、平滑化スプライン曲面は厳密曲面をよく近似していることがわかる。計算で使用した計算機のメモリ容量は５０４ＭＢ、ＭＰＵはＰｅｎｔｉｕｍ（登録商標）Ｍである。
式（３８）の等価パラメータ数による最適な平滑化パラメータの計算には、３３．３秒かかったのに対し、式（６６）の等価パラメータ数近似式による最適な平滑化パラメータの計算は８．７秒で済み、大幅に計算時間を短縮することができた。
【００３９】
次に、本発明の第三の実施形態について説明する。
ＦＡＲＯ社製レーザースキャナー装置を用いて、板の表面形状を計測した。
図８に、計測データのｘｙｚ座標系プロットを示すが、２２２４２個のデータ点が計測された。レーザースキャナーから得られるデータは、レーザー発信源を中心とした球面座標系（ｒ，θ，φ）である。距離ｒの計測データには、±３ｍｍの計測誤差がある。
図９に、（θ，φ）投影面上にデータ点を投影したプロットを示す。計測データには、板形状以外のデータも含まれているので、板形状のデータ部分のみを抽出する。さらに、データ数が多いと、計算負荷が高くなるので、計測データの間引きを行い、曲面推定のためのデータを作成した。
得られたデータを、図１０に示す。データ点数は４７１個である。
そして、図１０のデータを平滑化処理して推定された、平滑化スプライン曲面を図１１に示す。計測データには計測誤差が含まれているが、推定曲面の推定誤差は、最小二乗法に則って、極小化される。
図１１に示したように、本発明に従う推定曲面は滑らかであり、微分演算を行うことによって、曲面の勾配や曲率等の計量データを計算することが可能である。
【符号の説明】
【００４０】
１  データ値計測装置
２  計算機
３  出力装置
４  サンプリング測度演算器
５  平滑化演算器
６  画像処理器
７  記憶装置
【図１】
0005540977.tif 000080
【図２】
0005540977.tif 000081
【図３】
0005540977.tif 000082
【図４】
0005540977.tif 000083
【図５】
0005540977.tif 000084
【図６】
0005540977.tif 000085
【図７】
0005540977.tif 000086
【図８】
0005540977.tif 000087
【図９】
0005540977.tif 000088
【図１０】
0005540977.tif 000089
【図１１】
0005540977.tif 000090