(57)【特許請求の範囲】
【請求項１】
第１時点と当該第１時点から微小時間経過後の第２時点とで同一のカメラを用いて粒子を撮像して得られる第１画像および第２画像から、計測対象となる前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出するステップと、
前記第１粒子像および前記第２粒子像をフーリエ変換するステップと、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出するステップと、
既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データを参照し、算出した評価値に対応する奥行き変位量を特定するステップと、を含む流体計測方法。
【請求項２】
前記評価値は式（８）により算出する請求項１に記載の流体計測方法。
【数１】
0006399517.tif 000008
ただし、ΔＤ（ξ，η，Δｔ）は上記評価値であり、ξは、前記第１粒子像および粒子像のｘ方向における空間周波数であり、ηは、前記第１粒子像および第２粒子像のｙ方向における空間周波数であり、Δｔは前記微小時間であり、ｔは前記第１時点であり、（ｔ＋Δｔ）は前記第２時点であり、Ｇ（ξ，η，ｔ）はフーリエ変換した第１粒子像であり、Ｇ（ξ，η，ｔ＋Δｔ）はフーリエ変換した第２粒子像である。
【請求項３】
フーリエ変換で得られる複数の周波数成分に対してローパスフィルタを介して得られる値の平均値を前記演算式の各入力値とする請求項１又は２に記載の流体計測方法。
【請求項４】
フーリエ変換で得られる複数の周波数成分について中央値を算出して前記演算式の各入力値とする請求項１又は２に記載の流体計測方法。
【請求項５】
前記第１画像における前記第１粒子像のｘｙ座標、前記第２画像における前記第２粒子像のｘｙ座標、前記取得した奥行き変位量、及び前記微小時間に基づいて、計測対象の粒子の速度ベクトルを算出するステップ、を含む請求項１〜４のいずれかに記載の流体計測方法。
【請求項６】
第１時点から微小時間経過後の第２時点までの奥行き変位量が既知である粒子を同一のカメラを用いて撮像して得られる前記第１時点での第１画像および前記第２時点での第２画像から、前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出するステップと、
前記第１粒子像および前記第２粒子像をフーリエ変換するステップと、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出するステップと、
前記既知の奥行き変位量と前記算出した評価値とを関連づけて流体計測用データを生成するステップと、を含む流体計測用データ生成方法。
【請求項７】
請求項１〜６のいずれかに記載の方法をコンピュータに実行させるコンピュータプログラム。
【請求項８】
第１時点と当該第１時点から微小時間経過後の第２時点とで同一のカメラを用いて粒子を撮像して得られる第１画像および第２画像から、計測対象となる前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出する抽出部と、
前記第１粒子像および前記第２粒子像をフーリエ変換するフーリエ変換部と、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出する評価値算出部と、
既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データを参照し、算出した評価値に対応する奥行き変位量を特定する奥行き変位量特定部と、を備える流体計測装置。
【請求項９】
前記評価値算出部は、前記評価値を式（８）により算出する請求項８に記載の流体計測装置。
【数２】
0006399517.tif 000009
ただし、ΔＤ（ξ，η，Δｔ）は上記評価値であり、ξは、前記第１粒子像および粒子像のｘ方向における空間周波数であり、ηは、前記第１粒子像および第２粒子像のｙ方向における空間周波数であり、Δｔは前記微小時間であり、ｔは前記第１時点であり、（ｔ＋Δｔ）は前記第２時点であり、Ｇ（ξ，η，ｔ）はフーリエ変換した第１粒子像であり、Ｇ（ξ，η，ｔ＋Δｔ）はフーリエ変換した第２粒子像である。
【請求項１０】
前記評価値算出部は、フーリエ変換で得られる複数の周波数成分に対してローパスフィルタを介して得られる値の平均値を前記演算式の各入力値とする請求項８又は９に記載の流体計測装置。
【請求項１１】
前記評価値算出部は、フーリエ変換で得られる複数の周波数成分について中央値を算出して前記演算式の各入力値とする請求項８又は９に記載の流体計測装置。
【請求項１２】
前記第１画像における前記第１粒子像のｘｙ座標、前記第２画像における前記第２粒子像のｘｙ座標、前記取得した奥行き変位量、及び前記微小時間に基づいて、計測対象の粒子の速度ベクトルを算出する速度ベクトル算出部を更に備える流体計測装置。
【請求項１３】
第１時点から微小時間経過後の第２時点までの奥行き変位量が既知である粒子を同一のカメラを用いて撮像して得られる前記第１時点での第１画像および前記第２時点での第２画像から、前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出する抽出部と、
前記第１粒子像および前記第２粒子像をフーリエ変換するフーリエ変換部と、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出する評価値算出部と、
前記既知の奥行き変位量と前記算出した評価値とを関連づけて流体計測用データを生成する計測用データ生成部と、を備える流体計測用データ生成装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流体計測方法、流体計測装置、流体計測用データ生成方法、流体計測用データ生成装置、及びコンピュータプログラムに関する。
【背景技術】
【０００２】
流体の動きを可視化するために、三次元空間での流体の動きを精度良く計測できる方法が求められている。特許文献１には、二次元空間における流体の計測に関する記載がある。具体的には、粒子を撮像した二次元画像であって時間的に連続する画像を複数撮像しておき、粒子追跡法によって、各画像における粒子を対応付けて粒子の速度ベクトルを獲得するとの記載がある。
【０００３】
しかしながら、カメラから見た奥行き方向を含む三次元空間における粒子の変位量の計測方法については、特許文献１にはなんら記載がない。
【０００４】
粒子を撮像した二次元画像から粒子の奥行き位置を得るための一つの方法として、撮像した画像ぼけの大きさから粒子の奥行き位置を求めるDFD（Depth-from-defocus）が知られている。粒子のある１点からでた光はスクリーン上の所定の結像位置に集まるところ、結像位置からスクリーンがずれるとそのずれ量に対応する領域に光が拡がり、画像ぼけとなる。DFDは、画像ぼけの大きさを画像処理で求め、画像ぼけの大きさに基づき結像位置からのスクリーン位置のずれ量を算出し、このずれ量とスクリーン位置とから結像位置がわかるので、レンズと粒子間の距離がわかるという原理を利用した手法である。DFDを利用した計測では、一般的に、焦点位置の異なる２台のカメラを用いて同じ対象物（粒子）を同時に撮像し、照明強度のバラツキや粒子の散乱条件の不均一等の影響を相殺して奥行き位置を求める手法がとられる。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開平１０−２２１３５７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上記の通り、ＤＦＤでは流体（対象物）の奥行き位置を算出することができるが、流体の奥行き変位量を直接算出することができない。
【０００７】
ＤＦＤを用いて流体の単位時間あたりの奥行き変位量を算出するためには、焦点位置の異なる２台のカメラで同じ粒子を或る時点ｔで同時に撮像し、或る時点ｔにおける２枚の画像に基づき粒子の奥行き位置Ｚ（ｔ）を算出するとともに、或る時点（ｔ＋Δｔ）においても２台のカメラで同じ粒子を同時に撮像し、或る時点（ｔ＋Δｔ）における２枚の画像に基づき粒子の奥行き位置Ｚ（ｔ＋Δｔ）を算出し、２つの奥行き位置から奥行き変位量を得ることが考えられる。
【０００８】
しかしながら、上記手法で得られる奥行き変位量には、奥行き位置を算出する際に含まれる誤差や粒子の対応付け誤差などのノイズが含まれることが多く、奥行き位置の時間微分により奥行き方向の粒子速度を算出すれば、微分（厳密には差分）による極小評価でノイズが算出結果に多大な影響を与えてしまい、適切な計測ができないことが考えられる。
【０００９】
本発明は、このような課題に着目してなされたものであって、その目的は、奥行き変化量を直接算出可能であり且つ精度を向上させた流体計測方法、流体計測装置、流体計測用データ生成方法、流体計測用データ生成装置、及びコンピュータプログラムを提供することである。
【課題を解決するための手段】
【００１０】
本発明は、上記目的を達成するために、次のような手段を講じている。
【００１１】
すなわち、本発明の流体計測方法は、第１時点と当該第１時点から微小時間経過後の第２時点とで同一のカメラを用いて粒子を撮像して得られる第１画像および第２画像から、計測対象となる前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出するステップと、前記第１粒子像および前記第２粒子像をフーリエ変換するステップと、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出するステップと、既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データを参照し、算出した評価値に対応する奥行き変位量を特定するステップと、を含む。
【００１２】
本発明の流体計測装置は、第１時点と当該第１時点から微小時間経過後の第２時点とで同一のカメラを用いて粒子を撮像して得られる第１画像および第２画像から、計測対象となる前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出する抽出部と、前記第１粒子像および前記第２粒子像をフーリエ変換するフーリエ変換部と、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出する評価値算出部と、既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データを参照し、算出した評価値に対応する奥行き変位量を特定する奥行き変位量特定部と、を備える。
【００１３】
このように、計測対象となる第１時点での第１粒子像および第２時点での第２粒子像を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすれば、各々の粒子像に含まれる合焦点画像成分同士が同じであるとして相殺され、画像ぼけ成分同士の除算値に対応する評価値が算出できる。画像ぼけ成分は、対象物までの距離に対応するので、画像ぼけ成分の除算値に対応する評価値を既知の奥行き変位量と予め関連づけた流体計測用データを用いて、評価値から奥行き変位量を直接算出可能となる。
【００１４】
本発明の流体計測用データ生成方法は、第１時点から微小時間経過後の第２時点までの奥行き変位量が既知である粒子を同一のカメラを用いて撮像して得られる前記第１時点での第１画像および前記第２時点での第２画像から、前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出するステップと、前記第１粒子像および前記第２粒子像をフーリエ変換するステップと、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出するステップと、前記既知の奥行き変位量と前記算出した評価値とを関連づけて流体計測用データを生成するステップと、を含む。
【００１５】
本発明の流体計測用データ生成装置は、第１時点から微小時間経過後の第２時点までの奥行き変位量が既知である粒子を同一のカメラを用いて撮像して得られる前記第１時点での第１画像および前記第２時点での第２画像から、前記第１時点での第１粒子像および前記第２時点での第２粒子像を抽出する抽出部と、前記第１粒子像および前記第２粒子像をフーリエ変換するフーリエ変換部と、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出する評価値算出部と、前記既知の奥行き変位量と前記算出した評価値とを関連づけて流体計測用データを生成する計測用データ生成部と、を備える。
【００１６】
このように、奥行き変位量が既知である第１時点での第１粒子像および第２時点での第２粒子像を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすれば、各々の粒子像に含まれる合焦点画像成分同士が同じであるとして相殺され、画像ぼけ成分同士の除算値に対応する評価値が算出できる。画像ぼけ成分は、対象物までの距離に対応する。画像ぼけ成分の除算値に対応する評価値を既知の奥行き変位量と関連づけて流体計測用データを生成するので、このデータを用いれば、奥行き変位量が未知である粒子の撮像結果から得られた評価値に基づき奥行き変位量を知ることが可能となる。
【図面の簡単な説明】
【００１７】
【図１】本発明の流体計測装置を示すブロック図。
【図２】第１画像から第１粒子像の抽出、第２画像から第２粒子像の抽出に関する説明図。
【図３】流体計測用データに関する説明図。
【図４】流体計測装置が実行するフローチャート。
【図５】本発明の流体計測用データ生成装置を示すブロック図。
【図６】流体計測用データ生成装置が実行するフローチャート。
【発明を実施するための形態】
【００１８】
以下、本発明の一実施形態について、図面を参照して説明する。
【００１９】
［流体計測装置］
本実施形態の流体計測装置１は、例えば粒子の三次元流速分布などを得るために、計測対象となる粒子を撮像した二次元画像を用いることで三次元空間における粒子の流動（変位量、速度など）を計測する装置である。具体的には、或る第１時点ｔにおいて単眼の固定カメラを用いて粒子を撮像した第１画像と、第１時点から微小時間（Δｔ）経過後の第２時点（ｔ＋Δｔ）において同一の単眼の固定カメラを用いて同じ粒子を撮像した第２画像と、を用いる。流体計測装置１は、同一粒子像について第１画像におけるＸＹ座標と第２画像におけるＸＹ座標とに基づき、画像の面内方向（Ｘ方向、Ｙ方向）の変位量を算出し、ＸＹ方向に直交する奥行き方向（Ｚ方向）の変位量を後述する処理を実行して算出し、ＸＹＺ方向それぞれの変位量と微小時間とに基づいて粒子の速度ベクトルを算出する。
【００２０】
図１に示すように、流体計測装置１は、入力部１０と、抽出部１１と、フーリエ変換部１２と、評価値算出部１３と、流体計測用データ１４ａを記憶する記憶部１４と、奥行き変位量特定部１５と、速度ベクトル算出部１６と、を有する。これら各部１０〜１６は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置においてＣＰＵが予め記憶されている図示しない処理ルーチンを実行することによりソフトウェア及びハードウェアが協働して実現される。
【００２１】
図１に示す入力部１０は、第１画像および第２画像を装置１（コンピュータ）に取り込み、記憶部１４に記憶させる。入力部１０は、例えばカメラと装置１とを接続して、カメラで撮像した第１画像および第２画像を自動的に記憶部１４に記憶するように構成してもよいし、ユーザからの操作をトリガーとして、装置１に接続した記憶媒体や通信ネットワークを介して第１画像および第２画像を取り込み、記憶部１４に記憶するように構成してもよい。
【００２２】
上述の通り、第１画像は、第１時点ｔにカメラで粒子を撮像して得られる画像である。第２画像は、第１時点ｔから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）に同一のカメラで同一の粒子を撮像して得られる画像である。図２に第１画像および第２画像の例を示している。図２の例では、第１画像および第２画像に複数の粒子が写っているので、注目する粒子に矢印を示している。図２の例の撮像環境や条件は次の通りであるが、これに限定されるものではない。
照明：室内灯
撮像粒子：オルガソール（粒径４８〜５２ミクロン、密度１．０３ｇ／ｃｃ）
カメラレンズ：５０ｍｍ標準（焦点距離５０ｍｍ）
レンズ径：５０ｍｍ（フィルタ径４９ｍｍ）
レンズ絞り：２．８
シャッタースピード：１／１５秒
構成画素数：２４５６×２０５８画素
セルサイズ：３．４５×３．４５ミクロン
【００２３】
図１に示す抽出部１１は、図２に示すように、第１画像および第２画像から、計測対象となる第１時点ｔでの第１粒子像および第２時点（ｔ＋Δｔ）での第２粒子像を抽出する。抽出方法は、既知の種々の方法を利用できる。本実施形態では、第１時点ｔの第１画像中の全ての粒子像を粒子マスク法で抽出し、抽出した位置を測定点として、第２時点（ｔ＋Δｔ）の対応する粒子を直接相互相関法で対応付けて抽出した。粒子マスク法は、標準的な粒子像を予め用意しておき、その画像との相関関数の極大値位置を粒子位置とする方法である。直接相互相関法は、輝度値パターンの類似度を相互相関で評価する方法である。これ以外の一般的な動画像解析方法を利用可能である。なお、抽出部１１が第１画像から第１粒子像を抽出することによって、第１画像における第１粒子像のｘｙ座標が求まる。同様に、抽出部１１が第２画像から第２粒子像を抽出することによって、第２画像における第２粒子像のｘｙ座標が求まる。
【００２４】
図１に示すフーリエ変換部１２は、第１粒子像および第２粒子像をフーリエ変換する。
【００２５】
図１に示す評価値算出部１３は、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した前記第１粒子像の値および前記第２粒子像の値のいずれか一方に対する他方の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値を算出する。
【００２６】
上記評価値の算出方法について原理を説明する。
焦点ずれの粒子像をｇ（ｘ，ｙ）とし、粒子の合焦点画像をｆ（ｘ，ｙ）とし、画像ぼけ関数をｂ（ｘ，ｙ）とすると、式（１）が成立する。
ｇ（ｘ，ｙ）＝ｆ（ｘ，ｙ）＊ｂ（ｘ，ｙ）    …（１）
＊は畳み込み積分を表すコンボリューション演算を示す。
式（１）の両辺をフーリエ変換すると、式（２）となる。
Ｇ（ξ，η）＝Ｆ（ξ，η）Ｂ（ξ，η）    …（２）
画像ぼけ関数ｂ（ｘ，ｙ）は、式（３）のように一般的にガウス分布で表現できる。
【数１】
0006399517.tif 000002
なお、σ２が分散である。
式（３）をフーリエ変換すると、式（４）となる。
【数２】
0006399517.tif 000003
式（４）を分散σ２について解くと、式（５）となる。
【数３】
0006399517.tif 000004
第１粒子像Ｇ１及び第２粒子像Ｇ２について、微小時間が極めて小さいために各々の合焦点画像Ｆ１，Ｆ２が同じであるとすれば、各々の粒子像をフーリエ変換した値を除算することで合焦点画像成分同士を相殺でき、画像ぼけ成分同士の除算値を得ることができる。式（２）を用いて式（６）が導ける。
【数４】
0006399517.tif 000005
第１粒子像Ｇ１及び第２粒子像Ｇ２に対し、式（５）の分散の差を求め、式（６）を適用すると、式（７）となる。
【数５】
0006399517.tif 000006
ここで、評価値をΔＤ（ξ，η，Δｔ）と表し、第１時点ｔでの第１粒子像の表記をＧ１（ξ，η）からＧ（ξ，η，ｔ）とし、第２時点（ｔ＋Δｔ）での第２粒子像の表記をＧ２（ξ，η）からＧ（ξ，η，ｔ＋Δｔ）とすれば、式（７）は式（８）のように表せる。
【数６】
0006399517.tif 000007
ξ：第１粒子像および第２粒子像のｘ方向における空間周波数
η：第１粒子像および第２粒子像のｙ方向における空間周波数
Δｔ：微小時間
ｔ：第１時点
ｔ＋Δｔ：第２時点
Ｇ（ξ，η，ｔ）：フーリエ変換した第１粒子像
Ｇ（ξ，η，ｔ＋Δｔ）：フーリエ変換した第２粒子像
【００２７】
第１粒子像および第２粒子像をフーリエ変換すれば、画素数分の複数の周波数成分が画素の数得られる。これらは理論的にはいずれも同じ値を取るので、複数の周波数のうち少なくとも１つの周波数成分を選択して、式（８）の演算式に入力すればよい。
【００２８】
ただし、現実にはノイズによるバラツキが考えられるので、精度を向上させるために、例えば次の２通りの方法が考えられる。一つは、フーリエ変換で得られる複数の周波数に対してローパスフィルタを介して得られる値の平均値を演算式の入力値とすることが考えられる。このようにすれば、ノイズが多い高周波数成分を除いた平均値を用いるので、精度を向上させることが可能となる。
【００２９】
２つ目は、フーリエ変換で得られる複数の周波数成分について中央値を算出して演算式の入力値とすることが考えられる。理論的には全ての周波数成分の各値が同じであるところ、中央値を入力値とすることでバラツキが大きな特異値を除外でき、精度を向上させることが可能となる。
【００３０】
なお、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値「Ｇ（ξ，η，ｔ＋Δｔ）／Ｇ（ξ，η，ｔ）」を入力値の一つとしているが、これに限定されない。評価値をマイナスの値とするのであれば、フーリエ変換した第１粒子像の値に対する第２粒子像の値の除算値「Ｇ（ξ，η，ｔ）／Ｇ（ξ，η，ｔ＋Δｔ）」を入力値としてもよい。
【００３１】
図１に示す記憶部１４は流体計測用データ１４ａを記憶する。流体計測用データ１４ａは、図３に示すように、既知の奥行き変位量Δｚと評価値ΔＤとを予め関連づけたデータである。このデータは、後述する流体計測用データ生成装置又は方法によって生成しておく。本実施形態では、Δｚ＝ａ・ΔＤ＋ｂといった直線を表す関数（ａ、ｂは定数）で表現されているが、これに限定されない。評価値と奥行き変位量を対応付けたテーブル形式のデータであってもよい。
【００３２】
図１に示す変位量特定部１５は、流体計測用データ１４ａを参照し、評価値算出部１３が算出した評価値ΔＤに対応する奥行き変位量Δｚを特定する。
【００３３】
図１に示す速度ベクトル算出部１６は、第１画像における第１粒子像のｘｙ座標、第２画像における第２粒子像のｘｙ座標、取得した奥行き変位量Δｚ、及び微小時間Δｔに基づいて、計測対象の粒子の速度ベクトルを算出する。具体的には、抽出部１１の処理で判明した上記ｘｙ座標から、ｘ方向の変位量Δｘおよびｙ方向の変位量Δｙを算出し、流速ベクトル（ｕ，ｖ，ｗ）＝（Δｘ／Δｔ，Δｙ／Δｔ，Δｚ／Δｔ）を算出する。
【００３４】
＜流体計測方法＞
流体計測方法を、図４を用いて説明する。
【００３５】
まず、ステップＳＴ１において、第１時点ｔと当該第１時点ｔから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）とで同一のカメラを用いて粒子を撮像し、第１画像および第２画像を取得する。このステップは、装置で自動的に実行してもよいし、人の手で行ってもよい。
【００３６】
次のステップＳＴ２において、流体計測装置１の入力部１０は、第１画像および第２画像を記憶部１４に記憶する。
【００３７】
次のステップＳＴ３において、抽出部１１は、図２に示すように、第１画像および第２画像から、計測対象となる第１時点での第１粒子像および第２時点での第２粒子像を抽出する。
【００３８】
次のステップＳＴ４において、フーリエ変換部１２は、第１粒子像および第２粒子像をフーリエ変換する。
【００３９】
次のステップＳＴ５において、評価値算出部１３は、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値ΔＤを算出する。具体的には、式（８）を演算する。本実施形態では、評価値算出部１３は、フーリエ変換で得られる複数の周波数成分（全周波数成分）に対してローパスフィルタを介して得られる値の平均値を演算式（８）の入力値としている。勿論、これに限られない。例えば、評価値算出部１３が、フーリエ変換で得られる複数の周波数成分（全周波数成分）について中央値を算出して演算式（８）の入力値とするようにしてもよい。また、例えば、評価値算出部１３が、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を所定ルールに従って選択し、選択した周波数成分の値を演算式（８）の入力値にしてもよい。
【００４０】
次のステップＳＴ６において、変位量特定部１５は、既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データ１４ａを参照し、算出した評価値ΔＤに対応する奥行き変位量Δｚを特定する。
【００４１】
次のステップＳＴ７において、速度ベクトル算出部１６は、第１画像における第１粒子像のｘｙ座標、第２画像における第２粒子像のｘｙ座標、取得した奥行き変位量Δｚ、及び微小時間Δｔに基づいて、計測対象の粒子の速度ベクトル（ｕ，ｖ，ｗ）を算出する。
【００４２】
ステップＳＴ３〜７は、画像中の一つの粒子についての処理であるので、画像中に複数の粒子がある場合には、粒子毎にステップＳＴ３〜７を繰り返し実行すればよい。
【００４３】
＜流体計測用データ生成装置＞
図５に示すように、流体計測用データ生成装置２は、上記装置１が使用する流体計測用データ１４ａを生成する装置である。具体的には、奥行き変位量が既知である粒子を撮像した画像を用いて、式（８）により評価値ΔＤを算出し、算出した評価値ΔＤと既知の奥行き変位量Δｚを関連づけてデータとして保存する。流体計測用データ生成装置２について、流体計測装置１と同一部分には同じ符号を付して説明を省略する。
【００４４】
図５に示すように、流体計測用データ生成装置２は、入力部１０と、抽出部１１と、フーリエ変換部１２と、評価値算出部１３と、奥行き変位量特定部１５と、計測用データ生成部２０と、を有する。流体計測装置１との違いについて、流体計測装置１の入力部１０と抽出部１１は、変位量が未知である粒子が撮像された第１画像と第２画像を取り合うのに対し、流体計測用データ生成装置２の入力部１０と抽出部１１は、奥行き方向の変位量が既知である粒子が撮像された第１画像と第２画像を取り合う点で異なる。それ以外の各部１２〜１５は同一である。また、計測用データ生成部２０が設けられている。
【００４５】
図５に示す計測用データ生成部２０は、既知の奥行き変位量Δｚと算出した評価値ΔＤとを関連づけて流体計測用データ１４ａを生成する。具体的には、例えば図３に示すように、奥行き変位量が既知である複数の粒子を撮像して、粒子毎に評価値ΔＤを算出し、プロットする。図３では、一つの粒子について点で示しており、奥行き変位量が１．０ｍｍであっても粒子によって評価値ΔＤにバラツキがでる。また、奥行き変位量を異ならせた複数の粒子を撮像し、複数の奥行き変位量について各々の粒子毎に評価値ΔＤを算出しておく。図３は、奥行き変位量Δｚ＝０．１ｍｍ、０．２ｍｍ、０．３ｍｍ、０．４ｍｍ、…、０．９ｍｍ、１．０ｍｍにおいて、複数の粒子の評価値ΔＤを算出し、プロットした図である。計測用データ生成部２０は、これら複数の評価値ΔＤを用いて最小二乗法などで回帰直線を求め、流体計測用データ１４ａを生成する。本実施形態では、流体計測用データ１４ａが直線を表す関数としているが、テーブル形式であってもよい。
【００４６】
上記評価値算出部１３は、式（８）の演算式を用いるにあたり、第１粒子像および第２粒子像をフーリエ変換して得られる複数の周波数成分のうち少なくとも１つの周波数成分を選択する手法として２通り紹介したが、更に次の手法が可能である。すなわち、空間的に近い位置にある複数の粒子の平均値を用いる方法である。具体的には、同じ奥行き位置にある粒子のみに照明を当てて画像を撮像し、撮像画像の対応する粒子（部分）毎に粒子像を纏めてその平均値で表す手法である。この方法によってもノイズによるバラツキを減少して、精度を向上させることが可能となる。
【００４７】
＜流体計測用データ生成方法＞
上記流体計測用データ生成装置２を用いた流体計測用データ生成方法について図６を参照して説明する。
【００４８】
まず、ステップＳＴ１０１において、第１時点ｔとから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）までの奥行き変位量Δｚが既知である粒子を同一のカメラを用いて撮像し、第１画像および第２画像を取得する。このステップは、装置で自動的に実行してもよいし、人の手で行ってもよい。
【００４９】
次のステップＳＴ１０２において、入力部１０は、第１画像および第２画像を記憶部１４に記憶する。上記ステップＳＴ２と同様である。
【００５０】
次のステップＳＴ１０３において、抽出部１１は、図２に示すように、第１画像および第２画像から、第１時点での第１粒子像および第２時点での第２粒子像を抽出する。上記ステップＳＴ３と同様である。
【００５１】
次のステップＳＴ１０４において、フーリエ変換部１２は、第１粒子像および第２粒子像をフーリエ変換する。上記ステップＳＴ４と同様である。
【００５２】
次のステップＳＴ１０５において、評価値算出部１３は、フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値ΔＤを算出する。上記ステップＳＴ５と同様である。
【００５３】
上記ステップＳＴ１０１〜１０５は、図３に示すように、複数の奥行き変位量について所定数の評価値ΔＤが得られるまで、撮像する粒子や奥行き変位量を変えて、繰り返し実行する。
【００５４】
複数の奥行き変位量について所定数の評価値ΔＤが得られれば（ステップＳＴ１０６：ＹＥＳ）、次のステップＳＴ１０７において、計測用データ生成部２０は、既知の奥行き変位量Δｚと算出した評価値ΔＤとを関連づけて流体計測用データ１４ａを生成する。この流体計測用データ１４ａは、上記流体計測装置１や流体計測方法に用いる。
【００５５】
以上のように、本実施形態の流体計測方法は、
第１時点ｔと第１時点ｔから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）とで同一のカメラを用いて粒子を撮像して得られる第１画像および第２画像から、計測対象となる第１時点での第１粒子像および第２時点での第２粒子像を抽出するステップ（ＳＴ３）と、
第１粒子像および第２粒子像をフーリエ変換するステップ（ＳＴ４）と、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値ΔＤを算出するステップ（ＳＴ５）と、
既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データ１４ａを参照し、算出した評価値ΔＤに対応する奥行き変位量を特定するステップ（ＳＴ６）と、を含む。
【００５６】
本実施形態の流体計測装置１は、
第１時点ｔと第１時点ｔから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）とで同一のカメラを用いて粒子を撮像して得られる第１画像および第２画像から、計測対象となる第１時点での第１粒子像および第２時点での第２粒子像を抽出する抽出部１１と、
第１粒子像および第２粒子像をフーリエ変換するフーリエ変換部１２と、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値ΔＤを算出する評価値算出部１３と、
既知の奥行き変位量と評価値とが予め関連づけられた流体計測用データ１４ａを参照し、算出した評価値ΔＤに対応する奥行き変位量Δｚを特定する奥行き変位量特定部１５と、を備える。
【００５７】
このように、計測対象となる第１時点での第１粒子像および第２時点での第２粒子像を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすれば、各々の粒子像に含まれる合焦点画像成分同士が同じであるとして相殺され、画像ぼけ成分同士の除算値に対応する評価値が算出できる。画像ぼけ成分は、対象物までの距離に対応するので、画像ぼけ成分の除算値に対応する評価値を既知の奥行き変位量と予め関連づけた流体計測用データ１４ａを用いて、評価値ΔＤから奥行き変位量Δｚを直接算出可能となる。
【００５８】
具体的に、本実施形態では、評価値ΔＤを式（８）により算出している。
【００５９】
本実施形態では、フーリエ変換で得られる複数の周波数成分に対してローパスフィルタを介して得られる値の平均値を演算式（８）の各入力値としている。このようにすれば、ノイズが多い高周波数成分を除いた平均値を用いるので、精度を向上させることが可能となる。
【００６０】
本実施形態では、フーリエ変換で得られる複数の周波数成分について中央値を算出して演算式（８）の各入力値としている。このように、理論的には全ての周波数成分の各値が同じであるところ、中央値を入力値とすることでバラツキが大きな特異値を除外でき、精度を向上させることが可能となる。
【００６１】
本実施形態の方法では、第１画像における第１粒子像のｘｙ座標、第２画像における第２粒子像のｘｙ座標、取得した奥行き変位量Δｚ、及び微小時間Δｔに基づいて、計測対象の粒子の速度ベクトルを算出するステップ（ＳＴ７）を含む。
【００６２】
本実施形態の装置では、第１画像における第１粒子像のｘｙ座標、第２画像における第２粒子像のｘｙ座標、取得した奥行き変位量Δｚ、及び微小時間Δｔに基づいて、計測対象の粒子の速度ベクトルを算出する速度ベクトル算出部１６を更に備える。
【００６３】
上記方法及び装置によれば、奥行き方向を含んだ三次元空間における粒子の速度ベクトルを算出可能となる。
【００６４】
本実施形態の流体計測用データ生成方法は、
第１時点ｔから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）までの奥行き変位量Δｚが既知である粒子を同一のカメラを用いて撮像して得られる第１時点ｔでの第１画像および第２時点（ｔ＋Δｔ）での第２画像から、第１時点ｔでの第１粒子像および第２時点（ｔ＋Δｔ）での第２粒子像を抽出するステップ（ＳＴ１０３）と、
第１粒子像および第２粒子像をフーリエ変換するステップ（ＳＴ１０４）と、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値ΔＤを算出するステップ（ＳＴ１０５）と、
既知の奥行き変位量Δｚと算出した評価値ΔＤとを関連づけて流体計測用データ１４ａを生成するステップ（ＳＴ１０７）と、を含む。
【００６５】
本実施形態の流体計測用データ生成装置２は、
第１時点ｔから微小時間Δｔ経過後の第２時点（ｔ＋Δｔ）までの奥行き変位量Δｚが既知である粒子を同一のカメラを用いて撮像して得られる第１時点ｔでの第１画像および第２時点（ｔ＋Δｔ）での第２画像から、第１時点ｔでの第１粒子像および第２時点（ｔ＋Δｔ）での第２粒子像を抽出する抽出部１１と、
第１粒子像および第２粒子像をフーリエ変換するフーリエ変換部１２と、
フーリエ変換で得られる複数の周波数成分のうち少なくとも１つの周波数成分を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすることによって、各々の粒子像に含まれる画像ぼけ成分同士の除算値に対応する評価値ΔＤを算出する評価値算出部１３と、
既知の奥行き変位量Δｚと算出した評価値ΔＤとを関連づけて流体計測用データ１４ａを生成する計測用データ生成部２０と、を備える。
【００６６】
このように、奥行き変位量が既知である第１時点での第１粒子像および第２時点での第２粒子像を用い、フーリエ変換した第２粒子像の値に対する第１粒子像の値の除算値を演算式の入力値の一つとすれば、各々の粒子像に含まれる合焦点画像成分同士が同じであるとして相殺され、画像ぼけ成分同士の除算値に対応する評価値が算出できる。画像ぼけ成分は、対象物までの距離に対応する。画像ぼけ成分の除算値に対応する評価値を既知の奥行き変位量と関連づけて流体計測用データを生成するので、このデータを用いれば、奥行き変位量が未知である粒子の撮像結果から得られた評価値に基づき奥行き変位量を知ることが可能となる。
【００６７】
本実施形態のコンピュータプログラムは、上記方法をコンピュータに実行させる。また、本実施形態では、コンピュータに読取り可能な記録媒体であって、上記コンピュータプログラムを記憶する記録媒体も有用である。
【００６８】
以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【００６９】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【００７０】
１…流体計測装置
１１…抽出部
１２…フーリエ変換部
１３…評価値算出部
１４ａ…流体計測用データ
１５…奥行き変位量特定部
１６…速度ベクトル算出部
２…流体計測用データ生成装置
２０…計測用データ生成部
【図１】
0006399517.tif 000010
【図２】
0006399517.tif 000011
【図３】
0006399517.tif 000012
【図４】
0006399517.tif 000013
【図５】
0006399517.tif 000014
【図６】
0006399517.tif 000015