"円筒物の計測点から円筒の中心点、中心軸、半径を算出する方法において、当該計測点と球表面の距離の二乗の総和が最小になる球の中心を算出し、かつ、当該中心を仮の円筒中心とし、当該計測点を平面に投影するための2つの角度を変化させ、投影した計測点と投影した仮の円筒中心と円周との距離の二乗の総和を算出し、当該総和が最小になる2つの角度の組を少なくとも1組以上抽出し、抽出した角度毎に、当該角度で仮の円筒中心と計測点を投影し、投影した仮の円筒中心を中心とする円の円周と、投影した計測点の距離の二乗の総和が最小になるよう抽出した角度を修正し、修正した角度と投影した仮の円筒中心を、円周と投影した計測点の距離の二乗の総和が最小になるよう、当該角度と当該中心を修正して円筒の中心点、中心軸、半径を算出することを特徴とする計測点からの円筒情報の抽出方法。請求項1に記載の計測点からの円筒情報の抽出方法において、当該中心を仮の円筒中心とし、計測点を平面に投影するための2つの角度を変化させ、投影した計測点と投影した仮の円筒中心と、円周との距離の二乗の総和を算出し、当該総和が最小になる2つの角度の組を少なくとも1組以上抽出する際に、2つの角度を離散化するための分割数を3、当該分割数で離散化した二つの角度における円周との距離の二乗の総和の最大を格納するための距離Mを0、繰り返し終了を判定するためのしきい値を設定し、更に、2つの角度を離散化した格子リストを作成し、当該格子リストの2つの角度毎に、二乗距離と当該二乗距離の平均の差の二乗和を算出し、算出した二乗和と、当該二乗和を用いて近似曲面を定め、当該二乗和と当該近似曲面の距離の二乗の最大を求め、当該最大が距離Mより大きいかを判定し、当該最大が距離Mより大きい場合に、距離Mを当該最大に更新し、分割数を当該分割数×2-1に更新し、2つの角度を離散化した格子リストを作成する処理から繰り返し、当該最大が距離Mより小さい場合に、当該最大を距離Mで除した値がしきい値より大きいかを判定し、当該最大を距離Mで除した値がしきい値より大きい場合に、分割数を当該分割数×2-1に更新し、2つの角度を離散化した格子リストを作成する処理から繰り返し、当該最大が距離M以下の場合に、当該格子リストの2つの角度毎に算出した当該二乗和が最低となる角度を少なくとも1組以上抽出することを特徴とする計測点からの円筒情報の抽出方法。請求項1又は2に記載の計測点からの円筒情報の抽出方法を用いたプログラム。請求項1又は2に記載の計測点からの円筒情報の抽出方法を用いたプログラムを格納した記憶媒体。請求項3に記載のプログラムを組み込んで実行する計算機を用いた装置。しかしながら、非特許文献1では、3点に限定していること及び3点で作成した平面上の円しか作成できないことから、複数の計測点に最も近い円の中心を求める場合には、適用が難しい。また、非特許文献2では、3次元空間で円筒の中心軸を定めるには、非特許文献3を用いて、3次元空間の複数の点を回転変換し、点から2つの成分を取り出すことにより平面に投影した後で適用する必要があり、この方法では、回転変換時の角度および投影後の中心点が変数になり、回転変換は三角関数を使用するため、角度で無限回の微分が可能である。即ち、極値を複数持つ複雑な関数となるので、複数の計測点に最も近い円の中心を求める場合には適用が難しい。更に、3の複数の点に近い球では、球ではなく円筒を定める場合には利用できない。本発明は上述の点に鑑みなされたもので、その目的とするところは、円筒表面の計測点から円筒データを確実に作成できると共に、従来方式に比べて、高精度な近似が可能になる計測点からの円筒情報の抽出方法を提供することにある。本発明の計測点からの円筒情報の抽出方法は、上記目的を達成するために、計測点が円筒物の表面を計測していることを条件として、全計測点から球の中心を抽出する手段と、全計測点および球の中心から角度の初期値を作成する手段と、当該角度の初期値の精度向上を行う手段と、球中心を精度向上した角度で投影変換し、初期円筒中心を定める手段と、初期円筒中心および精度向上した角度を初期値とし、円筒中心および角度を定める手段と、を用意し計測点から円筒の半径、中心軸からなる円筒データを抽出することを特徴とするものである。更に、全計測点および球の中心から角度の初期値を作成する場合に、角度を離散化するための分割数、当該分割数と当該分割数における近似曲面との差の二乗の最大Bの初期値、当該分割数における近似曲面との差の二乗の最大Aを最大Bで除した値で繰り返し終了を判定するための初期値を設定する手段と、2つの角度を当該分割数で離散化した格子リストを作成する手段と、当該格子リストの角度毎に、計測点と球の中心を座標変換して平面上に投影し、投影した計測点と投影した中心の二乗距離を作成し、当該二乗距離の平均を作成し、当該二乗距離と当該平均の差の二乗和を作成する手段と、当該格子リストに対応した二乗和から近似曲面を作成する手段と、当該格子リストに対応した二乗和と近似曲面との距離の二乗の最大Aを求める手段と、最大Aが最大Bより大きいかを判定する手段と、最大Bに最大Aを代入する手段と、最大Aを最大Bで除した値がしきい値より大きいかを判定する手段と、分割数を増加させる手段と、を用意し、全計測点および球の中心から角度の初期値を作成することを特徴とするものである。本発明の計測点からの円筒情報の抽出方法によれば、円筒表面の計測点から円筒データを確実に作成できると共に、従来方式に比べて、高精度な近似が可能になる。以下、本発明の計測点からの円筒情報の抽出方法の実施例を、図1乃至図10を用いて説明する。図1は、本発明の計測点からの円筒情報の抽出方法を実施する装置後述する計算機で実行するプログラムでもあるを示すものであり、該図において、101は円筒物の表面を3次元座標X、Y、Zで計測する計測機器、102は計測機器101で計測した計測点、収束判定に用いる0とみなす値、初期角度算出部106で使用する角度のピッチを保存する記憶装置、103は計算機で実行する円筒データ抽出部、104は記憶装置102から3次元の計測点Pi、収束判定に用いる0とみなす値、初期角度算出部106で使用する角度のピッチを読み取る計測点収束判定値角度ピッチ読取部、105は計測点Piから球の中心を算出する球中心算出部、106は記憶装置に格納されている角度ピッチで分割した角度で、計測点および105で算出した球中心を回転変換し、回転変換した計測点と球中心のベクトルの2成分を取り出すことにより、投影した球中心と投影した計測点を定め、投影した球中心と投影した計測点の距離を求めると共に、当該距離の平均を求め、当該距離と当該平均の差の二乗和を求め、当該二乗和が最小になる初期角度を算出する初期角度算出部、107は106で作成した初期角度の精度向上を行う角度精度向上部、108は、106で作成した投影後の球中心を初期円筒中心、107で作成した角度を初期角度とし、円筒中心、角度、半径を算出する円筒中心角度算出部、109は108で作成した円筒中心、角度、半径を記憶装置102に格納する円筒データ格納部である。図2は、図1の円筒データ抽出処理を実施する計算機を示すものであり、200は、円筒データ抽出部103を実行する計算機である。計算機200において、201は通信網で接続された他の計算機とデータのやり取りをおこなうための通信装置、202はマウスやキイボードなどのデータの入力をおこなう入力装置、203はグラフィックディスプレイなどの計算結果を文字やグラフなどで表示するための表示装置、204は計算結果を蓄積するFD等の記憶装置記憶媒体、205はメモリ、206は103の円筒データ抽出処理をおこなう中央処理装置である。図3は、本発明の円筒データ抽出処理の流れを示すフローチャートである。以下、図3のフローチャートと数式を交えながら、本発明の一実施例を述べる。尚、図3は、本発明のプログラムに相当する。最初に計測点収束判定値角度ピッチ読取部104が、記憶装置102から計測点データ、収束判定値、角度ピッチを読み取る。計測点データは、数5に示すようにPiでi番目の計測点を表し、x成分、y成分、z成分を持つ3次元空間上の点である。ここで、iは1から計測点数Nまで変化する。収束判定値は、例えば1.0E-20のように更新量を0とみなす値である。角度ピッチは、0から360°までの角度を離散化するときの幅であり、例えば20°とすれば、0°、20°、、340°、360°の角度に離散化する。次に、球中心算出部105がN個の計測点と球表面との距離の二乗の総和が最小になる球以下、最小二乗球というの中心点を算出する。以下、その算出法について述べる。最小二乗球の中心Osは、例えば数3で算出することができるが、本実施例では、初期角度算出部106、角度精度向上部107、円筒中心角度算出部108を説明する都合上、数6に示す方法を用いる。数6において、Osは球の中心、giは点Piと中心Osの距離の二乗半径の二乗、gavrはgiの平均である。即ち、各点Piと中心Osの半径の二乗と、当該半径の二乗の平均との差の二乗和が最小になるよう中心Osを決定する。数6は厳密に言えば、計測点と球表面の距離ではなく、距離の二乗を最小化する。しかしながら、距離の場合は、平方根を用いた関数とせねばならず、微分時に複雑さが増してしまう。そのため、距離の二乗を最小化している。この方法は一般的に最小二乗法を用いるときの常套手段である。尚、数3と数6の解が一致することは実験で確認済である。また、数6のようにgiとgavrの差を作成すると、gi中に同一値を持つ変数成分及び定数成分がある場合に0となるので、数値が大きくなることがないため、計算機処理する場合の丸め誤差を少なくする上で都合が良い。さらに、gavrが加算と定数による除算なので、giとgavrの次数が同一となり、次数が増加しないので関数の複雑さが増すことがない。数6を解くために、数6を球の中心Osの各成分Osx、Osy、Oszで偏微分した関数が0になる式を作成し、当該式を1次までTaylor展開して連立一次方程式を作成する。当該連立一次方程式を、初期中心Osを与えて解くことにより、中心更新量を算出する。本方法は、Newton-Raphson法と呼ばれている。当該連立一次方程式は、2回微分が0または定数となる比較的単純な関数であるため、任意の初期中心でも収束する。数7は当該連立一次方程式である。数3が4行4列であるのに対し、数7は3行3列である。このように数6で示すようにgiとgavrの差を求める方法を用いると、行数および列数を1行1列減少させることができるので、計算時間の短縮が可能なため、本実施例ではこの方法を採用している。数7において各行列成分Msは数8である。数8においてgi及びgavrは数6、giの中心Osの成分による微分及び平均は、数9及び数10で表せる。数10に示すように、giの中心Osの成分による2回偏微分は定数2か0である、従って、平均も定数2か0である。そのため、数8は数11のように簡単化される。また、数11で使用するgiの中心Osの成分による2回偏微分及び平均は、数9のみとなる。以上、数式により数6の解き方を説明したが、次に、図3に示す球中心算出部105のフローチャートを用いて球中心の算出過程を説明する。図3の301から310は、数6を解く過程の一実施例を示すフローチャートである。最初に球の中心点の初期値を例えば0,0,0のように与えるステップ301。当該初期値は0,0,0に限定されるものではなく、任意の値、例えば、100,100,100や-1,500,0などでもかまわないが、数値計算上は予想できる中心点を与えることが望ましい。本実施例では、予想できないこと、及び更新量がそのまま球の中心点となることから0,0,0を与えている。次に、計測点Piと中心点の差ベクトルを求めるステップ302。次に、数6のgiに示す差ベクトルの内積を求め、計測点数分の内積値を作成するステップ303。当該内積値は中心点と計測点Piの距離の二乗であることから、内積値を二乗距離と呼ぶ。次に、数9のgi,os,x、gi,os,y、gi,os,zに示す二乗距離の中心点による偏微分を作成するステップ304。次に、数6のgavr、数9のgavr,os,x、gavr,os,y、gavr,oszに示す二乗距離および二乗距離の中心点による偏微分の平均を作成するステップ305。次に、ステップ304で求めた偏微分およびステップ305で求めた偏微分の平均から、数11のMs,0,0、Ms,2,0、Ms,0,1、Ms,2,1、Ms,0,2、Ms,2,2に示す二乗距離の総和の2回偏微分を作成するとともに、ステップ303で求めた二乗距離とステップ304で求めた偏微分およびステップ305で求めた二乗距離の平均と二乗距離の偏微分の平均から数11のMs,3,0、Ms,3,1、Ms,3,2に示す二乗距離の総和の偏微分を作成するステップ306。次に、ステップ306で求めた2回偏微分と偏微分から数7に示す連立一次方程式の係数を作成するステップ307。次に、連立一次方程式を解き、中心点の更新量を求めるステップ308。連立一次方程式の解法は、例えば、掃き出し法やLU分解で解くことができる。また、3行3列なので解析的に解いてもかまわない。次に、中心点の更新量が104で読み込んだ0とみなせる値になっているかを判定するステップ309。判定時には、例えば、中心点の更新量δOsx,δOsy,δOszの絶対値を作成して,当該絶対値がそれぞれ0とみなせる値になるかどうかで確認する。そのほかに、更新量の内積δOsx×δOsxδOsy×δOsyδOsz×δOszを作成して0となるかどうかを判定してもかまわない。本実施例では、更新量の内積で判定している。ステップ309で0に近くないと判定された場合、中心点の更新を行うステップ310。中心点の更新は、現在の中心点Osx,Osy,Oszに、更新量δOsx,δOsy,δOszを加算すればよい。この加算した値を現在の中心点Osx,Osy,Oszとし、ステップ302から繰り返す。ステップ309で0に近いと判定された場合、球中心算出部105の処理を終了する。以上のステップ302からステップ309を繰り返すことで、球の中心点を定めることができる。実際には、球中心算出部105の処理は、殆ど1回で繰り返し計算が終了する。これは、giの中心Osの成分による偏微分が、1回偏微分が変数、2回の偏微分および相互偏微分が定数または0、3回以降の偏微分および相互偏微分が0であることから、数7に示す関数が線形関数であることに関連している。但し、数値計算誤差が含まれることも想定されるので、本実施例に示したように、ステップ302~ステップ309、ステップ310の処理を行うことが精度を高める上で望ましい。3次元で定義される点に球をあてはめる場合、数7に示すように、中心Osのみを求めることとなり回転角度には依存しない。このため、特異な場合を除き確実に中心点Osを定めることができるので、あてはめ時間を最小限に抑制することができる。ここで、特異な場合とは、例えば、3次元で定義された点が直線上に配置された場合である。本状態は、計測物が円筒面の場合に限定すると、円筒の中心軸に平行となるよう計測したときに生じる。このとき、数7で繰り返し計算をおこなうと無限に繰り返すか、または0除算が発生してしまう。そのため、繰り返し回数に上限を設け、上限を超えた場合、及び数7の連立一次方程式を解くときに0除算が発生するか確認し、エラー出力をおこなうことが望ましい。当該エラーが出力された場合、計測位置を変更する必要がある。そのため、数7の処理を行う計算機、及び計算の結果をモニタできる装置、例えばブザーや小型の表示装置を設け、計測結果をリアルタイムにモニタリングすることが望ましい。本構成を採用する場合、数7の計算処理が単純であり、繰り返し回数も少ないことから計算能力が劣っている計算機でも採用することができるとともに、計測位置の誤りがリアルタイムで手元確認できるため、すぐに計測位置を修正することができるので計測間違いが少なくなる。次に、初期角度算出部106で、円筒の中心軸となる角度の初期値を算出する。初期角度を算出するために数12を用いる。数12において、Qiは、計測点PiをY軸で角度Bだけ回転してからX軸で角度Aだけ回転した後のx、y成分である。Oは、球の中心OsをY軸で角度Bだけ回転してからX軸で角度Aだけ回転した後のx、y成分である。数12では、投影平面をx、y平面に限定し、回転をY軸回転後にX軸回転しているが、投影平面はx、y平面に限定されるものではなく、y、z平面、z、x平面でもかまわないし、回転は、Y軸回転後のZ軸回転、X軸回転後のY軸回転、X軸回転後のZ軸回転、Z軸回転後のX軸回転、Z軸回転後のY軸回転のいずれでもかまわない。重要なことは、初期角度算出部106、角度精度向上部107、円筒中心角度算出部108で回転順序を一致させておくことである。これを誤ると角度の高精度化、円筒の中心や角度の算出ができなくなるので注意が必要である。数12は三角関数が入るため、角度AおよびBで無限に偏微分、相互偏微分することが可能であり極値を複数持つ複雑な関数である。そのため、本願では数13で示すように、数12の角度A,Bを0~360°の中で離散化して数12の値を求め、最小値を持つ角度A,Bを探す構成とした。以上、数式により初期角度算出部106の処理を説明したが、次に図3に示す初期角度算出部106のフローチャートで初期角度の算出課程を説明する。図3において、311から316は、初期角度算出部106の処理の一実施例を示すフローチャートである。最初に数13に示す角度A,Bを、計測点収束判定値角度ピッチ読取104で読み込んだ角度ピッチに従って離散化するステップ311。本実施例では、角度ピッチを20°としている。20°という値は、テストデータを用いたとき角度精度向上部106で収束するピッチを実験して決定したものである。本結果については後述する。尚、初期角度算出部106において、円筒中心も変数とした場合、円筒中心が3成分毎に-∞~∞まで変化するため設定が難しい。また、変数が5個となり組合せ数が一気に増大するため計算量が増大してしまう。更に、変数5個それぞれに極値を持つため大変複雑になってしまう。そのため、球中心算出部105を用意して角度だけで数13を算出する構成としている。次に、計測点Pi、球中心算出手段105で作成した球中心を離散化した角度毎に座標変換するステップ312。座標変換の方法は、上述した43次元から2次元への座標変換で説明した処理でかまわない。次に、数12のfiに示すように座標変換した計測点と球の中心のx、y成分を取り出して差ベクトルを作成し、差ベクトルの内積を作成してx、y成分の二乗距離を作成するステップ313。次に、数12のfavrに示すように当該二乗距離の平均を作成するステップ314。次に、二乗距離と平均の差を作成し、当該差の二乗和を作成するステップ315。最後に、当該二乗和の最低値を取出し、最低値と一致する離散化した角度A,Bを抽出するステップ316。最低値を作成する方法としては、当該二乗和を小さい値から昇順に並べ替え最初の値を抽出する方法や、最初の値を最低値の初期値とし2番目から計測点数番目までの間で当該最低値より小さい値があらわれたとき置き換えるなどの方法を用いてもかまわない。並べ替えの方法としては、バブルソートやクイックソートなどのソーティングアルゴリズムを利用してもかまわない。ここで、最低値を持つ角度A,Bが複数存在することを考慮し、ステップ315の二乗和は小数点以下の有効桁数を考慮して整数化しておくことが望ましい。整数化のかわりに四捨五入した実数を使うことも考えられるが、実数値に丸め誤差が累積して微小に異なる場合があるため、整数化しておく方が良い。本実施例では、小数点以下1桁までを有効とし、実数値を10倍して切り捨てる方法を使用している。図4aは、計測点8点を与え、角度Aおよび角度Bをそれぞれ0から360°まで20°ピッチで分割して数13のFA,Bをプロットした3次元グラフの上面図、図4bは、図4aの斜視図、図4cは、角度Aを横軸、角度Bを縦軸としFA,Bの値を記入した図である。図4cでは、最低値を四角で囲んで示している。図4a、図4b、図4cより、FA,Bが極値を多数持つ複雑な関数であり、最低値が4個程度存在することがわかる。図5aの1~4は、球の中心点を複数求めた角度で座標変換した位置のx、y座標である。図5aに示すように、中心点が角度毎に象限を分けて算出されている。このように、中心点が分れて検出されることから、角度と中心を変数として繰り返し計算すると、中心の更新量に従い角度が大きく変化するため解が大きく変動することに繋がる。そのため、初期角度算出部106では球中心算出部105で求めた中心を円筒の中心とみなし、角度の初期値を定めているので、解が大きく変動するのを防ぐことができる。次に初期角度算出部106を実現する別の処理について、図7及び図8a、図8b、図8cを用いて説明する。図7は、初期角度算出部106の別の実施方法を示すフローチャートである。該図に示す如く、最初に分割数、最大B、しきい値を設定するステップ901。分割数は2つの角度A,Bを離散化するときの点数、最大Bは分割数を変数としたときの当該分割数における関数FA,Bと当該関数の近似曲面の距離の二乗の最大を示す関数の最大の初期値、しきい値は分割数の決定処理が終了したかどうかを表す値である。尚、図7は、本発明の別の実施方法のプログラムに相当する。本実施例では、分割数は3を設定しているが、関数FA,Bの複雑さに応じ、3より大きい値を設定してもかまわない。また、最大Bは0を設定しているが、関数FA,Bの最小が0となるからであり、当該最小が0以上となることが確実な場合、0より大きい値を設定してもかまわない。また、しきい値は0.1を設定しているが、角度精度向上部107の処理が発散する場合0~0.1の値を設定する必要がある。更に、角度精度向上部107の処理が確実に収束するなら、0.1~1の範囲に設定してもかまわない。本実施例でいろいろな計測点を与え、実験したところ、発散しなかったことから0.1に設定している。次に、角度Aの範囲0~360°、角度Bの範囲0~360°を当該分割数個に離散化するステップ902。本実施例では、関数FA,Bの傾向が不明なため、角度が等感覚になるような離散化をおこなっている。例えば、分割数が3なら、0,0、180,0、360,0、0,180、180,180、360,180、0,360、180,360、360,360の格子リストとなる。次に、格子リストの角度毎に関数FA,Bを作成するステップ903。例えば、分割数が3なら、F0,0、F180,0、F360,0、F0,180、F180,180、F360,180、F0,360、F180,360、F360,360を作成する。903の処理過程は、図2のステップ312、313、314、315と同様である。次に、格子リストの角度毎に関数FAi-1,Bj-1、FAi1,Bj-1、FAi-1,Bj1、FAi1,Bj1から近似曲面Fai,ju,vを作成し、FAi,BjとFai,j05,05の距離を作成し、当該距離の最大である最大Aを求めるステップ904。近似曲面Fai,ju,vは、m×n次の多項式曲面、Bezier曲面、B-Spline曲面、Coonz曲面など任意の形式でかまわない。また、これらの曲面を有理式化したものでもかまわない。本実施例では、処理を簡単にするため、双1次パッチを採用し、双1次パッチのu,v方向それぞれの中央で評価する構成としている。数14に双1次パッチによる近似曲面Fai,ju,vを示す。数14において、Fai,ju,vは、パラメータu,vを変数とする双1次曲面双1次Bezier曲面でもある。パラメータu,vはそれぞれ、0~1まで変化する変数。FAi-1,Bj-1、FAi1,Bj-1、FAi-1,Bj1、FAi1,Bj1は双1次曲面の4隅の点双1次Bezier曲面の制御点であり、関数FAi,Bjにおける1つおきの値である。近似評価に用いる点は数15で表される。数15において、gi,jは5成分を持つ行列で、FAi-1,Bj-1、FAi1,Bj-1、FAi-1,Bj1、FAi1,Bj1と、パラメータ05,0、0,05、05,05、1,05、05,1の位置を持つ双1次曲面上の点との差を成分に持つ。行列gi,jの内積を数16で表す。数15は、厳密にいえば,FAi,Bj-1とFai,j05,0、FAi-1,BjとFai,j0,05、FAi,BjとFai,j05,05、FAi1,BjとFai,j10,05、FAi,Bj1とFai,j05,10が垂点距離ではない。垂点は、例えば、FAi,Bjと、Fai,ju,v、Fai,ju,vのパラメータu,vによる偏微分で定義される値で、非線形の式であり、解くには初期値を与えて収束計算することが必要である。そのため、複雑な計算が必要である。本発明では、分割数増加による近似評価のために使用することを前提として数16を距離の二乗としている。数16のgi,jの内積i,jを0,0から分割数―1,分割数―1まで1つおきに変化させ、当該内積の最大を最大Aとしている。次に、最大Aが最大Bより大きいか判定するステップ905。ステップ905の判定結果が大きければ、最大Bに最大Aを代入する更新をおこなうステップ907。次に、分割数×2-1の値を分割数に代入しステップ908、ステップ902から繰り返す。本処理は、最大Aの関数値が分割数を増大するのに従い、一旦上昇してから下がっていく傾向を示すため、当該関数値の最大を基準とした近似度合いの判定をおこなうために用いている。このようにすることで、最初の近似度がよい場合に誤って処理が終了することを防ぐことができる。ステップ905の判定結果が小さければ、最大Aを最大Bで除した値がしきい値より大きいかを判定するステップ906。ステップ906の判定結果が大きいときは、ステップ905で分割数を増加させ、ステップ902から繰り返す。ステップ906の判定結果が小さければ、分割数が決定したので、当該分割数でFA,Bが最低となる角度を抽出するステップ316。ステップ316は図3のステップ316と同様の処理をおこなう。以上のように、ステップ901からステップ316の処理をおこなうことにより、角度A,Bの初期値を決定する。図8aは、分割数が増加する様子を示したイメージ図である。最初に3×3のFA,Bを定め、次に5×5、次に9×9のように間に1個づつFA,Bを増加させていく。従って、i番目の分割数は2Ni-1一つ前の分割数-1となる。このようにFA,Bを増加させていくため、分割数は、3、5、9、17、33、65、129、257、513、1025、と増えていく。角度の範囲は0~360なので、10回目で約0.35°の間隔となる。分割数が増加すると計算時間が増大するため、角度間隔と角度精度向上部107の収束度合いとの関連を実験したところ、5から6回繰り返したときに安定したので、繰り返しの回数は余裕をみて10回に抑えている。図8bは、FA,BとFai,ju,vの評価位置の関係を示している。FA,BとFai,ju,vは、数14の定義からもわかるように、4隅の点が一致する。そのため、評価位置を中間位置にしている。図8cは、分割数の増加に従い、最大A/最大Bがどのように変化するかを示した特性図である。図8cからわかるように、評価値は分割数が5か9で最大に達し、以後0に近づいていく。この結果より、0.110をしきい値の初期値としている。しきい値は、最大値からの分散を定義し、その分散範囲以下になる値をしきい値としてもかまわない。この場合、3σや6σなどの位置をしきい値としてもかまわない。次に、角度精度向上部107で複数の初期角度をそれぞれ高精度化する。初期角度を高精度化するために数12を角度A,Bで偏微分した値が0になる位置、すなわち数12の極値を算出する。しかしながら、数12は三角関数で記述された式であることから、角度A,Bでの2回以上の偏微分、相互偏微分が無限に作成できる。そのため、図4a、図4b、図4cに示したように、極値を複数持つ複雑な関数となる。そこで、数12を角度A,Bで偏微分した値が0になる式を1次Taylor展開して0とする式に近似し、数17に示す連立一次方程式を作成する。数17は初期角度A,Bを与えて解くことにより角度更新量δA,δBを算出する、所謂Newton-Raphson法である。数17は複雑な関数であるため、良い初期角度を設定しないと解が発散してしまう。そのため。初期角度算出部106で初期角度を定めている。数17において各行列成分Maは数18である。数18においてfi及びfiの平均、fi及びfiの角度A,Bによる偏微分と当該偏微分の平均は、数19で表せる。数19において、投影した計測点Qi及び投影した球の中心O、投影した計測点Qiの角度A,Bによる偏微分は数20で表せる。数20において、計測点Piを投影するためのx軸方向方向ベクトル、当該ベクトルの角度A,Bによる偏微分は数21で表せる。数20において、計測点Piを投影するためのy軸方向方向ベクトル、当該ベクトルの角度A,Bによる偏微分は数22で表せる。数21に示すように、計測点Piを投影するためのx軸方向方向ベクトルの角度Aによる偏微分、相互偏微分、2回以上の偏微分はすべて0ベクトルである。そのため、数20は数23に示すように簡単化される。また、数21は数24のように簡単化される。以上、数式により数17の解き方を説明したが、次に、図3に示す角度精度向上部107のフローチャートを用いて投影したときに円となる角度A,Bを高精度化する過程を説明する。図3のステップ321からステップ330は、数17を解く過程の一実施例を示すフローチャートである。尚、初期角度算出部106で求めた角度は複数存在するので、角度精度向上部107を角度の個数分繰り返して、当該角度毎に高精度化する。最初に初期角度算出部106で求めた角度を角度初期値に設定するステップ321。次に、角度初期値に従い数24で投影面のx軸方向のベクトルt、当該ベクトルtの角度Bによる偏微分値tbと2回偏微分値tbb、数22で投影面のy軸方向のベクトルu、当該ベクトルuの角度Aによる偏微分値uaと2回偏微分値uaa、角度Aで微分してから角度Bで微分する相互偏微分値uab、当該ベクトルuの角度Bによる偏微分値ub、2回偏微分値ubbを求める。次に、数23で計測点Piを投影したベクトルQi、当該ベクトルの角度Aによる偏微分Qi,a、2回偏微分Qi,aa、角度Aで微分後に角度Bで微分する相互偏微分Qi,ab、当該ベクトルの角度Bによる偏微分値Qi,b、2回偏微分値Qi,bb、を求めるステップ322。次に、数19のfiにより投影した計測点と中心点の差ベクトルの内積で二乗距離を作成するステップ323。次に、数19のfi,aにより当該二乗距離の角度Aによる偏微分、fi,aaにより角度Aによる2回偏微分、fi,abにより角度Aで偏微分してから角度Bで偏微分する相互偏微分、fi,bにより当該二乗距離の角度Bによる偏微分、fi,bbにより角度Bによる2回偏微分を作成するステップ324。次に、数19のfavrにより二乗距離の平均、favr,aにより二乗距離の平均の角度Aによる偏微分、favr,aaにより二乗距離の平均の角度Aによる2回偏微分、favr,abにより二乗距離の平均の角度Aで偏微分してから角度Bで偏微分する相互偏微分、favr,bにより当該二乗距離の平均の角度Bによる偏微分、favr,bbにより二乗距離の平均の角度Bによる2回偏微分を作成するステップ325。次に、数18により二乗距離と二乗距離の平均の差の二乗和、当該二乗和の角度Aによる偏微分、角度Aによる2回偏微分、角度Aで微分してから角度Bで微分する相互偏微分、角度Bによる偏微分、角度Bによる2回偏微分を作成するステップ326。次に、ステップ326で求めた値を連立一次方程式の係数に代入するステップ327。次に、当該連立一次方程式を解き、更新量δA,δBを算出するステップ328。当該連立一次方程式を解くには一般的な連立一次方程式の解法、例えば、掃き出し法やLU分解を使えば良い。また、2行2列なので解析的に解いてもかまわない。次に、更新量δA,δBが0とみなせるか判定するステップ329。判定方法は、更新量の成分δA,δBの絶対値がそれぞれ0とみなせるかを判定する方法のほかに、更新量の内積δA×δAδB×δBが0とみなせるかを判定してもかまわない。本実施例では、更新量の内積が0とみなせるかを判定している。更新量が0とみなせない場合、角度更新をおこなうステップ330。角度更新は、角度初期値A,Bに更新量δA,δBを加算することでおこなう。更新した角度初期値A,Bでステップ322から繰り返す。更新量が0とみなせる場合、処理を終了する。以上のステップ321からステップ330を繰り返すことで、球の中心を円筒中心とした場合の座標変換角度A,Bを定めることができる。実際に計測点を与え実験してみたところ、初期角度として初期角度算出部106で算出した角度を与えたとき、5から6回で収束した。しかしながら、初期角度として任意の値、例えば0°、0°などを指定した場合、図4で示すグラフの0°、0°近傍の極値が定まってしまい、極大値が解になってしまう現象が確認されており、初期角度算出部106が有効であることが確認されている。尚、特異な場合は球中心算出部104で除外すること、計測物が円筒面に限定されていることから角度精度向上部107が有効に動作することを確認している。次に、円筒中心角度算出部108で精度向上した角度と円筒中心を、更に高精度化する手法を説明する。当該角度と円筒中心を高精度化するために数12を角度A,Bおよび中心Ox,Oyで偏微分した値が0になる位置、即ち、数12の極値を算出する。しかしながら、数12は、三角関数で記述された式であることから、極値を複数持つ複雑な関数となる。そこで、数12を角度A,B及び円筒中心Ox,Oyで偏微分した値が0になる式を1次Taylor展開して0とする式に近似し、数25に示す連立一次方程式を作成する。数25は、初期角度A,Bおよび初期円筒中心Ox,Oyを与えて解くことにより更新量δA,δB,δOx,δOyを算出する。いわゆるNewton-Raphson法である。数25は複雑な関数であるため、良い初期角度を設定しないと解が発散してしまう。そのため、角度精度向上部で高精度化した初期角度、当該初期角度で球の中心を座標変換した初期円筒中心を定めている。数25において各行列成分Mbは数26である。数26においてfi及びfiの平均、fi及びfiの角度A,Bによる偏微分、fi及びfiの円筒中心Ox,Oyによる偏微分、と当該偏微分の平均は、数27で表せる。数27において、投影した計測点Qiおよび投影した球の中心O、投影した計測点Qiの角度A,Bは数23で表せる。数23において偏微分した値が数27に反映されているもの、及び偏微分または相互偏微分して0ベクトルになるものは記述していない。数27において、計測点Piを投影するためのx軸方向方向ベクトル、当該ベクトルの角度A,Bによる偏微分は数24で表せる。数24において偏微分または相互偏微分して0ベクトルになるものは記述していない。また、数27において、計測点Piを投影するためのy軸方向方向ベクトル、当該ベクトルの角度A,Bによる偏微分は数22で表せる。数22において偏微分または相互偏微分して0ベクトルになるものは記述していない。数27に示すように、円筒中心Ox、Oyによる相互偏微分、2回上の偏微分はすべて0か2である。そのため、数26は数28、数27は数29に示すように簡単化される。以上、数式により数25の解き方を説明したが、次に、図3に示す円筒中心角度算出部108のフローチャートを用いて投影したときに円となる角度A,B、中心点Ox,Oyを高精度化する過程を説明する。図3のステップ331からステップ340は数25を解く過程の一実施例を示すフローチャートである。なお、角度精度向上部107で精度向上した角度は複数存在するので、円筒中心角度算出部108を角度の個数分繰り返して、当該角度毎に角度と中心点を高精度化する。最初に、角度精度向上部107で求めた角度を角度初期値に設定するとともに、球の中心点を当該角度初期値で座標変換した値のx、y成分を中心点初期値に設定するステップ331。次に、角度初期値と中心点初期値に従い数24で投影面のx軸方向のベクトルt、当該ベクトルtの角度Bによる偏微分値tbと2回偏微分値tbb、数22で投影面のy軸方向のベクトルu、当該ベクトルuの角度Aによる偏微分値uaと2回偏微分値uaa、角度Aで微分してから角度Bで微分する相互偏微分値uab、当該ベクトルuの角度Bによる偏微分値ub、2回偏微分値ubbを求める。次に、数23で計測点Piを投影したベクトルQi、当該ベクトルの角度Aによる偏微分Qi,a、2回偏微分Qi,aa、角度Aで微分後に角度Bで微分する相互偏微分Qi,ab、当該ベクトルの角度Bによる偏微分値Qi,b、2回偏微分値Qi,bb、を求めるステップ332。次に、数29のfiにより計測点と中心点の差ベクトルの内積で二乗距離を作成するステップ333。次に、数29のfi,aにより当該二乗距離の角度Aによる偏微分、fi,aaにより角度Aによる2回偏微分、fi,abにより角度Aで偏微分してから角度Bで偏微分する相互偏微分、fi,aoxにより角度Aで偏微分してから中心Oxで偏微分する相互偏微分、fi,aoyにより角度Aで偏微分してから中心Oyで偏微分する相互偏微分、fi,bにより当該二乗距離の角度Bによる偏微分、fi,bbにより角度Bによる2回偏微分、fi,boxにより角度Bで偏微分してから中心Oxで偏微分する相互偏微分、fi,boyにより角度Bで偏微分してから中心Oyで偏微分する相互偏微分、fi,oxにより中心Oxでの偏微分、fi,oyにより中心Oyでの偏微分を作成するステップ334。次に、数29のfavr,aにより当該二乗距離の角度Aによる偏微分の平均、favr,aaにより角度Aによる2回偏微分の平均、favr,abにより角度Aで偏微分してから角度Bで偏微分する相互偏微分の平均、favr,aoxにより角度Aで偏微分してから中心Oxで偏微分する相互偏微分の平均、favr,aoyにより角度Aで偏微分してから中心Oyで偏微分する相互偏微分の平均、favr,bにより当該二乗距離の角度Bによる偏微分の平均、favr,bbにより角度Bによる2回偏微分の平均、favr,boxにより角度Bで偏微分してから中心Oxで偏微分する相互偏微分の平均、favr,boyにより角度Bで偏微分してから中心Oyで偏微分する相互偏微分の平均、favr,oxにより中心Oxでの偏微分の平均、favr,oyにより中心Oyでの偏微分の平均を作成するステップ335。次に、数28により二乗距離と二乗距離の平均の差の二乗和、当該二乗和の角度Aによる偏微分、角度Aによる2回偏微分、角度Aで微分してから角度Bで微分する相互偏微分、角度Aで微分してから中心Oxで微分する相互偏微分、角度Aで微分してから中心Oyで微分する相互偏微分、角度Bによる偏微分、角度Bによる2回偏微分、角度Bで微分してから中心Oxで微分する相互偏微分、角度Bで微分してから中心Oyで微分する相互偏微分、中心Oxでの偏微分、中心Oyでの偏微分を作成するステップ336。次に、ステップ336で求めた値を連立一次方程式の係数に代入するステップ337。次に、当該連立一次方程式を解き、更新量δA,δBを算出するステップ338。当該連立一次方程式を解くには一般的な連立一次方程式の解法、例えば、掃き出し法やLU分解を使えば良い。次に、更新量δA,δB,δOx,δOyが0とみなせるか判定するステップ339。判定方法は、更新量の成分δA,δB,δOx,δOyの絶対値がそれぞれ0とみなせるかを判定する方法のほかに、更新量の内積δA×δAδB×δBδOx×δOxδOy×δOyが0とみなせるかを判定してもかまわない。本実施例では、更新量の内積が0とみなせるかを判定している。更新量が0とみなせない場合、角度中心点更新をおこなうステップ340。更新は、初期値A,B,Ox,Oyに更新量δA,δB,δOx,δOyを加算することでおこなう。更新した初期値A,B,Ox,Oyでステップ332から繰り返す。更新量が0とみなせる場合、処理を終了する。以上のステップ331からステップ340を繰り返すことで、円筒の中心Ox、Oyと角度A,Bを定めることができる。実際に計測点を与え実験してみたところ、初期角度として角度精度向上部107で算出した角度、初期中心として球中心算出部105で作成した球中心を当該角度で座標変換した中心を与えたとき、5から6回で収束した。尚、初期角度算出部106、角度精度向上部107の効果を確認するため、実際に計測点を与え、分割ピッチを360°/N、Nは1から19まで変化させ、角度精度向上部107をスキップした場合としない場合で円筒半径の違いを確認した。角度精度向上部107をスキップした場合、Nが15より大きくなった場合に円筒半径が一定値となった。これは、24°ピッチに相当する。これに対し、角度精度向上部107をスキップしない場合、Nが7より大きい場合に一定値となった。これは約51.42°ピッチに相当する。よって、角度精度向上部107は角度ピッチを2倍まで大きくできることがわかった。このことから、現在、角度ピッチは余裕をみて20°ピッチとしている。しかしながら、角度ピッチは与える計測点によって変化する可能性があるので記憶装置204にあらかじめ格納しておき、適宜変更可能な構成としている。また、初期角度算出部106の計算量は角度ピッチを2倍にすると約4倍に増加するため、角度精度向上部107が計算量を減少させることがわかった。計算量の減少は計算時間の減少に繋がるため、その効果は大きいことがわかった。図3のフローチャートの説明に戻る。円筒中心角度算出部108の処理が終了すると、円筒データ格納部109の処理を実行する。円筒データ格納部109は、円筒中心角度算出部108で作成した円の中心O、角度を記憶装置102に格納する。ここで、円の中心Oは2次元座標であるため、Z座標を0とし、数3に示す座標変換の逆変換をおこなえば、円筒の中心点が得られる。また、方向ベクトル0,0,1を同様に逆変換すれば、円筒の中心軸が得られる。半径は、数29で作成したfavrの平方根である。円筒データが計測点を含むようにするには、中心点と計測点の差ベクトルと、中心軸の内積を作成して最低値と最大値を作成し、最低値が円筒の底面になるよう中心点を中心軸方向に移動すれば良い。そのとき、円筒の高さは最大値-最小値である。以上、図1~図5を用いて本発明を説明したが、最後に図6を用いて本発明の方式を説明する。図6は、本発明の各処理に応じて作成される形状のイメージ図である。図6では、最初に計測点群が与えられる。次に、計測点群に球をあてはめ、球中心を作成する。次に、球中心を円筒の中心と仮定し、投影したときに円との分散が最小となる角度を少なくとも一つ作成する。次に当該角度を高精度化し、最後に、当該角度および円筒の中心点を高精度化して円筒あてはめを終了する。以上により、従来の問題点である、あてはめ時間増大、別解に収束する問題を解決する。図9は、本発明で説明した球、円筒の3次元モデルを別々、及び重ねて表示した結果であり、最初に計測点から球をあてはめ、仮円筒と円筒をあてはめていった結果を記している。また、最下段は計測点、球、円筒を重ねて表示した結果である。計測点と球、円筒の位置関係が判別できる。図10は、何回目の計測がどの程度の時間なのかを示す図である。各処理は1000回繰り返しているので、左端の目盛の単位はmsecである。円弧溝の計算時間が多いのは、円筒3、円筒4の計測点数が8であるのと比較して、円弧溝の計測点数が150点になっているためである。これにより、十分高速な処理ができていることが理解できる。"