計測点からの円筒情報の抽出方法
【課題】
本発明の課題は、円筒表面の計測点から円筒データを確実に作成できると共に、高精度な近似が可能な計測点からの円筒情報の抽出方法を提供することにある。
【解決手段】
上記課題を解決するために、本発明では、最初に計測点から球の中心を定め、次に、球の中心を円筒の仮中心点とし２つの角度を変化させて投影方向を定めて計測点を投影し、円周と投影した計測点の距離の二乗の総和が最小になる角度を選択し、次に、当該角度を収束計算で精度向上させ、最後に、精度向上させた角度、仮中心点を当該角度で投影した値を円筒中心の初期値とし、円周と投影した距離の二乗の総和が最小になる角度および円の中心点を定めることによって、短時間で精度の良い円筒の半径、中心、中心軸を得ることを特徴とする。
本発明の課題は、円筒表面の計測点から円筒データを確実に作成できると共に、高精度な近似が可能な計測点からの円筒情報の抽出方法を提供することにある。
【解決手段】
上記課題を解決するために、本発明では、最初に計測点から球の中心を定め、次に、球の中心を円筒の仮中心点とし２つの角度を変化させて投影方向を定めて計測点を投影し、円周と投影した計測点の距離の二乗の総和が最小になる角度を選択し、次に、当該角度を収束計算で精度向上させ、最後に、精度向上させた角度、仮中心点を当該角度で投影した値を円筒中心の初期値とし、円周と投影した距離の二乗の総和が最小になる角度および円の中心点を定めることによって、短時間で精度の良い円筒の半径、中心、中心軸を得ることを特徴とする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、円筒形状の対象物を計測して得た計測点から円筒の中心点、中心軸、半径（以下、円筒データ）を抽出する計測点からの円筒情報の抽出方法に関するものである。
【背景技術】
【０００２】
計測点などの測定誤差を持つ点群から円筒データを抽出する従来の方法の一例を、以下に説明する。
【０００３】
従来の方法は、次のようなステップで行われている。即ち、ステップ１）：円筒データの初期値を設定する。ステップ２）：当該円筒表面と計測点群の距離を作る。ステップ３）：距離の二乗和が最低になるよう円筒データの更新量を定める。ステップ４）：当該更新量が収束判定値以下なら終了。そして、当該二乗和が、収束判定値以上なら当該更新量により初期値を更新してステップ２）から繰り返す。
【０００４】
ステップ１）〜４）の方法は、解に十分近い初期値を与えれば収束させることができる。しかし、円筒の表面と計測点群との距離の二乗和は、極値を複数持つ複雑な関数であるため、解から離れた初期値が与えられた場合、収束が遅くなったり、別の解に収束してしまう問題がある。
【０００５】
この問題を解決する方法として、あらかじめ円筒データが存在する範囲を限定し、当該範囲を細かく分割して良好な初期値を定める方法が知られているが、範囲の設定方法が不明確なこと、どの程度の細かさまで分割すれば良いかが不明確なこと、分割数が多くなった場合の処理時間の増大が問題になっている。
【０００６】
また、関連する技術として、非特許文献１、２及び３に記載されている以下の（１）、（２）、（４）、特に、先行技術文献は提示しないが（３）の技術が知られている。
（１）３点を通過する円
非特許文献１に記載されているのは、３次元空間内に３点Ｐ０、Ｐ１、Ｐ２を与え、３点を通過する円の中心点を作成する方法である。
【０００７】
これは、先ず３点中の２点を２組取出し、２点を結ぶ単位ベクトルを法線とし、２点の中間点を通過する平面を２枚作成する。更に、３点を通過する平面を１枚作成する。次に、当該３枚の平面が交差する１点を求める。これにより、３点を通過する円の中心点が作成できる。
【０００８】
（数１）は、３点Ｐ０、Ｐ１、Ｐ２と、円の中心点Ｏの関係を示す連立一次方程式である。（数１）において、添え字０、１、２は点を識別する番号、ｘ、ｙ、ｚは点の座標成分であり、符号×は、ベクトルの外積、符号・はベクトルの内積である。
【０００９】
【数１】
【００１０】
（数１）を連立一次方程式の解を求める方法、例えば、掃き出し法やＬＵ分解などで解けば、３点を通過する円の中心点を３次元空間内に定めることができる。
（２）複数の点に近い円
非特許文献２に記載されているのは、２次元空間内の複数の点Ｑｉを与え、複数の点と、円との距離の二乗和が最小になる円の中心点を作成する方法である。
【００１１】
（数２）は、当該方法の関係を整理した結果であり、（数２）を点数分作成した連立１次方程式を一般化逆行列を用いて解くことにより中心点を求める。
【００１２】
【数２】
【００１３】
（数２）は平面上の複数の点に対し、円の中心点を作成するために効率の良い方法である。
（３）複数の点に近い球
これは、３次元空間内の複数の点Ｐｉを与え、複数の点と球との距離の二乗和が最小になる球の中心点を作成する方法である。
【００１４】
（数３）は当該方法の関係を整理した結果であり、（数３）を点数分作成した連立１次方程式を一般化逆行列を用いて解くことにより中心点を求める。
【００１５】
【数３】
【００１６】
（数３）は、空間内の複数の点に対し、球の中心点を作成するために効率の良い方法である。
（４）３次元から２次元への座標変換
非特許文献３に記載されているのは、３次元の点Ｐｉを、１番目の軸まわりに角度Ｂだけ回転した後、２番目の軸まわりに角度Ａだけ回転する例である。ここで、３次元座標Ｐｉから２次元座標Ｑｉへの変換は、３次元の点Ｐｉを、１番目の軸まわりに角度Ｂだけ回転した後、２番目の軸まわりに角度Ａだけ回転した座標（ｘ、ｙ）成分を取り出したものである。例えば、１番目をｙ軸、２番目をｘ軸とした場合、（数４）となる。
【００１７】
【数４】
【００１８】
（数４）において、Ｏｉ，ｘとＱｉ，ｙを成分とするベクトルＱｉを角度Ａ、角度Ｂでｎ回偏微分または相互偏微分した値を求めることは公知である。また、（数４）は、Ｙ軸回転してからＸ軸回転して（ｘ、ｙ）成分を取り出しているが、Ｘ軸回転、Ｙ軸回転、Ｚ軸回転を任意に組み合わせることも公知である。その場合、（数４）の回転行列成分が変化する。
【先行技術文献】
【非特許文献】
【００１９】
【非特許文献１】コンピュータディスプレイによる形状処理工学［Ｉ］、ｐ.１９２ ３点を通る円、日刊興行新聞社、昭和５７年１０月３０日発行
【非特許文献２】Coope, I.D., Circle fitting by linear and nonlinear least squares, Journal of Optimization Theory and Applications Volume 76, Issue 2, New York: Plenum Press, February 1993
【非特許文献３】参考文献３：コンピュータディスプレイによる図形処理工学、ｐ.１１２ ｘ／ｙ／ｚ軸まわりの回転、日刊興行新聞社、１９８１年６月３０日発行
【発明の概要】
【発明が解決しようとする課題】
【００２０】
しかしながら、非特許文献１では、３点に限定していること及び３点で作成した平面上の円しか作成できないことから、複数の計測点に最も近い円の中心を求める場合には、適用が難しい。また、非特許文献２では、３次元空間で円筒の中心軸を定めるには、非特許文献３を用いて、３次元空間の複数の点を回転変換し、点から２つの成分を取り出すことにより平面に投影した後で適用する必要があり、この方法では、回転変換時の角度および投影後の中心点が変数になり、回転変換は三角関数を使用するため、角度で無限回の微分が可能である。即ち、極値を複数持つ複雑な関数となるので、複数の計測点に最も近い円の中心を求める場合には適用が難しい。更に、（３）の複数の点に近い球では、球ではなく円筒を定める場合には利用できない。
【００２１】
本発明は上述の点に鑑みなされたもので、その目的とするところは、円筒表面の計測点から円筒データを確実に作成できると共に、従来方式に比べて、高精度な近似が可能になる計測点からの円筒情報の抽出方法を提供することにある。
【課題を解決するための手段】
【００２２】
本発明の計測点からの円筒情報の抽出方法は、上記目的を達成するために、計測点が円筒物の表面を計測していることを条件として、全計測点から球の中心を抽出する手段と、全計測点および球の中心から角度の初期値を作成する手段と、当該角度の初期値の精度向上を行う手段と、球中心を精度向上した角度で投影変換し、初期円筒中心を定める手段と、初期円筒中心および精度向上した角度を初期値とし、円筒中心および角度を定める手段と、を用意し計測点から円筒の半径、中心軸からなる円筒データを抽出することを特徴とするものである。
【００２３】
更に、全計測点および球の中心から角度の初期値を作成する場合に、角度を離散化するための分割数、当該分割数と当該分割数における近似曲面との差の二乗の最大Ｂの初期値、当該分割数における近似曲面との差の二乗の最大Ａを最大Ｂで除した値で繰り返し終了を判定するための初期値を設定する手段と、２つの角度を当該分割数で離散化した格子リストを作成する手段と、当該格子リストの角度毎に、計測点と球の中心を座標変換して平面上に投影し、投影した計測点と投影した中心の二乗距離を作成し、当該二乗距離の平均を作成し、当該二乗距離と当該平均の差の二乗和を作成する手段と、当該格子リストに対応した二乗和から近似曲面を作成する手段と、当該格子リストに対応した二乗和と近似曲面との距離の二乗の最大Ａを求める手段と、最大Ａが最大Ｂより大きいかを判定する手段と、最大Ｂに最大Ａを代入する手段と、最大Ａを最大Ｂで除した値がしきい値より大きいかを判定する手段と、分割数を増加させる手段と、を用意し、全計測点および球の中心から角度の初期値を作成することを特徴とするものである。
【発明の効果】
【００２４】
本発明の計測点からの円筒情報の抽出方法によれば、円筒表面の計測点から円筒データを確実に作成できると共に、従来方式に比べて、高精度な近似が可能になる。
【図面の簡単な説明】
【００２５】
【図１】本発明の計測点からの円筒情報の抽出方法を実施する装置を示す構成図である。
【図２】図１の円筒データ抽出処理を実施する計算機の構成図である。
【図３】図１の円筒データ抽出処理の処理を示すフローチャートである。
【図４（ａ）】投影するための角度を変数とし、投影した計測点と円との距離の二乗の総和を縦軸とした３次元グラフの上面図である。
【図４（ｂ）】図４（ａ）の斜視図である。
【図４（ｃ）】図４（ａ）及び図４（ｂ）の基となる値を記入した図である。
【図５（ａ）】球中心を、投影した計測点と円との距離の二乗の総和が最小になる４組の角度で投影したときの円中心の図である。
【図５（ｂ）】４組の角度の角度値と円中心を記述した図である。
【図６】本発明の処理の内容を図示した説明図である。
【図７】初期角度算出部の別の実施方法を示すフローチャートである。
【図８（ａ）】角度Ａ，Ｂを離散化するときの説明図である。
【図８（ｂ）】二乗和と二乗和を近似した曲面との差を説明する説明図である。
【図８（ｃ）】分割数と当該分割数における当該差の二乗和の最大の変化を示す特性図である。
【図９】本発明で、球あてはめ、仮円筒および円筒あてはめをおこなった結果の３次元モデル、当該計測点、球、円筒を重ね合わせた３次元モデルを表示した例を示す図である。
【図１０】本発明を実施したときの球中心算出処理、初期角度算出処理、角度精度向上処理、円筒中心・角度算出処理の計測順番と計算時間を示す図である。
【発明を実施するための形態】
【００２６】
以下、本発明の計測点からの円筒情報の抽出方法の実施例を、図１乃至図１０を用いて説明する。
【実施例】
【００２７】
図１は、本発明の計測点からの円筒情報の抽出方法を実施する装置（後述する計算機で実行するプログラムでもある）を示すものであり、該図において、１０１は円筒物の表面を３次元座標（Ｘ、Ｙ、Ｚ）で計測する計測機器、１０２は計測機器１０１で計測した計測点、収束判定に用いる０とみなす値、初期角度算出部１０６で使用する角度のピッチを保存する記憶装置、１０３は計算機で実行する円筒データ抽出部、１０４は記憶装置１０２から３次元の計測点Ｐｉ、収束判定に用いる０とみなす値、初期角度算出部１０６で使用する角度のピッチを読み取る計測点・収束判定値・角度ピッチ読取部、１０５は計測点Ｐｉから球の中心を算出する球中心算出部、１０６は記憶装置に格納されている角度ピッチで分割した角度で、計測点および１０５で算出した球中心を回転変換し、回転変換した計測点と球中心のベクトルの２成分を取り出すことにより、投影した球中心と投影した計測点を定め、投影した球中心と投影した計測点の距離を求めると共に、当該距離の平均を求め、当該距離と当該平均の差の二乗和を求め、当該二乗和が最小になる初期角度を算出する初期角度算出部、１０７は１０６で作成した初期角度の精度向上を行う角度精度向上部、１０８は、１０６で作成した投影後の球中心を初期円筒中心、１０７で作成した角度を初期角度とし、円筒中心、角度、半径を算出する円筒中心・角度算出部、１０９は１０８で作成した円筒中心、角度、半径を記憶装置１０２に格納する円筒データ格納部である。
【００２８】
図２は、図１の円筒データ抽出処理を実施する計算機を示すものであり、２００は、円筒データ抽出部１０３を実行する計算機である。計算機２００において、２０１は通信網で接続された他の計算機とデータのやり取りをおこなうための通信装置、２０２はマウスやキイボードなどのデータの入力をおこなう入力装置、２０３はグラフィックディスプレイなどの計算結果を文字やグラフなどで表示するための表示装置、２０４は計算結果を蓄積するＦＤ等の記憶装置（記憶媒体）、２０５はメモリ、２０６は１０３の円筒データ抽出処理をおこなう中央処理装置である。
【００２９】
図３は、本発明の円筒データ抽出処理の流れを示すフローチャートである。以下、図３のフローチャートと数式を交えながら、本発明の一実施例を述べる。尚、図３は、本発明のプログラムに相当する。
【００３０】
最初に計測点・収束判定値・角度ピッチ読取部１０４が、記憶装置１０２から計測点データ、収束判定値、角度ピッチを読み取る。計測点データは、（数５）に示すようにＰｉでｉ番目の計測点を表し、（ｘ成分、ｙ成分、ｚ成分）を持つ３次元空間上の点である。ここで、ｉは１から計測点数Ｎまで変化する。
【００３１】
【数５】
【００３２】
収束判定値は、例えば１．０Ｅ−２０のように更新量を０とみなす値である。角度ピッチは、０から３６０°までの角度を離散化するときの幅であり、例えば２０°とすれば、０°、２０°、・・・、３４０°、３６０°の角度に離散化する。
【００３３】
次に、球中心算出部１０５がＮ個の計測点と球表面との距離の二乗の総和が最小になる球（以下、最小二乗球という）の中心点を算出する。以下、その算出法について述べる。
【００３４】
最小二乗球の中心Ｏｓは、例えば（数３）で算出することができるが、本実施例では、初期角度算出部１０６、角度精度向上部１０７、円筒中心・角度算出部１０８を説明する都合上、（数６）に示す方法を用いる。
【００３５】
【数６】
【００３６】
（数６）において、Ｏｓは球の中心、ｇｉは点Ｐｉと中心Ｏｓの距離の二乗（半径の二乗）、ｇａｖｒはｇｉの平均である。即ち、各点Ｐｉと中心Ｏｓの半径の二乗と、当該半径の二乗の平均との差の二乗和が最小になるよう中心Ｏｓを決定する。
【００３７】
（数６）は厳密に言えば、計測点と球表面の距離ではなく、距離の二乗を最小化する。しかしながら、距離の場合は、平方根を用いた関数とせねばならず、微分時に複雑さが増してしまう。そのため、距離の二乗を最小化している。この方法は一般的に最小二乗法を用いるときの常套手段である。尚、（数３）と（数６）の解が一致することは実験で確認済である。
【００３８】
また、（数６）のようにｇｉとｇａｖｒの差を作成すると、ｇｉ中に同一値を持つ変数成分及び定数成分がある場合に０となるので、数値が大きくなることがないため、計算機処理する場合の丸め誤差を少なくする上で都合が良い。
【００３９】
さらに、ｇａｖｒが加算と定数による除算なので、ｇｉとｇａｖｒの次数が同一となり、次数が増加しないので関数の複雑さが増すことがない。
【００４０】
（数６）を解くために、（数６）を球の中心Ｏｓの各成分（Ｏｓｘ、Ｏｓｙ、Ｏｓｚ）で偏微分した関数が０になる式を作成し、当該式を１次までＴａｙｌｏｒ展開して連立一次方程式を作成する。当該連立一次方程式を、初期中心Ｏｓを与えて解くことにより、中心更新量を算出する。本方法は、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法と呼ばれている。
【００４１】
当該連立一次方程式は、２回微分が０または定数となる比較的単純な関数であるため、任意の初期中心でも収束する。
【００４２】
（数７）は当該連立一次方程式である。
【００４３】
【数７】
【００４４】
（数３）が４行４列であるのに対し、（数７）は３行３列である。このように（数６）で示すようにｇｉとｇａｖｒの差を求める方法を用いると、行数および列数を１行１列減少させることができるので、計算時間の短縮が可能なため、本実施例ではこの方法を採用している。
【００４５】
（数７）において各行列成分Ｍｓは（数８）である。
【００４６】
【数８】
【００４７】
（数８）においてｇｉ及びｇａｖｒは（数６）、ｇｉの中心Ｏｓの成分による微分及び平均は、（数９）及び（数１０）で表せる。
【００４８】
【数９】
【００４９】
【数１０】
【００５０】
（数１０）に示すように、ｇｉの中心Ｏｓの成分による２回偏微分は定数２か０である、従って、平均も定数２か０である。そのため、（数８）は（数１１）のように簡単化される。
【００５１】
【数１１】
【００５２】
また、（数１１）で使用するｇｉの中心Ｏｓの成分による２回偏微分及び平均は、（数９）のみとなる。
【００５３】
以上、数式により（数６）の解き方を説明したが、次に、図３に示す球中心算出部１０５のフローチャートを用いて球中心の算出過程を説明する。
【００５４】
図３の３０１から３１０は、（数６）を解く過程の一実施例を示すフローチャートである。最初に球の中心点の初期値を例えば（０，０，０）のように与える（ステップ３０１）。当該初期値は（０，０，０）に限定されるものではなく、任意の値、例えば、（１００，１００，１００）や（−１，５００，０）などでもかまわないが、数値計算上は予想できる中心点を与えることが望ましい。本実施例では、予想できないこと、及び更新量がそのまま球の中心点となることから（０，０，０）を与えている。
【００５５】
次に、計測点Ｐｉと中心点の差ベクトルを求める（ステップ３０２）。
【００５６】
次に、（数６）のｇｉに示す差ベクトルの内積を求め、計測点数分の内積値を作成する（ステップ３０３）。当該内積値は中心点と計測点Ｐｉの距離の二乗であることから、内積値を二乗距離と呼ぶ。
【００５７】
次に、（数９）の（ｇｉ，ｏｓ，ｘ）、（ｇｉ，ｏｓ，ｙ）、（ｇｉ，ｏｓ，ｚ）に示す二乗距離の中心点による偏微分を作成する（ステップ３０４）。
【００５８】
次に、（数６）のｇａｖｒ、（数９）の（ｇａｖｒ，ｏｓ，ｘ）、（ｇａｖｒ，ｏｓ，ｙ）、（ｇａｖｒ，ｏｓ,ｚ）に示す二乗距離および二乗距離の中心点による偏微分の平均を作成する（ステップ３０５）。
【００５９】
次に、ステップ３０４で求めた偏微分およびステップ３０５で求めた偏微分の平均から、（数１１）の（Ｍｓ，０，０）、（Ｍｓ，２，０）、（Ｍｓ，０，１）、（Ｍｓ，２，１）、（Ｍｓ，０，２）、（Ｍｓ，２，２）に示す二乗距離の総和の２回偏微分を作成するとともに、ステップ３０３で求めた二乗距離とステップ３０４で求めた偏微分およびステップ３０５で求めた二乗距離の平均と二乗距離の偏微分の平均から（数１１）の（Ｍｓ，３，０）、（Ｍｓ，３，１）、（Ｍｓ，３，２）に示す二乗距離の総和の偏微分を作成する（ステップ３０６）。
【００６０】
次に、ステップ３０６で求めた２回偏微分と偏微分から（数７）に示す連立一次方程式の係数を作成する（ステップ３０７）。
【００６１】
次に、連立一次方程式を解き、中心点の更新量を求める（ステップ３０８）。連立一次方程式の解法は、例えば、掃き出し法やＬＵ分解で解くことができる。また、３行３列なので解析的に解いてもかまわない。
【００６２】
次に、中心点の更新量が１０４で読み込んだ０とみなせる値になっているかを判定する（ステップ３０９）。判定時には、例えば、中心点の更新量（δＯｓｘ，δＯｓｙ，δＯｓｚ）の絶対値を作成して，当該絶対値がそれぞれ０とみなせる値になるかどうかで確認する。そのほかに、更新量の内積（δＯｓｘ×δＯｓｘ＋δＯｓｙ×δＯｓｙ＋δＯｓｚ×δＯｓｚ）を作成して０となるかどうかを判定してもかまわない。本実施例では、更新量の内積で判定している。
【００６３】
ステップ３０９で０に近くないと判定された場合、中心点の更新を行う（ステップ３１０）。中心点の更新は、現在の中心点（Ｏｓｘ，Ｏｓｙ，Ｏｓｚ）に、更新量（δＯｓｘ，δＯｓｙ，δＯｓｚ）を加算すればよい。この加算した値を現在の中心点（Ｏｓｘ，Ｏｓｙ，Ｏｓｚ）とし、ステップ３０２から繰り返す。ステップ３０９で０に近いと判定された場合、球中心算出部１０５の処理を終了する。
【００６４】
以上のステップ３０２からステップ３０９を繰り返すことで、球の中心点を定めることができる。
【００６５】
実際には、球中心算出部１０５の処理は、殆ど１回で繰り返し計算が終了する。これは、ｇｉの中心Ｏｓの成分による偏微分が、１回偏微分が変数、２回の偏微分および相互偏微分が定数または０、３回以降の偏微分および相互偏微分が０であることから、（数７）に示す関数が線形関数であることに関連している。但し、数値計算誤差が含まれることも想定されるので、本実施例に示したように、ステップ３０２〜ステップ３０９、ステップ３１０の処理を行うことが精度を高める上で望ましい。
【００６６】
３次元で定義される点に球をあてはめる場合、（数７）に示すように、中心Ｏｓのみを求めることとなり回転角度には依存しない。このため、特異な場合を除き確実に中心点Ｏｓを定めることができるので、あてはめ時間を最小限に抑制することができる。
【００６７】
ここで、特異な場合とは、例えば、３次元で定義された点が直線上に配置された場合である。本状態は、計測物が円筒面の場合に限定すると、円筒の中心軸に平行となるよう計測したときに生じる。このとき、（数７）で繰り返し計算をおこなうと無限に繰り返すか、または０除算が発生してしまう。そのため、繰り返し回数に上限を設け、上限を超えた場合、及び（数７）の連立一次方程式を解くときに０除算が発生するか確認し、エラー出力をおこなうことが望ましい。
【００６８】
当該エラーが出力された場合、計測位置を変更する必要がある。そのため、（数７）の処理を行う計算機、及び計算の結果をモニタできる装置、例えばブザーや小型の表示装置を設け、計測結果をリアルタイムにモニタリングすることが望ましい。本構成を採用する場合、（数７）の計算処理が単純であり、繰り返し回数も少ないことから計算能力が劣っている計算機でも採用することができるとともに、計測位置の誤りがリアルタイムで手元確認できるため、すぐに計測位置を修正することができるので計測間違いが少なくなる。
【００６９】
次に、初期角度算出部１０６で、円筒の中心軸となる角度の初期値を算出する。初期角度を算出するために（数１２）を用いる。
【００７０】
【数１２】
【００７１】
（数１２）において、Ｑｉは、計測点ＰｉをＹ軸で角度Ｂだけ回転してからＸ軸で角度Ａだけ回転した後の（ｘ、ｙ）成分である。Ｏは、球の中心ＯｓをＹ軸で角度Ｂだけ回転してからＸ軸で角度Ａだけ回転した後の（ｘ、ｙ）成分である。
【００７２】
（数１２）では、投影平面を（ｘ、ｙ）平面に限定し、回転をＹ軸回転後にＸ軸回転しているが、投影平面は（ｘ、ｙ）平面に限定されるものではなく、（ｙ、ｚ）平面、（ｚ、ｘ）平面でもかまわないし、回転は、Ｙ軸回転後のＺ軸回転、Ｘ軸回転後のＹ軸回転、Ｘ軸回転後のＺ軸回転、Ｚ軸回転後のＸ軸回転、Ｚ軸回転後のＹ軸回転のいずれでもかまわない。
【００７３】
重要なことは、初期角度算出部１０６、角度精度向上部１０７、円筒中心・角度算出部１０８で回転順序を一致させておくことである。これを誤ると角度の高精度化、円筒の中心や角度の算出ができなくなるので注意が必要である。
【００７４】
（数１２）は三角関数が入るため、角度ＡおよびＢで無限に偏微分、相互偏微分することが可能であり極値を複数持つ複雑な関数である。そのため、本願では（数１３）で示すように、（数１２）の角度（Ａ，Ｂ）を０〜３６０°の中で離散化して（数１２）の値を求め、最小値を持つ角度（Ａ，Ｂ）を探す構成とした。
【００７５】
【数１３】
【００７６】
以上、数式により初期角度算出部１０６の処理を説明したが、次に図３に示す初期角度算出部１０６のフローチャートで初期角度の算出課程を説明する。図３において、３１１から３１６は、初期角度算出部１０６の処理の一実施例を示すフローチャートである。
【００７７】
最初に（数１３）に示す角度（Ａ，Ｂ）を、計測点・収束判定値・角度ピッチ読取１０４で読み込んだ角度ピッチに従って離散化する（ステップ３１１）。本実施例では、角度ピッチを２０°としている。２０°という値は、テストデータを用いたとき角度精度向上部１０６で収束するピッチを実験して決定したものである。本結果については後述する。
【００７８】
尚、初期角度算出部１０６において、円筒中心も変数とした場合、円筒中心が３成分毎に−∞〜∞まで変化するため設定が難しい。また、変数が５個となり組合せ数が一気に増大するため計算量が増大してしまう。更に、変数５個それぞれに極値を持つため大変複雑になってしまう。そのため、球中心算出部１０５を用意して角度だけで（数１３）を算出する構成としている。
【００７９】
次に、計測点Ｐｉ、球中心算出手段１０５で作成した球中心を離散化した角度毎に座標変換する（ステップ３１２）。座標変換の方法は、上述した“（４）３次元から２次元への座標変換”で説明した処理でかまわない。
【００８０】
次に、（数１２）のｆｉに示すように座標変換した計測点と球の中心の（ｘ、ｙ）成分を取り出して差ベクトルを作成し、差ベクトルの内積を作成して（ｘ、ｙ）成分の二乗距離を作成する（ステップ３１３）。
【００８１】
次に、（数１２）のｆａｖｒに示すように当該二乗距離の平均を作成する（ステップ３１４）。
【００８２】
次に、二乗距離と平均の差を作成し、当該差の二乗和を作成する（ステップ３１５）。
【００８３】
最後に、当該二乗和の最低値を取出し、最低値と一致する離散化した角度（Ａ，Ｂ）を抽出する（ステップ３１６）。
【００８４】
最低値を作成する方法としては、当該二乗和を小さい値から昇順に並べ替え最初の値を抽出する方法や、最初の値を最低値の初期値とし２番目から計測点数番目までの間で当該最低値より小さい値があらわれたとき置き換えるなどの方法を用いてもかまわない。並べ替えの方法としては、バブルソートやクイックソートなどのソーティングアルゴリズムを利用してもかまわない。
【００８５】
ここで、最低値を持つ角度（Ａ，Ｂ）が複数存在することを考慮し、ステップ３１５の二乗和は小数点以下の有効桁数を考慮して整数化しておくことが望ましい。整数化のかわりに四捨五入した実数を使うことも考えられるが、実数値に丸め誤差が累積して微小に異なる場合があるため、整数化しておく方が良い。
【００８６】
本実施例では、小数点以下１桁までを有効とし、実数値を１０倍して切り捨てる方法を使用している。
【００８７】
図４（ａ）は、計測点８点を与え、角度Ａおよび角度Ｂをそれぞれ０から３６０°まで２０°ピッチで分割して（数１３）のＦ（Ａ，Ｂ）をプロットした３次元グラフの上面図、図４（ｂ）は、図４（ａ）の斜視図、図４（ｃ）は、角度Ａを横軸、角度Ｂを縦軸としＦ（Ａ，Ｂ）の値を記入した図である。図４（ｃ）では、最低値を四角で囲んで示している。
【００８８】
図４（ａ）、図４（ｂ）、図４（ｃ）より、Ｆ（Ａ，Ｂ）が極値を多数持つ複雑な関数であり、最低値が４個程度存在することがわかる。
【００８９】
図５（ａ）の１〜４は、球の中心点を複数求めた角度で座標変換した位置の（ｘ、ｙ）座標である。図５（ａ）に示すように、中心点が角度毎に象限を分けて算出されている。
【００９０】
このように、中心点が分れて検出されることから、角度と中心を変数として繰り返し計算すると、中心の更新量に従い角度が大きく変化するため解が大きく変動することに繋がる。そのため、初期角度算出部１０６では球中心算出部１０５で求めた中心を円筒の中心とみなし、角度の初期値を定めているので、解が大きく変動するのを防ぐことができる。
【００９１】
次に初期角度算出部１０６を実現する別の処理について、図７及び図８（ａ）、図８（ｂ）、図８（ｃ）を用いて説明する。
【００９２】
図７は、初期角度算出部１０６の別の実施方法を示すフローチャートである。該図に示す如く、最初に分割数、最大Ｂ、しきい値を設定する（ステップ９０１）。分割数は２つの角度（Ａ，Ｂ）を離散化するときの点数、最大Ｂは分割数を変数としたときの当該分割数における関数Ｆ（Ａ，Ｂ）と当該関数の近似曲面の距離の二乗の最大を示す関数の最大の初期値、しきい値は分割数の決定処理が終了したかどうかを表す値である。尚、図７は、本発明の別の実施方法のプログラムに相当する。
【００９３】
本実施例では、分割数は３を設定しているが、関数Ｆ（Ａ，Ｂ）の複雑さに応じ、３より大きい値を設定してもかまわない。また、最大Ｂは０を設定しているが、関数Ｆ（Ａ，Ｂ）の最小が０となるからであり、当該最小が０以上となることが確実な場合、０より大きい値を設定してもかまわない。また、しきい値は０．１を設定しているが、角度精度向上部１０７の処理が発散する場合０〜０．１の値を設定する必要がある。
【００９４】
更に、角度精度向上部１０７の処理が確実に収束するなら、０．１〜１の範囲に設定してもかまわない。本実施例でいろいろな計測点を与え、実験したところ、発散しなかったことから０．１に設定している。
【００９５】
次に、角度Ａの範囲０〜３６０°、角度Ｂの範囲０〜３６０°を当該分割数個に離散化する（ステップ９０２）。本実施例では、関数Ｆ（Ａ，Ｂ）の傾向が不明なため、角度が等感覚になるような離散化をおこなっている。例えば、分割数が３なら、（０，０）、（１８０，０）、（３６０，０）、（０，１８０）、（１８０，１８０）、（３６０，１８０）、（０，３６０）、（１８０，３６０）、（３６０，３６０）の格子リストとなる。
【００９６】
次に、格子リストの角度毎に関数Ｆ（Ａ，Ｂ）を作成する（ステップ９０３）。例えば、分割数が３なら、Ｆ（０，０）、Ｆ（１８０，０）、Ｆ（３６０，０）、Ｆ（０，１８０）、Ｆ（１８０，１８０）、Ｆ（３６０，１８０）、Ｆ（０，３６０）、Ｆ（１８０，３６０）、（Ｆ３６０，３６０）を作成する。９０３の処理過程は、図２のステップ３１２、３１３、３１４、３１５と同様である。
【００９７】
次に、格子リストの角度毎に関数Ｆ（Ａｉ−１，Ｂｊ−１）、Ｆ（Ａｉ＋１，Ｂｊ−１）、Ｆ（Ａｉ−１，Ｂｊ＋１）、Ｆ（Ａｉ＋１，Ｂｊ＋１）から近似曲面Ｆａｉ，ｊ（ｕ，ｖ）を作成し、Ｆ（Ａｉ，Ｂｊ）とＦａｉ，ｊ（０.５，０.５）の距離を作成し、当該距離の最大である最大Ａを求める（ステップ９０４）。近似曲面Ｆａｉ，ｊ（ｕ，ｖ）は、ｍ×ｎ次の多項式曲面、Ｂｅｚｉｅｒ曲面、Ｂ−Ｓｐｌｉｎｅ曲面、Ｃｏｏｎｚ曲面など任意の形式でかまわない。また、これらの曲面を有理式化したものでもかまわない。
【００９８】
本実施例では、処理を簡単にするため、双１次パッチを採用し、双１次パッチの（ｕ，ｖ）方向それぞれの中央で評価する構成としている。
【００９９】
（数１４）に双１次パッチによる近似曲面Ｆａｉ，ｊ（ｕ，ｖ）を示す。
【０１００】
【数１４】
【０１０１】
（数１４）において、Ｆａｉ，ｊ（ｕ，ｖ）は、パラメータ（ｕ，ｖ）を変数とする双１次曲面（双１次Ｂｅｚｉｅｒ曲面でもある）。パラメータ（ｕ，ｖ）はそれぞれ、０〜１まで変化する変数。Ｆ（Ａｉ−１，Ｂｊ−１）、Ｆ（Ａｉ＋１，Ｂｊ−１）、Ｆ（Ａｉ−１，Ｂｊ＋１）、Ｆ（Ａｉ＋１，Ｂｊ＋１）は双１次曲面の４隅の点（双１次Ｂｅｚｉｅｒ曲面の制御点）であり、関数Ｆ（Ａｉ，Ｂｊ）における１つおきの値である。
【０１０２】
近似評価に用いる点は（数１５）で表される。
【０１０３】
【数１５】
【０１０４】
（数１５）において、ｇｉ，ｊは５成分を持つ行列で、Ｆ（Ａｉ−１，Ｂｊ−１）、Ｆ（Ａｉ＋１，Ｂｊ−１）、Ｆ（Ａｉ−１，Ｂｊ＋１）、Ｆ（Ａｉ＋１，Ｂｊ＋１）と、パラメータ（０.５，０）、（０，０.５）、（０.５，０.５）、（１，０.５）、（０.５，１）の位置を持つ双１次曲面上の点との差を成分に持つ。
【０１０５】
行列ｇｉ，ｊの内積を（数１６）で表す。
【０１０６】
【数１６】
【０１０７】
（数１５）は、厳密にいえば，Ｆ（Ａｉ，Ｂｊ−１）とＦａｉ，ｊ（０.５，０）、Ｆ（Ａｉ−１，Ｂｊ）とＦａｉ，ｊ（０，０.５）、Ｆ（Ａｉ，Ｂｊ）とＦａｉ，ｊ（０.５，０.５）、Ｆ（Ａｉ＋１，Ｂｊ）とＦａｉ，ｊ（１.０，０.５）、Ｆ（Ａｉ，Ｂｊ＋１）とＦａｉ，ｊ（０.５，１.０）が垂点（距離）ではない。垂点は、例えば、Ｆ（Ａｉ，Ｂｊ）と、Ｆａｉ，ｊ（ｕ，ｖ）、Ｆａｉ，ｊ（ｕ，ｖ）のパラメータ（ｕ，ｖ）による偏微分で定義される値で、非線形の式であり、解くには初期値を与えて収束計算することが必要である。そのため、複雑な計算が必要である。
【０１０８】
本発明では、分割数増加による近似評価のために使用することを前提として（数１６）を距離の二乗としている。（数１６）のｇｉ，ｊの内積（ｉ，ｊ）を（０，０）から（分割数―１，分割数―１）まで１つおきに変化させ、当該内積の最大を最大Ａとしている。
【０１０９】
次に、最大Ａが最大Ｂより大きいか判定する（ステップ９０５）。ステップ９０５の判定結果が大きければ、最大Ｂに最大Ａを代入する更新をおこなう（ステップ９０７）。
【０１１０】
次に、分割数×２−１の値を分割数に代入し（ステップ９０８）、ステップ９０２から繰り返す。本処理は、最大Ａの関数値が分割数を増大するのに従い、一旦上昇してから下がっていく傾向を示すため、当該関数値の最大を基準とした近似度合いの判定をおこなうために用いている。
【０１１１】
このようにすることで、最初の近似度がよい場合に誤って処理が終了することを防ぐことができる。ステップ９０５の判定結果が小さければ、最大Ａを最大Ｂで除した値がしきい値より大きいかを判定する（ステップ９０６）。ステップ９０６の判定結果が大きいときは、ステップ９０５で分割数を増加させ、ステップ９０２から繰り返す。ステップ９０６の判定結果が小さければ、分割数が決定したので、当該分割数でＦ（Ａ，Ｂ）が最低となる角度を抽出する（ステップ３１６）。ステップ３１６は図３のステップ３１６と同様の処理をおこなう。
【０１１２】
以上のように、ステップ９０１からステップ３１６の処理をおこなうことにより、角度（Ａ，Ｂ）の初期値を決定する。
【０１１３】
図８（ａ）は、分割数が増加する様子を示したイメージ図である。最初に３×３のＦ（Ａ，Ｂ）を定め、次に５×５、次に９×９のように間に１個づつＦ（Ａ，Ｂ）を増加させていく。従って、ｉ番目の分割数は２Ｎｉ−１（一つ前の分割数）−１となる。
【０１１４】
このようにＦ（Ａ，Ｂ）を増加させていくため、分割数は、３、５、９、１７、３３、６５、１２９、２５７、５１３、１０２５、・・・と増えていく。角度の範囲は０〜３６０なので、１０回目で約０．３５°の間隔となる。
【０１１５】
分割数が増加すると計算時間が増大するため、角度間隔と角度精度向上部１０７の収束度合いとの関連を実験したところ、５から６回繰り返したときに安定したので、繰り返しの回数は余裕をみて１０回に抑えている。
【０１１６】
図８（ｂ）は、Ｆ（Ａ，Ｂ）とＦａｉ，ｊ（ｕ，ｖ）の評価位置の関係を示している。Ｆ（Ａ，Ｂ）とＦａｉ，ｊ（ｕ，ｖ）は、（数１４）の定義からもわかるように、４隅の点が一致する。そのため、評価位置を中間位置にしている。
【０１１７】
図８（ｃ）は、分割数の増加に従い、最大Ａ／最大Ｂがどのように変化するかを示した特性図である。
【０１１８】
図８（ｃ）からわかるように、評価値は分割数が５か９で最大に達し、以後０に近づいていく。この結果より、０．１（１０％）をしきい値の初期値としている。しきい値は、最大値からの分散を定義し、その分散範囲以下になる値をしきい値としてもかまわない。この場合、３σや６σなどの位置をしきい値としてもかまわない。
【０１１９】
次に、角度精度向上部１０７で複数の初期角度をそれぞれ高精度化する。初期角度を高精度化するために（数１２）を角度（Ａ，Ｂ）で偏微分した値が０になる位置、すなわち（数１２）の極値を算出する。しかしながら、（数１２）は三角関数で記述された式であることから、角度（Ａ，Ｂ）での２回以上の偏微分、相互偏微分が無限に作成できる。そのため、図４（ａ）、図４（ｂ）、図４（ｃ）に示したように、極値を複数持つ複雑な関数となる。
【０１２０】
そこで、（数１２）を角度（Ａ，Ｂ）で偏微分した値が０になる式を１次Ｔａｙｌｏｒ展開して０とする式に近似し、（数１７）に示す連立一次方程式を作成する。
【０１２１】
（数１７）は初期角度（Ａ，Ｂ）を与えて解くことにより角度更新量（δＡ，δＢ）を算出する、所謂Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法である。（数１７）は複雑な関数であるため、良い初期角度を設定しないと解が発散してしまう。そのため。初期角度算出部１０６で初期角度を定めている。
【０１２２】
【数１７】
【０１２３】
（数１７）において各行列成分Ｍａは（数１８）である。
【０１２４】
【数１８】
【０１２５】
（数１８）においてｆｉ及びｆｉの平均、ｆｉ及びｆｉの角度（Ａ，Ｂ）による偏微分と当該偏微分の平均は、（数１９）で表せる。
【０１２６】
【数１９】
【０１２７】
（数１９）において、投影した計測点Ｑｉ及び投影した球の中心Ｏ、投影した計測点Ｑｉの角度（Ａ，Ｂ）による偏微分は（数２０）で表せる。
【０１２８】
【数２０】
【０１２９】
（数２０）において、計測点Ｐｉを投影するためのｘ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２１）で表せる。
【０１３０】
【数２１】
【０１３１】
（数２０）において、計測点Ｐｉを投影するためのｙ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２２）で表せる。
【０１３２】
【数２２】
【０１３３】
（数２１）に示すように、計測点Ｐｉを投影するためのｘ軸方向方向ベクトルの角度（Ａ）による偏微分、相互偏微分、２回以上の偏微分はすべて０ベクトルである。そのため、（数２０）は（数２３）に示すように簡単化される。また、（数２１）は（数２４）のように簡単化される。
【０１３４】
【数２３】
【０１３５】
【数２４】
【０１３６】
以上、数式により（数１７）の解き方を説明したが、次に、図３に示す角度精度向上部１０７のフローチャートを用いて投影したときに円となる角度（Ａ，Ｂ）を高精度化する過程を説明する。
【０１３７】
図３のステップ３２１からステップ３３０は、（数１７）を解く過程の一実施例を示すフローチャートである。尚、初期角度算出部１０６で求めた角度は複数存在するので、角度精度向上部１０７を角度の個数分繰り返して、当該角度毎に高精度化する。
【０１３８】
最初に初期角度算出部１０６で求めた角度を角度初期値に設定する（ステップ３２１）。
【０１３９】
次に、角度初期値に従い（数２４）で投影面のｘ軸方向のベクトルｔ、当該ベクトルｔの角度Ｂによる偏微分値ｔｂと２回偏微分値ｔｂｂ、（数２２）で投影面のｙ軸方向のベクトルｕ、当該ベクトルｕの角度Ａによる偏微分値ｕａと２回偏微分値ｕａａ、角度Ａで微分してから角度Ｂで微分する相互偏微分値ｕａｂ、当該ベクトルｕの角度Ｂによる偏微分値ｕｂ、２回偏微分値ｕｂｂを求める。次に、（数２３）で計測点Ｐｉを投影したベクトルＱｉ、当該ベクトルの角度Ａによる偏微分Ｑｉ，ａ、２回偏微分Ｑｉ，ａａ、角度Ａで微分後に角度Ｂで微分する相互偏微分Ｑｉ，ａｂ、当該ベクトルの角度Ｂによる偏微分値Ｑｉ，ｂ、２回偏微分値Ｑｉ，ｂｂ、を求める（ステップ３２２）。
【０１４０】
次に、（数１９）のｆｉにより投影した計測点と中心点の差ベクトルの内積で二乗距離を作成する（ステップ３２３）。
【０１４１】
次に、（数１９）のｆｉ，ａにより当該二乗距離の角度Ａによる偏微分、ｆｉ，ａａにより角度Ａによる２回偏微分、ｆｉ，ａｂにより角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分、ｆｉ，ｂにより当該二乗距離の角度Ｂによる偏微分、ｆｉ，ｂｂにより角度Ｂによる２回偏微分を作成する（ステップ３２４）。
【０１４２】
次に、（数１９）のｆａｖｒにより二乗距離の平均、ｆａｖｒ，ａにより二乗距離の平均の角度Ａによる偏微分、ｆａｖｒ，ａａにより二乗距離の平均の角度Ａによる２回偏微分、ｆａｖｒ，ａｂにより二乗距離の平均の角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分、ｆａｖｒ，ｂにより当該二乗距離の平均の角度Ｂによる偏微分、ｆａｖｒ，ｂｂにより二乗距離の平均の角度Ｂによる２回偏微分を作成する（ステップ３２５）。
【０１４３】
次に、（数１８）により二乗距離と二乗距離の平均の差の二乗和、当該二乗和の角度Ａによる偏微分、角度Ａによる２回偏微分、角度Ａで微分してから角度Ｂで微分する相互偏微分、角度Ｂによる偏微分、角度Ｂによる２回偏微分を作成する（ステップ３２６）。
【０１４４】
次に、ステップ３２６で求めた値を連立一次方程式の係数に代入する（ステップ３２７）。
【０１４５】
次に、当該連立一次方程式を解き、更新量（δＡ，δＢ）を算出する（ステップ３２８）。当該連立一次方程式を解くには一般的な連立一次方程式の解法、例えば、掃き出し法やＬＵ分解を使えば良い。また、２行２列なので解析的に解いてもかまわない。
【０１４６】
次に、更新量（δＡ，δＢ）が０とみなせるか判定する（ステップ３２９）。判定方法は、更新量の成分（δＡ，δＢ）の絶対値がそれぞれ０とみなせるかを判定する方法のほかに、更新量の内積（δＡ×δＡ＋δＢ×δＢ）が０とみなせるかを判定してもかまわない。本実施例では、更新量の内積が０とみなせるかを判定している。
【０１４７】
更新量が０とみなせない場合、角度更新をおこなう（ステップ３３０）。角度更新は、角度初期値（Ａ，Ｂ）に更新量（δＡ，δＢ）を加算することでおこなう。更新した角度初期値（Ａ，Ｂ）でステップ３２２から繰り返す。更新量が０とみなせる場合、処理を終了する。
【０１４８】
以上のステップ３２１からステップ３３０を繰り返すことで、球の中心を円筒中心とした場合の座標変換角度（Ａ，Ｂ）を定めることができる。実際に計測点を与え実験してみたところ、初期角度として初期角度算出部１０６で算出した角度を与えたとき、５から６回で収束した。
【０１４９】
しかしながら、初期角度として任意の値、例えば（０°、０°）などを指定した場合、図４で示すグラフの（０°、０°）近傍の極値が定まってしまい、極大値が解になってしまう現象が確認されており、初期角度算出部１０６が有効であることが確認されている。尚、特異な場合は球中心算出部１０４で除外すること、計測物が円筒面に限定されていることから角度精度向上部１０７が有効に動作することを確認している。
【０１５０】
次に、円筒中心・角度算出部１０８で精度向上した角度と円筒中心を、更に高精度化する手法を説明する。
【０１５１】
当該角度と円筒中心を高精度化するために（数１２）を角度（Ａ，Ｂ）および中心（Ｏｘ，Ｏｙ）で偏微分した値が０になる位置、即ち、（数１２）の極値を算出する。しかしながら、（数１２）は、三角関数で記述された式であることから、極値を複数持つ複雑な関数となる。
【０１５２】
そこで、（数１２）を角度（Ａ，Ｂ）及び円筒中心（Ｏｘ，Ｏｙ）で偏微分した値が０になる式を１次Ｔａｙｌｏｒ展開して０とする式に近似し、（数２５）に示す連立一次方程式を作成する。
【０１５３】
（数２５）は、初期角度（Ａ，Ｂ）および初期円筒中心（Ｏｘ，Ｏｙ）を与えて解くことにより更新量（δＡ，δＢ，δＯｘ，δＯｙ）を算出する。いわゆるＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法である。（数２５）は複雑な関数であるため、良い初期角度を設定しないと解が発散してしまう。そのため、角度精度向上部で高精度化した初期角度、当該初期角度で球の中心を座標変換した初期円筒中心を定めている。
【０１５４】
【数２５】
【０１５５】
（数２５）において各行列成分Ｍｂは（数２６）である。
【０１５６】
【数２６】
【０１５７】
（数２６）においてｆｉ及びｆｉの平均、ｆｉ及びｆｉの角度（Ａ，Ｂ）による偏微分、ｆｉ及びｆｉの円筒中心（Ｏｘ，Ｏｙ）による偏微分、と当該偏微分の平均は、（数２７）で表せる。
【０１５８】
【数２７】
【０１５９】
（数２７）において、投影した計測点Ｑｉおよび投影した球の中心Ｏ、投影した計測点Ｑｉの角度（Ａ，Ｂ）は（数２３）で表せる。（数２３）において偏微分した値が（数２７）に反映されているもの、及び偏微分または相互偏微分して０ベクトルになるものは記述していない。
【０１６０】
（数２７）において、計測点Ｐｉを投影するためのｘ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２４）で表せる。（数２４）において偏微分または相互偏微分して０ベクトルになるものは記述していない。また、（数２７）において、計測点Ｐｉを投影するためのｙ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２２）で表せる。（数２２）において偏微分または相互偏微分して０ベクトルになるものは記述していない。
【０１６１】
（数２７）に示すように、円筒中心（Ｏｘ、Ｏｙ）による相互偏微分、２回上の偏微分はすべて０か２である。そのため、（数２６）は（数２８）、（数２７）は（数２９）に示すように簡単化される。
【０１６２】
【数２８】
【０１６３】
【数２９】
【０１６４】
以上、数式により（数２５）の解き方を説明したが、次に、図３に示す円筒中心・角度算出部１０８のフローチャートを用いて投影したときに円となる角度（Ａ，Ｂ）、中心点（Ｏｘ，Ｏｙ）を高精度化する過程を説明する。
【０１６５】
図３のステップ３３１からステップ３４０は（数２５）を解く過程の一実施例を示すフローチャートである。なお、角度精度向上部１０７で精度向上した角度は複数存在するので、円筒中心・角度算出部１０８を角度の個数分繰り返して、当該角度毎に角度と中心点を高精度化する。
【０１６６】
最初に、角度精度向上部１０７で求めた角度を角度初期値に設定するとともに、球の中心点を当該角度初期値で座標変換した値の（ｘ、ｙ）成分を中心点初期値に設定する（ステップ３３１）。
【０１６７】
次に、角度初期値と中心点初期値に従い（数２４）で投影面のｘ軸方向のベクトルｔ、当該ベクトルｔの角度Ｂによる偏微分値ｔｂと２回偏微分値ｔｂｂ、（数２２）で投影面のｙ軸方向のベクトルｕ、当該ベクトルｕの角度Ａによる偏微分値ｕａと２回偏微分値ｕａａ、角度Ａで微分してから角度Ｂで微分する相互偏微分値ｕａｂ、当該ベクトルｕの角度Ｂによる偏微分値ｕｂ、２回偏微分値ｕｂｂを求める。次に、（数２３）で計測点Ｐｉを投影したベクトルＱｉ、当該ベクトルの角度Ａによる偏微分Ｑｉ，ａ、２回偏微分Ｑｉ，ａａ、角度Ａで微分後に角度Ｂで微分する相互偏微分Ｑｉ，ａｂ、当該ベクトルの角度Ｂによる偏微分値Ｑｉ，ｂ、２回偏微分値Ｑｉ，ｂｂ、を求める（ステップ３３２）。
【０１６８】
次に、（数２９）のｆｉにより計測点と中心点の差ベクトルの内積で二乗距離を作成する（ステップ３３３）。
【０１６９】
次に、（数２９）のｆｉ，ａにより当該二乗距離の角度Ａによる偏微分、ｆｉ，ａａにより角度Ａによる２回偏微分、ｆｉ，ａｂにより角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分、ｆｉ，ａｏｘにより角度Ａで偏微分してから中心Ｏｘで偏微分する相互偏微分、ｆｉ，ａｏｙにより角度Ａで偏微分してから中心Ｏｙで偏微分する相互偏微分、ｆｉ，ｂにより当該二乗距離の角度Ｂによる偏微分、ｆｉ，ｂｂにより角度Ｂによる２回偏微分、ｆｉ，ｂｏｘにより角度Ｂで偏微分してから中心Ｏｘで偏微分する相互偏微分、ｆｉ，ｂｏｙにより角度Ｂで偏微分してから中心Ｏｙで偏微分する相互偏微分、ｆｉ，ｏｘにより中心Ｏｘでの偏微分、ｆｉ，ｏｙにより中心Ｏｙでの偏微分を作成する（ステップ３３４）。
【０１７０】
次に、（数２９）のｆａｖｒ，ａにより当該二乗距離の角度Ａによる偏微分の平均、ｆａｖｒ，ａａにより角度Ａによる２回偏微分の平均、ｆａｖｒ，ａｂにより角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分の平均、ｆａｖｒ，ａｏｘにより角度Ａで偏微分してから中心Ｏｘで偏微分する相互偏微分の平均、ｆａｖｒ，ａｏｙにより角度Ａで偏微分してから中心Ｏｙで偏微分する相互偏微分の平均、ｆａｖｒ，ｂにより当該二乗距離の角度Ｂによる偏微分の平均、ｆａｖｒ，ｂｂにより角度Ｂによる２回偏微分の平均、ｆａｖｒ，ｂｏｘにより角度Ｂで偏微分してから中心Ｏｘで偏微分する相互偏微分の平均、ｆａｖｒ，ｂｏｙにより角度Ｂで偏微分してから中心Ｏｙで偏微分する相互偏微分の平均、ｆａｖｒ，ｏｘにより中心Ｏｘでの偏微分の平均、ｆａｖｒ，ｏｙにより中心Ｏｙでの偏微分の平均を作成する（ステップ３３５）。
【０１７１】
次に、（数２８）により二乗距離と二乗距離の平均の差の二乗和、当該二乗和の角度Ａによる偏微分、角度Ａによる２回偏微分、角度Ａで微分してから角度Ｂで微分する相互偏微分、角度Ａで微分してから中心Ｏｘで微分する相互偏微分、角度Ａで微分してから中心Ｏｙで微分する相互偏微分、角度Ｂによる偏微分、角度Ｂによる２回偏微分、角度Ｂで微分してから中心Ｏｘで微分する相互偏微分、角度Ｂで微分してから中心Ｏｙで微分する相互偏微分、中心Ｏｘでの偏微分、中心Ｏｙでの偏微分を作成する（ステップ３３６）。
【０１７２】
次に、ステップ３３６で求めた値を連立一次方程式の係数に代入する（ステップ３３７）。
【０１７３】
次に、当該連立一次方程式を解き、更新量（δＡ，δＢ）を算出する（ステップ３３８）。当該連立一次方程式を解くには一般的な連立一次方程式の解法、例えば、掃き出し法やＬＵ分解を使えば良い。
【０１７４】
次に、更新量（δＡ，δＢ，δＯｘ，δＯｙ）が０とみなせるか判定する（ステップ３３９）。判定方法は、更新量の成分（δＡ，δＢ，δＯｘ，δＯｙ）の絶対値がそれぞれ０とみなせるかを判定する方法のほかに、更新量の内積（δＡ×δＡ＋δＢ×δＢ＋δＯｘ×δＯｘ＋δＯｙ×δＯｙ）が０とみなせるかを判定してもかまわない。本実施例では、更新量の内積が０とみなせるかを判定している。
【０１７５】
更新量が０とみなせない場合、角度・中心点更新をおこなう（ステップ３４０）。更新は、初期値（Ａ，Ｂ，Ｏｘ，Ｏｙ）に更新量（δＡ，δＢ，δＯｘ，δＯｙ）を加算することでおこなう。更新した初期値（Ａ，Ｂ，Ｏｘ，Ｏｙ）でステップ３３２から繰り返す。更新量が０とみなせる場合、処理を終了する。
【０１７６】
以上のステップ３３１からステップ３４０を繰り返すことで、円筒の中心（Ｏｘ、Ｏｙ）と角度（Ａ，Ｂ）を定めることができる。
【０１７７】
実際に計測点を与え実験してみたところ、初期角度として角度精度向上部１０７で算出した角度、初期中心として球中心算出部１０５で作成した球中心を当該角度で座標変換した中心を与えたとき、５から６回で収束した。
【０１７８】
尚、初期角度算出部１０６、角度精度向上部１０７の効果を確認するため、実際に計測点を与え、分割ピッチを３６０°／Ｎ、Ｎは１から１９まで変化させ、角度精度向上部１０７をスキップした場合としない場合で円筒半径の違いを確認した。
【０１７９】
角度精度向上部１０７をスキップした場合、Ｎが１５より大きくなった場合に円筒半径が一定値となった。これは、２４°ピッチに相当する。これに対し、角度精度向上部１０７をスキップしない場合、Ｎが７より大きい場合に一定値となった。これは約５１．４２°ピッチに相当する。よって、角度精度向上部１０７は角度ピッチを２倍まで大きくできることがわかった。このことから、現在、角度ピッチは余裕をみて２０°ピッチとしている。
【０１８０】
しかしながら、角度ピッチは与える計測点によって変化する可能性があるので記憶装置２０４にあらかじめ格納しておき、適宜変更可能な構成としている。また、初期角度算出部１０６の計算量は角度ピッチを２倍にすると約４倍に増加するため、角度精度向上部１０７が計算量を減少させることがわかった。計算量の減少は計算時間の減少に繋がるため、その効果は大きいことがわかった。
【０１８１】
図３のフローチャートの説明に戻る。円筒中心・角度算出部１０８の処理が終了すると、円筒データ格納部１０９の処理を実行する。円筒データ格納部１０９は、円筒中心・角度算出部１０８で作成した円の中心Ｏ、角度を記憶装置１０２に格納する。
【０１８２】
ここで、円の中心Ｏは２次元座標であるため、Ｚ座標を０とし、（数３）に示す座標変換の逆変換をおこなえば、円筒の中心点が得られる。また、方向ベクトル（０，０，１）を同様に逆変換すれば、円筒の中心軸が得られる。半径は、（数２９）で作成したｆａｖｒの平方根である。
【０１８３】
円筒データが計測点を含むようにするには、中心点と計測点の差ベクトルと、中心軸の内積を作成して最低値と最大値を作成し、最低値が円筒の底面になるよう中心点を中心軸方向に移動すれば良い。そのとき、円筒の高さは最大値−最小値である。
【０１８４】
以上、図１〜図５を用いて本発明を説明したが、最後に図６を用いて本発明の方式を説明する。図６は、本発明の各処理に応じて作成される形状のイメージ図である。
【０１８５】
図６では、最初に計測点群が与えられる。次に、計測点群に球をあてはめ、球中心を作成する。次に、球中心を円筒の中心と仮定し、投影したときに円との分散が最小となる角度を少なくとも一つ作成する。次に当該角度を高精度化し、最後に、当該角度および円筒の中心点を高精度化して円筒あてはめを終了する。
【０１８６】
以上により、従来の問題点である、あてはめ時間増大、別解に収束する問題を解決する。
【０１８７】
図９は、本発明で説明した球、円筒の３次元モデルを別々、及び重ねて表示した結果であり、最初に計測点から球をあてはめ、仮円筒と円筒をあてはめていった結果を記している。また、最下段は計測点、球、円筒を重ねて表示した結果である。計測点と球、円筒の位置関係が判別できる。
【０１８８】
図１０は、何回目の計測がどの程度の時間なのかを示す図である。
【０１８９】
各処理は１０００回繰り返しているので、左端の目盛の単位はｍｓｅｃである。円弧溝の計算時間が多いのは、円筒３、円筒４の計測点数が８であるのと比較して、円弧溝の計測点数が１５０点になっているためである。これにより、十分高速な処理ができていることが理解できる。
【符号の説明】
【０１９０】
１０１…計測機器、１０２…記憶装置、１０３…円筒データ抽出部、１０４…計測点・収束判定値・角度ピッチ読取部、１０５…球中心算出部、１０６…初期角度算出部、１０７…角度精度向上部、１０８…円筒中心・角度算出部、１０９…円筒データ格納部、２００…計算機、２０１…通信装置、２０２…入力装置、２０３…表示装置、２０４…記憶装置、２０５…メモリ、２０６…中央処理装置。
【技術分野】
【０００１】
本発明は、円筒形状の対象物を計測して得た計測点から円筒の中心点、中心軸、半径（以下、円筒データ）を抽出する計測点からの円筒情報の抽出方法に関するものである。
【背景技術】
【０００２】
計測点などの測定誤差を持つ点群から円筒データを抽出する従来の方法の一例を、以下に説明する。
【０００３】
従来の方法は、次のようなステップで行われている。即ち、ステップ１）：円筒データの初期値を設定する。ステップ２）：当該円筒表面と計測点群の距離を作る。ステップ３）：距離の二乗和が最低になるよう円筒データの更新量を定める。ステップ４）：当該更新量が収束判定値以下なら終了。そして、当該二乗和が、収束判定値以上なら当該更新量により初期値を更新してステップ２）から繰り返す。
【０００４】
ステップ１）〜４）の方法は、解に十分近い初期値を与えれば収束させることができる。しかし、円筒の表面と計測点群との距離の二乗和は、極値を複数持つ複雑な関数であるため、解から離れた初期値が与えられた場合、収束が遅くなったり、別の解に収束してしまう問題がある。
【０００５】
この問題を解決する方法として、あらかじめ円筒データが存在する範囲を限定し、当該範囲を細かく分割して良好な初期値を定める方法が知られているが、範囲の設定方法が不明確なこと、どの程度の細かさまで分割すれば良いかが不明確なこと、分割数が多くなった場合の処理時間の増大が問題になっている。
【０００６】
また、関連する技術として、非特許文献１、２及び３に記載されている以下の（１）、（２）、（４）、特に、先行技術文献は提示しないが（３）の技術が知られている。
（１）３点を通過する円
非特許文献１に記載されているのは、３次元空間内に３点Ｐ０、Ｐ１、Ｐ２を与え、３点を通過する円の中心点を作成する方法である。
【０００７】
これは、先ず３点中の２点を２組取出し、２点を結ぶ単位ベクトルを法線とし、２点の中間点を通過する平面を２枚作成する。更に、３点を通過する平面を１枚作成する。次に、当該３枚の平面が交差する１点を求める。これにより、３点を通過する円の中心点が作成できる。
【０００８】
（数１）は、３点Ｐ０、Ｐ１、Ｐ２と、円の中心点Ｏの関係を示す連立一次方程式である。（数１）において、添え字０、１、２は点を識別する番号、ｘ、ｙ、ｚは点の座標成分であり、符号×は、ベクトルの外積、符号・はベクトルの内積である。
【０００９】
【数１】
【００１０】
（数１）を連立一次方程式の解を求める方法、例えば、掃き出し法やＬＵ分解などで解けば、３点を通過する円の中心点を３次元空間内に定めることができる。
（２）複数の点に近い円
非特許文献２に記載されているのは、２次元空間内の複数の点Ｑｉを与え、複数の点と、円との距離の二乗和が最小になる円の中心点を作成する方法である。
【００１１】
（数２）は、当該方法の関係を整理した結果であり、（数２）を点数分作成した連立１次方程式を一般化逆行列を用いて解くことにより中心点を求める。
【００１２】
【数２】
【００１３】
（数２）は平面上の複数の点に対し、円の中心点を作成するために効率の良い方法である。
（３）複数の点に近い球
これは、３次元空間内の複数の点Ｐｉを与え、複数の点と球との距離の二乗和が最小になる球の中心点を作成する方法である。
【００１４】
（数３）は当該方法の関係を整理した結果であり、（数３）を点数分作成した連立１次方程式を一般化逆行列を用いて解くことにより中心点を求める。
【００１５】
【数３】
【００１６】
（数３）は、空間内の複数の点に対し、球の中心点を作成するために効率の良い方法である。
（４）３次元から２次元への座標変換
非特許文献３に記載されているのは、３次元の点Ｐｉを、１番目の軸まわりに角度Ｂだけ回転した後、２番目の軸まわりに角度Ａだけ回転する例である。ここで、３次元座標Ｐｉから２次元座標Ｑｉへの変換は、３次元の点Ｐｉを、１番目の軸まわりに角度Ｂだけ回転した後、２番目の軸まわりに角度Ａだけ回転した座標（ｘ、ｙ）成分を取り出したものである。例えば、１番目をｙ軸、２番目をｘ軸とした場合、（数４）となる。
【００１７】
【数４】
【００１８】
（数４）において、Ｏｉ，ｘとＱｉ，ｙを成分とするベクトルＱｉを角度Ａ、角度Ｂでｎ回偏微分または相互偏微分した値を求めることは公知である。また、（数４）は、Ｙ軸回転してからＸ軸回転して（ｘ、ｙ）成分を取り出しているが、Ｘ軸回転、Ｙ軸回転、Ｚ軸回転を任意に組み合わせることも公知である。その場合、（数４）の回転行列成分が変化する。
【先行技術文献】
【非特許文献】
【００１９】
【非特許文献１】コンピュータディスプレイによる形状処理工学［Ｉ］、ｐ.１９２ ３点を通る円、日刊興行新聞社、昭和５７年１０月３０日発行
【非特許文献２】Coope, I.D., Circle fitting by linear and nonlinear least squares, Journal of Optimization Theory and Applications Volume 76, Issue 2, New York: Plenum Press, February 1993
【非特許文献３】参考文献３：コンピュータディスプレイによる図形処理工学、ｐ.１１２ ｘ／ｙ／ｚ軸まわりの回転、日刊興行新聞社、１９８１年６月３０日発行
【発明の概要】
【発明が解決しようとする課題】
【００２０】
しかしながら、非特許文献１では、３点に限定していること及び３点で作成した平面上の円しか作成できないことから、複数の計測点に最も近い円の中心を求める場合には、適用が難しい。また、非特許文献２では、３次元空間で円筒の中心軸を定めるには、非特許文献３を用いて、３次元空間の複数の点を回転変換し、点から２つの成分を取り出すことにより平面に投影した後で適用する必要があり、この方法では、回転変換時の角度および投影後の中心点が変数になり、回転変換は三角関数を使用するため、角度で無限回の微分が可能である。即ち、極値を複数持つ複雑な関数となるので、複数の計測点に最も近い円の中心を求める場合には適用が難しい。更に、（３）の複数の点に近い球では、球ではなく円筒を定める場合には利用できない。
【００２１】
本発明は上述の点に鑑みなされたもので、その目的とするところは、円筒表面の計測点から円筒データを確実に作成できると共に、従来方式に比べて、高精度な近似が可能になる計測点からの円筒情報の抽出方法を提供することにある。
【課題を解決するための手段】
【００２２】
本発明の計測点からの円筒情報の抽出方法は、上記目的を達成するために、計測点が円筒物の表面を計測していることを条件として、全計測点から球の中心を抽出する手段と、全計測点および球の中心から角度の初期値を作成する手段と、当該角度の初期値の精度向上を行う手段と、球中心を精度向上した角度で投影変換し、初期円筒中心を定める手段と、初期円筒中心および精度向上した角度を初期値とし、円筒中心および角度を定める手段と、を用意し計測点から円筒の半径、中心軸からなる円筒データを抽出することを特徴とするものである。
【００２３】
更に、全計測点および球の中心から角度の初期値を作成する場合に、角度を離散化するための分割数、当該分割数と当該分割数における近似曲面との差の二乗の最大Ｂの初期値、当該分割数における近似曲面との差の二乗の最大Ａを最大Ｂで除した値で繰り返し終了を判定するための初期値を設定する手段と、２つの角度を当該分割数で離散化した格子リストを作成する手段と、当該格子リストの角度毎に、計測点と球の中心を座標変換して平面上に投影し、投影した計測点と投影した中心の二乗距離を作成し、当該二乗距離の平均を作成し、当該二乗距離と当該平均の差の二乗和を作成する手段と、当該格子リストに対応した二乗和から近似曲面を作成する手段と、当該格子リストに対応した二乗和と近似曲面との距離の二乗の最大Ａを求める手段と、最大Ａが最大Ｂより大きいかを判定する手段と、最大Ｂに最大Ａを代入する手段と、最大Ａを最大Ｂで除した値がしきい値より大きいかを判定する手段と、分割数を増加させる手段と、を用意し、全計測点および球の中心から角度の初期値を作成することを特徴とするものである。
【発明の効果】
【００２４】
本発明の計測点からの円筒情報の抽出方法によれば、円筒表面の計測点から円筒データを確実に作成できると共に、従来方式に比べて、高精度な近似が可能になる。
【図面の簡単な説明】
【００２５】
【図１】本発明の計測点からの円筒情報の抽出方法を実施する装置を示す構成図である。
【図２】図１の円筒データ抽出処理を実施する計算機の構成図である。
【図３】図１の円筒データ抽出処理の処理を示すフローチャートである。
【図４（ａ）】投影するための角度を変数とし、投影した計測点と円との距離の二乗の総和を縦軸とした３次元グラフの上面図である。
【図４（ｂ）】図４（ａ）の斜視図である。
【図４（ｃ）】図４（ａ）及び図４（ｂ）の基となる値を記入した図である。
【図５（ａ）】球中心を、投影した計測点と円との距離の二乗の総和が最小になる４組の角度で投影したときの円中心の図である。
【図５（ｂ）】４組の角度の角度値と円中心を記述した図である。
【図６】本発明の処理の内容を図示した説明図である。
【図７】初期角度算出部の別の実施方法を示すフローチャートである。
【図８（ａ）】角度Ａ，Ｂを離散化するときの説明図である。
【図８（ｂ）】二乗和と二乗和を近似した曲面との差を説明する説明図である。
【図８（ｃ）】分割数と当該分割数における当該差の二乗和の最大の変化を示す特性図である。
【図９】本発明で、球あてはめ、仮円筒および円筒あてはめをおこなった結果の３次元モデル、当該計測点、球、円筒を重ね合わせた３次元モデルを表示した例を示す図である。
【図１０】本発明を実施したときの球中心算出処理、初期角度算出処理、角度精度向上処理、円筒中心・角度算出処理の計測順番と計算時間を示す図である。
【発明を実施するための形態】
【００２６】
以下、本発明の計測点からの円筒情報の抽出方法の実施例を、図１乃至図１０を用いて説明する。
【実施例】
【００２７】
図１は、本発明の計測点からの円筒情報の抽出方法を実施する装置（後述する計算機で実行するプログラムでもある）を示すものであり、該図において、１０１は円筒物の表面を３次元座標（Ｘ、Ｙ、Ｚ）で計測する計測機器、１０２は計測機器１０１で計測した計測点、収束判定に用いる０とみなす値、初期角度算出部１０６で使用する角度のピッチを保存する記憶装置、１０３は計算機で実行する円筒データ抽出部、１０４は記憶装置１０２から３次元の計測点Ｐｉ、収束判定に用いる０とみなす値、初期角度算出部１０６で使用する角度のピッチを読み取る計測点・収束判定値・角度ピッチ読取部、１０５は計測点Ｐｉから球の中心を算出する球中心算出部、１０６は記憶装置に格納されている角度ピッチで分割した角度で、計測点および１０５で算出した球中心を回転変換し、回転変換した計測点と球中心のベクトルの２成分を取り出すことにより、投影した球中心と投影した計測点を定め、投影した球中心と投影した計測点の距離を求めると共に、当該距離の平均を求め、当該距離と当該平均の差の二乗和を求め、当該二乗和が最小になる初期角度を算出する初期角度算出部、１０７は１０６で作成した初期角度の精度向上を行う角度精度向上部、１０８は、１０６で作成した投影後の球中心を初期円筒中心、１０７で作成した角度を初期角度とし、円筒中心、角度、半径を算出する円筒中心・角度算出部、１０９は１０８で作成した円筒中心、角度、半径を記憶装置１０２に格納する円筒データ格納部である。
【００２８】
図２は、図１の円筒データ抽出処理を実施する計算機を示すものであり、２００は、円筒データ抽出部１０３を実行する計算機である。計算機２００において、２０１は通信網で接続された他の計算機とデータのやり取りをおこなうための通信装置、２０２はマウスやキイボードなどのデータの入力をおこなう入力装置、２０３はグラフィックディスプレイなどの計算結果を文字やグラフなどで表示するための表示装置、２０４は計算結果を蓄積するＦＤ等の記憶装置（記憶媒体）、２０５はメモリ、２０６は１０３の円筒データ抽出処理をおこなう中央処理装置である。
【００２９】
図３は、本発明の円筒データ抽出処理の流れを示すフローチャートである。以下、図３のフローチャートと数式を交えながら、本発明の一実施例を述べる。尚、図３は、本発明のプログラムに相当する。
【００３０】
最初に計測点・収束判定値・角度ピッチ読取部１０４が、記憶装置１０２から計測点データ、収束判定値、角度ピッチを読み取る。計測点データは、（数５）に示すようにＰｉでｉ番目の計測点を表し、（ｘ成分、ｙ成分、ｚ成分）を持つ３次元空間上の点である。ここで、ｉは１から計測点数Ｎまで変化する。
【００３１】
【数５】
【００３２】
収束判定値は、例えば１．０Ｅ−２０のように更新量を０とみなす値である。角度ピッチは、０から３６０°までの角度を離散化するときの幅であり、例えば２０°とすれば、０°、２０°、・・・、３４０°、３６０°の角度に離散化する。
【００３３】
次に、球中心算出部１０５がＮ個の計測点と球表面との距離の二乗の総和が最小になる球（以下、最小二乗球という）の中心点を算出する。以下、その算出法について述べる。
【００３４】
最小二乗球の中心Ｏｓは、例えば（数３）で算出することができるが、本実施例では、初期角度算出部１０６、角度精度向上部１０７、円筒中心・角度算出部１０８を説明する都合上、（数６）に示す方法を用いる。
【００３５】
【数６】
【００３６】
（数６）において、Ｏｓは球の中心、ｇｉは点Ｐｉと中心Ｏｓの距離の二乗（半径の二乗）、ｇａｖｒはｇｉの平均である。即ち、各点Ｐｉと中心Ｏｓの半径の二乗と、当該半径の二乗の平均との差の二乗和が最小になるよう中心Ｏｓを決定する。
【００３７】
（数６）は厳密に言えば、計測点と球表面の距離ではなく、距離の二乗を最小化する。しかしながら、距離の場合は、平方根を用いた関数とせねばならず、微分時に複雑さが増してしまう。そのため、距離の二乗を最小化している。この方法は一般的に最小二乗法を用いるときの常套手段である。尚、（数３）と（数６）の解が一致することは実験で確認済である。
【００３８】
また、（数６）のようにｇｉとｇａｖｒの差を作成すると、ｇｉ中に同一値を持つ変数成分及び定数成分がある場合に０となるので、数値が大きくなることがないため、計算機処理する場合の丸め誤差を少なくする上で都合が良い。
【００３９】
さらに、ｇａｖｒが加算と定数による除算なので、ｇｉとｇａｖｒの次数が同一となり、次数が増加しないので関数の複雑さが増すことがない。
【００４０】
（数６）を解くために、（数６）を球の中心Ｏｓの各成分（Ｏｓｘ、Ｏｓｙ、Ｏｓｚ）で偏微分した関数が０になる式を作成し、当該式を１次までＴａｙｌｏｒ展開して連立一次方程式を作成する。当該連立一次方程式を、初期中心Ｏｓを与えて解くことにより、中心更新量を算出する。本方法は、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法と呼ばれている。
【００４１】
当該連立一次方程式は、２回微分が０または定数となる比較的単純な関数であるため、任意の初期中心でも収束する。
【００４２】
（数７）は当該連立一次方程式である。
【００４３】
【数７】
【００４４】
（数３）が４行４列であるのに対し、（数７）は３行３列である。このように（数６）で示すようにｇｉとｇａｖｒの差を求める方法を用いると、行数および列数を１行１列減少させることができるので、計算時間の短縮が可能なため、本実施例ではこの方法を採用している。
【００４５】
（数７）において各行列成分Ｍｓは（数８）である。
【００４６】
【数８】
【００４７】
（数８）においてｇｉ及びｇａｖｒは（数６）、ｇｉの中心Ｏｓの成分による微分及び平均は、（数９）及び（数１０）で表せる。
【００４８】
【数９】
【００４９】
【数１０】
【００５０】
（数１０）に示すように、ｇｉの中心Ｏｓの成分による２回偏微分は定数２か０である、従って、平均も定数２か０である。そのため、（数８）は（数１１）のように簡単化される。
【００５１】
【数１１】
【００５２】
また、（数１１）で使用するｇｉの中心Ｏｓの成分による２回偏微分及び平均は、（数９）のみとなる。
【００５３】
以上、数式により（数６）の解き方を説明したが、次に、図３に示す球中心算出部１０５のフローチャートを用いて球中心の算出過程を説明する。
【００５４】
図３の３０１から３１０は、（数６）を解く過程の一実施例を示すフローチャートである。最初に球の中心点の初期値を例えば（０，０，０）のように与える（ステップ３０１）。当該初期値は（０，０，０）に限定されるものではなく、任意の値、例えば、（１００，１００，１００）や（−１，５００，０）などでもかまわないが、数値計算上は予想できる中心点を与えることが望ましい。本実施例では、予想できないこと、及び更新量がそのまま球の中心点となることから（０，０，０）を与えている。
【００５５】
次に、計測点Ｐｉと中心点の差ベクトルを求める（ステップ３０２）。
【００５６】
次に、（数６）のｇｉに示す差ベクトルの内積を求め、計測点数分の内積値を作成する（ステップ３０３）。当該内積値は中心点と計測点Ｐｉの距離の二乗であることから、内積値を二乗距離と呼ぶ。
【００５７】
次に、（数９）の（ｇｉ，ｏｓ，ｘ）、（ｇｉ，ｏｓ，ｙ）、（ｇｉ，ｏｓ，ｚ）に示す二乗距離の中心点による偏微分を作成する（ステップ３０４）。
【００５８】
次に、（数６）のｇａｖｒ、（数９）の（ｇａｖｒ，ｏｓ，ｘ）、（ｇａｖｒ，ｏｓ，ｙ）、（ｇａｖｒ，ｏｓ,ｚ）に示す二乗距離および二乗距離の中心点による偏微分の平均を作成する（ステップ３０５）。
【００５９】
次に、ステップ３０４で求めた偏微分およびステップ３０５で求めた偏微分の平均から、（数１１）の（Ｍｓ，０，０）、（Ｍｓ，２，０）、（Ｍｓ，０，１）、（Ｍｓ，２，１）、（Ｍｓ，０，２）、（Ｍｓ，２，２）に示す二乗距離の総和の２回偏微分を作成するとともに、ステップ３０３で求めた二乗距離とステップ３０４で求めた偏微分およびステップ３０５で求めた二乗距離の平均と二乗距離の偏微分の平均から（数１１）の（Ｍｓ，３，０）、（Ｍｓ，３，１）、（Ｍｓ，３，２）に示す二乗距離の総和の偏微分を作成する（ステップ３０６）。
【００６０】
次に、ステップ３０６で求めた２回偏微分と偏微分から（数７）に示す連立一次方程式の係数を作成する（ステップ３０７）。
【００６１】
次に、連立一次方程式を解き、中心点の更新量を求める（ステップ３０８）。連立一次方程式の解法は、例えば、掃き出し法やＬＵ分解で解くことができる。また、３行３列なので解析的に解いてもかまわない。
【００６２】
次に、中心点の更新量が１０４で読み込んだ０とみなせる値になっているかを判定する（ステップ３０９）。判定時には、例えば、中心点の更新量（δＯｓｘ，δＯｓｙ，δＯｓｚ）の絶対値を作成して，当該絶対値がそれぞれ０とみなせる値になるかどうかで確認する。そのほかに、更新量の内積（δＯｓｘ×δＯｓｘ＋δＯｓｙ×δＯｓｙ＋δＯｓｚ×δＯｓｚ）を作成して０となるかどうかを判定してもかまわない。本実施例では、更新量の内積で判定している。
【００６３】
ステップ３０９で０に近くないと判定された場合、中心点の更新を行う（ステップ３１０）。中心点の更新は、現在の中心点（Ｏｓｘ，Ｏｓｙ，Ｏｓｚ）に、更新量（δＯｓｘ，δＯｓｙ，δＯｓｚ）を加算すればよい。この加算した値を現在の中心点（Ｏｓｘ，Ｏｓｙ，Ｏｓｚ）とし、ステップ３０２から繰り返す。ステップ３０９で０に近いと判定された場合、球中心算出部１０５の処理を終了する。
【００６４】
以上のステップ３０２からステップ３０９を繰り返すことで、球の中心点を定めることができる。
【００６５】
実際には、球中心算出部１０５の処理は、殆ど１回で繰り返し計算が終了する。これは、ｇｉの中心Ｏｓの成分による偏微分が、１回偏微分が変数、２回の偏微分および相互偏微分が定数または０、３回以降の偏微分および相互偏微分が０であることから、（数７）に示す関数が線形関数であることに関連している。但し、数値計算誤差が含まれることも想定されるので、本実施例に示したように、ステップ３０２〜ステップ３０９、ステップ３１０の処理を行うことが精度を高める上で望ましい。
【００６６】
３次元で定義される点に球をあてはめる場合、（数７）に示すように、中心Ｏｓのみを求めることとなり回転角度には依存しない。このため、特異な場合を除き確実に中心点Ｏｓを定めることができるので、あてはめ時間を最小限に抑制することができる。
【００６７】
ここで、特異な場合とは、例えば、３次元で定義された点が直線上に配置された場合である。本状態は、計測物が円筒面の場合に限定すると、円筒の中心軸に平行となるよう計測したときに生じる。このとき、（数７）で繰り返し計算をおこなうと無限に繰り返すか、または０除算が発生してしまう。そのため、繰り返し回数に上限を設け、上限を超えた場合、及び（数７）の連立一次方程式を解くときに０除算が発生するか確認し、エラー出力をおこなうことが望ましい。
【００６８】
当該エラーが出力された場合、計測位置を変更する必要がある。そのため、（数７）の処理を行う計算機、及び計算の結果をモニタできる装置、例えばブザーや小型の表示装置を設け、計測結果をリアルタイムにモニタリングすることが望ましい。本構成を採用する場合、（数７）の計算処理が単純であり、繰り返し回数も少ないことから計算能力が劣っている計算機でも採用することができるとともに、計測位置の誤りがリアルタイムで手元確認できるため、すぐに計測位置を修正することができるので計測間違いが少なくなる。
【００６９】
次に、初期角度算出部１０６で、円筒の中心軸となる角度の初期値を算出する。初期角度を算出するために（数１２）を用いる。
【００７０】
【数１２】
【００７１】
（数１２）において、Ｑｉは、計測点ＰｉをＹ軸で角度Ｂだけ回転してからＸ軸で角度Ａだけ回転した後の（ｘ、ｙ）成分である。Ｏは、球の中心ＯｓをＹ軸で角度Ｂだけ回転してからＸ軸で角度Ａだけ回転した後の（ｘ、ｙ）成分である。
【００７２】
（数１２）では、投影平面を（ｘ、ｙ）平面に限定し、回転をＹ軸回転後にＸ軸回転しているが、投影平面は（ｘ、ｙ）平面に限定されるものではなく、（ｙ、ｚ）平面、（ｚ、ｘ）平面でもかまわないし、回転は、Ｙ軸回転後のＺ軸回転、Ｘ軸回転後のＹ軸回転、Ｘ軸回転後のＺ軸回転、Ｚ軸回転後のＸ軸回転、Ｚ軸回転後のＹ軸回転のいずれでもかまわない。
【００７３】
重要なことは、初期角度算出部１０６、角度精度向上部１０７、円筒中心・角度算出部１０８で回転順序を一致させておくことである。これを誤ると角度の高精度化、円筒の中心や角度の算出ができなくなるので注意が必要である。
【００７４】
（数１２）は三角関数が入るため、角度ＡおよびＢで無限に偏微分、相互偏微分することが可能であり極値を複数持つ複雑な関数である。そのため、本願では（数１３）で示すように、（数１２）の角度（Ａ，Ｂ）を０〜３６０°の中で離散化して（数１２）の値を求め、最小値を持つ角度（Ａ，Ｂ）を探す構成とした。
【００７５】
【数１３】
【００７６】
以上、数式により初期角度算出部１０６の処理を説明したが、次に図３に示す初期角度算出部１０６のフローチャートで初期角度の算出課程を説明する。図３において、３１１から３１６は、初期角度算出部１０６の処理の一実施例を示すフローチャートである。
【００７７】
最初に（数１３）に示す角度（Ａ，Ｂ）を、計測点・収束判定値・角度ピッチ読取１０４で読み込んだ角度ピッチに従って離散化する（ステップ３１１）。本実施例では、角度ピッチを２０°としている。２０°という値は、テストデータを用いたとき角度精度向上部１０６で収束するピッチを実験して決定したものである。本結果については後述する。
【００７８】
尚、初期角度算出部１０６において、円筒中心も変数とした場合、円筒中心が３成分毎に−∞〜∞まで変化するため設定が難しい。また、変数が５個となり組合せ数が一気に増大するため計算量が増大してしまう。更に、変数５個それぞれに極値を持つため大変複雑になってしまう。そのため、球中心算出部１０５を用意して角度だけで（数１３）を算出する構成としている。
【００７９】
次に、計測点Ｐｉ、球中心算出手段１０５で作成した球中心を離散化した角度毎に座標変換する（ステップ３１２）。座標変換の方法は、上述した“（４）３次元から２次元への座標変換”で説明した処理でかまわない。
【００８０】
次に、（数１２）のｆｉに示すように座標変換した計測点と球の中心の（ｘ、ｙ）成分を取り出して差ベクトルを作成し、差ベクトルの内積を作成して（ｘ、ｙ）成分の二乗距離を作成する（ステップ３１３）。
【００８１】
次に、（数１２）のｆａｖｒに示すように当該二乗距離の平均を作成する（ステップ３１４）。
【００８２】
次に、二乗距離と平均の差を作成し、当該差の二乗和を作成する（ステップ３１５）。
【００８３】
最後に、当該二乗和の最低値を取出し、最低値と一致する離散化した角度（Ａ，Ｂ）を抽出する（ステップ３１６）。
【００８４】
最低値を作成する方法としては、当該二乗和を小さい値から昇順に並べ替え最初の値を抽出する方法や、最初の値を最低値の初期値とし２番目から計測点数番目までの間で当該最低値より小さい値があらわれたとき置き換えるなどの方法を用いてもかまわない。並べ替えの方法としては、バブルソートやクイックソートなどのソーティングアルゴリズムを利用してもかまわない。
【００８５】
ここで、最低値を持つ角度（Ａ，Ｂ）が複数存在することを考慮し、ステップ３１５の二乗和は小数点以下の有効桁数を考慮して整数化しておくことが望ましい。整数化のかわりに四捨五入した実数を使うことも考えられるが、実数値に丸め誤差が累積して微小に異なる場合があるため、整数化しておく方が良い。
【００８６】
本実施例では、小数点以下１桁までを有効とし、実数値を１０倍して切り捨てる方法を使用している。
【００８７】
図４（ａ）は、計測点８点を与え、角度Ａおよび角度Ｂをそれぞれ０から３６０°まで２０°ピッチで分割して（数１３）のＦ（Ａ，Ｂ）をプロットした３次元グラフの上面図、図４（ｂ）は、図４（ａ）の斜視図、図４（ｃ）は、角度Ａを横軸、角度Ｂを縦軸としＦ（Ａ，Ｂ）の値を記入した図である。図４（ｃ）では、最低値を四角で囲んで示している。
【００８８】
図４（ａ）、図４（ｂ）、図４（ｃ）より、Ｆ（Ａ，Ｂ）が極値を多数持つ複雑な関数であり、最低値が４個程度存在することがわかる。
【００８９】
図５（ａ）の１〜４は、球の中心点を複数求めた角度で座標変換した位置の（ｘ、ｙ）座標である。図５（ａ）に示すように、中心点が角度毎に象限を分けて算出されている。
【００９０】
このように、中心点が分れて検出されることから、角度と中心を変数として繰り返し計算すると、中心の更新量に従い角度が大きく変化するため解が大きく変動することに繋がる。そのため、初期角度算出部１０６では球中心算出部１０５で求めた中心を円筒の中心とみなし、角度の初期値を定めているので、解が大きく変動するのを防ぐことができる。
【００９１】
次に初期角度算出部１０６を実現する別の処理について、図７及び図８（ａ）、図８（ｂ）、図８（ｃ）を用いて説明する。
【００９２】
図７は、初期角度算出部１０６の別の実施方法を示すフローチャートである。該図に示す如く、最初に分割数、最大Ｂ、しきい値を設定する（ステップ９０１）。分割数は２つの角度（Ａ，Ｂ）を離散化するときの点数、最大Ｂは分割数を変数としたときの当該分割数における関数Ｆ（Ａ，Ｂ）と当該関数の近似曲面の距離の二乗の最大を示す関数の最大の初期値、しきい値は分割数の決定処理が終了したかどうかを表す値である。尚、図７は、本発明の別の実施方法のプログラムに相当する。
【００９３】
本実施例では、分割数は３を設定しているが、関数Ｆ（Ａ，Ｂ）の複雑さに応じ、３より大きい値を設定してもかまわない。また、最大Ｂは０を設定しているが、関数Ｆ（Ａ，Ｂ）の最小が０となるからであり、当該最小が０以上となることが確実な場合、０より大きい値を設定してもかまわない。また、しきい値は０．１を設定しているが、角度精度向上部１０７の処理が発散する場合０〜０．１の値を設定する必要がある。
【００９４】
更に、角度精度向上部１０７の処理が確実に収束するなら、０．１〜１の範囲に設定してもかまわない。本実施例でいろいろな計測点を与え、実験したところ、発散しなかったことから０．１に設定している。
【００９５】
次に、角度Ａの範囲０〜３６０°、角度Ｂの範囲０〜３６０°を当該分割数個に離散化する（ステップ９０２）。本実施例では、関数Ｆ（Ａ，Ｂ）の傾向が不明なため、角度が等感覚になるような離散化をおこなっている。例えば、分割数が３なら、（０，０）、（１８０，０）、（３６０，０）、（０，１８０）、（１８０，１８０）、（３６０，１８０）、（０，３６０）、（１８０，３６０）、（３６０，３６０）の格子リストとなる。
【００９６】
次に、格子リストの角度毎に関数Ｆ（Ａ，Ｂ）を作成する（ステップ９０３）。例えば、分割数が３なら、Ｆ（０，０）、Ｆ（１８０，０）、Ｆ（３６０，０）、Ｆ（０，１８０）、Ｆ（１８０，１８０）、Ｆ（３６０，１８０）、Ｆ（０，３６０）、Ｆ（１８０，３６０）、（Ｆ３６０，３６０）を作成する。９０３の処理過程は、図２のステップ３１２、３１３、３１４、３１５と同様である。
【００９７】
次に、格子リストの角度毎に関数Ｆ（Ａｉ−１，Ｂｊ−１）、Ｆ（Ａｉ＋１，Ｂｊ−１）、Ｆ（Ａｉ−１，Ｂｊ＋１）、Ｆ（Ａｉ＋１，Ｂｊ＋１）から近似曲面Ｆａｉ，ｊ（ｕ，ｖ）を作成し、Ｆ（Ａｉ，Ｂｊ）とＦａｉ，ｊ（０.５，０.５）の距離を作成し、当該距離の最大である最大Ａを求める（ステップ９０４）。近似曲面Ｆａｉ，ｊ（ｕ，ｖ）は、ｍ×ｎ次の多項式曲面、Ｂｅｚｉｅｒ曲面、Ｂ−Ｓｐｌｉｎｅ曲面、Ｃｏｏｎｚ曲面など任意の形式でかまわない。また、これらの曲面を有理式化したものでもかまわない。
【００９８】
本実施例では、処理を簡単にするため、双１次パッチを採用し、双１次パッチの（ｕ，ｖ）方向それぞれの中央で評価する構成としている。
【００９９】
（数１４）に双１次パッチによる近似曲面Ｆａｉ，ｊ（ｕ，ｖ）を示す。
【０１００】
【数１４】
【０１０１】
（数１４）において、Ｆａｉ，ｊ（ｕ，ｖ）は、パラメータ（ｕ，ｖ）を変数とする双１次曲面（双１次Ｂｅｚｉｅｒ曲面でもある）。パラメータ（ｕ，ｖ）はそれぞれ、０〜１まで変化する変数。Ｆ（Ａｉ−１，Ｂｊ−１）、Ｆ（Ａｉ＋１，Ｂｊ−１）、Ｆ（Ａｉ−１，Ｂｊ＋１）、Ｆ（Ａｉ＋１，Ｂｊ＋１）は双１次曲面の４隅の点（双１次Ｂｅｚｉｅｒ曲面の制御点）であり、関数Ｆ（Ａｉ，Ｂｊ）における１つおきの値である。
【０１０２】
近似評価に用いる点は（数１５）で表される。
【０１０３】
【数１５】
【０１０４】
（数１５）において、ｇｉ，ｊは５成分を持つ行列で、Ｆ（Ａｉ−１，Ｂｊ−１）、Ｆ（Ａｉ＋１，Ｂｊ−１）、Ｆ（Ａｉ−１，Ｂｊ＋１）、Ｆ（Ａｉ＋１，Ｂｊ＋１）と、パラメータ（０.５，０）、（０，０.５）、（０.５，０.５）、（１，０.５）、（０.５，１）の位置を持つ双１次曲面上の点との差を成分に持つ。
【０１０５】
行列ｇｉ，ｊの内積を（数１６）で表す。
【０１０６】
【数１６】
【０１０７】
（数１５）は、厳密にいえば，Ｆ（Ａｉ，Ｂｊ−１）とＦａｉ，ｊ（０.５，０）、Ｆ（Ａｉ−１，Ｂｊ）とＦａｉ，ｊ（０，０.５）、Ｆ（Ａｉ，Ｂｊ）とＦａｉ，ｊ（０.５，０.５）、Ｆ（Ａｉ＋１，Ｂｊ）とＦａｉ，ｊ（１.０，０.５）、Ｆ（Ａｉ，Ｂｊ＋１）とＦａｉ，ｊ（０.５，１.０）が垂点（距離）ではない。垂点は、例えば、Ｆ（Ａｉ，Ｂｊ）と、Ｆａｉ，ｊ（ｕ，ｖ）、Ｆａｉ，ｊ（ｕ，ｖ）のパラメータ（ｕ，ｖ）による偏微分で定義される値で、非線形の式であり、解くには初期値を与えて収束計算することが必要である。そのため、複雑な計算が必要である。
【０１０８】
本発明では、分割数増加による近似評価のために使用することを前提として（数１６）を距離の二乗としている。（数１６）のｇｉ，ｊの内積（ｉ，ｊ）を（０，０）から（分割数―１，分割数―１）まで１つおきに変化させ、当該内積の最大を最大Ａとしている。
【０１０９】
次に、最大Ａが最大Ｂより大きいか判定する（ステップ９０５）。ステップ９０５の判定結果が大きければ、最大Ｂに最大Ａを代入する更新をおこなう（ステップ９０７）。
【０１１０】
次に、分割数×２−１の値を分割数に代入し（ステップ９０８）、ステップ９０２から繰り返す。本処理は、最大Ａの関数値が分割数を増大するのに従い、一旦上昇してから下がっていく傾向を示すため、当該関数値の最大を基準とした近似度合いの判定をおこなうために用いている。
【０１１１】
このようにすることで、最初の近似度がよい場合に誤って処理が終了することを防ぐことができる。ステップ９０５の判定結果が小さければ、最大Ａを最大Ｂで除した値がしきい値より大きいかを判定する（ステップ９０６）。ステップ９０６の判定結果が大きいときは、ステップ９０５で分割数を増加させ、ステップ９０２から繰り返す。ステップ９０６の判定結果が小さければ、分割数が決定したので、当該分割数でＦ（Ａ，Ｂ）が最低となる角度を抽出する（ステップ３１６）。ステップ３１６は図３のステップ３１６と同様の処理をおこなう。
【０１１２】
以上のように、ステップ９０１からステップ３１６の処理をおこなうことにより、角度（Ａ，Ｂ）の初期値を決定する。
【０１１３】
図８（ａ）は、分割数が増加する様子を示したイメージ図である。最初に３×３のＦ（Ａ，Ｂ）を定め、次に５×５、次に９×９のように間に１個づつＦ（Ａ，Ｂ）を増加させていく。従って、ｉ番目の分割数は２Ｎｉ−１（一つ前の分割数）−１となる。
【０１１４】
このようにＦ（Ａ，Ｂ）を増加させていくため、分割数は、３、５、９、１７、３３、６５、１２９、２５７、５１３、１０２５、・・・と増えていく。角度の範囲は０〜３６０なので、１０回目で約０．３５°の間隔となる。
【０１１５】
分割数が増加すると計算時間が増大するため、角度間隔と角度精度向上部１０７の収束度合いとの関連を実験したところ、５から６回繰り返したときに安定したので、繰り返しの回数は余裕をみて１０回に抑えている。
【０１１６】
図８（ｂ）は、Ｆ（Ａ，Ｂ）とＦａｉ，ｊ（ｕ，ｖ）の評価位置の関係を示している。Ｆ（Ａ，Ｂ）とＦａｉ，ｊ（ｕ，ｖ）は、（数１４）の定義からもわかるように、４隅の点が一致する。そのため、評価位置を中間位置にしている。
【０１１７】
図８（ｃ）は、分割数の増加に従い、最大Ａ／最大Ｂがどのように変化するかを示した特性図である。
【０１１８】
図８（ｃ）からわかるように、評価値は分割数が５か９で最大に達し、以後０に近づいていく。この結果より、０．１（１０％）をしきい値の初期値としている。しきい値は、最大値からの分散を定義し、その分散範囲以下になる値をしきい値としてもかまわない。この場合、３σや６σなどの位置をしきい値としてもかまわない。
【０１１９】
次に、角度精度向上部１０７で複数の初期角度をそれぞれ高精度化する。初期角度を高精度化するために（数１２）を角度（Ａ，Ｂ）で偏微分した値が０になる位置、すなわち（数１２）の極値を算出する。しかしながら、（数１２）は三角関数で記述された式であることから、角度（Ａ，Ｂ）での２回以上の偏微分、相互偏微分が無限に作成できる。そのため、図４（ａ）、図４（ｂ）、図４（ｃ）に示したように、極値を複数持つ複雑な関数となる。
【０１２０】
そこで、（数１２）を角度（Ａ，Ｂ）で偏微分した値が０になる式を１次Ｔａｙｌｏｒ展開して０とする式に近似し、（数１７）に示す連立一次方程式を作成する。
【０１２１】
（数１７）は初期角度（Ａ，Ｂ）を与えて解くことにより角度更新量（δＡ，δＢ）を算出する、所謂Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法である。（数１７）は複雑な関数であるため、良い初期角度を設定しないと解が発散してしまう。そのため。初期角度算出部１０６で初期角度を定めている。
【０１２２】
【数１７】
【０１２３】
（数１７）において各行列成分Ｍａは（数１８）である。
【０１２４】
【数１８】
【０１２５】
（数１８）においてｆｉ及びｆｉの平均、ｆｉ及びｆｉの角度（Ａ，Ｂ）による偏微分と当該偏微分の平均は、（数１９）で表せる。
【０１２６】
【数１９】
【０１２７】
（数１９）において、投影した計測点Ｑｉ及び投影した球の中心Ｏ、投影した計測点Ｑｉの角度（Ａ，Ｂ）による偏微分は（数２０）で表せる。
【０１２８】
【数２０】
【０１２９】
（数２０）において、計測点Ｐｉを投影するためのｘ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２１）で表せる。
【０１３０】
【数２１】
【０１３１】
（数２０）において、計測点Ｐｉを投影するためのｙ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２２）で表せる。
【０１３２】
【数２２】
【０１３３】
（数２１）に示すように、計測点Ｐｉを投影するためのｘ軸方向方向ベクトルの角度（Ａ）による偏微分、相互偏微分、２回以上の偏微分はすべて０ベクトルである。そのため、（数２０）は（数２３）に示すように簡単化される。また、（数２１）は（数２４）のように簡単化される。
【０１３４】
【数２３】
【０１３５】
【数２４】
【０１３６】
以上、数式により（数１７）の解き方を説明したが、次に、図３に示す角度精度向上部１０７のフローチャートを用いて投影したときに円となる角度（Ａ，Ｂ）を高精度化する過程を説明する。
【０１３７】
図３のステップ３２１からステップ３３０は、（数１７）を解く過程の一実施例を示すフローチャートである。尚、初期角度算出部１０６で求めた角度は複数存在するので、角度精度向上部１０７を角度の個数分繰り返して、当該角度毎に高精度化する。
【０１３８】
最初に初期角度算出部１０６で求めた角度を角度初期値に設定する（ステップ３２１）。
【０１３９】
次に、角度初期値に従い（数２４）で投影面のｘ軸方向のベクトルｔ、当該ベクトルｔの角度Ｂによる偏微分値ｔｂと２回偏微分値ｔｂｂ、（数２２）で投影面のｙ軸方向のベクトルｕ、当該ベクトルｕの角度Ａによる偏微分値ｕａと２回偏微分値ｕａａ、角度Ａで微分してから角度Ｂで微分する相互偏微分値ｕａｂ、当該ベクトルｕの角度Ｂによる偏微分値ｕｂ、２回偏微分値ｕｂｂを求める。次に、（数２３）で計測点Ｐｉを投影したベクトルＱｉ、当該ベクトルの角度Ａによる偏微分Ｑｉ，ａ、２回偏微分Ｑｉ，ａａ、角度Ａで微分後に角度Ｂで微分する相互偏微分Ｑｉ，ａｂ、当該ベクトルの角度Ｂによる偏微分値Ｑｉ，ｂ、２回偏微分値Ｑｉ，ｂｂ、を求める（ステップ３２２）。
【０１４０】
次に、（数１９）のｆｉにより投影した計測点と中心点の差ベクトルの内積で二乗距離を作成する（ステップ３２３）。
【０１４１】
次に、（数１９）のｆｉ，ａにより当該二乗距離の角度Ａによる偏微分、ｆｉ，ａａにより角度Ａによる２回偏微分、ｆｉ，ａｂにより角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分、ｆｉ，ｂにより当該二乗距離の角度Ｂによる偏微分、ｆｉ，ｂｂにより角度Ｂによる２回偏微分を作成する（ステップ３２４）。
【０１４２】
次に、（数１９）のｆａｖｒにより二乗距離の平均、ｆａｖｒ，ａにより二乗距離の平均の角度Ａによる偏微分、ｆａｖｒ，ａａにより二乗距離の平均の角度Ａによる２回偏微分、ｆａｖｒ，ａｂにより二乗距離の平均の角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分、ｆａｖｒ，ｂにより当該二乗距離の平均の角度Ｂによる偏微分、ｆａｖｒ，ｂｂにより二乗距離の平均の角度Ｂによる２回偏微分を作成する（ステップ３２５）。
【０１４３】
次に、（数１８）により二乗距離と二乗距離の平均の差の二乗和、当該二乗和の角度Ａによる偏微分、角度Ａによる２回偏微分、角度Ａで微分してから角度Ｂで微分する相互偏微分、角度Ｂによる偏微分、角度Ｂによる２回偏微分を作成する（ステップ３２６）。
【０１４４】
次に、ステップ３２６で求めた値を連立一次方程式の係数に代入する（ステップ３２７）。
【０１４５】
次に、当該連立一次方程式を解き、更新量（δＡ，δＢ）を算出する（ステップ３２８）。当該連立一次方程式を解くには一般的な連立一次方程式の解法、例えば、掃き出し法やＬＵ分解を使えば良い。また、２行２列なので解析的に解いてもかまわない。
【０１４６】
次に、更新量（δＡ，δＢ）が０とみなせるか判定する（ステップ３２９）。判定方法は、更新量の成分（δＡ，δＢ）の絶対値がそれぞれ０とみなせるかを判定する方法のほかに、更新量の内積（δＡ×δＡ＋δＢ×δＢ）が０とみなせるかを判定してもかまわない。本実施例では、更新量の内積が０とみなせるかを判定している。
【０１４７】
更新量が０とみなせない場合、角度更新をおこなう（ステップ３３０）。角度更新は、角度初期値（Ａ，Ｂ）に更新量（δＡ，δＢ）を加算することでおこなう。更新した角度初期値（Ａ，Ｂ）でステップ３２２から繰り返す。更新量が０とみなせる場合、処理を終了する。
【０１４８】
以上のステップ３２１からステップ３３０を繰り返すことで、球の中心を円筒中心とした場合の座標変換角度（Ａ，Ｂ）を定めることができる。実際に計測点を与え実験してみたところ、初期角度として初期角度算出部１０６で算出した角度を与えたとき、５から６回で収束した。
【０１４９】
しかしながら、初期角度として任意の値、例えば（０°、０°）などを指定した場合、図４で示すグラフの（０°、０°）近傍の極値が定まってしまい、極大値が解になってしまう現象が確認されており、初期角度算出部１０６が有効であることが確認されている。尚、特異な場合は球中心算出部１０４で除外すること、計測物が円筒面に限定されていることから角度精度向上部１０７が有効に動作することを確認している。
【０１５０】
次に、円筒中心・角度算出部１０８で精度向上した角度と円筒中心を、更に高精度化する手法を説明する。
【０１５１】
当該角度と円筒中心を高精度化するために（数１２）を角度（Ａ，Ｂ）および中心（Ｏｘ，Ｏｙ）で偏微分した値が０になる位置、即ち、（数１２）の極値を算出する。しかしながら、（数１２）は、三角関数で記述された式であることから、極値を複数持つ複雑な関数となる。
【０１５２】
そこで、（数１２）を角度（Ａ，Ｂ）及び円筒中心（Ｏｘ，Ｏｙ）で偏微分した値が０になる式を１次Ｔａｙｌｏｒ展開して０とする式に近似し、（数２５）に示す連立一次方程式を作成する。
【０１５３】
（数２５）は、初期角度（Ａ，Ｂ）および初期円筒中心（Ｏｘ，Ｏｙ）を与えて解くことにより更新量（δＡ，δＢ，δＯｘ，δＯｙ）を算出する。いわゆるＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法である。（数２５）は複雑な関数であるため、良い初期角度を設定しないと解が発散してしまう。そのため、角度精度向上部で高精度化した初期角度、当該初期角度で球の中心を座標変換した初期円筒中心を定めている。
【０１５４】
【数２５】
【０１５５】
（数２５）において各行列成分Ｍｂは（数２６）である。
【０１５６】
【数２６】
【０１５７】
（数２６）においてｆｉ及びｆｉの平均、ｆｉ及びｆｉの角度（Ａ，Ｂ）による偏微分、ｆｉ及びｆｉの円筒中心（Ｏｘ，Ｏｙ）による偏微分、と当該偏微分の平均は、（数２７）で表せる。
【０１５８】
【数２７】
【０１５９】
（数２７）において、投影した計測点Ｑｉおよび投影した球の中心Ｏ、投影した計測点Ｑｉの角度（Ａ，Ｂ）は（数２３）で表せる。（数２３）において偏微分した値が（数２７）に反映されているもの、及び偏微分または相互偏微分して０ベクトルになるものは記述していない。
【０１６０】
（数２７）において、計測点Ｐｉを投影するためのｘ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２４）で表せる。（数２４）において偏微分または相互偏微分して０ベクトルになるものは記述していない。また、（数２７）において、計測点Ｐｉを投影するためのｙ軸方向方向ベクトル、当該ベクトルの角度（Ａ，Ｂ）による偏微分は（数２２）で表せる。（数２２）において偏微分または相互偏微分して０ベクトルになるものは記述していない。
【０１６１】
（数２７）に示すように、円筒中心（Ｏｘ、Ｏｙ）による相互偏微分、２回上の偏微分はすべて０か２である。そのため、（数２６）は（数２８）、（数２７）は（数２９）に示すように簡単化される。
【０１６２】
【数２８】
【０１６３】
【数２９】
【０１６４】
以上、数式により（数２５）の解き方を説明したが、次に、図３に示す円筒中心・角度算出部１０８のフローチャートを用いて投影したときに円となる角度（Ａ，Ｂ）、中心点（Ｏｘ，Ｏｙ）を高精度化する過程を説明する。
【０１６５】
図３のステップ３３１からステップ３４０は（数２５）を解く過程の一実施例を示すフローチャートである。なお、角度精度向上部１０７で精度向上した角度は複数存在するので、円筒中心・角度算出部１０８を角度の個数分繰り返して、当該角度毎に角度と中心点を高精度化する。
【０１６６】
最初に、角度精度向上部１０７で求めた角度を角度初期値に設定するとともに、球の中心点を当該角度初期値で座標変換した値の（ｘ、ｙ）成分を中心点初期値に設定する（ステップ３３１）。
【０１６７】
次に、角度初期値と中心点初期値に従い（数２４）で投影面のｘ軸方向のベクトルｔ、当該ベクトルｔの角度Ｂによる偏微分値ｔｂと２回偏微分値ｔｂｂ、（数２２）で投影面のｙ軸方向のベクトルｕ、当該ベクトルｕの角度Ａによる偏微分値ｕａと２回偏微分値ｕａａ、角度Ａで微分してから角度Ｂで微分する相互偏微分値ｕａｂ、当該ベクトルｕの角度Ｂによる偏微分値ｕｂ、２回偏微分値ｕｂｂを求める。次に、（数２３）で計測点Ｐｉを投影したベクトルＱｉ、当該ベクトルの角度Ａによる偏微分Ｑｉ，ａ、２回偏微分Ｑｉ，ａａ、角度Ａで微分後に角度Ｂで微分する相互偏微分Ｑｉ，ａｂ、当該ベクトルの角度Ｂによる偏微分値Ｑｉ，ｂ、２回偏微分値Ｑｉ，ｂｂ、を求める（ステップ３３２）。
【０１６８】
次に、（数２９）のｆｉにより計測点と中心点の差ベクトルの内積で二乗距離を作成する（ステップ３３３）。
【０１６９】
次に、（数２９）のｆｉ，ａにより当該二乗距離の角度Ａによる偏微分、ｆｉ，ａａにより角度Ａによる２回偏微分、ｆｉ，ａｂにより角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分、ｆｉ，ａｏｘにより角度Ａで偏微分してから中心Ｏｘで偏微分する相互偏微分、ｆｉ，ａｏｙにより角度Ａで偏微分してから中心Ｏｙで偏微分する相互偏微分、ｆｉ，ｂにより当該二乗距離の角度Ｂによる偏微分、ｆｉ，ｂｂにより角度Ｂによる２回偏微分、ｆｉ，ｂｏｘにより角度Ｂで偏微分してから中心Ｏｘで偏微分する相互偏微分、ｆｉ，ｂｏｙにより角度Ｂで偏微分してから中心Ｏｙで偏微分する相互偏微分、ｆｉ，ｏｘにより中心Ｏｘでの偏微分、ｆｉ，ｏｙにより中心Ｏｙでの偏微分を作成する（ステップ３３４）。
【０１７０】
次に、（数２９）のｆａｖｒ，ａにより当該二乗距離の角度Ａによる偏微分の平均、ｆａｖｒ，ａａにより角度Ａによる２回偏微分の平均、ｆａｖｒ，ａｂにより角度Ａで偏微分してから角度Ｂで偏微分する相互偏微分の平均、ｆａｖｒ，ａｏｘにより角度Ａで偏微分してから中心Ｏｘで偏微分する相互偏微分の平均、ｆａｖｒ，ａｏｙにより角度Ａで偏微分してから中心Ｏｙで偏微分する相互偏微分の平均、ｆａｖｒ，ｂにより当該二乗距離の角度Ｂによる偏微分の平均、ｆａｖｒ，ｂｂにより角度Ｂによる２回偏微分の平均、ｆａｖｒ，ｂｏｘにより角度Ｂで偏微分してから中心Ｏｘで偏微分する相互偏微分の平均、ｆａｖｒ，ｂｏｙにより角度Ｂで偏微分してから中心Ｏｙで偏微分する相互偏微分の平均、ｆａｖｒ，ｏｘにより中心Ｏｘでの偏微分の平均、ｆａｖｒ，ｏｙにより中心Ｏｙでの偏微分の平均を作成する（ステップ３３５）。
【０１７１】
次に、（数２８）により二乗距離と二乗距離の平均の差の二乗和、当該二乗和の角度Ａによる偏微分、角度Ａによる２回偏微分、角度Ａで微分してから角度Ｂで微分する相互偏微分、角度Ａで微分してから中心Ｏｘで微分する相互偏微分、角度Ａで微分してから中心Ｏｙで微分する相互偏微分、角度Ｂによる偏微分、角度Ｂによる２回偏微分、角度Ｂで微分してから中心Ｏｘで微分する相互偏微分、角度Ｂで微分してから中心Ｏｙで微分する相互偏微分、中心Ｏｘでの偏微分、中心Ｏｙでの偏微分を作成する（ステップ３３６）。
【０１７２】
次に、ステップ３３６で求めた値を連立一次方程式の係数に代入する（ステップ３３７）。
【０１７３】
次に、当該連立一次方程式を解き、更新量（δＡ，δＢ）を算出する（ステップ３３８）。当該連立一次方程式を解くには一般的な連立一次方程式の解法、例えば、掃き出し法やＬＵ分解を使えば良い。
【０１７４】
次に、更新量（δＡ，δＢ，δＯｘ，δＯｙ）が０とみなせるか判定する（ステップ３３９）。判定方法は、更新量の成分（δＡ，δＢ，δＯｘ，δＯｙ）の絶対値がそれぞれ０とみなせるかを判定する方法のほかに、更新量の内積（δＡ×δＡ＋δＢ×δＢ＋δＯｘ×δＯｘ＋δＯｙ×δＯｙ）が０とみなせるかを判定してもかまわない。本実施例では、更新量の内積が０とみなせるかを判定している。
【０１７５】
更新量が０とみなせない場合、角度・中心点更新をおこなう（ステップ３４０）。更新は、初期値（Ａ，Ｂ，Ｏｘ，Ｏｙ）に更新量（δＡ，δＢ，δＯｘ，δＯｙ）を加算することでおこなう。更新した初期値（Ａ，Ｂ，Ｏｘ，Ｏｙ）でステップ３３２から繰り返す。更新量が０とみなせる場合、処理を終了する。
【０１７６】
以上のステップ３３１からステップ３４０を繰り返すことで、円筒の中心（Ｏｘ、Ｏｙ）と角度（Ａ，Ｂ）を定めることができる。
【０１７７】
実際に計測点を与え実験してみたところ、初期角度として角度精度向上部１０７で算出した角度、初期中心として球中心算出部１０５で作成した球中心を当該角度で座標変換した中心を与えたとき、５から６回で収束した。
【０１７８】
尚、初期角度算出部１０６、角度精度向上部１０７の効果を確認するため、実際に計測点を与え、分割ピッチを３６０°／Ｎ、Ｎは１から１９まで変化させ、角度精度向上部１０７をスキップした場合としない場合で円筒半径の違いを確認した。
【０１７９】
角度精度向上部１０７をスキップした場合、Ｎが１５より大きくなった場合に円筒半径が一定値となった。これは、２４°ピッチに相当する。これに対し、角度精度向上部１０７をスキップしない場合、Ｎが７より大きい場合に一定値となった。これは約５１．４２°ピッチに相当する。よって、角度精度向上部１０７は角度ピッチを２倍まで大きくできることがわかった。このことから、現在、角度ピッチは余裕をみて２０°ピッチとしている。
【０１８０】
しかしながら、角度ピッチは与える計測点によって変化する可能性があるので記憶装置２０４にあらかじめ格納しておき、適宜変更可能な構成としている。また、初期角度算出部１０６の計算量は角度ピッチを２倍にすると約４倍に増加するため、角度精度向上部１０７が計算量を減少させることがわかった。計算量の減少は計算時間の減少に繋がるため、その効果は大きいことがわかった。
【０１８１】
図３のフローチャートの説明に戻る。円筒中心・角度算出部１０８の処理が終了すると、円筒データ格納部１０９の処理を実行する。円筒データ格納部１０９は、円筒中心・角度算出部１０８で作成した円の中心Ｏ、角度を記憶装置１０２に格納する。
【０１８２】
ここで、円の中心Ｏは２次元座標であるため、Ｚ座標を０とし、（数３）に示す座標変換の逆変換をおこなえば、円筒の中心点が得られる。また、方向ベクトル（０，０，１）を同様に逆変換すれば、円筒の中心軸が得られる。半径は、（数２９）で作成したｆａｖｒの平方根である。
【０１８３】
円筒データが計測点を含むようにするには、中心点と計測点の差ベクトルと、中心軸の内積を作成して最低値と最大値を作成し、最低値が円筒の底面になるよう中心点を中心軸方向に移動すれば良い。そのとき、円筒の高さは最大値−最小値である。
【０１８４】
以上、図１〜図５を用いて本発明を説明したが、最後に図６を用いて本発明の方式を説明する。図６は、本発明の各処理に応じて作成される形状のイメージ図である。
【０１８５】
図６では、最初に計測点群が与えられる。次に、計測点群に球をあてはめ、球中心を作成する。次に、球中心を円筒の中心と仮定し、投影したときに円との分散が最小となる角度を少なくとも一つ作成する。次に当該角度を高精度化し、最後に、当該角度および円筒の中心点を高精度化して円筒あてはめを終了する。
【０１８６】
以上により、従来の問題点である、あてはめ時間増大、別解に収束する問題を解決する。
【０１８７】
図９は、本発明で説明した球、円筒の３次元モデルを別々、及び重ねて表示した結果であり、最初に計測点から球をあてはめ、仮円筒と円筒をあてはめていった結果を記している。また、最下段は計測点、球、円筒を重ねて表示した結果である。計測点と球、円筒の位置関係が判別できる。
【０１８８】
図１０は、何回目の計測がどの程度の時間なのかを示す図である。
【０１８９】
各処理は１０００回繰り返しているので、左端の目盛の単位はｍｓｅｃである。円弧溝の計算時間が多いのは、円筒３、円筒４の計測点数が８であるのと比較して、円弧溝の計測点数が１５０点になっているためである。これにより、十分高速な処理ができていることが理解できる。
【符号の説明】
【０１９０】
１０１…計測機器、１０２…記憶装置、１０３…円筒データ抽出部、１０４…計測点・収束判定値・角度ピッチ読取部、１０５…球中心算出部、１０６…初期角度算出部、１０７…角度精度向上部、１０８…円筒中心・角度算出部、１０９…円筒データ格納部、２００…計算機、２０１…通信装置、２０２…入力装置、２０３…表示装置、２０４…記憶装置、２０５…メモリ、２０６…中央処理装置。
【特許請求の範囲】
【請求項１】
円筒物の計測点から円筒の中心点、中心軸、半径を算出する方法において、
当該計測点と球表面の距離の二乗の総和が最小になる球の中心を算出し、かつ、当該中心を仮の円筒中心とし、当該計測点を平面に投影するための２つの角度を変化させ、投影した計測点と投影した仮の円筒中心と円周との距離の二乗の総和を算出し、当該総和が最小になる２つの角度の組を少なくとも１組以上抽出し、抽出した角度毎に、当該角度で仮の円筒中心と計測点を投影し、投影した仮の円筒中心を中心とする円の円周と、投影した計測点の距離の二乗の総和が最小になるよう抽出した角度を修正し、修正した角度と投影した仮の円筒中心を、円周と投影した計測点の距離の二乗の総和が最小になるよう、当該角度と当該中心を修正して円筒の中心点、中心軸、半径を算出することを特徴とする計測点からの円筒情報の抽出方法。
【請求項２】
請求項１に記載の計測点からの円筒情報の抽出方法において、
当該中心を仮の円筒中心とし、計測点を平面に投影するための２つの角度を変化させ、投影した計測点と投影した仮の円筒中心と、円周との距離の二乗の総和を算出し、当該総和が最小になる２つの角度の組を少なくとも１組以上抽出する際に、２つの角度を離散化するための分割数を３、当該分割数で離散化した二つの角度における円周との距離の二乗の総和の最大を格納するための距離Ｍを０、繰り返し終了を判定するためのしきい値を設定し、更に、２つの角度を離散化した格子リストを作成し、当該格子リストの２つの角度毎に、二乗距離と当該二乗距離の平均の差の二乗和を算出し、算出した二乗和と、当該二乗和を用いて近似曲面を定め、当該二乗和と当該近似曲面の距離の二乗の最大を求め、当該最大が距離Ｍより大きいかを判定し、当該最大が距離Ｍより大きい場合に、距離Ｍを当該最大に更新し、分割数を当該分割数×２−１に更新し、２つの角度を離散化した格子リストを作成する処理から繰り返し、当該最大が距離Ｍより小さい場合に、当該最大を距離Ｍで除した値がしきい値より大きいかを判定し、当該最大を距離Ｍで除した値がしきい値より大きい場合に、分割数を当該分割数×２−１に更新し、２つの角度を離散化した格子リストを作成する処理から繰り返し、当該最大が距離Ｍ以下の場合に、当該格子リストの２つの角度毎に算出した当該二乗和が最低となる角度を少なくとも１組以上抽出することを特徴とする計測点からの円筒情報の抽出方法。
【請求項３】
請求項１又は２に記載の計測点からの円筒情報の抽出方法を用いたプログラム。
【請求項４】
請求項１又は２に記載の計測点からの円筒情報の抽出方法を用いたプログラムを格納した記憶媒体。
【請求項５】
請求項３に記載のプログラムを組み込んで実行する計算機を用いた装置。
【請求項１】
円筒物の計測点から円筒の中心点、中心軸、半径を算出する方法において、
当該計測点と球表面の距離の二乗の総和が最小になる球の中心を算出し、かつ、当該中心を仮の円筒中心とし、当該計測点を平面に投影するための２つの角度を変化させ、投影した計測点と投影した仮の円筒中心と円周との距離の二乗の総和を算出し、当該総和が最小になる２つの角度の組を少なくとも１組以上抽出し、抽出した角度毎に、当該角度で仮の円筒中心と計測点を投影し、投影した仮の円筒中心を中心とする円の円周と、投影した計測点の距離の二乗の総和が最小になるよう抽出した角度を修正し、修正した角度と投影した仮の円筒中心を、円周と投影した計測点の距離の二乗の総和が最小になるよう、当該角度と当該中心を修正して円筒の中心点、中心軸、半径を算出することを特徴とする計測点からの円筒情報の抽出方法。
【請求項２】
請求項１に記載の計測点からの円筒情報の抽出方法において、
当該中心を仮の円筒中心とし、計測点を平面に投影するための２つの角度を変化させ、投影した計測点と投影した仮の円筒中心と、円周との距離の二乗の総和を算出し、当該総和が最小になる２つの角度の組を少なくとも１組以上抽出する際に、２つの角度を離散化するための分割数を３、当該分割数で離散化した二つの角度における円周との距離の二乗の総和の最大を格納するための距離Ｍを０、繰り返し終了を判定するためのしきい値を設定し、更に、２つの角度を離散化した格子リストを作成し、当該格子リストの２つの角度毎に、二乗距離と当該二乗距離の平均の差の二乗和を算出し、算出した二乗和と、当該二乗和を用いて近似曲面を定め、当該二乗和と当該近似曲面の距離の二乗の最大を求め、当該最大が距離Ｍより大きいかを判定し、当該最大が距離Ｍより大きい場合に、距離Ｍを当該最大に更新し、分割数を当該分割数×２−１に更新し、２つの角度を離散化した格子リストを作成する処理から繰り返し、当該最大が距離Ｍより小さい場合に、当該最大を距離Ｍで除した値がしきい値より大きいかを判定し、当該最大を距離Ｍで除した値がしきい値より大きい場合に、分割数を当該分割数×２−１に更新し、２つの角度を離散化した格子リストを作成する処理から繰り返し、当該最大が距離Ｍ以下の場合に、当該格子リストの２つの角度毎に算出した当該二乗和が最低となる角度を少なくとも１組以上抽出することを特徴とする計測点からの円筒情報の抽出方法。
【請求項３】
請求項１又は２に記載の計測点からの円筒情報の抽出方法を用いたプログラム。
【請求項４】
請求項１又は２に記載の計測点からの円筒情報の抽出方法を用いたプログラムを格納した記憶媒体。
【請求項５】
請求項３に記載のプログラムを組み込んで実行する計算機を用いた装置。
【図１】
【図２】
【図３】
【図４（ａ）】
【図４（ｂ）】
【図４（ｃ）】
【図５（ａ）】
【図５（ｂ）】
【図６】
【図７】
【図８（ａ）】
【図８（ｂ）】
【図８（ｃ）】
【図９】
【図１０】
【図２】
【図３】
【図４（ａ）】
【図４（ｂ）】
【図４（ｃ）】
【図５（ａ）】
【図５（ｂ）】
【図６】
【図７】
【図８（ａ）】
【図８（ｂ）】
【図８（ｃ）】
【図９】
【図１０】
【公開番号】特開２０１２−１６８０８２（Ｐ２０１２−１６８０８２Ａ）
【公開日】平成２４年９月６日（２０１２．９．６）
【国際特許分類】
【出願番号】特願２０１１−３０６４４（Ｐ２０１１−３０６４４）
【出願日】平成２３年２月１６日（２０１１．２．１６）
【出願人】（０００００５１０８）株式会社日立製作所 (27,607)
【Ｆターム（参考）】
【公開日】平成２４年９月６日（２０１２．９．６）
【国際特許分類】
【出願日】平成２３年２月１６日（２０１１．２．１６）
【出願人】（０００００５１０８）株式会社日立製作所 (27,607)
【Ｆターム（参考）】
[ Back to top ]