JP2015152535A - Phase analysis method of lattice image using weighting - Google Patents
Phase analysis method of lattice image using weighting Download PDFInfo
- Publication number
- JP2015152535A JP2015152535A JP2014028830A JP2014028830A JP2015152535A JP 2015152535 A JP2015152535 A JP 2015152535A JP 2014028830 A JP2014028830 A JP 2014028830A JP 2014028830 A JP2014028830 A JP 2014028830A JP 2015152535 A JP2015152535 A JP 2015152535A
- Authority
- JP
- Japan
- Prior art keywords
- phase
- weighting
- phase analysis
- analysis
- grating
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Length Measuring Devices By Optical Means (AREA)
Abstract
Description
本発明は、２次元格子画像の位相解析方法に関する。 The present invention relates to a phase analysis method for a two-dimensional lattice image.
従来、２台のカメラを用いて２次元格子模様のゴムシートを貼り付けた物体を撮影し、その撮影された２次元格子シートの画像を位相解析して、物体表面の三次元形状とひずみ分布の変化を高精度に測定する測定システムの開発されている。２次元格子の画像を解析する位相解析手法として、高精度かつ高速な位相解析手法であるサンプリングモアレ法（Ｓａｍｐｌｉｎｇ ｍｏｉｒｅ ｍｅｔｈｏｄ）が用いられてきた。 Conventionally, an object with a two-dimensional lattice pattern rubber sheet attached is photographed using two cameras, the image of the photographed two-dimensional lattice sheet is phase-analyzed, and the three-dimensional shape and strain distribution of the object surface. A measurement system has been developed that measures changes in the air with high accuracy. As a phase analysis method for analyzing an image of a two-dimensional grating, a sampling moire method (Sampling moire method), which is a high-precision and high-speed phase analysis method, has been used.
特許文献１には、物体の表面に存在する格子模様の所定の領域を光学式カメラで撮影し、撮影された画像に対して、等間隔の画素ごとのサンプリングを、起点の画素を変えながら３回以上の複数回実行し、このサンプリング処理によって得られた間引き画像を補間処理することによってモアレ縞画像を生成し、位相シフト法によって得られるモアレ縞の位相分布を求める格子画像の位相解析方法が開示されている。 In Patent Document 1, a predetermined area of a lattice pattern existing on the surface of an object is photographed with an optical camera, and sampling for each pixel at equal intervals is performed on the photographed image while changing the starting pixel. A phase analysis method for a lattice image that is executed a plurality of times and interpolates the thinned image obtained by the sampling process to generate a moire fringe image and obtains the phase distribution of the moire fringes obtained by the phase shift method. It is disclosed.
特開平２００９−２６４８５２号公報 JP 2009-264852 A
しかし、サンプリングモアレ法は、等間隔な格子を対象とした位相解析手法であるため、間引き間隔と格子ピッチが大きく異なる場合、位相解析精度が低下するという問題がある。つまり、カメラの配置や格子を貼付する表面が曲面である場合が要因で、撮影された２次元格子が変形して格子のピッチが変化してしまうため、サンプリングモアレ法では計測精度が大幅に低下してしまう問題が発生する。 However, since the sampling moire method is a phase analysis method for equally spaced gratings, there is a problem that the phase analysis accuracy is lowered when the thinning interval and the grating pitch are greatly different. In other words, due to the fact that the surface of the camera and the surface to which the grid is attached are curved surfaces, the captured two-dimensional grid is deformed and the pitch of the grid changes, so the sampling moire method significantly reduces the measurement accuracy. The problem that will occur.
そこで、本発明の目的は、新たな２次元格子の位相解析手法として、重み付け位相解析法（Ｗｅｉｇｈｔｉｎｇ ｐｈａｓｅ ａｎａｌｙｓｉｓ ｍｅｔｈｏｄ）を提供することである。 Therefore, an object of the present invention is to provide a weighting phase analysis method as a new two-dimensional grating phase analysis method.
本発明は、求めたい画素の周辺画素に重み付けを行い、位相シフト法の原理を用いて位相解析することで、高精度に格子の位相分布を求めることができる方法である。この位相解析方法の一番の特徴は、容易に重み付けを変化させることができることであり、重み付けを適切に行うことで高精度に２次元格子画像の位相分布を求めることができる。
よって、格子ピッチが大きく変化するような大変形している格子に対しても、重み付け位相解析法を用いた位相解析は、サンプリングモアレ法に比べて位相解析誤差が小さくなるので、精度よく位相解析が可能となる。
そして、本願の請求項１に係る発明は、注目点を起点にして、注目点の前後の微小領域の画素値から求めた位相値を微小領域の位置ごとに複数個求め、微小領域の位置に応じた重みをかけて平均化することによって注目点の位相値を求める方法である。
請求項２に係る発明は、２次元格子画像において、縦方向と横方向にそれぞれ平均化して得られた２つの画像に対して請求項１に係る位相値を求める方法をそれぞれ適用し、２次元格子の位相分布を求める方法である。
The present invention is a method that can obtain a phase distribution of a grating with high accuracy by weighting peripheral pixels of a pixel to be obtained and performing phase analysis using the principle of the phase shift method. The most important feature of this phase analysis method is that the weighting can be easily changed, and the phase distribution of the two-dimensional lattice image can be obtained with high accuracy by appropriately performing the weighting.
Therefore, phase analysis using the weighted phase analysis method has a smaller phase analysis error than the sampling moire method even for large deformed lattices whose lattice pitch changes greatly. Is possible.
In the invention according to claim 1 of the present application, a plurality of phase values obtained from the pixel values of the minute area before and after the attention point are obtained for each position of the minute area, and the position of the minute area is determined. In this method, the phase value of the attention point is obtained by averaging with a corresponding weight.
The invention according to claim 2 applies the method for obtaining the phase value according to claim 1 to two images obtained by averaging in the vertical direction and the horizontal direction in a two-dimensional lattice image, respectively. This is a method for obtaining the phase distribution of the grating.
本発明により、格子画像の位相解析に要する時間が短縮される、格子画像の位相解析の解析精度が向上し、格子画像を用いた三次元計測、変位計測、ひずみ計測の計測精度が向上する。 According to the present invention, the time required for the phase analysis of the lattice image is reduced, the analysis accuracy of the phase analysis of the lattice image is improved, and the measurement accuracy of three-dimensional measurement, displacement measurement, and strain measurement using the lattice image is improved.
以下、本発明の実施形態を図面とともに説明する。
１．位相シフト法による格子の位相解析
位相シフト法とは、格子１周期の輝度情報から、その格子の位相を求めることができる手法である。格子画像が撮影された場合、位相シフト法を用いて位相解析することで、その格子の位相分布を求めることができる。
Hereinafter, embodiments of the present invention will be described with reference to the drawings.
1. Phase analysis of grating by phase shift method The phase shift method is a technique that can determine the phase of the grating from the luminance information of one period of the grating. When a lattice image is taken, the phase distribution of the lattice can be obtained by performing phase analysis using the phase shift method.
ここで、ｘ軸方向に周期的に輝度が変化する格子の位相解析について説明する。図１にｘ軸方向に周期的に輝度が変化する格子が撮影された場合の、撮影画像の座標（ｘ，ｙ）付近における格子の輝度と位相の関係について示す。ｘ軸方向に周期的に輝度が変化する格子が撮影された場合、ｘ軸方向にその格子の位相が変化していると考えることができる。よって、撮影された格子画像の輝度を位相が変化している方向に１画素ごとにサンプリングした場合、格子１周期の輝度情報を得ることができ、位相シフト法を適用することが可能となる。 Here, a phase analysis of a grating whose luminance periodically changes in the x-axis direction will be described. FIG. 1 shows the relationship between the luminance of the grating and the phase in the vicinity of the coordinates (x, y) of the photographed image when a grating whose luminance periodically changes in the x-axis direction is photographed. When a grating whose luminance periodically changes in the x-axis direction is photographed, it can be considered that the phase of the grating changes in the x-axis direction. Therefore, when the luminance of the captured lattice image is sampled for each pixel in the direction in which the phase changes, luminance information for one period of the lattice can be obtained, and the phase shift method can be applied.
ここで、撮影された格子１周期の画素数がＮ画素であった場合について考える。位相シフト法を適用するためには、求めたい格子のピッチＮをあらかじめ決定する必要がある。数１式に、離散フーリエ変換を用いた位相シフト法による座標（ｘ，ｙ）のｘ方向の格子の位相算出式を示す。 Here, consider a case where the number of pixels in one captured period of the grid is N pixels. In order to apply the phase shift method, it is necessary to determine in advance the pitch N of the grating to be obtained. Equation (1) shows the phase calculation formula of the lattice in the x direction of the coordinates (x, y) by the phase shift method using discrete Fourier transform.
数１式は、位相シフト法において、座標（ｘ，ｙ）における位相、つまり初期位相を求める位相算出式である。Ｉ（ｘ，ｙ）は座標（ｘ，ｙ）における輝度値であり、Ｎは求めたいｘ方向格子の１周期の画素数である。
ここで、ｘ軸方向に１画素ごと変化すると位相が１回位相シフトすると考えると、ｓはｘ軸方向の座標の変化と対応して変化しているため、ｓは座標（ｘ，ｙ）からのｘ方向への位相シフト回数と考えることができる。つまり、座標（ｘ，ｙ）からｘ軸方向に位相シフトしていると考えたとき、座標（ｘ，ｙ）の位相シフト回数ｎは０回となり、座標（ｘ＋ｓ，ｙ）における位相シフト回数ｎはｓ回となる。
Formula 1 is a phase calculation formula for obtaining a phase at coordinates (x, y), that is, an initial phase, in the phase shift method. I (x, y) is a luminance value at the coordinates (x, y), and N is the number of pixels in one period of the x-direction lattice to be obtained.
Here, assuming that the phase shifts once when the pixel changes in the x-axis direction, s changes corresponding to the change in the coordinate in the x-axis direction, so s is determined from the coordinate (x, y). The number of phase shifts in the x direction can be considered. That is, when the phase shift from the coordinate (x, y) to the x-axis direction is considered, the number of phase shifts n of the coordinates (x, y) is 0, and the number of phase shifts n at the coordinates (x + s, y). Becomes s times.
次に、座標（ｘ，ｙ）における位相シフト回数ｎが０≦ｎ≦Ｎ−１である場合の位相シフト法による位相解析について考える。そのときの座標（ｘ，ｙ）の位相は、数２式で求めることができる。また、数２式に０≦ｎ≦３，Ｎ＝４の場合の数１式による位相解析について示す。数２式における＋２・ｓ／Ｎは、ｎ回の位相シフト回数における位相の変化量を表している。 Next, phase analysis by the phase shift method when the number n of phase shifts at coordinates (x, y) is 0 ≦ n ≦ N−1 will be considered. The phase of the coordinates (x, y) at that time can be obtained by Equation 2. In addition, the phase analysis according to Formula 1 in the case of 0 ≦ n ≦ 3 and N = 4 is shown in Formula 2. In Expression 2, + 2 · s / N represents the amount of change in phase at the number of n phase shifts.
また、数２式の左辺は位相であり、右辺の＋２πｎ／Ｎも位相であるため、数２式においては、右辺の＋２πｎ／Ｎは、数３式に示す複素平面での回転の座標変換と考えることができる。 Also, since the left side of Equation 2 is the phase and + 2πn / N on the right side is also the phase, in Equation 2, + 2πn / N on the right side is the coordinate conversion of the rotation on the complex plane shown in Equation 3. Can think.
ゆえに、数３式によって数２式は数４式，数５式と変形することができる。 Therefore, Formula 2 can be transformed into Formula 4 and Formula 5 by Formula 3.
ここで、ｋ＝ｓ−ｎと定義する。ｋを用いることで、数５式は、さらに数６式に変形することができる。また、ｋは座標（ｘ，ｙ）の位相シフト回数を０と考えた場合のそれぞれの座標の位相シフト回数と考えることができる。 Here, it is defined as k = s−n. By using k, Formula 5 can be further transformed into Formula 6. Further, k can be considered as the number of phase shifts of each coordinate when the number of phase shifts of the coordinate (x, y) is considered to be zero.
図３に０≦ｎ≦３，Ｎ＝４の場合における数６式を用いた位相解析について示す。ｎが変化した場合、位相シフト法の計算式が変化する。つまり、数６式は、ｎの数の分だけ、位相解析パターンが考えられることを示している。 FIG. 3 shows the phase analysis using Equation 6 in the case of 0 ≦ n ≦ 3 and N = 4. When n changes, the calculation formula of the phase shift method changes. That is, Equation 6 shows that the phase analysis pattern can be considered by the number of n.
２．重み付け位相解析法による格子の位相解析
重み付け位相解析法とは、注目画素の周辺画素に対して重み付けを行い、位相シフト法の原理を利用した計算式から、高精度に位相解析することができる手法である。
2. Phase analysis of grating by weighted phase analysis method Weighted phase analysis method is a method that weights the surrounding pixels of the pixel of interest and performs phase analysis with high accuracy from a formula using the principle of the phase shift method. It is.
まず、重み付け位相解析法の原理について述べる。重み付け位相解析法の原理は、数６式で得られる複数の位相算出式に対して重み付けを行い、位相解析することである。位相シフト法による位相解析は、誤差やばらつきの影響を受けやすく、位相解析精度はあまり良くないが、それぞれの位相解析結果に対して適切な重み付けを行うことで、誤差やばらつきの影響を低減させることができる。 First, the principle of the weighted phase analysis method will be described. The principle of the weighted phase analysis method is to perform phase analysis by weighting a plurality of phase calculation formulas obtained by Formula 6. Phase analysis using the phase shift method is easily affected by errors and variations, and the phase analysis accuracy is not very good. However, by appropriately weighting each phase analysis result, the effects of errors and variations are reduced. be able to.
次に、数６式の位相シフト法を用いた重み付け位相解析法によるｘ方向格子の位相解析について説明する。数７式に、数６式の位相シフト法に対して重み付けを行う重み付け関数Ｗ．（ｎ，ｋ）について示す。ｗ（ｎ）は０以上の任意の実数の重みである。Ｗ．（ｎ，ｋ）を用いて数６式の位相シフト法の位相解析結果に対して重み付けを行うことで、重み付け位相解析法の位相算出式が得られる。数８式に、ｎが任意の整数の範囲内の値である場合の重み付け位相解析法による位相算出式について示す。また、数８式のｎとｋは、数８式の位相シフト法のｎとｋとそれぞれ対応している。 Next, the phase analysis of the x-direction grating by the weighted phase analysis method using the phase shift method of Formula 6 will be described. Weighting function W. for weighting Equation 7 to the phase shift method of Equation 6. (N, k) will be described. w (n) is an arbitrary real number weight of 0 or more. W. By weighting the phase analysis result of the phase shift method of Formula 6 using (n, k), the phase calculation formula of the weighted phase analysis method is obtained. Equation 8 shows a phase calculation formula based on the weighted phase analysis method when n is a value within an arbitrary integer range. Further, n and k in Expression 8 correspond to n and k in the phase shift method of Expression 8, respectively.
数８式において、Ｗ．（ｎ，ｋ）は座標（ｘ，ｙ）の輝度Ｉ（ｘ，ｙ）にかかる重みである。数８式における、それぞれの位相シフト法が座標に掛ける重み付けの総和をＷx（ｋ）とすると、Ｗx（ｋ）は数９式で表される。また、ｋは注目座標（ｘ，ｙ）からの距離となる。 In Formula 8, W. (N, k) is a weight applied to the luminance I (x, y) at the coordinates (x, y). In Equation 8 wherein the total sum of the weighting respective phase shift method subjected to coordinate and W x (k), W x (k) is expressed by equation (9). K is a distance from the coordinate of interest (x, y).
つまり、数９式のＷx（ｋ）を用いると、数８式は数１０式で表すことができる。ただし、数１０式のＷx（ｋ）は数９式で表される重み付け関数でなければならない。また、ｋは任意の整数の範囲内の値である。 In other words, when W x (k) in Expression 9 is used, Expression 8 can be expressed by Expression 10. However, W x (k) in Formula 10 must be a weighting function expressed by Formula 9. K is a value within an arbitrary integer range.
２次元格子の位相解析の場合、ｙ方向の画素に対しても重み付けを行うことで２次元格子のｘ方向の格子の位相解析が可能となる。ｙ方向の画素に対しても重み付けを行ったｘ方向の格子の位相解析は、ｙ方向への重み付け関数をＷy（ｌ）と定義すると、数８式より数１１式で表される。また、ｌは任意の整数の範囲内の値である。また、数１１式は数９式より数１２式として表される。 In the case of the phase analysis of the two-dimensional grating, the phase analysis of the grating in the x direction of the two-dimensional grating can be performed by weighting the pixels in the y direction. The phase analysis of the lattice in the x direction in which the pixels in the y direction are also weighted is expressed by the following equation (11) from the equation (8) when the weighting function in the y direction is defined as W y (l). Further, l is a value within an arbitrary integer range. In addition, Expression 11 is expressed as Expression 12 from Expression 9.
ｘ方向への重み付け関数Ｗx（ｋ）とｙ方向への重み付け関数Ｗy（ｌ）は、２次元方向への重み付け関数Ｗxy（ｋ，ｌ）として表すことができる。Ｗxy（ｋ，ｌ）を数１３式として定義した場合、数１２式は数１４式として表される。 The weighting function W x (k) in the x direction and the weighting function W y (l) in the y direction can be expressed as a weighting function W xy (k, l) in the two-dimensional direction. When W xy (k, l) is defined as Equation 13, Equation 12 is expressed as Equation 14.
また、ｙ方向格子の位相解析は前述の位相算出式のｘとｙを入れ替えると位相解析できる。よって、ｙ方向格子の位相解析における２次元方向への重み付け関数Ｗyx（ｋ，ｌ）を数１５式として定義すると、数１６式を用いることで位相解析することができる。 The phase analysis of the y-direction grating can be performed by exchanging x and y in the above-described phase calculation formula. Therefore, when the weighting function W yx (k, l) in the two-dimensional direction in the phase analysis of the y-direction grating is defined as Equation 15, the phase analysis can be performed using Equation 16.
３．重み付け位相解析法における重み付け関数について
数８式のｎの範囲を０≦ｎ≦Ｎ−１と定義した場合における位相解析について説明する。ｎの範囲が０≦ｎ≦Ｎ−１の場合は、数８式の位相算出パターンはＮ通り考えることができ、その場合の数８式は数１７式で表される。
3. Regarding Weighting Function in Weighted Phase Analysis Method Phase analysis when the range of n in Equation 8 is defined as 0 ≦ n ≦ N−1 will be described. When the range of n is 0 ≦ n ≦ N−1, the phase calculation pattern of Formula 8 can be considered in N ways, and Formula 8 in that case is expressed by Formula 17.
まず、複数の位相算出式に均等の重み付けを行う場合の重み付け位相解析法による格子の位相解析について述べる。均等に重み付けを行う場合の重み付け関数Ｗ．（ｎ，ｋ）は、数８式より数１８式で表すことができる。 First, the phase analysis of the grating by the weighted phase analysis method when equal weighting is applied to a plurality of phase calculation formulas will be described. Weighting function when weighting equally W. (N, k) can be expressed by Equation 18 from Equation 8.
また、数１８式の重み付け関数を用いた場合の数１７式は数１９式で表される。 Further, Formula 17 when the weighting function of Formula 18 is used is expressed by Formula 19.
数１９式を用いた場合の、座標（ｘ，ｙ）のｘ方向格子の位相の位相解析について説明する。位相シフト法に掛ける重み付け関数を数１８式とした場合の、それぞれの位相シフト法が座標に掛ける重み付けについて図に説明したものを図４に示す。また、位相シフト法に均等に重み付けを行った場合の重み付け位相解析法の画素にかける重み付け関数Ｗx（ｋ）は、１画素における最大の重みの比率を１とした場合、数２０式となる。また、数２０式の重み付け関数Ｗx（ｋ）は、三角波状の重み付け関数となる。図５にＮ＝１０の場合の重み付け関数について示す。 The phase analysis of the phase of the x-direction grating at the coordinates (x, y) when Equation 19 is used will be described. FIG. 4 illustrates the weighting that each phase shift method applies to the coordinates when the weighting function to be applied to the phase shift method is expressed by equation (18). Further, the weighting function W x (k) applied to the pixels of the weighted phase analysis method when weighting is equally performed in the phase shift method is expressed by Equation 20 when the maximum weight ratio in one pixel is 1. . In addition, the weighting function W x (k) of Equation 20 is a triangular wave weighting function. FIG. 5 shows the weighting function when N = 10.
数２１式に数２０式のＷx（ｋ）を用いた場合の位相算出式について示す。数２１式は、線形補間を利用したサンプリングモアレ法による格子の位相算出式と同じ計算式となる。 A phase calculation formula when W x (k) of Formula 20 is used as Formula 21 will be described. Equation (21) is the same as the equation for calculating the phase of the grating by the sampling moire method using linear interpolation.
次に、位相シフト法に掛ける重み付けを複数の位相シフト法の中心画素の重みが大きくなるように三角波状の重み付け関数、余弦波状の重み付け関数を用いた場合を考える。三角波状の重み付け関数を数２２式に、余弦波状の重み付け関数を数２３式に示す。 Next, consider a case where a weighting function having a triangular wave shape or a weighting function having a cosine wave shape is used so that the weighting applied to the phase shift method increases the weight of the center pixel of the plurality of phase shift methods. The triangular wave-like weighting function is shown in Equation 22, and the cosine wave-like weighting function is shown in Equation 23.
数２２式の三角波状の重み付けを用いた場合の重み付け位相解析法の説明図を図６に示し、最大の重みの比率を１とした場合の数２２式の重み付け関数Ｗx（ｋ），Ｎ＝１０を図７に示す。また、数２３式の余弦波状の重み付けを用いた場合の重み付け位相解析法の説明図を図８に示し、最大の重みの比率を１とした場合の数２３式の重み付け関数Ｗx（ｋ），Ｎ＝１０を図９に示す。
ｎの範囲が０≦ｎ≦Ｎ−１の場合における、図７と図９の画素への重み付け関数は、余弦波状に似た分布になることから、重み付け位相解析法では一般的に、数２４式に示す余弦波の重み付け関数を用いて位相解析を行う。また、数２５式に数２４式の重み付け関数を用いた重み付け位相解析法による位相算出式を示す。数２５式はハニング窓を用いた窓関数付フーリエ変換の周波数２の位相算出式と同じ計算式となる。
FIG. 6 shows an explanatory diagram of the weighting phase analysis method using the triangular wave weighting of Formula 22, and the weighting function W x (k), N of Formula 22 when the maximum weight ratio is 1. = 10 is shown in FIG. Also, FIG. 8 shows an explanatory diagram of the weighting phase analysis method in the case where the cosine wave-like weighting of Expression 23 is used, and the weighting function W x (k) of Expression 23 when the maximum weight ratio is 1. , N = 10 is shown in FIG.
In the case where the range of n is 0 ≦ n ≦ N−1, the weighting function for the pixels in FIGS. 7 and 9 has a distribution similar to a cosine wave. Phase analysis is performed using the cosine wave weighting function shown in the equation. In addition, Equation (25) shows a phase calculation formula by a weighted phase analysis method using the weighting function of Equation (24). Equation 25 is the same calculation formula as the phase calculation formula of the frequency 2 of the Fourier transform with window function using the Hanning window.
また、ｎの範囲が０≦ｎ≦Ｎ−１の範囲における重み付け関数の例について述べる。放物線を利用した重み付け関数の例を数２６式と数２７式に、三角関数を利用した重み付け関数の例を数２８式，数２９式に示す。数２６式〜数２９式は、任意のパラメータｐ（＞０）を変化させることで用意に重みを変化させることができ、これらの式で得られる重みは組み合わせて使用することが可能である。また、ｐ＝１のとき、数２６式と数２７式は三角波状の重み付け関数に、数２８式と数２９式は余弦波状の重み付け関数となる。 An example of a weighting function in the range where n is 0 ≦ n ≦ N−1 will be described. Examples of weighting functions using parabolas are shown in Equations 26 and 27, and examples of weighting functions using trigonometric functions are shown in Equations 28 and 29. In Expressions 26 to 29, weights can be changed easily by changing an arbitrary parameter p (> 0), and the weights obtained by these expressions can be used in combination. Further, when p = 1, Equations 26 and 27 are triangular wave weighting functions, and Equations 28 and 29 are cosine wave weighting functions.
４．シミュレーションによる位相解析精度の数値解析
数２５式を用いた重み付け位相解析法の位相解析精度のシミュレーションをおこなった。位相解析精度の精度評価の方法について説明する。まず、図１０（ａ）に示すコンピュータで作成した２００×２００［ｐｉｘｅｌ］、格子間ピッチ１０［ｐｉｘｅｌ］の理想的な２次元矩形波格子画像に対して、縦方向に２０［ｐｉｘｅｌ］の領域の平滑化処理を行い、図１０（ｂ）に示す、ｘ方向格子画像の作成を行った。その後、ｘ方向格子画像に対して、重み付け位相解析法（数２５式）と線形補間を用いたサンプリングモアレ法（数２１式）それぞれで位相解析を行った。位相解析は、解析パラメータである格子１ピッチの間隔Ｎ（サンプリングモアレ法の場合は間引き間隔）を６〜１９［ｐｉｘｅｌ］に変化させて、ｘ方向格子の位相解析を重み付け位相解析法とサンプリングモアレ法を用いて行った。
4). Numerical analysis of phase analysis accuracy by simulation A simulation of the phase analysis accuracy of the weighted phase analysis method using Equation 25 was performed. A method for evaluating the accuracy of the phase analysis accuracy will be described. First, an area of 20 [pixel] in the vertical direction with respect to an ideal two-dimensional rectangular wave lattice image having a size of 200 × 200 [pixel] and an inter-lattice pitch of 10 [pixel] created by the computer shown in FIG. The smoothing process was performed, and an x-direction lattice image shown in FIG. 10B was created. Thereafter, phase analysis was performed on the x-direction lattice image by a weighted phase analysis method (Equation 25) and a sampling moire method (Equation 21) using linear interpolation. In the phase analysis, the pitch N pitch of the grating, which is the analysis parameter (in the case of the sampling moire method) is changed from 6 to 19 [pixel], and the phase analysis of the x-direction grating is changed to the weighted phase analysis method and the sampling moire. The method was used.
その後、図１０（ｃ）に示す、ｘ方向格子の位相分布に対して位相接続を行い、Ｆｉｇ．４．４に示す、位相接続されたｘ方向格子の位相を１周期分ラインで抜き出し、理論値との誤差を計算し、位相解析精度の評価を行った。
図１０（ｄ）の座標（９５，１００）〜（１０４，１００）のサンプリングモアレ法と重み付け位相解析法の位相解析結果の位相解析誤差と解析格子ピッチの関係を表にまとめたものを図１１に示し、位相解析誤差の標準偏差と解析格子ピッチの関係をグラフにしたものを図１２に示す。また、解析格子ピッチＮが８〜１２［ｐｉｘｅｌ］であった場合のサンプリングモアレ法の位相解析誤差のグラフを図１３に示し、重み付け位相解析法の位相解析誤差のグラフを図１４に示す。
Thereafter, phase connection is made to the phase distribution of the x-direction grating shown in FIG. The phase of the phase-connected x-direction grating shown in 4.4 was extracted by a line for one period, the error from the theoretical value was calculated, and the phase analysis accuracy was evaluated.
FIG. 11 is a table summarizing the relationship between the phase analysis error and the analysis lattice pitch of the phase analysis results of the sampling moire method and the weighted phase analysis method of the coordinates (95, 100) to (104, 100) in FIG. FIG. 12 is a graph showing the relationship between the standard deviation of the phase analysis error and the analysis lattice pitch. FIG. 13 shows a phase analysis error graph of the sampling moire method when the analysis lattice pitch N is 8 to 12 [pixel], and FIG. 14 shows a phase analysis error graph of the weighted phase analysis method.
５．変位量計測実験による計測精度評価
２次元格子シートを貼付した移動ステージを用いた変位量計測実験をおこなった。実験環境の風景と装置の位置関係を図１５に示す。移動ステージに変位を与える際に、Ｍｉｔｕｔｏｙｏ社製のマイクロメーターを、カメラは、画像サイズ６４０×４８０［ｐｉｘｅｌ］のＤｒａｇｏｎｆｌｙＥｘｐｒｅｓｓを使用した。また、２次元格子シートには４．０［ｍｍ］ピッチの２次元矩形波格子を使用した。
5. Measurement accuracy evaluation by displacement measurement experiment A displacement measurement experiment using a moving stage with a two-dimensional lattice sheet was performed. FIG. 15 shows the positional relationship between the scenery of the experimental environment and the apparatus. When giving displacement to the moving stage, a Mitutoyo micrometer and a Dragonfly Express with an image size of 640 × 480 [pixel] were used as the camera. A two-dimensional rectangular wave grating having a pitch of 4.0 [mm] was used as the two-dimensional grating sheet.
変位量計測実験の実験方法について説明する。図１５（ｂ）に示すように、カメラの光軸を基準面の表面に対してｚ軸方向に３０度傾けて配置し、撮影画像の格子のピッチが変化するようにした。その後、移動ステージを用いて、ｘ方向に０．１［ｍｍ］ずつ、最大０．５［ｍｍ］の変位を与え、変位を与えた位置でそれぞれ撮影を行い、その後、重み付け位相解析法とサンプリングモアレ法それぞれで位相解析を行い、得られた格子のｘ方向の位相分布から変位量の計測結果を算出し、その計測結果の精度評価を行った。 The experiment method of the displacement measurement experiment will be described. As shown in FIG. 15B, the optical axis of the camera is arranged to be inclined by 30 degrees in the z-axis direction with respect to the surface of the reference plane so that the pitch of the grid of the captured image changes. After that, using a moving stage, a displacement of 0.1 [mm] in the x direction and a maximum of 0.5 [mm] is given, and images are taken at the positions where the displacement is given. The phase analysis was performed by each of the moire methods, the measurement result of the displacement amount was calculated from the phase distribution in the x direction of the obtained grating, and the accuracy of the measurement result was evaluated.
また、画像に対する縦方向の平滑化画素数を７０［ｐｉｘｅｌ］、解析格子ピッチＮを２０［ｐｉｘｅｌ］に設定して位相解析を行った。撮影におけるランダムノイズは今回の精度評価には不必要であるため、１００回撮影を行って平均化した画像を撮影画像として用いることで、ランダムノイズの低減を行った。２次元格子を撮影した画像を図１６（ａ）に、図１６（ａ）の撮影画像のｘ方向格子の位相解析を行った画像を図１６（ｂ）にそれぞれ示す。 Further, the phase analysis was performed with the number of smoothed pixels in the vertical direction of the image set to 70 [pixel] and the analysis grid pitch N set to 20 [pixel]. Since random noise in shooting is unnecessary for this accuracy evaluation, random noise was reduced by using an image averaged after shooting 100 times as a shot image. FIG. 16A shows an image obtained by photographing a two-dimensional lattice, and FIG. 16B shows an image obtained by performing phase analysis of the x-direction lattice of the photographed image shown in FIG.
変位量の計測は、変位前と変位後の位相と２次元格子のピッチｐ［ｍｍ］を用いて計測することができる。変位量をｄ［ｍｍ］、格子ピッチをｐ［ｍｍ］、変位前後の位相の変化量を・・、・・・≦・・・≦・とすると、変位量は数３０式で求められる。 The displacement amount can be measured using the phase before and after displacement and the pitch p [mm] of the two-dimensional grating. If the displacement amount is d [mm], the grating pitch is p [mm], and the phase change amount before and after the displacement is...,.
サンプリングモアレ法と重み付け位相解析法、それぞれの手法で解析した変位量計測結果の座標（１００，２４０）〜（６００，２４０）の１ラインの計測結果を抜き出して、計測精度の評価を行った。
計測結果をグラフにしたもので、重み付け位相解析法を用いた場合を図１７に、サンプリングモアレ法を用いた場合を図１８に示す。また、図１９に、移動ステージに０．５［ｍｍ］の変位を与えた場合の変位量計測結果に対して、計測結果のｘ方向２１［ｐｉｘｅｌ］の標準偏差の計算を行い、グラフにしたものを示す。
また、格子ピッチは位相分布から数３１式を用いて計算することができる。撮影されたｘ方向格子のピッチを位相分布から計算した結果をグラフにしたものを図２０に示す。
A measurement result of one line of the coordinates (100, 240) to (600, 240) of the displacement measurement result analyzed by the sampling moire method and the weighted phase analysis method was extracted, and the measurement accuracy was evaluated.
FIG. 17 is a graph showing the measurement results, and FIG. 17 shows a case where the weighted phase analysis method is used, and FIG. 18 shows a case where the sampling moire method is used. Further, in FIG. 19, the standard deviation in the x direction 21 [pixel] of the measurement result is calculated with respect to the displacement amount measurement result when the displacement of 0.5 [mm] is given to the moving stage, and the result is plotted in a graph. Show things.
Further, the grating pitch can be calculated from the phase distribution using Equation 31. FIG. 20 is a graph showing the result of calculating the pitch of the photographed x-direction grating from the phase distribution.
６．まとめ
重み付け位相解析法は従来手法であるサンプリングモアレ法に比べ、解析格子ピッチと実際のピッチが大きく異なる場合でも、位相解析誤差が少ない手法であることが、シミュレーションによる解析結果と変位量計測実験の計測結果から確認できた。
6). Conclusion Compared to the conventional sampling moire method, the weighted phase analysis method is a method with less phase analysis error even when the analysis grid pitch and the actual pitch are significantly different. It was confirmed from the measurement results.
シミュレーションによる解析結果から、解析格子ピッチが実際の格子ピッチより大きい場合と小さい場合では、大きい場合のほうが、計測精度の低下が少ないことが確認できた。つまり、格子ピッチが急激に変化する格子に対して位相解析を行う場合、解析格子ピッチを実際のピッチより大きく設定するほうがよいと考えられる。 From the analysis results by simulation, it has been confirmed that the measurement accuracy is less degraded when the analysis lattice pitch is larger or smaller than the actual lattice pitch. That is, when performing phase analysis on a grating whose grating pitch changes rapidly, it is considered better to set the analysis grating pitch to be larger than the actual pitch.
変位量計測実験の結果から、シミュレーションの解析結果の傾向が実計測においても同様の傾向が確認できた。また、移動ステージに与える変位量が大きくなるにつれ、変位計測の計測誤差が大きくなることが確認できた。これは、位相解析によって得られる位相に解析格子ピッチと実際の格子ピッチの違いから発生するランダムではない系統的な位相解析誤差が原因だと考えられる。重み付け位相解析法を用いて解析を行うことで、サンプリングモアレ法に比べ、計測結果のばらつきが低減されていることから、実計測における重み付け位相解析法の有用性が確認できた。
上述したように、本願発明により、格子画像の位相解析に要する時間が短縮される、格子画像の位相解析の解析精度が向上する、ＦＰＧＡで解析する場合に、容量が小さくて済むために、安価なものを利用することができる。格子画像を用いた三次元計測、変位計測、ひずみ計測の計測精度が向上する。
From the results of the displacement measurement experiment, it was confirmed that the trend of the simulation analysis result was the same in the actual measurement. It was also confirmed that the measurement error of the displacement measurement increases as the amount of displacement applied to the moving stage increases. This is considered to be caused by a non-random systematic phase analysis error generated from the difference between the analysis grating pitch and the actual grating pitch in the phase obtained by the phase analysis. The analysis using the weighted phase analysis method has reduced the variation in measurement results compared to the sampling moire method, so the usefulness of the weighted phase analysis method in actual measurement has been confirmed.
As described above, according to the present invention, the time required for the phase analysis of the lattice image can be shortened, the analysis accuracy of the phase analysis of the lattice image can be improved, and the capacity can be reduced when the analysis is performed by the FPGA. You can use anything. Measurement accuracy of three-dimensional measurement, displacement measurement, and strain measurement using lattice images is improved.
本発明は、求めたい画素の周辺画素に重み付けを行い、位相シフト法の原理を用いて位相解析することで、高精度に格子の位相分布を求めることができる方法である。この位相解析方法の一番の特徴は、容易に重み付けを変化させることができることであり、重み付けを適切に行うことで高精度に２次元格子画像の位相分布を求めることができる。
よって、格子ピッチが大きく変化するような大変形している格子に対しても、重み付け位相解析法を用いた位相解析は、サンプリングモアレ法に比べて位相解析誤差が小さくなるので、精度よく位相解析が可能となる。
そして、本願の請求項１に係る発明は、注目点を起点にして、注目点とその周辺の微小領域の画素値から求められる注目点の位相値を、解析領域を微小領域内で複数回変化させることにより、注目点の複素数値を複数個の解析領域でそれぞれ求め、得られた複数個の複素数値に対し、微小領域の位置に応じた重みを前記複数個の複素数値の実部と虚部にかけて積算し、さらに得られた複素数値の偏角を求めることで注目点の位相値を求める方法である。
請求項２に係る発明は、注目点を起点として、注目点の周辺の微小領域の画素値から求めた注目点とその周辺画素の画素値に、その画素の位置に応じた重みをかけて合わせ、それを積算した値を用いて注目点の位相値を求める方法である。
請求項３に係る発明は、２次元格子画像において、縦方向と横方向の位相値を複数個求め、請求項１または請求項２に係る位相値を求める方法を適用して、注目点の縦方向と横方向の位相値を求める、２次元格子の位相分布を求める方法である。
The present invention is a method that can obtain a phase distribution of a grating with high accuracy by weighting peripheral pixels of a pixel to be obtained and performing phase analysis using the principle of the phase shift method. The most important feature of this phase analysis method is that the weighting can be easily changed, and the phase distribution of the two-dimensional lattice image can be obtained with high accuracy by appropriately performing the weighting.
Therefore, phase analysis using the weighted phase analysis method has a smaller phase analysis error than the sampling moire method even for large deformed lattices whose lattice pitch changes greatly. Is possible.
In the invention according to claim 1 of the present application, the phase value of the target point obtained from the pixel value of the target point and the surrounding micro region is changed a plurality of times within the micro region from the target point. by, respectively determined complex values of the target point by a plurality of analysis regions, with respect to a plurality of complex values obtained, the real and imaginary of the plurality of complex-valued weights corresponding to the position of the minute region In other words, the phase value of the target point is obtained by calculating the deviation angle of the obtained complex value .
The invention according to claim 2 is based on the attention point as a starting point, and the attention point obtained from the pixel values of the minute region around the attention point and the pixel values of the surrounding pixels are multiplied by a weight corresponding to the position of the pixel. In this method, the phase value of the point of interest is obtained using a value obtained by integrating the values .
In the invention according to claim 3, in the two-dimensional lattice image, a plurality of phase values in the vertical direction and the horizontal direction are obtained, and the method for obtaining the phase value according to claim 1 or claim 2 is applied to This is a method of obtaining a phase distribution of a two-dimensional grating for obtaining a direction value in a direction and a transverse direction.
本発明は、求めたい画素の周辺画素に重み付けを行い、位相シフト法の原理を用いて位相解析することで、高精度に格子の位相分布を求めることができる方法である。この位相解析方法の一番の特徴は、容易に重み付けを変化させることができることであり、重み付けを適切に行うことで高精度に２次元格子画像の位相分布を求めることができる。
よって、格子ピッチが大きく変化するような大変形している格子に対しても、重み付け位相解析法を用いた位相解析は、サンプリングモアレ法に比べて位相解析誤差が小さくなるので、精度よく位相解析が可能となる。
本願の請求項１に係る発明は、注目点を起点にして、注目点とその周辺の微小領域に含まれる各画素の画素値から求められる注目点の位相値を、解析領域を微小領域内で複数回変化させることにより、注目点の複素数値を複数個の解析領域でそれぞれ求め、得られた複数個の複素数値に対し、微小領域に含まれる各画素の位置に応じた重みを前記複数個の複素数値の実部と虚部にかけて積算し、さらに得られた複素数値の偏角を求めることで注目点の位相値を求める方法である。
請求項２に係る発明は、注目点を起点として、注目点の周辺の微小領域に含まれる各画素の画素値に、それぞれの画素の位置に応じた重みと、それぞれの画素の位置に応じた複素数値を掛け合わせ、それらを積算して得られた複素数値の偏角として注目点の位相値を求める方法である。
請求項３に係る発明は、２次元格子画像において、縦方向と横方向の位相値を複数個求め、請求項１または請求項２に係る位相値を求める方法を適用して、注目点の縦方向と横方向の位相値を求める、２次元格子の位相分布を求める方法である。
The present invention is a method that can obtain a phase distribution of a grating with high accuracy by weighting peripheral pixels of a pixel to be obtained and performing phase analysis using the principle of the phase shift method. The most important feature of this phase analysis method is that the weighting can be easily changed, and the phase distribution of the two-dimensional lattice image can be obtained with high accuracy by appropriately performing the weighting.
Therefore, phase analysis using the weighted phase analysis method has a smaller phase analysis error than the sampling moire method even for large deformed lattices whose lattice pitch changes greatly. Is possible.
The invention according to claim 1 of the present gun is to the target point as the starting point, the phase value of the point of interest obtained from the pixel value of each pixel included target point and the minute region around the micro area analysis region The complex value of the attention point is obtained in each of the plurality of analysis regions by changing the number of times in the plurality of analysis values, and the plurality of complex values obtained are weighted according to the position of each pixel included in the minute region. This is a method of obtaining the phase value of the attention point by accumulating over the real part and the imaginary part of each complex value and further obtaining the declination of the obtained complex value.
According to the second aspect of the present invention, the pixel value of each pixel included in the minute area around the target point is set to the weight corresponding to the position of each pixel and the position of each pixel. This is a method of obtaining a phase value of an attention point as a declination of a complex value obtained by multiplying complex values and integrating them .
In the invention according to claim 3, in the two-dimensional lattice image, a plurality of phase values in the vertical direction and the horizontal direction are obtained, and the method for obtaining the phase value according to claim 1 or claim 2 is applied to This is a method of obtaining a phase distribution of a two-dimensional grating for obtaining a direction value in a direction and a transverse direction.
Claims (2)
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2014028830A JP5795095B2 (en) | 2014-02-18 | 2014-02-18 | Phase analysis method of lattice image using weighting |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2014028830A JP5795095B2 (en) | 2014-02-18 | 2014-02-18 | Phase analysis method of lattice image using weighting |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2015152535A true JP2015152535A (en) | 2015-08-24 |
| JP5795095B2 JP5795095B2 (en) | 2015-10-14 |
Family
ID=53894932
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2014028830A Active JP5795095B2 (en) | 2014-02-18 | 2014-02-18 | Phase analysis method of lattice image using weighting |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP5795095B2 (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2018155115A1 (en) * | 2017-02-23 | 2018-08-30 | 国立研究開発法人産業技術総合研究所 | Deformation measurement method, deformation measurement device, and program therefor |
| JP2018194738A (en) * | 2017-05-19 | 2018-12-06 | キヤノン株式会社 | Position measurement device, lithography device, and article production method |
| JP2021121786A (en) * | 2020-01-31 | 2021-08-26 | 株式会社島津製作所 | Displacement distribution measuring device, displacement distribution measuring method, and control program of displacement distribution measuring device |
| CN113358063A (en) * | 2021-06-04 | 2021-09-07 | 华中科技大学 | Surface structured light three-dimensional measurement method and system based on phase weighted fusion |
-
2014
- 2014-02-18 JP JP2014028830A patent/JP5795095B2/en active Active
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2018155115A1 (en) * | 2017-02-23 | 2018-08-30 | 国立研究開発法人産業技術総合研究所 | Deformation measurement method, deformation measurement device, and program therefor |
| JPWO2018155115A1 (en) * | 2017-02-23 | 2019-11-07 | 国立研究開発法人産業技術総合研究所 | Deformation measuring method, deformation measuring apparatus, and program thereof |
| JP2018194738A (en) * | 2017-05-19 | 2018-12-06 | キヤノン株式会社 | Position measurement device, lithography device, and article production method |
| JP2021121786A (en) * | 2020-01-31 | 2021-08-26 | 株式会社島津製作所 | Displacement distribution measuring device, displacement distribution measuring method, and control program of displacement distribution measuring device |
| JP7338495B2 (en) | 2020-01-31 | 2023-09-05 | 株式会社島津製作所 | Displacement distribution measuring device, displacement distribution measuring method, and control program for displacement distribution measuring device |
| CN113358063A (en) * | 2021-06-04 | 2021-09-07 | 华中科技大学 | Surface structured light three-dimensional measurement method and system based on phase weighted fusion |
| CN113358063B (en) * | 2021-06-04 | 2022-03-18 | 华中科技大学 | Surface structured light three-dimensional measurement method and system based on phase weighted fusion |
Also Published As
| Publication number | Publication date |
|---|---|
| JP5795095B2 (en) | 2015-10-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP6583761B2 (en) | Three-dimensional shape / displacement / strain measuring device, method and program using periodic pattern | |
| JP4873485B2 (en) | Shape measuring method and shape measuring apparatus using a number of reference surfaces | |
| WO2013136620A1 (en) | Phase distribution analysis method and device for fringe image using high-dimensional brightness information, and program therefor | |
| TWI573984B (en) | System and method of adjusting matching image | |
| KR101796129B1 (en) | Method and device for measuring displacement distribution using regular pattern, and program for such method and device | |
| JP5795095B2 (en) | Phase analysis method of lattice image using weighting | |
| CN105579809B (en) | Measuring method, measuring apparatus, and computer-readable recording medium | |
| WO2017175341A1 (en) | Measurement method, measurement device, measurement program, and computer-readable recording medium having measurement program recorded thereon | |
| US11512946B2 (en) | Method and system for automatic focusing for high-resolution structured light 3D imaging | |
| KR101627950B1 (en) | Image processing device, distortion correction map creation device, and semiconductor measurement device | |
| JP6035031B2 (en) | Three-dimensional shape measuring device using multiple grids | |
| JP2011133341A (en) | Displacement measuring device, displacement measuring method, and displacement measuring program | |
| Malesa et al. | Multi-camera DIC system with spatial data stitching procedure for measurements of engineering objects | |
| JP2019139030A (en) | Method and device for projecting information related to measurement result on surface of three-dimensional measurement target object | |
| CN109242905B (en) | Method for determining object displacement and direction based on digital image diffraction | |
| CN104897065A (en) | Measurement system for surface displacement field of shell structure | |
| JP5466325B1 (en) | Method to measure physical quantity of object from image of grid attached to object | |
| KR20190103271A (en) | Cross reflection detection apparatus, cross reflection detection method and program | |
| JP6533914B2 (en) | Computer readable recording medium recording measurement method, measurement device, measurement program and measurement program | |
| JP2018179577A (en) | Position measuring device | |
| JP6977205B2 (en) | Displacement acquisition device, displacement acquisition method, and program | |
| JP6666670B2 (en) | 3D shape measurement method using curved surface as reference surface | |
| Rodríguez | Online self-calibration for mobile vision based on laser imaging and computer algorithms | |
| KR20160119412A (en) | Method and system for estimating 6-dof displacement of large member | |
| CN107588741B (en) | Method and system for measuring camera depth change based on moire fringes |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20150331 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20150622 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20150622 |
|
| A911 | Transfer to examiner for re-examination before appeal (zenchi) |
Free format text: JAPANESE INTERMEDIATE CODE: A911 Effective date: 20150713 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20150804 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20150811 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 5795095 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |