ホーム > 特許ランキング > 国立大学法人室蘭工業大学
特開2023-19924背景指向シュリーレン法における測定分解能の向上を目的とした流れ場の密度勾配可視化方法、プログラム、システム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023019924
(43)【公開日】2023-02-09
(54)【発明の名称】背景指向シュリーレン法における測定分解能の向上を目的とした流れ場の密度勾配可視化方法、プログラム、システム
(51)【国際特許分類】
G06T 7/254 20170101AFI20230202BHJP
G01N 21/17 20060101ALI20230202BHJP
G01N 21/41 20060101ALI20230202BHJP
【ＦＩ】
G06T7/254 B
G01N21/17 A
G01N21/41 Z
【審査請求】未請求
【請求項の数】15
【出願形態】ＯＬ
(21)【出願番号】P 2021125006
(22)【出願日】2021-07-30
【新規性喪失の例外の表示】特許法第３０条第２項適用申請有り （１）２０２０年度衝撃波シンポジウムの講演予稿集にて公開 ▲１▼ウェブサイトの掲載日 令和３年３月２日 ▲２▼ウェブサイトのアドレス ｈｔｔｐｓ：／／ｃｕｆｔｒａｎｓ．ｃｈｉｂａ－ｕ．ｊｐ／ｐｕｂｌｉｃ／ＵＥＳｗＱＡ１ｈｎＥ７Ａｖ－ＵＢｗ６ｘ３Ｗ５ＤｕＴｄＯａｎｗ－ｗｒＹｄＢｒＫ－ｋ１ｏｊ９ （２）２０２０年度衝撃波シンポジウムにて発表 ▲１▼開催日 令和３年３月３日 ▲２▼集会名 ２０２０年度衝撃波シンポジウム（ウェブ開催） ▲３▼ウェブサイトのアドレス ｈｔｔｐｓ：／／ｕｓ０２ｗｅｂ．ｚｏｏｍ．ｕｓ／ｊ／８７９９６５６３８５６？ｐｗｄ＝ＵＥ９ＳＶＴｄＹＵｚＦＺＶＵｈＮｃｌＮＨＱＷＮｖＬ３ＢｏＵＴ０９
(71)【出願人】
【識別番号】504193837
【氏名又は名称】国立大学法人室蘭工業大学
(74)【代理人】
【識別番号】230117846
【弁護士】
【氏名又は名称】長友 隆典
(74)【代理人】
【識別番号】100217032
【弁理士】
【氏名又は名称】常本 俊幸
(72)【発明者】
【氏名】畠中 和明
(72)【発明者】
【氏名】廣田 光智
(72)【発明者】
【氏名】伊藤 瑞基
【テーマコード（参考）】
2G059
5L096
【Ｆターム（参考）】
2G059AA05
2G059BB01
2G059EE04
2G059FF01
2G059KK04
2G059MM01
5L096BA08
5L096DA01
5L096FA06
5L096FA39
5L096GA02
5L096GA08
5L096GA17
5L096GA55
5L096HA02
(57)【要約】 （修正有）
【課題】複雑な光学系の装置を必要せず、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の密度勾配を適切に評価し、現象をよく再現する密度勾配可視化方法及びシステムを提供する。
【解決手段】密度勾配可視化システムは、背景画像を用いて密度分布のない状態で計測領域を撮影して参照画像を取得する参照画像取得手段と、背景画像を用いて密度分布が発生した状態での計測領域を撮影して計測画像を取得する計測画像取得手段と、取得した参照画像及び計測画像に対し輪郭抽出処理を行い、参照画像から参照輪郭画像、計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出手段と、輪郭抽出手段によって生成した参照輪郭画像及び計測輪郭画像のペアに基づいて密度勾配を算出する密度勾配算出手段と、算出した密度勾配に基づいて計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力手段と、を備える。
【選択図】図２４
【解決手段】密度勾配可視化システムは、背景画像を用いて密度分布のない状態で計測領域を撮影して参照画像を取得する参照画像取得手段と、背景画像を用いて密度分布が発生した状態での計測領域を撮影して計測画像を取得する計測画像取得手段と、取得した参照画像及び計測画像に対し輪郭抽出処理を行い、参照画像から参照輪郭画像、計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出手段と、輪郭抽出手段によって生成した参照輪郭画像及び計測輪郭画像のペアに基づいて密度勾配を算出する密度勾配算出手段と、算出した密度勾配に基づいて計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力手段と、を備える。
【選択図】図２４
【特許請求の範囲】
【請求項1】
計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって、観測対象となる流れ場の密度勾配を可視化する方法であって、
前記背景画像を用いて、密度分布が発生していない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得ステップと、
前記背景画像を用いて、密度分布が発生した状態で計測領域を撮影することにより、前記計測画像を取得する計測画像取得ステップと、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記輪郭抽出ステップによって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出ステップと、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力ステップと、
を備えたことを特徴とする密度勾配可視化方法。
前記背景画像を用いて、密度分布が発生していない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得ステップと、
前記背景画像を用いて、密度分布が発生した状態で計測領域を撮影することにより、前記計測画像を取得する計測画像取得ステップと、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記輪郭抽出ステップによって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出ステップと、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力ステップと、
を備えたことを特徴とする密度勾配可視化方法。
【請求項2】
前記背景画像として、所定の空間周波数を持つ波で構成され周期的な輝度分布を有する縞画像を用いること、
を特徴とする請求項１に記載の密度勾配可視化方法。
を特徴とする請求項１に記載の密度勾配可視化方法。
【請求項3】
前記の参照輪郭画像および前記計測輪郭画像のペアに基づいて密度勾配を算出する方法として、
前記参照輪郭画像、前記参照輪郭画像を微分処理した微分画像、および前記計測輪郭画像の３種の画像から算出した、前記参照輪郭画像と前記計測輪郭画像との位相差、に基づいて密度勾配を算出すること、
を特徴とする請求項２に記載の密度勾配可視化方法。
前記参照輪郭画像、前記参照輪郭画像を微分処理した微分画像、および前記計測輪郭画像の３種の画像から算出した、前記参照輪郭画像と前記計測輪郭画像との位相差、に基づいて密度勾配を算出すること、
を特徴とする請求項２に記載の密度勾配可視化方法。
【請求項4】
前記輪郭抽出処理として、処理対象となる画素とその周囲の画素を含めた領域内の画素値を用いて、前記処理対象となる画素についての新たな画素値を算出する空間フィルタリング処理を利用すること、
を特徴とする請求項３に記載の密度勾配可視化方法。
を特徴とする請求項３に記載の密度勾配可視化方法。
【請求項5】
前記空間フィルタリング処理として、線形フィルタの畳み込みフィルタを用いる場合に、ＳＯＢＥＬフィルタ、又はプレヴィットフィルタを含む微分フィルタを用いること、
を特徴とする請求項４に記載の密度勾配可視化方法。
を特徴とする請求項４に記載の密度勾配可視化方法。
【請求項6】
前記計測輪郭画像における位相差を算出し、当該算出した位相差から推定される複数の変位候補を出力するステップと、
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力するステップと、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定するステップ、
を備えたことを特徴とする請求項３～５のいずれか１つに記載の密度勾配可視化方法。
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力するステップと、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定するステップ、
を備えたことを特徴とする請求項３～５のいずれか１つに記載の密度勾配可視化方法。
【請求項7】
前記参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成するステップと、
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出ステップと、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出ステップと、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出ステップと、を備え、
前記理論変位量算出ステップで算出した変位量の理論値と、前記元画像変位量算出ステップで算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出ステップで算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証ステップを備えたこと、
を特徴とする請求項３～６のいずれか一つに記載の密度勾配可視化方法。
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出ステップと、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出ステップと、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出ステップと、を備え、
前記理論変位量算出ステップで算出した変位量の理論値と、前記元画像変位量算出ステップで算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出ステップで算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証ステップを備えたこと、
を特徴とする請求項３～６のいずれか一つに記載の密度勾配可視化方法。
【請求項8】
コンピュータシステムにおいて、請求項１～７のいずれか１つに記載の情報処理方法を実行する密度勾配可視化プログラム。
【請求項9】
計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって、観測対象となる流れ場の密度勾配を可視化するシステムであって、
前記背景画像を用いて、密度分布のない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得手段と、
前記背景画像を用いて、密度分布が発生した状態での計測領域を撮影することにより、前記計測画像を取得する計測画像取得手段と、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記輪郭抽出手段によって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出手段と、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力手段と、
を備えたことを特徴とする密度勾配可視化システム。
前記背景画像を用いて、密度分布のない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得手段と、
前記背景画像を用いて、密度分布が発生した状態での計測領域を撮影することにより、前記計測画像を取得する計測画像取得手段と、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記輪郭抽出手段によって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出手段と、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力手段と、
を備えたことを特徴とする密度勾配可視化システム。
【請求項10】
前記背景画像として、所定の空間周波数を持つ波で構成され周期的な輝度分布を有する縞画像を用いること、
を特徴とする請求項９に記載の密度勾配可視化システム。
を特徴とする請求項９に記載の密度勾配可視化システム。
【請求項11】
前記の参照輪郭画像および前記計測輪郭画像のペアに基づいて密度勾配を算出する手段として、
前記参照輪郭画像、前記参照輪郭画像を微分処理した微分画像、および前記計測輪郭画像の３種の画像から算出した、前記参照輪郭画像と前記計測輪郭画像との位相差、に基づいて密度勾配を算出すること、
を特徴とする請求項１０に記載の密度勾配可視化システム。
前記参照輪郭画像、前記参照輪郭画像を微分処理した微分画像、および前記計測輪郭画像の３種の画像から算出した、前記参照輪郭画像と前記計測輪郭画像との位相差、に基づいて密度勾配を算出すること、
を特徴とする請求項１０に記載の密度勾配可視化システム。
【請求項12】
前記輪郭抽出処理として、処理対象となる画素とその周囲の画素を含めた領域内の画素値を用いて、前記処理対象となる画素についての新たな画素値を算出する空間フィルタリング処理を利用すること、
を特徴とする請求項１１に記載の密度勾配可視化システム。
を特徴とする請求項１１に記載の密度勾配可視化システム。
【請求項13】
前記空間フィルタリング処理として、線形フィルタの畳み込みフィルタを用いる場合に、ＳＯＢＥＬフィルタ、又はプレヴィットフィルタを含む微分フィルタを用いること、
を特徴とする請求項１２に記載の密度勾配可視化システム。
を特徴とする請求項１２に記載の密度勾配可視化システム。
【請求項14】
前記計測輪郭画像における位相差を算出し、当該算出した位相差から推定される複数の変位候補を出力する手段と、
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力する手段と、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定する手段、
を備えたことを特徴とする請求項１１～１３のいずれか１つに記載の密度勾配可視化システム。
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力する手段と、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定する手段、
を備えたことを特徴とする請求項１１～１３のいずれか１つに記載の密度勾配可視化システム。
【請求項15】
前記参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成する手段と、
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出手段と、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出手段と、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出手段と、を備え、
前記理論変位量算出手段で算出した変位量の理論値と、前記元画像変位量算出手段で算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出手段で算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証手段を備えたこと、
を特徴とする請求項１１～１４のいずれか一つに記載の密度勾配可視化システム。
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出手段と、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出手段と、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出手段と、を備え、
前記理論変位量算出手段で算出した変位量の理論値と、前記元画像変位量算出手段で算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出手段で算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証手段を備えたこと、
を特徴とする請求項１１～１４のいずれか一つに記載の密度勾配可視化システム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって密度勾配を可視化する方法、プログラム、システムに関する。
特に、従来の背景指向シュリーレン法における流れ場の密度勾配可視化方法を改良し、大幅な測定分解能の向上を可能とする新たな密度勾配可視化方法、プログラム、システムを提供することに関する。
特に、従来の背景指向シュリーレン法における流れ場の密度勾配可視化方法を改良し、大幅な測定分解能の向上を可能とする新たな密度勾配可視化方法、プログラム、システムを提供することに関する。
【背景技術】
【0002】
ロケットエンジンの噴射、衝撃波等、目には見えにくい流れ場の現象を把握するために、流れ場の可視化が必要となる。
光の速度は媒質の屈折率に依存し、屈折率は媒質の密度とともに増加するので、擾乱が発生した部分など、密度の異なる部分を通過する光の速度は異なり、密度勾配のある部分を通過する光は屈折する。
この原理を応用して密度勾配のある流れ場を光学的に可視化する方法を光学的可視化法という。
光の速度は媒質の屈折率に依存し、屈折率は媒質の密度とともに増加するので、擾乱が発生した部分など、密度の異なる部分を通過する光の速度は異なり、密度勾配のある部分を通過する光は屈折する。
この原理を応用して密度勾配のある流れ場を光学的に可視化する方法を光学的可視化法という。
【0003】
従来、光学的可視化法としては、流れ場の密度勾配を把握して可視化する手法が用いられ、例えば、シャドウグラフ法、シュリーレン法、およびマッハツェンダ干渉計による干渉法などがあるが、主に、高いＳＮ比で明瞭な映像が得られるシュリーレン法が用いられていた。
しかし、シュリーレン法はかなり大掛かりで複雑な光学装置を必要とするため、設置するための場所、撮影方向の制約、光軸調整にも経験と技術を要する等の問題があった。
しかし、シュリーレン法はかなり大掛かりで複雑な光学装置を必要とするため、設置するための場所、撮影方向の制約、光軸調整にも経験と技術を要する等の問題があった。
【0004】
このため、流れ場の密度勾配を把握して可視化する手法として、背景指向シュリーレン法（英語表記で「Background oriented schlieren法」（ 以下「ＢＯＳ法」と呼ぶ。））という手法が用いられるようになった。
ＢＯＳ法では、当初、背景画像としてランダムドットを用いて撮影した画像をもとに、相互相関法を用いて密度勾配を定量化する方法が採用されていたが、計算負荷が大きく、エラー(誤ベクトル)の除去を行う処理が不可欠であるといった課題があった。
ＢＯＳ法では、当初、背景画像としてランダムドットを用いて撮影した画像をもとに、相互相関法を用いて密度勾配を定量化する方法が採用されていたが、計算負荷が大きく、エラー(誤ベクトル)の除去を行う処理が不可欠であるといった課題があった。
【0005】
そこで、従来のＢＯＳ法で用いるランダムドット画像と相互相関法の組み合わせに代えて、特許文献１に記載の発明および、これに対応する非特許文献１に記載のＳ－ＢＯＳ法（以下、両者を合わせて「従来のＳ－ＢＯＳ法」という）のように、周期的な輝度変化を持つ背景画像を用いて、相互相関計算を経ずに密度勾配を可視化する方法が提案された。
【0006】
具体的には、特許文献１に記載の発明および、これに対応する非特許文献１に記載の従来のＳ－ＢＯＳ法では、下記のイ）とロ）の理論に基づく定式化により、密度勾配を算出している。
イ）密度勾配を検出する方向に周期的な輝度変化を持つ背景を用意し、擾乱が発生した場合の撮影画像との差分をとると、位相差に応じた輝度値の変化を算出することができる。
ロ）この時、振動の中心に対して正負の両方の輝度が現れるため、差分のみでは、画像の変異方向に応じた出力を得ることができない。そのため、輝度値の差分と、輝度の勾配との積を取ることにより、画像の変異方向に応じて、所定のオフセット値を中心として、正負に輝度値を振り分けるようにして、ある座標の輝度値を算出することとすると、最終的な出力結果の輝度値は、参照画像と測定画像との位相差φのみが変数となるsin関数で得ることができる。
イ）密度勾配を検出する方向に周期的な輝度変化を持つ背景を用意し、擾乱が発生した場合の撮影画像との差分をとると、位相差に応じた輝度値の変化を算出することができる。
ロ）この時、振動の中心に対して正負の両方の輝度が現れるため、差分のみでは、画像の変異方向に応じた出力を得ることができない。そのため、輝度値の差分と、輝度の勾配との積を取ることにより、画像の変異方向に応じて、所定のオフセット値を中心として、正負に輝度値を振り分けるようにして、ある座標の輝度値を算出することとすると、最終的な出力結果の輝度値は、参照画像と測定画像との位相差φのみが変数となるsin関数で得ることができる。
【先行技術文献】
【特許文献】
【0007】
【特許文献1】特開２０１２－１４５４３０号公報
【非特許文献1】「Background Oriented Schlieren 法に基づく密度勾配の可視化法の改良」日本機械学会論文集（B編）77巻784号（2011-12）
【発明の概要】
【発明が解決しようとする課題】
【0008】
しかし、従来のＳ－ＢＯＳ法によると、後述の位相反転の問題を避けるため、空間分解能をある程度制限した背景画像を用いて計測するしかなく（図１参照）、その結果、急激な密度変化を伴う衝撃波面について、流れ場の擾乱を可視化した結果画像の密度勾配の再現度が低いレベルに留まっていた。すなわち、空間分解能の不足によって、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の密度勾配を過小評価し、現象をよく再現できないという課題である。
【0009】
これを図示したのが、図１および図２である。
図１（ａ）～（ｃ）は、背景画像として、所定の空間周波数を持つ波（例えば正弦波）で構成され、周期的な輝度分布を有する縞画像を用いる場合を図示したものである。ここで、図１（ｂ）に示した縞画像では、縞が白と黒の領域で明瞭に区別されているが、実際には、図１（ａ）に示すように、所定の空間周波数を持つ波の輝度値のレベルに応じて白から灰色、そして黒へと濃度が変化するものであり、図１（ｂ）は、これを白黒で簡略化して表現したものである。
これを図１（ｃ）との関係でみてみると、所定の空間周波数を持つ波（例えば正弦波）の輝度値のレベルが高い個所では白に近付き、輝度値のレベルが下がるにつれて次第に灰色、ないし黒に近づくことがわかる。
図１によれば、後述の図３と比較して、解像度を抑えた背景画像が示されており、その周期もやや長めとなっていることが分かる。
なお、図では、背景画像として横縞を使用した場合を例示したが、背景画像として縦縞を利用する場合も同様である。
図１（ａ）～（ｃ）は、背景画像として、所定の空間周波数を持つ波（例えば正弦波）で構成され、周期的な輝度分布を有する縞画像を用いる場合を図示したものである。ここで、図１（ｂ）に示した縞画像では、縞が白と黒の領域で明瞭に区別されているが、実際には、図１（ａ）に示すように、所定の空間周波数を持つ波の輝度値のレベルに応じて白から灰色、そして黒へと濃度が変化するものであり、図１（ｂ）は、これを白黒で簡略化して表現したものである。
これを図１（ｃ）との関係でみてみると、所定の空間周波数を持つ波（例えば正弦波）の輝度値のレベルが高い個所では白に近付き、輝度値のレベルが下がるにつれて次第に灰色、ないし黒に近づくことがわかる。
図１によれば、後述の図３と比較して、解像度を抑えた背景画像が示されており、その周期もやや長めとなっていることが分かる。
なお、図では、背景画像として横縞を使用した場合を例示したが、背景画像として縦縞を利用する場合も同様である。
【0010】
図２（ａ）は、衝撃波と混合層を含む超音速流れ場を従来のＳ－ＢＯＳ法で可視化した結果画像を示している。図の各画素は、密度勾配の正負に応じて、正（赤）、負（青）などと色分けがされるが、モノクロ画像であるため、密度勾配が正の領域を点線で示し、それ以外の領域は負の領域であることを示すにとどめた。なお、実験装置の構成については後述するのでここでは説明を省略する。
図２（ｂ）は、図２（ａ）の密度勾配分布において、Ｙ軸に平行なＸ＝１０ｍｍのＡとＢを結ぶ線上の密度勾配値の数値計算による理論値と、従来のＳ－ＢＯＳ法による出力結果を抜き出したグラフである。
図２（ｂ）は、図２（ａ）の密度勾配分布において、Ｙ軸に平行なＸ＝１０ｍｍのＡとＢを結ぶ線上の密度勾配値の数値計算による理論値と、従来のＳ－ＢＯＳ法による出力結果を抜き出したグラフである。
【0011】
ここで、図２（ａ）によれば、Ｙ＝－３ｍｍと８ｍｍ付近に、衝撃波による大きな密度勾配があることが分かる。このうち例えば、Ｙ＝－３ｍｍ付近に対応する密度勾配は、図２（ｂ）によれば、従来のＳ－ＢＯＳ法による出力結果では、衝撃波による大きな密度勾配がある箇所では、理論値に対して５０％を下回る程度の密度勾配の値を示すに留まり、衝撃波などの急峻な擾乱を十分に再現できていないことが分かる。
このように、従来のＳ－ＢＯＳ法では、空間分解能の不足によって、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の密度勾配を過小評価し、現象をよく再現できないという課題があった。
このように、従来のＳ－ＢＯＳ法では、空間分解能の不足によって、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の密度勾配を過小評価し、現象をよく再現できないという課題があった。
【0012】
また、最終的に得られる結果画像の解像度を上げようとすると、背景画像の縞画像の間隔を狭くする必要がある（図３参照）。
しかし、単純に輝度パターンの周期を小さくして分解能を上げたとしても、今度は背景画像の変位が検出可能範囲（所定の位相差）を超えるような大きな密度勾配が生じる現象を正確に測定できなくなる問題が発生する。検出可能範囲を超えた位相差は、正負の符号の反転を起こすという問題である（以下「位相反転」という）。観測対象となる現象が、大きな密度勾配と小さな密度勾配を同時に可視化する必要がある場合には、この問題が顕著になる。
このため、空間分解能を向上させつつ位相反転を防止するといった、相反する課題を同時に解決する技術の開発が望まれていた。
しかし、単純に輝度パターンの周期を小さくして分解能を上げたとしても、今度は背景画像の変位が検出可能範囲（所定の位相差）を超えるような大きな密度勾配が生じる現象を正確に測定できなくなる問題が発生する。検出可能範囲を超えた位相差は、正負の符号の反転を起こすという問題である（以下「位相反転」という）。観測対象となる現象が、大きな密度勾配と小さな密度勾配を同時に可視化する必要がある場合には、この問題が顕著になる。
このため、空間分解能を向上させつつ位相反転を防止するといった、相反する課題を同時に解決する技術の開発が望まれていた。
【0013】
なお、特許文献１に記載の発明および非特許文献１に記載の従来のＳ－ＢＯＳ法では、位相差φが±π／２の範囲にあるときしか、密度勾配を算出することができないという不都合があり、急激な密度変化を伴う擾乱を可視化する際に、位相反転が発生しやすいという課題があった。これは、非特許文献１の（７）式において位相差φが±π／２の範囲の時に単調増加となることを利用して、移動平均を背景画像の１／２周期の範囲で行って（８）式を導出していることによるものである。
このため、計算アルゴリズムを改善して参照画像と測定画像との位相差φが±πの範囲まで密度勾配を検出可能とするＳ－ＢＯＳ法の処理技術（以下、「新たなＳ－ＢＯＳ法」という）、の開発が望まれていた。
このため、計算アルゴリズムを改善して参照画像と測定画像との位相差φが±πの範囲まで密度勾配を検出可能とするＳ－ＢＯＳ法の処理技術（以下、「新たなＳ－ＢＯＳ法」という）、の開発が望まれていた。
【0014】
このような位相反転の問題について図を用いて詳述するとともに、同じ現象を観測した場合でも、位相反転が発生する程度、および位相反転を検出して補正できるかどうかに関して、本発明と従来技術との間に、大きな差異があるので、以下、図４を用いて説明する。
図４は、図３に示すような空間分解能を上げた背景画像を使用した場合の測定画像を模式的に示した図であって、密度変化の大きさが、（イ）（ロ）（ハ）の順に大きくなることが示されている。なお、図では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
図４は、図３に示すような空間分解能を上げた背景画像を使用した場合の測定画像を模式的に示した図であって、密度変化の大きさが、（イ）（ロ）（ハ）の順に大きくなることが示されている。なお、図では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
【0015】
ここで、（イ）の部位では、位相の変異が±１／２πの範囲内に収まっており、従来のＳ－ＢＯＳ法においても位相反転は発生しないことが分かる。
他方、（ロ）の部位では、位相の変異が±πの範囲内であるので、本発明で提案する新たなＳ－ＢＯＳ法の計算手法では位相反転が発生しないが、±１／２πの範囲を超えているので、従来のＳ－ＢＯＳ法においては位相反転が発生する。
また、（ハ）の部位では、大きな変位があり、位相の変異が±πの範囲を超えているので、新たなＳ－ＢＯＳ法の計算手法でも位相反転が生じる。
もっとも、後述のように、（ハ）の部位のような変位が、元画像ではなく生成した輪郭画像において発生した場合には、位相反転を防止することができる。
他方、（ロ）の部位では、位相の変異が±πの範囲内であるので、本発明で提案する新たなＳ－ＢＯＳ法の計算手法では位相反転が発生しないが、±１／２πの範囲を超えているので、従来のＳ－ＢＯＳ法においては位相反転が発生する。
また、（ハ）の部位では、大きな変位があり、位相の変異が±πの範囲を超えているので、新たなＳ－ＢＯＳ法の計算手法でも位相反転が生じる。
もっとも、後述のように、（ハ）の部位のような変位が、元画像ではなく生成した輪郭画像において発生した場合には、位相反転を防止することができる。
【0016】
すなわち、本発明では、後述のように、単純に縞画像の間隔を狭めて解像度を上げる手法を用いずに、元画像と輪郭画像を用いて位相反転の検出と補正を行うので、輪郭画像において、±πの範囲を超える位相反転が発生した場合でも、輪郭抽出処理を行う前の元画像においては、位相反転が生じていないので、元画像と輪郭画像を用いて位相反転を防止することができる（（図３１～図３３の説明箇所参照）。
【0017】
このように、複雑な光学系の装置を必要せず、急激な密度変化を伴う擾乱の密度変化を精度よく再現できる技術、位相反転を検出して補正する技術、及び、参照画像と測定画像との位相差φが±πの範囲まで密度勾配を検出可能とする技術の開発が望まれていた。
【0018】
そこで、本発明では、図５に示すように、複雑な光学系の装置を必要せず、参照画像と計測画像に対し、輪郭抽出処理を行うことにより、参照輪郭画像と計測輪郭画像を生成し、この輪郭画像ペアを用いて新たなＳ－ＢＯＳ法の適用を行うことで、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の密度勾配を適切に評価し、現象をよく再現する手法を提供すること、を目的とする。
【0019】
また、従来のＳ－ＢＯＳ法では、縞の湾曲（参照画像に対して縞がずれた平均的な変位）を算出するのにとどまり、太さの変化には着目せず、流れ場の密度勾配を可視化するのに重要な特徴量が見過ごされてきたのに対し、本発明では、図６に示すように、輪郭抽出処理により、縞の両端の変化傾向を把握するような演算を行う関係で、揺らぎよって縞の太さが変化すること（縞の太さの変化）に着目することにより、流れ場の密度勾配を可視化するのに重要な特徴量をもれなく考慮することを目的とする。
【0020】
また、本発明では、位相反転を検出して補正するという課題に対しては、図５に示すように、元画像ペアにおける変位量と、輪郭画像ペアにおける変位量を併せて考慮することで解消することを目的とする。
【0021】
また、本発明では、従来のＳ－ＢＯＳ法において位相差φが±π／２の範囲にあるときしか密度勾配を算出することができないという不都合（位相差の検出可能範囲が±π／２の範囲に限定されること）についても、参照画像と、参照画像に対して微分処理を行った微分画像、および計測画像の３種の画像を用いた、新たなＳ－ＢＯＳ法を提案することで解決することを目的とする（図１０～図１３参照）。
【課題を解決するための手段】
【0022】
上記の目的を達成するために、第１の発明は、
計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって、観測対象となる流れ場の密度勾配を可視化する方法であって、
前記背景画像を用いて、密度分布が発生していない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得ステップと、
前記背景画像を用いて、密度分布が発生した状態で計測領域を撮影することにより、前記計測画像を取得する計測画像取得ステップと、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記輪郭抽出ステップによって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出ステップと、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力ステップと、
を備えたことを特徴とする密度勾配可視化方法である。
計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって、観測対象となる流れ場の密度勾配を可視化する方法であって、
前記背景画像を用いて、密度分布が発生していない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得ステップと、
前記背景画像を用いて、密度分布が発生した状態で計測領域を撮影することにより、前記計測画像を取得する計測画像取得ステップと、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記輪郭抽出ステップによって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出ステップと、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力ステップと、
を備えたことを特徴とする密度勾配可視化方法である。
【0023】
第２の発明は、第１の発明において、前記背景画像として、所定の空間周波数を持つ波で構成され周期的な輝度分布を有する縞画像を用いること、を特徴とする。
【0024】
第３の発明は、第２の発明において、前記の参照輪郭画像および前記計測輪郭画像のペアに基づいて密度勾配を算出する方法として、前記参照輪郭画像、前記参照輪郭画像を微分処理した微分画像、および前記計測輪郭画像の３種の画像から算出した、前記参照輪郭画像と前記計測輪郭画像との位相差、に基づいて密度勾配を算出すること、を特徴とする。
【0025】
第４の発明は、第３の発明において、前記輪郭抽出処理として、処理対象となる画素とその周囲の画素を含めた領域内の画素値を用いて、前記処理対象となる画素についての新たな画素値を算出する空間フィルタリング処理を利用すること、を特徴とする。
【0026】
第５の発明は、前記空間フィルタリング処理として、線形フィルタの畳み込みフィルタを用いる場合に、ＳＯＢＥＬフィルタ、又はプレヴィットフィルタを含む微分フィルタを用いること、を特徴とする。
【0027】
第６の発明は、第３～５の発明のいずれか１つに記載の発明において、前記計測輪郭画像における位相差を算出し、当該算出した位相差から推定される複数の変位候補を出力するステップと、
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力するステップと、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定するステップ、
を備えたことを特徴とする。
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力するステップと、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定するステップ、
を備えたことを特徴とする。
【0028】
第７の発明は、第３～６の発明のいずれか１つに記載の発明において、
前記参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成するステップと、
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出ステップと、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出ステップと、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出ステップと、を備え、
前記理論変位量算出ステップで算出した変位量の理論値と、前記元画像変位量算出ステップで算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出ステップで算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証ステップを備えたこと、
を特徴とする。
前記参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成するステップと、
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出ステップと、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出ステップと、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出ステップと、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出ステップと、を備え、
前記理論変位量算出ステップで算出した変位量の理論値と、前記元画像変位量算出ステップで算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出ステップで算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証ステップを備えたこと、
を特徴とする。
【0029】
第８の発明は、コンピュータシステムにおいて、請求項１～７のいずれか１つに記載の情報処理方法を実行する密度勾配可視化プログラムである。
【0030】
第９の発明は、
計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって、観測対象となる流れ場の密度勾配を可視化するシステムであって、
前記背景画像を用いて、密度分布のない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得手段と、
前記背景画像を用いて、密度分布が発生した状態での計測領域を撮影することにより、前記計測画像を取得する計測画像取得手段と、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記輪郭抽出手段によって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出手段と、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力手段と、
を備えたことを特徴とする密度勾配可視化システムである。
この場合において、各手段は一つの装置内に設けてもよいし、複数の装置、又はサーバー、又はクラウドに分散して設けてもよい。
計測領域の後方に背景画像を、前方に撮影手段を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した参照画像と、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した計測画像の２つの画像から密度勾配を算出することによって、観測対象となる流れ場の密度勾配を可視化するシステムであって、
前記背景画像を用いて、密度分布のない状態で計測領域を撮影することにより、前記参照画像を取得する参照画像取得手段と、
前記背景画像を用いて、密度分布が発生した状態での計測領域を撮影することにより、前記計測画像を取得する計測画像取得手段と、
取得した前記参照画像および前記計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記計測画像から計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記輪郭抽出手段によって生成した、前記参照輪郭画像および前記計測輪郭画像のペアに基づいて、密度勾配を算出する密度勾配算出手段と、
前記算出した密度勾配に基づいて、計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する結果画像出力手段と、
を備えたことを特徴とする密度勾配可視化システムである。
この場合において、各手段は一つの装置内に設けてもよいし、複数の装置、又はサーバー、又はクラウドに分散して設けてもよい。
【0031】
第１０の発明は、第９の発明において、前記背景画像として、所定の空間周波数を持つ波で構成され周期的な輝度分布を有する縞画像を用いること、を特徴とする。
【0032】
第１１の発明は、第１０の発明において、前記の参照輪郭画像および前記計測輪郭画像のペアに基づいて密度勾配を算出する手段として、前記参照輪郭画像、前記参照輪郭画像を微分処理した微分画像、および前記計測輪郭画像の３種の画像から算出した、前記参照輪郭画像と前記計測輪郭画像との位相差、に基づいて密度勾配を算出すること、を特徴とする。
【0033】
第１２の発明は、第１１の発明において、前記輪郭抽出処理として、処理対象となる画素とその周囲の画素を含めた領域内の画素値を用いて、前記処理対象となる画素についての新たな画素値を算出する空間フィルタリング処理を利用すること、を特徴とする。
【0034】
第１３の発明は、第１２の発明において、前記空間フィルタリング処理として、線形フィルタの畳み込みフィルタを用いる場合に、ＳＯＢＥＬフィルタ、又はプレヴィットフィルタを含む微分フィルタを用いること、を特徴とする。
【0035】
第１４の発明は、第１１～１３のいずれか一つに記載の発明において、前記計測輪郭画像における位相差を算出し、当該算出した位相差から推定される複数の変位候補を出力する手段と、
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力する手段と、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定する手段、
を備えたことを特徴とする。
輪郭抽出処理を行う前の元画像である前記計測画像における位相差を算出し、当該算出した元画像である前記計測画像における位相差から判断される変位量を出力する手段と、
当該出力した元画像である前記計測画像における変位量を基礎として、前記の複数の変位候補から真の変位量を特定する手段、
を備えたことを特徴とする。
【0036】
第１５の発明は、第１１～１４のいずれか一つに記載の発明において、
前記参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成する手段と、
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出手段と、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出手段と、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出手段と、を備え、
前記理論変位量算出手段で算出した変位量の理論値と、前記元画像変位量算出手段で算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出手段で算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証手段を備えたこと、を特徴とする。
前記参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成する手段と、
前記参照画像と前記精度検証用の計測画像との位相差に基づいて、元画像である前記精度検証用の計測画像の変位量を算出する元画像変位量算出手段と、
取得した前記参照画像および生成した前記精度検証用の計測画像に対し、輪郭抽出処理を行って、前記参照画像から参照輪郭画像、前記精度検証用の計測画像から精度検証用の計測輪郭画像をそれぞれ生成する輪郭抽出手段と、
前記参照輪郭画像および前記精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像である前記精度検証用の計測輪郭画像の変位量を算出する輪郭画像変位量算出手段と、
前記所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する理論変位量算出手段と、を備え、
前記理論変位量算出手段で算出した変位量の理論値と、前記元画像変位量算出手段で算出した元画像である前記精度検証用の計測画像の変位量と、前記輪郭画像変位量算出手段で算出した輪郭画像である前記精度検証用の計測輪郭画像の変位量とを比較することで、密度勾配可視化の精度を検証する精度検証手段を備えたこと、を特徴とする。
【発明の効果】
【0037】
本発明によれば、複雑な光学系の装置を必要せず、参照画像と計測画像に対し、輪郭抽出処理を行うことにより、参照輪郭画像と計測輪郭画像を生成し、この輪郭画像ペアを用いて、Ｓ－ＢＯＳ法の適用を行うことで、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の密度勾配を適切に評価し、現象をよく再現する手法を提供することができる。
【0038】
また、従来のＳ－ＢＯＳ法による出力結果では、衝撃波による大きな密度勾配がある箇所では、理論値に対して５０～７２％程度の密度勾配の値を示すに留まり、衝撃波などの急峻な擾乱を十分に再現できていなかったのに対し、本発明によれば、衝撃波による大きな密度勾配がある箇所においても、理論値の６５～９４％程度にまで再現性を高めることができ、衝撃波などの急峻な擾乱を十分に再現できるようになった。
【0039】
また、従来のＳ－ＢＯＳにおいて位相差φが±π／２の範囲にあるときしか密度勾配を算出することができないのに対し、本発明の新たなＳ－ＢＯＳ法をでは、参照画像と、参照画像に対して微分処理を行った微分画像、および計測画像の３種の画像を用いた算出手法を提案することで位相差φが±πの範囲まで位相反転を起こさない手法を提供することができる。
【0040】
また、本発明では、位相反転が発生した場合でも、元画像ペアにおける変位量と、輪郭画像ペアにおける変位量を併せて考慮することで、位相反転を検出して補正することができるという効果を奏する。
【0041】
また、本発明では、新たな評価手法を用いることで、実際の現象の確認前に、密度勾配を可視化する手法の精度を検証することができるようになった。
【0042】
さらに、参照輪郭画像と計測輪郭画像のペアを用いる本発明を活用すれば、既存の各種のＳ－ＢＯＳ法を利用した可搬型シュリーレン装置やデジタルシュリーレン装置などの既存の測定システムに対し、本発明の手法を適用したソフトウェアを搭載することで、測定分解能の向上が可能になるため、多様な活用が考えられる。
【図面の簡単な説明】
【0043】
【図1】従来技術の課題の一例を示す図であって、背景画像として、所定の空間周波数を持つ波（例えばＳＩＮ波）で構成され、周期的な輝度分布を有する縞画像を用いる場合を図示した図である。
【図2】従来技術の課題の一例を示す図であって、図２（ａ）は、衝撃波と混合層を含む超音速流れ場をＳ－ＢＯＳ法で可視化した結果画像を示している。図２（ｂ）は、図２（ａ）の密度勾配分布において、Ｙ軸に平行なＸ＝１０ｍｍのＡとＢを結ぶ線上の密度勾配値の数値計算による理論値と、従来のＳ－ＢＯＳ法による出力結果を抜き出したグラフである。
【図3】背景画像として、所定の空間周波数を持つ波（例えばＳＩＮ波）で構成され、周期的な輝度分布を有する縞画像を用いる場合を図示した図であって、図１に対し、短い周期の波で構成され、解像度の高い背景画像を用いた場合を示している。
【図4】従来技術の課題の一例を示す図であって、空間分解能を上げた背景画像を使用した場合の測定画像を模式的に示した図である。図４（イ）は１／２πの範囲内に収まる変位が発生した場合の一例を示し、図４（ロ）は１／２πを超え、πの範囲に収まる変位が発生した場合の一例を示し、図４（ハ）はπを超えた変位が発生した場合の一例を示す図である。
【図5】本発明の概要を示す図であって、参照画像と計測画像に対し、輪郭抽出処理を行うことにより、参照輪郭画像と計測輪郭画像を生成し、この輪郭画像ペアを用いてＢＯＳ法の適用を行うことを示す図である。
【図6】本発明と従来技術との着眼点の相違を示す図であって、従来のＳ－ＢＯＳ法では、縞の湾曲を算出するのにとどまり、太さの変化には着目せず、流れ場の密度勾配を可視化するのに重要な特徴量が見過ごされてきたのに対し、本発明では、縞の両端の変化傾向を把握するような演算を行う関係で、揺らぎよって縞の太さが変化すること（縞の太さの変化）に着目することにより、流れ場の密度勾配を可視化するのに重要な特徴量をできる限り考慮すること示した図である。
【図7】本発明の測定システムの構成の一例を示す正面図である。
【図8】本発明の測定システムの構成の一例を示す平面図である。
【図9】図７、図８の測定システムの構成例において、風洞２２０の内部の観測部位３００を抜き出して記載するとともに、撮影装置２１０を構成するレンズ２１２及び撮像素子２１４を抜き出して記載した、光学系からみた測定システムの構成例を示す図である。
【図10】本発明の位相差検出型ＢＯＳ法の原理を示す図である。
【図11】本発明の位相差検出型ＢＯＳ法の原理を示す図である。
【図12】本発明の位相差検出型ＢＯＳ法の原理を示す図である。
【図13】本発明の位相差検出型ＢＯＳ法の原理を示す図である。
【図14】本発明の測定システムの風洞の内部の様子と衝撃波と混合層を含む超音速流れ場を可視化した様子を示す図である。
【図15】本発明において輪郭画像を生成するための輪郭抽出処理の原理を示す図である。
【図16】本発明において元画像である参照画像と参照画像に対して輪郭抽出処理を行って生成した参照輪郭画像の様子を示す図である。
【図17】本発明の輪郭抽出処理を行い作成した輪郭画像（参照輪郭画像）とその元画像（参照画像）の一部、および両者の縞画像の輝度値の分布（所定の空間周波数を持つ波）を比較した図である。
【図18】本発明の輪郭抽出処理を行い作成した輪郭画像（計測輪郭画像）とその元画像（計測画像）の一部、および両者の縞画像の輝度値の分布（所定の空間周波数を持つ波）を比較した図である。
【図19】衝撃波と混合層を含む超音速流れ場において、得られた画像に対し、Ｓ－ＢＯＳ法を適用して密度勾配を算出する過程で得た位相差をもとに結果画像を出力したものであって、図１９（ａ）は、元画像（参照画像および計測画像のペア）に対してＳ－ＢＯＳ法を適用して密度勾配を算出する過程で算出した位相差をプロットした画像である。図１９（ｂ）は、輪郭抽出処理を行った輪郭画像（参照輪郭画像および計測輪郭画像のペア）に対してＳ－ＢＯＳ法を適用して密度勾配を算出する過程で算出した位相差をプロットした画像である。
【図20】図２０（ａ）は、衝撃波と混合層を含む超音速流れ場において、Ｓ－ＢＯＳ法を適用して可視化した様子を示す図であって、密度勾配の変化が急峻な衝撃波の存在箇所を大まかに示すとともに、衝撃波を含む断面ＡＢの位置を示す図である。図２０（ｂ）（ｃ）は、密度勾配の変化が急峻な衝撃波を含む断面ＡＢについて、元画像に基づいてＳ－ＢＯＳ法を適用した場合の密度勾配と、輪郭画像に基づいてＳ－ＢＯＳ法を適用した場合の密度勾配の再現率の違いを、グラフで定量的に示した図である。
【図21】本発明で提案する新たな精度検証方法の原理を示す図である。
【図22】本発明の改良されたＢＯＳ法に対して、新たな精度検証方法を適用して得た、精度検証の結果（本発明の効果）を定量的に示す図である。
【図23】本発明の改良されたＢＯＳ法に対して、新たな精度検証方法を適用して得た、精度検証の結果（本発明の効果）を定量的に示す図である。
【図24】本発明の密度勾配可視化システムの各構成手段の一例を示した機能ブロック図である。
【図25】本発明の密度勾配可視化システムのシステム構成の一例を示す図であって、画像処理等を行う情報処理システムがスタンドアローン構成の場合のシステム構成を示す図である。
【図26】本発明の密度勾配可視化システムのシステム構成の一例を示す図であって、画像処理等を行う情報処理システムの各構成手段がネットワークなどを介してクライアントサーバ、クラウド、ピアツーピアなど各装置に分散して配置されている場合のシステム構成を示す図である。
【図27】本発明の密度勾配可視化処理の全体の一例を示すフローチャートである。
【図28】本発明の輪郭抽出処理の一例を示すフローチャートである。
【図29】本発明の密度勾配算出処理の一例を示すフローチャートである。
【図30】本発明において密度勾配から、流れ場の擾乱を可視化した結果画像を出力する処理の一例を示すフローチャートである。
【図31】本発明において、輪郭抽出処理を行う前の元画像である計測画像における変位量と、輪郭抽出処理を行った計測輪郭画像における変位量を模式的に示した図であって、周期が約半分になった計測輪郭画像では±πの範囲を超えて位相反転が発生しているが、輪郭抽出処理を行う前の元画像である計測画像では位相反転が発生していないことを示す図である。
【図32】本発明において、輪郭抽出処理を行った計測輪郭画像における位相差から推定される複数の変位候補について、輪郭抽出処理を行う前の元画像である計測画像における位相差から判断される変位量を基礎として、複数の変位候補から真の変位量を特定する処理の仕組みを示した図である。
【図33】本発明の位相反転検出および位相反転防止処理の一例を示すフローチャートである。
【図34】本発明の精度検証処理の一例を示すフローチャートである。
【図35】本発明において、輪郭抽出処理に用いるフィルタリング処理として、線形フィルタである畳み込みフィルタを用いる場合のフィルタの種別、サイズおよびフィルタの重みとなるフィルタの値（以下「パラメータ」ともいう）を格納したフィルタテーブルの一例を示す図である。
【0044】
＜用語の説明＞
◇流れ場：液体、気体又は粒子などが流れている場をいい、３次元の流れ場、２次元の流れ場などがある。本発明で取り扱うのは、主に２次元の流れ場であるが、これに限定されない。
◇擾乱：一般的には定常状態からの乱れのことをいうが、本発明では、流れ場において時間とともに刻々と変化する気体や微粒子の乱れをいう。
◇衝撃波：ロケットやジェットエンジンの噴射又は気体中で火薬の爆発などにより局部的に急激な強い圧縮を生じたとき、この強い圧力変化は音速以上の速度で気体中を伝播する。 これを衝撃波といい、波の前方と後方で気体の圧力や密度が不連続的に変化する性質を有する。
◇流れ場：液体、気体又は粒子などが流れている場をいい、３次元の流れ場、２次元の流れ場などがある。本発明で取り扱うのは、主に２次元の流れ場であるが、これに限定されない。
◇擾乱：一般的には定常状態からの乱れのことをいうが、本発明では、流れ場において時間とともに刻々と変化する気体や微粒子の乱れをいう。
◇衝撃波：ロケットやジェットエンジンの噴射又は気体中で火薬の爆発などにより局部的に急激な強い圧縮を生じたとき、この強い圧力変化は音速以上の速度で気体中を伝播する。 これを衝撃波といい、波の前方と後方で気体の圧力や密度が不連続的に変化する性質を有する。
【0045】
◇混合層：異なる一様な速度を持つ二つの流れが混じり合う領域をいう。
◇密度分布：高速気流では気体の圧縮性のために、圧縮された部分の密度が高くなり、場所によって密度の変化が発生するが、このような空間的な密度変化の分布のこと密度分布という。
◇密度勾配とは、密度分布の１階微分をいい、観測者に対して直交する面における直交座標系のＸ方向の流れにおいて密度ρがＹ方向の勾配をもつとすると、背景画像として横縞を使用してＹ軸方向の変位を検出しようとする場合には、密度勾配はδρ／δｙで表される(図９のＸ、Ｙ方向の矢印を参照)。
◇流れ場の可視化：観測者に向かう方向であるz 方向に入射させた光は屈折率（すなわち密度）が増す方向に屈折し、その屈折角（歪み）は密度勾配と一定の関係があることから、流れ場の密度勾配によって生じた背景画像の歪みを測定することで流れ場を可視化することが可能となる(図９のＺ方向の矢印を参照)。
◇密度分布：高速気流では気体の圧縮性のために、圧縮された部分の密度が高くなり、場所によって密度の変化が発生するが、このような空間的な密度変化の分布のこと密度分布という。
◇密度勾配とは、密度分布の１階微分をいい、観測者に対して直交する面における直交座標系のＸ方向の流れにおいて密度ρがＹ方向の勾配をもつとすると、背景画像として横縞を使用してＹ軸方向の変位を検出しようとする場合には、密度勾配はδρ／δｙで表される(図９のＸ、Ｙ方向の矢印を参照)。
◇流れ場の可視化：観測者に向かう方向であるz 方向に入射させた光は屈折率（すなわち密度）が増す方向に屈折し、その屈折角（歪み）は密度勾配と一定の関係があることから、流れ場の密度勾配によって生じた背景画像の歪みを測定することで流れ場を可視化することが可能となる(図９のＺ方向の矢印を参照)。
【0046】
◇シュリーレン法：流れ場の中に屈折率がわずかに異なる部分があるとき、光線の進行方向の変化を利用してその部分が明確に見えるようにする光学的手法をいう。屈折率の変化を、平行度の高い光により大きな明暗の差に変えて観測する手法によって可視化する。凹面鏡、ミラー、レンズ、ナイフエッジなどを複数組み合わせた規模の大きな測定システムが必要となる。
◇背景指向シュリーレン法（ＢＯＳ法）：流れ場に擾乱が発生するときに流体の屈折率が変化する現象を利用して、流れ場の密度変化を光学的に可視化する測定方法の一つである。流れ場の密度勾配によって生じた背景画像の歪み（変位量）を測定することで流れ場を可視化する。このため、必要となる機材は、適当な光源と背景画像、及びデジタルカメラなどの撮影装置のみで足り、従来のシュリーレン法で必要とされる凹面鏡やミラー、レンズなどの光学部品を必要とせず小規模な測定システムで構成できる。
◇Ｇｒａｄｓｔｏｎｅ－Ｄａｌｅ則：物質の屈折率は温度や密度の影響を受けるので、気体の種類によって定まるグラッドストーン・デール定数Ｋを用いて密度と屈折率を関連づける法則をいう。屈折率をｎ、密度をρとしたとき、Ｋ＝（ｎ－１）／ρで表現される。
◇背景指向シュリーレン法（ＢＯＳ法）：流れ場に擾乱が発生するときに流体の屈折率が変化する現象を利用して、流れ場の密度変化を光学的に可視化する測定方法の一つである。流れ場の密度勾配によって生じた背景画像の歪み（変位量）を測定することで流れ場を可視化する。このため、必要となる機材は、適当な光源と背景画像、及びデジタルカメラなどの撮影装置のみで足り、従来のシュリーレン法で必要とされる凹面鏡やミラー、レンズなどの光学部品を必要とせず小規模な測定システムで構成できる。
◇Ｇｒａｄｓｔｏｎｅ－Ｄａｌｅ則：物質の屈折率は温度や密度の影響を受けるので、気体の種類によって定まるグラッドストーン・デール定数Ｋを用いて密度と屈折率を関連づける法則をいう。屈折率をｎ、密度をρとしたとき、Ｋ＝（ｎ－１）／ρで表現される。
【0047】
◇Ｓ－ＢＯＳ法：Ｓｉｍｐｌｉｆｉｅｄ ＢＯＳの略であり、ＢＯＳ法のうち、背景画像として周期的輝度パターンを有する画像を用いる手法をいう。
◇周期的輝度パターンを有する画像：所定の空間周波数を持つ波で構成され、周期的な輝度分布を有する縞画像をいう。周期的な輝度変化を有するものであれば足り、一つの正弦波や余弦波で構成される場合のほか複数の周期の正弦波が重畳して構成された合成波であってもかまわない。
◇位相反転：単純に輝度パターンの周期を小さくして分解能を上げた場合、背景画像の変位が検出可能範囲（所定の位相差）を超えるような大きな密度勾配が生じる現象を正確に測定できなくなる問題をいう。端的には検出可能範囲を超えた位相差は、正負の符号の反転を起こすという問題である。
◇周期的輝度パターンを有する画像：所定の空間周波数を持つ波で構成され、周期的な輝度分布を有する縞画像をいう。周期的な輝度変化を有するものであれば足り、一つの正弦波や余弦波で構成される場合のほか複数の周期の正弦波が重畳して構成された合成波であってもかまわない。
◇位相反転：単純に輝度パターンの周期を小さくして分解能を上げた場合、背景画像の変位が検出可能範囲（所定の位相差）を超えるような大きな密度勾配が生じる現象を正確に測定できなくなる問題をいう。端的には検出可能範囲を超えた位相差は、正負の符号の反転を起こすという問題である。
【0048】
◇背景画像：ＢＯＳ法で用いる画像であって、計測領域の後方に配置する画像である。背景画像は、紙や布に周期的輝度パターンを有する画像を印刷したものを使用するのが一般的であるが、液晶ディスプレイ、有機ＥＬディスプレイなどのディスプレイ機器で表示した画像、又はプロジェクターを用いてスクリーンに投影した画像であっても構わない。
◇横縞、縦縞：観測者に直交する座標でＸ軸方向に延びる縞で構成される縞を横縞、Ｙ軸方向に延びる縞で構成される縞を縦縞という。横縞の背景画像は主に縦方向（Ｙ軸方向）の変位を捉えるのに効果的であり、縦縞の背景画像は主に横方向（Ｘ軸方向）の変位を捉えるのに効果的である。
◇参照画像：計測領域の後方に背景画像を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した画像をいう。
◇計測画像：計測領域の後方に背景画像を配置し、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した画像をいう。
◇生成した計測画像：参照画像に基づいて生成した画像であって、参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させて生成した精度検証用の計測画像をいう。
◇撮影装置：デジタルカメラなど、レンズと撮像素子で構成され、計測領域を撮影する装置をいう。
◇横縞、縦縞：観測者に直交する座標でＸ軸方向に延びる縞で構成される縞を横縞、Ｙ軸方向に延びる縞で構成される縞を縦縞という。横縞の背景画像は主に縦方向（Ｙ軸方向）の変位を捉えるのに効果的であり、縦縞の背景画像は主に横方向（Ｘ軸方向）の変位を捉えるのに効果的である。
◇参照画像：計測領域の後方に背景画像を配置し、観測対象となる流れ場の密度分布が発生してない状態で計測領域を撮影した画像をいう。
◇計測画像：計測領域の後方に背景画像を配置し、観測対象となる流れ場の密度分布が発生した状態で計測領域を撮影した画像をいう。
◇生成した計測画像：参照画像に基づいて生成した画像であって、参照画像において所定の幅Ａの矩形領域を取得し、当該矩形領域を所定の距離Ｂだけ平行移動させて生成した精度検証用の計測画像をいう。
◇撮影装置：デジタルカメラなど、レンズと撮像素子で構成され、計測領域を撮影する装置をいう。
【0049】
◇元画像：輪郭抽出処理を行う前の画像をいい、参照画像、計測画像、又は生成した計測画像が含まれる。
◇輪郭画像：輪郭抽出処理を行った後の画像をいう。
◇参照輪郭画像：参照画像に対し輪郭抽出処理を行った画像をいう。
◇計測輪郭画像：計測画像に対し輪郭抽出処理を行った画像をいう。
◇精度検証用の輪郭計測画像：参照画像に基づいて生成した精度検証用の計測画像に対し、輪郭抽出処理を行った画像をいう。
◇結果画像：撮影した参照画像および計測画像、又は輪郭抽出処理を行った参照輪郭画像および計測輪郭画像に対し、Ｓ－ＢＯＳ法を適用して得られた、流れ場の擾乱を可視化した画像であって、主に、流れ場の密度勾配を可視化した画像をいうが、密度勾配を算出する過程で得られた位相差や変位量を可視化した図を含む。換言すると、撮影画像に対して画像処理を行って背景画像の歪み量を測定し、画素ごとに歪み量に応じた明るさを割り当てる処理を行った画像であるともいえる。結果画像は、背景画像として横縞を使用した場合、又は、背景画像として縦縞を使用した場合のいずれか一方の計算結果を利用したものであってもよいし、横縞と縦縞の両方の計算結果を二乗平均等で重畳した結果であってもよい。
◇輪郭画像：輪郭抽出処理を行った後の画像をいう。
◇参照輪郭画像：参照画像に対し輪郭抽出処理を行った画像をいう。
◇計測輪郭画像：計測画像に対し輪郭抽出処理を行った画像をいう。
◇精度検証用の輪郭計測画像：参照画像に基づいて生成した精度検証用の計測画像に対し、輪郭抽出処理を行った画像をいう。
◇結果画像：撮影した参照画像および計測画像、又は輪郭抽出処理を行った参照輪郭画像および計測輪郭画像に対し、Ｓ－ＢＯＳ法を適用して得られた、流れ場の擾乱を可視化した画像であって、主に、流れ場の密度勾配を可視化した画像をいうが、密度勾配を算出する過程で得られた位相差や変位量を可視化した図を含む。換言すると、撮影画像に対して画像処理を行って背景画像の歪み量を測定し、画素ごとに歪み量に応じた明るさを割り当てる処理を行った画像であるともいえる。結果画像は、背景画像として横縞を使用した場合、又は、背景画像として縦縞を使用した場合のいずれか一方の計算結果を利用したものであってもよいし、横縞と縦縞の両方の計算結果を二乗平均等で重畳した結果であってもよい。
【0050】
◇グレースケール画像：白黒の濃淡で表現した画像のことをいう。
◇画素・ピクセル：画素とは画像を構成する要素であって、画像がデジタルデータで構成される場合の最小単位をいう。ピクセルと呼ぶこともある。本明細書では、画素数をカウントするときの単位を示す用語として「ピクセル」を使用するものとする。
◇画素値・輝度値：画素値とは、画素の値のことをいう。カラー画像では、１つの画素の色を、Ｒ（赤）、Ｇ（緑）、Ｂ（青）の３原色で表し、１画素のＲＧＢの各要素をそれぞれ８ビットとした場合、合計２４ビットの画素値で表される。グレースケール画像の場合は、色情報は含まず明るさ情報のみを含んでいるので「輝度値」と呼ぶことにする。１画素を８ビットで表した場合、濃淡を２の８乗＝２５６階調まで表すことができ、画素値０が黒、画素値２５５が白となる。
◇画素値の変化率（又は輝度値の変化率）：空間フィルタリング処理を行う場合において、例えば、線形フィルタのＳＯＢＥＬフィルタなどの微分フィルタを利用して、周辺の画素も含めて畳み込み演算を行った場合、算出される新たな画素値（又は新たな輝度値）は、画素値の変化率（又は輝度値の変化率）を示すことになるので、新たな画素値（又は新たな輝度値）を画素値の変化率（又は輝度値の変化率）と定義することにする。
◇画素・ピクセル：画素とは画像を構成する要素であって、画像がデジタルデータで構成される場合の最小単位をいう。ピクセルと呼ぶこともある。本明細書では、画素数をカウントするときの単位を示す用語として「ピクセル」を使用するものとする。
◇画素値・輝度値：画素値とは、画素の値のことをいう。カラー画像では、１つの画素の色を、Ｒ（赤）、Ｇ（緑）、Ｂ（青）の３原色で表し、１画素のＲＧＢの各要素をそれぞれ８ビットとした場合、合計２４ビットの画素値で表される。グレースケール画像の場合は、色情報は含まず明るさ情報のみを含んでいるので「輝度値」と呼ぶことにする。１画素を８ビットで表した場合、濃淡を２の８乗＝２５６階調まで表すことができ、画素値０が黒、画素値２５５が白となる。
◇画素値の変化率（又は輝度値の変化率）：空間フィルタリング処理を行う場合において、例えば、線形フィルタのＳＯＢＥＬフィルタなどの微分フィルタを利用して、周辺の画素も含めて畳み込み演算を行った場合、算出される新たな画素値（又は新たな輝度値）は、画素値の変化率（又は輝度値の変化率）を示すことになるので、新たな画素値（又は新たな輝度値）を画素値の変化率（又は輝度値の変化率）と定義することにする。
【0051】
◇空間フィルタリング処理：画素値を操作する際に、当該画素の周辺に存在する画素も用いて計算を行うことをいう。大きく分けて、線形および非線形のフィルタリング処理がある。デジタルデータで構成される画像データを空間フィルタリング処理する場合は、積和演算（畳み込み演算）による線形のフィルタリング処理が用いられることが多いが、非線形のフィルタリング処理を用いてもよい。
◇エッジ検出用の空間フィルタリング処理：輪郭抽出処理を行うための空間フィルタリング処理に用いるフィルタとしては、線形フィルタに分類されるＳＯＢＥＬフィルタ、又はプレヴィットフィルタを含む微分フィルタが用いられることが多いがこれに限定されない。例えば、ノイズを除去しながらエッジを検出する非線形フィルタを用いてもよい。
◇エッジ検出用の空間フィルタリング処理：輪郭抽出処理を行うための空間フィルタリング処理に用いるフィルタとしては、線形フィルタに分類されるＳＯＢＥＬフィルタ、又はプレヴィットフィルタを含む微分フィルタが用いられることが多いがこれに限定されない。例えば、ノイズを除去しながらエッジを検出する非線形フィルタを用いてもよい。
【発明を実施するための形態】
【0052】
以下、本発明の実施例について説明する。
なお、説明中の方法や手段はあくまで例示であって、他の方法や手段にも適用できる。
なお、説明中の方法や手段はあくまで例示であって、他の方法や手段にも適用できる。
【0053】
１．本発明の測定システムの構成例
まず、図７、図８を用いて本発明の測定システムの構成の一例について説明する。
図７は、本発明の測定システムの構成の一例を示す正面図である。図７の測定システムは、左側から順に、背景２３０を照らす光源２４０、背景画像を設置した背景２３０、流れ場を観測するための風洞２２０、背景画像と風洞２２０内の流れ場の現象を併せて撮影するための撮影装置２１０で構成されている。
観測対象となる流れ場は、風洞２２０の手前から奥に向かう方向に形成される。
光源２４０は、背景２３０として印刷された背景画像を用いる場合に背景を照らすために用いられるが、背景２３０として液晶ディスプレイを用いる場合には不要となる。
まず、図７、図８を用いて本発明の測定システムの構成の一例について説明する。
図７は、本発明の測定システムの構成の一例を示す正面図である。図７の測定システムは、左側から順に、背景２３０を照らす光源２４０、背景画像を設置した背景２３０、流れ場を観測するための風洞２２０、背景画像と風洞２２０内の流れ場の現象を併せて撮影するための撮影装置２１０で構成されている。
観測対象となる流れ場は、風洞２２０の手前から奥に向かう方向に形成される。
光源２４０は、背景２３０として印刷された背景画像を用いる場合に背景を照らすために用いられるが、背景２３０として液晶ディスプレイを用いる場合には不要となる。
【0054】
ここで、背景画像を設置した背景２３０と風洞２２０の中心（観測対象の流れ場３００の中心）との距離ＺＤは、通常、３０～８０ｍｍ程度が採用されることが多いが、撮影装置２１０の焦点距離ｆ、現象の規模や密度勾配に応じて適宜調整される。例えば、密度勾配が小さい場合は、位相反転が微妙な場合の判定精度を上げる必要があることへの配慮の観点から、見かけ上の変位量を大きくとって、計測精度を上げることが望ましいので、ＺＤを大きくすることが望ましい。
本発明の実験では、特に断りがない場合、ＺＤ＝５０ｍｍで測定したデータをもとに説明する。
なお、背景画像の大きさについても、現象の規模や密度勾配に応じて適宜調整され、Ａ４サイズ程度の小規模なものからＡ０サイズを超える各種のサイズを利用可能であるが、本発明の実験では、特に断りがない場合、Ａ４サイズで測定したデータをもとに説明する。
本発明の実験では、特に断りがない場合、ＺＤ＝５０ｍｍで測定したデータをもとに説明する。
なお、背景画像の大きさについても、現象の規模や密度勾配に応じて適宜調整され、Ａ４サイズ程度の小規模なものからＡ０サイズを超える各種のサイズを利用可能であるが、本発明の実験では、特に断りがない場合、Ａ４サイズで測定したデータをもとに説明する。
【0055】
撮影装置２１０と風洞２２０の中心（観測対象の流れ場３００の中心）との距離ＺＡは、通常、１０００ｍｍ～２０００ｍｍ程度が採用されることが多く、撮影装置２１０の焦点距離ｆ、現象の規模や密度勾配に応じて適宜調整される。本発明の実験では、特に断りがない場合、ＺＡ＝１７３０ｍｍで測定したデータをもとに説明する。
風洞２２０の幅２ΔＺＤ（観測部位３００の幅に相当）も、現象の規模や密度勾配に応じて適宜調整されるが、本発明では、特に断りがない場合、２ΔＺＤ＝４０ｍｍで測定したデータをもとに説明する。
風洞２２０の幅２ΔＺＤ（観測部位３００の幅に相当）も、現象の規模や密度勾配に応じて適宜調整されるが、本発明では、特に断りがない場合、２ΔＺＤ＝４０ｍｍで測定したデータをもとに説明する。
【0056】
図８は、本発明の測定システムを上から見た平面図である。
図８において、大気２５０が図の矢印で示した方向で風洞２２０に導入され、真空タンク２６０で吸引されることで、風洞２２０内で超音速の流れを生成する。風洞２２０内に形成された流れ場は、観測窓２２２を介して、撮影装置２１０で撮影される。
真空タンク２６０の減圧レベルは、測定したい現象の規模や流速に応じて調整されるが、本発明の実験では、特に断りがない場合、減圧のレベルが５ｋPaabsに設定されていることを前提に説明する。
図８では図示を省略するが、風洞２２０内には、風洞センターボディ２２４が設けられ、大気２５０は、センターボディ２２４の上下で、マッハ２．０の気流Ｍ１（符号２５２）と、マッハ１．５の気流Ｍ２（符号２５４）に加速した、流速の異なる２つの流れを形成するように構成されている（図１４参照）。
図８において、大気２５０が図の矢印で示した方向で風洞２２０に導入され、真空タンク２６０で吸引されることで、風洞２２０内で超音速の流れを生成する。風洞２２０内に形成された流れ場は、観測窓２２２を介して、撮影装置２１０で撮影される。
真空タンク２６０の減圧レベルは、測定したい現象の規模や流速に応じて調整されるが、本発明の実験では、特に断りがない場合、減圧のレベルが５ｋPaabsに設定されていることを前提に説明する。
図８では図示を省略するが、風洞２２０内には、風洞センターボディ２２４が設けられ、大気２５０は、センターボディ２２４の上下で、マッハ２．０の気流Ｍ１（符号２５２）と、マッハ１．５の気流Ｍ２（符号２５４）に加速した、流速の異なる２つの流れを形成するように構成されている（図１４参照）。
【0057】
本発明の実験では、背景画像としては、平均波長０．４２ｍｍの縞模様をフォトマット紙に印刷したものを用いた。撮影装置２１０にはデジタルカメラ（焦点距離ｆが３００ｍｍのレンズを使用）を用いてシャッタースピード１／１０００ｓｅｃ、F値３２、ＩＳＯ感度５０００とした。
【0058】
２．本発明の位相差検出型ＢＯＳ法の原理
２－１．流れ場の擾乱の密度勾配の算出過程について
図９は、図７、図８の測定システムの構成例において、風洞２２０の内部の観測部位３００を抜き出して記載するとともに、撮影装置２１０を構成するレンズ２１２及び撮像素子２１４を抜き出して記載した、光学系からみた測定システムの構成例を示す図である。
図９では、観測部位３００の流れ場において、擾乱による密度変化が発生し、背景画像２３０からの光線が、流れ場内で、角度εで屈折する様子が示されている。
屈折した光線は、撮像素子２１４上で、Ｙ軸方向の所定の変位量δyとして観測される。
２－１．流れ場の擾乱の密度勾配の算出過程について
図９は、図７、図８の測定システムの構成例において、風洞２２０の内部の観測部位３００を抜き出して記載するとともに、撮影装置２１０を構成するレンズ２１２及び撮像素子２１４を抜き出して記載した、光学系からみた測定システムの構成例を示す図である。
図９では、観測部位３００の流れ場において、擾乱による密度変化が発生し、背景画像２３０からの光線が、流れ場内で、角度εで屈折する様子が示されている。
屈折した光線は、撮像素子２１４上で、Ｙ軸方向の所定の変位量δyとして観測される。
【0059】
次に、本発明で利用する位相差検出型ＢＯＳ法における新たなＳ－ＢＯＳ法の原理について、数式を用いて説明するとともに、数式が示す概念を図示した図１０～図１２を用いて説明する。
本発明の位相差検出型ＢＯＳ法では、擾乱のない参照画像と、流れ場を撮影した計測画像の２つの画像を用いるが、これを図示したのが図１０である。
図１０によれば、例えば、背景画像（参照画像）として正弦波を用いた場合、参照画像がsinθ、計測画像上の変位を位相差φとして表した撮影画像（計測画像）がsin(θ＋φ)で表現されることが示されている。
本発明の位相差検出型ＢＯＳ法では、擾乱のない参照画像と、流れ場を撮影した計測画像の２つの画像を用いるが、これを図示したのが図１０である。
図１０によれば、例えば、背景画像（参照画像）として正弦波を用いた場合、参照画像がsinθ、計測画像上の変位を位相差φとして表した撮影画像（計測画像）がsin(θ＋φ)で表現されることが示されている。
【0060】
そして、本発明の位相差検出型ＢＯＳ法の特徴としては、出発点として、図１１に示すように、例えば参照画像が正弦波の場合において、参照画像（sinθ）、及び参照画像を微分した微分画像である余弦波（cosθ）、及び計測画像であるsin(θ+φ)の３種の画像を利用することが挙げられる。
すなわち、参照画像（sinθ）、及び参照画像を微分した微分画像である余弦波（cosθ）、及び計測画像であるsin(θ+φ)の３種の画像をもとに、式（１）を立てて、位相差を算出する点に特徴がある。なお、微分計算には中心差分法を用いたが、これに限定されない。
すなわち、参照画像（sinθ）、及び参照画像を微分した微分画像である余弦波（cosθ）、及び計測画像であるsin(θ+φ)の３種の画像をもとに、式（１）を立てて、位相差を算出する点に特徴がある。なお、微分計算には中心差分法を用いたが、これに限定されない。
【数1】
【0061】
ここで、従来のＳ－ＢＯＳ法では、変位量が±１／２πの範囲を超えると位相反転が発生したが、本発明の新たなＳ－ＢＯＳ法では、式（１）を基礎とするので、±πの範囲まで位相反転が生じないことが分かる。
【0062】
式（１）にローパスフィルタをかけることで、右辺よりsinφおよびcosφが得られ逆正接関数より位相差を特定する。この時、位相差は(±π)の範囲で求められる。ローパスフィルタとしては、ＦＩＲフィルタを用いたが、これに限定されない。
式（１）から、位相差を求める過程を概念的に図示したのが、図１２である。
図１２によれば、参照画像（sinθ）、及び参照画像を微分した微分画像である余弦波（cosθ）、及び計測画像であるsin(θ+φ)の３種の画像を出発点として、ローパスフィルタ処理および式変形を経て、位相差を特定することが示されている。
式（１）から、位相差を求める過程を概念的に図示したのが、図１２である。
図１２によれば、参照画像（sinθ）、及び参照画像を微分した微分画像である余弦波（cosθ）、及び計測画像であるsin(θ+φ)の３種の画像を出発点として、ローパスフィルタ処理および式変形を経て、位相差を特定することが示されている。
【0063】
ここで、位相差は、背景画像を構成する所定の空間周波数を持つ波の波長λを用いて、式（２）に示すように、Ｙ軸方向の変位量δyへ変換できる。ここでｍは撮影倍率である。
【数2】
【0064】
背景画像からの光線が流れ場内で起こす屈折の角度εは、十分に小さな値であることを前提とすれば（tanε≒ε）、背景画像と流れ場中心との距離ＺＤ、流れ場の中心とカメラレンズとの距離ＺＡ、撮影装置のレンズの焦点距離ｆを用いて表すと、式（３）のようになる（図９参照）。
また、屈折角εは、流れ場の密度勾配∂ρ/∂r、グラッドストーン定数Ｇ、流れ場の幅ΔＺＤを用いて、式（４）のように表すことができる。ρ０は、観測部位の気体の基準密度である。
ここで、背景画像として横縞を使用する場合には、Ｙ軸方向の変位を検出しようとする場面なので、流れ場の密度勾配∂ρ/∂rは、∂ρ/∂ｙとして式を定立することとする。なお、背景画像として縦縞を使用する場合には、Ｘ軸方向の変位を検出しようとする場面なので、流れ場の密度勾配∂ρ/∂rは、∂ρ/∂ｘとして式を定立することとすればよく、後の説明は同様であるので、ここでは説明も省略する。
【数3】
ここで、背景画像として横縞を使用する場合には、Ｙ軸方向の変位を検出しようとする場面なので、流れ場の密度勾配∂ρ/∂rは、∂ρ/∂ｙとして式を定立することとする。なお、背景画像として縦縞を使用する場合には、Ｘ軸方向の変位を検出しようとする場面なので、流れ場の密度勾配∂ρ/∂rは、∂ρ/∂ｘとして式を定立することとすればよく、後の説明は同様であるので、ここでは説明も省略する。
【数4】
【0065】
ここで、観測対象が二次元流れであれば、式（４）の積分は、Ｚ軸方向に対し一定であることを利用して、式（３）（４）より、屈折角を消去し、式（４）を整理すると、最終的に流れ場の密度勾配を求めることができる。
以上の式（２）～（４）の導出過程を図示したのが図１３である。
図１３によれば、位相差を算出後、変位量を求め、測定システムの幾何関係から（図９参照）、屈折角を求め、屈折率勾配にＧｒａｄｓｔｏｎｅ－Ｄａｌｅ則を適用して、密度勾配を算出できることが分かる。
本発明では、以上の原理に基づくＳ－ＢＯＳ法を新たなＳ－ＢＯＳ法と呼ぶこととする。
なお、本発明では、位相反転をできるだけ発生させないようにするために、主に、この新たなＳ－ＢＯＳ法を用いて流れ場の擾乱の密度勾配を算出することとするが、従来のＳ－ＢＯＳ法を用いても、検出可能範囲が狭くなるというデメリットはあるものの、流れ場の擾乱の密度勾配を算出すること自体は可能であり、従来のＳ－ＢＯＳ法を適用しても構わない。
以上の式（２）～（４）の導出過程を図示したのが図１３である。
図１３によれば、位相差を算出後、変位量を求め、測定システムの幾何関係から（図９参照）、屈折角を求め、屈折率勾配にＧｒａｄｓｔｏｎｅ－Ｄａｌｅ則を適用して、密度勾配を算出できることが分かる。
本発明では、以上の原理に基づくＳ－ＢＯＳ法を新たなＳ－ＢＯＳ法と呼ぶこととする。
なお、本発明では、位相反転をできるだけ発生させないようにするために、主に、この新たなＳ－ＢＯＳ法を用いて流れ場の擾乱の密度勾配を算出することとするが、従来のＳ－ＢＯＳ法を用いても、検出可能範囲が狭くなるというデメリットはあるものの、流れ場の擾乱の密度勾配を算出すること自体は可能であり、従来のＳ－ＢＯＳ法を適用しても構わない。
【0066】
２－２．計算過程で得られた数値に基づく結果画像の生成について
以上の計算原理を利用して、画素ごとに密度勾配を算出し、算出した画素ごとの密度勾配の数値（アナログ）に対し階調を表現するビット数に応じてデジタル値に換算し、換算した値に応じて輝度を設定するとともに、各画素を正負で色分けを行って（正の場合は赤色、負の場合は青色等）、輝度や色分けに応じてプロットすると、図１４に示すような結果画像が得られる。
図１４によれば、急峻な密度勾配の変化がある衝撃波３１０の領域と、やや穏やかな密度勾配の変化がある混合層３２０の領域が再現されていることが分かる。
実際の結果画像では、密度勾配が正の領域を赤、密度勾配が負の領域を青で示すように作成されるが、本発明の図はグレースケール画像のため、密度勾配が正の領域を点線の領域で示し、点線の領域以外は密度勾配が負の領域であることを示している。
なお、結果画像としては、密度勾配を可視化したものだけでなく、密度勾配の算出過程で得られた位相差、変位量を可視化したものであっても構わない（図１９参照）。
以上の計算原理を利用して、画素ごとに密度勾配を算出し、算出した画素ごとの密度勾配の数値（アナログ）に対し階調を表現するビット数に応じてデジタル値に換算し、換算した値に応じて輝度を設定するとともに、各画素を正負で色分けを行って（正の場合は赤色、負の場合は青色等）、輝度や色分けに応じてプロットすると、図１４に示すような結果画像が得られる。
図１４によれば、急峻な密度勾配の変化がある衝撃波３１０の領域と、やや穏やかな密度勾配の変化がある混合層３２０の領域が再現されていることが分かる。
実際の結果画像では、密度勾配が正の領域を赤、密度勾配が負の領域を青で示すように作成されるが、本発明の図はグレースケール画像のため、密度勾配が正の領域を点線の領域で示し、点線の領域以外は密度勾配が負の領域であることを示している。
なお、結果画像としては、密度勾配を可視化したものだけでなく、密度勾配の算出過程で得られた位相差、変位量を可視化したものであっても構わない（図１９参照）。
【0067】
３．輪郭抽出処理および輪郭画像の生成
３－１．輪郭抽出処理
次に、図５～図６、および図１５～図１８を用いて、本発明における輪郭抽出処理および輪郭画像の生成について説明する。
図６に示すように、従来技術では、元画像をそのまま利用してＳ－ＢＯＳ法を適用して、縞の湾曲（背景画像に対して縞がずれた変位）を算出していたため、太さの変化には着目せず、流れ場の密度勾配を可視化するのに重要な特徴量が見過ごされてきた。
３－１．輪郭抽出処理
次に、図５～図６、および図１５～図１８を用いて、本発明における輪郭抽出処理および輪郭画像の生成について説明する。
図６に示すように、従来技術では、元画像をそのまま利用してＳ－ＢＯＳ法を適用して、縞の湾曲（背景画像に対して縞がずれた変位）を算出していたため、太さの変化には着目せず、流れ場の密度勾配を可視化するのに重要な特徴量が見過ごされてきた。
【0068】
そこで、本発明では、縞の湾曲（背景画像に対して縞がずれた変位）だけでなく、揺らぎよって縞の太さが変化すること（縞の太さの変化）にも着目し、流れ場の密度勾配を可視化するのに重要な特徴量をもれなく考慮する手法について研究を進めた。
その中で、縞の境目である輪郭付近の情報を抽出するとよい結果画像が得られることが判明し、輪郭付近の輝度値の変化傾向をうまく捉える輪郭抽出処理を利用することとした。
その中で、縞の境目である輪郭付近の情報を抽出するとよい結果画像が得られることが判明し、輪郭付近の輝度値の変化傾向をうまく捉える輪郭抽出処理を利用することとした。
【0069】
輪郭画像における輪郭線とは、その線上で輝度値が変化することを意味するが、そのような部分を特定するには微分演算を行えばよい。
微分演算としては、各種の線形および非線形の空間フィルタリング処理を用いることができるが、説明の便宜のため、以下、一例として、線形フィルタの微分フィルタに分類されるＳＯＢＥＬフィルタを用いて、積和演算（畳み込み演算）による線形のフィルタリング処理を行った場合について説明する。
微分演算としては、各種の線形および非線形の空間フィルタリング処理を用いることができるが、説明の便宜のため、以下、一例として、線形フィルタの微分フィルタに分類されるＳＯＢＥＬフィルタを用いて、積和演算（畳み込み演算）による線形のフィルタリング処理を行った場合について説明する。
【0070】
フィルタのフィルタ値（パラメータ）は、図３５に示すように、フィルタの種類ごとにフィルタテーブルに格納して、適宜、メモリ上にロードして利用することができる。
図３２によれば、例えば、フィルタサイズ３（３×３）のＳＯＢＥＬフィルタを選択した場合は（ＩＤ２０の列を参照）、横方向の輪郭検出用、縦方向の輪郭検出用のフィルタ値（パラメータ）を利用することができることが分かる。
図３２によれば、例えば、フィルタサイズ３（３×３）のＳＯＢＥＬフィルタを選択した場合は（ＩＤ２０の列を参照）、横方向の輪郭検出用、縦方向の輪郭検出用のフィルタ値（パラメータ）を利用することができることが分かる。
【0071】
以下、輪郭抽出処理の理論的な側面について、数式をもとに説明する。
まず、図５（ａ）（ｂ）に示したように、輪郭抽出処理の元になる参照画像および計測画像などの元画像（グレースケール画像）を用意する。
ここで、元画像の任意の画素をＩ(i,j)とし、その１つ上をＩ(i,j+1)、１つ右をＩ(i+1,j)と表現する（Ｘ方向、Ｙ方向の座標は図１５を参照）。
この場合、縦方向(ΔＹ)、横方向(ΔＸ)それぞれの輝度値の変化率は式（５）で表される。
まず、図５（ａ）（ｂ）に示したように、輪郭抽出処理の元になる参照画像および計測画像などの元画像（グレースケール画像）を用意する。
ここで、元画像の任意の画素をＩ(i,j)とし、その１つ上をＩ(i,j+1)、１つ右をＩ(i+1,j)と表現する（Ｘ方向、Ｙ方向の座標は図１５を参照）。
この場合、縦方向(ΔＹ)、横方向(ΔＸ)それぞれの輝度値の変化率は式（５）で表される。
【数5】
【0072】
ＳＯＢＥＬフィルタを選択した場合、隣り合う画素のほかに斜めに隣接する画素の情報も用いて畳み込み演算が行われる。
具体的には、元画像の任意の画素Ｉ(i,j)に対し、縦横斜めの周辺の画素も利用して、ＳＯＢＥＬフィルタを用いて、積和演算による畳み込み演算を行うことで、輪郭画像としての画素Ｉ(i,j)における縦方向(ΔＹ)、横方向(ΔＸ)それぞれの輝度値の変化率を得ることができる。
この輪郭抽出処理としての畳み込み演算は、式（６）で表される。
ここで、縦方向(ΔＹ)の数式は、背景画像として横縞を使用した場合の縦方向の変位を輪郭抽出するためのフィルタ処理を示すものであり、横方向(ΔＸ)の数式は、背景画像として縦縞を使用した場合の横方向の変位を輪郭抽出するためのフィルタ処理を示すものである。
具体的には、元画像の任意の画素Ｉ(i,j)に対し、縦横斜めの周辺の画素も利用して、ＳＯＢＥＬフィルタを用いて、積和演算による畳み込み演算を行うことで、輪郭画像としての画素Ｉ(i,j)における縦方向(ΔＹ)、横方向(ΔＸ)それぞれの輝度値の変化率を得ることができる。
この輪郭抽出処理としての畳み込み演算は、式（６）で表される。
ここで、縦方向(ΔＹ)の数式は、背景画像として横縞を使用した場合の縦方向の変位を輪郭抽出するためのフィルタ処理を示すものであり、横方向(ΔＸ)の数式は、背景画像として縦縞を使用した場合の横方向の変位を輪郭抽出するためのフィルタ処理を示すものである。
【数6】
【0073】
そして、こうして求めた縦横それぞれの輝度値の変化率に対して二乗平均を求める形で、輪郭画像上の各画素の輝度値Ｉ(i,j)contoursを定義すると、式（７）を得ることができる。
【数7】
【0074】
これを模式的に示したのが図１５である。
図１５によれば、ある画素について、ある画素を含めた周辺の３×３画素の範囲について、３×３のサイズの横方向と縦方向の輪郭抽出用のフィルタを用いて、積和演算による畳み込み演算を行うことが示されている。
より具体的には、画素Ｉ(i,j)における縦方向(ΔＹ)、横方向(ΔＸ)それぞれの輝度値の変化率を求め、そのうえで縦方向及び横方向の輝度値の変化率の二乗平均を求める形で、輪郭画像上の各画素の輝度値Ｉ(i,j)contoursを得ることが示されている。
以上のように、式（７）によって定義した輪郭画像の輝度値に基づいて算出した密度勾配をもとに生成するという基本原理を説明したが、実際の計算では、背景画像として横縞を使用した場合の輪郭画像（横縞）と、背景画像として縦縞を使用した場合の輪郭画像（縦縞）を生成し、それぞれの輪郭画像に対してＳ－ＢＯＳ法を適用して得られた密度勾配を二乗平均などを用いて重畳する形で、流れ場の擾乱を可視化した結果画像を生成することとする。
このため、以下の説明では、説明の便宜のため、主に、背景画像として横縞を使用した場合について実験結果や効果を説明するが、背景画像として縦縞を用いた場合も同様である。
図１５によれば、ある画素について、ある画素を含めた周辺の３×３画素の範囲について、３×３のサイズの横方向と縦方向の輪郭抽出用のフィルタを用いて、積和演算による畳み込み演算を行うことが示されている。
より具体的には、画素Ｉ(i,j)における縦方向(ΔＹ)、横方向(ΔＸ)それぞれの輝度値の変化率を求め、そのうえで縦方向及び横方向の輝度値の変化率の二乗平均を求める形で、輪郭画像上の各画素の輝度値Ｉ(i,j)contoursを得ることが示されている。
以上のように、式（７）によって定義した輪郭画像の輝度値に基づいて算出した密度勾配をもとに生成するという基本原理を説明したが、実際の計算では、背景画像として横縞を使用した場合の輪郭画像（横縞）と、背景画像として縦縞を使用した場合の輪郭画像（縦縞）を生成し、それぞれの輪郭画像に対してＳ－ＢＯＳ法を適用して得られた密度勾配を二乗平均などを用いて重畳する形で、流れ場の擾乱を可視化した結果画像を生成することとする。
このため、以下の説明では、説明の便宜のため、主に、背景画像として横縞を使用した場合について実験結果や効果を説明するが、背景画像として縦縞を用いた場合も同様である。
【0075】
なお、積和演算による畳み込み演算を行う場合、１列目の最初の画素Ｉ(１,１)など、周辺画素としてフィルタサイズと同じ３×３の周辺画素が得られない場合に、所定のデータを使って３×３の周辺画素を生成するパディング処理を行うのが通例であるが、観測対象となる現象がない個所からフィルタ処理を開始する必要はなく、２列目ないし、それ以降の列からフィルタ処理を行えば足りる。
【0076】
なお、ＳＯＢＥＬフィルタなどのように、周辺の画素を含めて畳み込み処理を行うことは、結果が平均化されノイズに強く輪郭線の途切れが生じにくいというメリットがある。
以上の処理では、ＳＯＢＥＬフィルタを用いたが、プレヴィットフィルタなどの他の微分フィルタを用いることもできるし、非線形フィルタを用いた処理を行うこともできる。
以上の処理では、ＳＯＢＥＬフィルタを用いたが、プレヴィットフィルタなどの他の微分フィルタを用いることもできるし、非線形フィルタを用いた処理を行うこともできる。
【0077】
３－２．輪郭画像の生成
以上の輪郭抽出処理を行う前の元画像の一部と、輪郭抽出処理を行った後の輪郭画像の一部を抜き出した様子を図５（ｃ）（ｄ）及び図１６に示す。
図１６によれば、本発明の測定システムのパラメータを用いた場合において、背景画像を適宜設定して、平均線幅０．２ｍｍ（８．５ピクセル）の縞画像（元画像）に対して輪郭抽出処理を行った場合、輪郭画像として、平均線幅０．１ｍｍ（４．３ピクセル）の縞画像が得られること分かる。なお、図では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
以上の輪郭抽出処理を行う前の元画像の一部と、輪郭抽出処理を行った後の輪郭画像の一部を抜き出した様子を図５（ｃ）（ｄ）及び図１６に示す。
図１６によれば、本発明の測定システムのパラメータを用いた場合において、背景画像を適宜設定して、平均線幅０．２ｍｍ（８．５ピクセル）の縞画像（元画像）に対して輪郭抽出処理を行った場合、輪郭画像として、平均線幅０．１ｍｍ（４．３ピクセル）の縞画像が得られること分かる。なお、図では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
【0078】
また、輪郭画像を構成する輝度値の分布（所定の空間周波数を持つ波）の波長についても、元画像の波長の約１／２になっているので以下説明する。
図１７（ａ）（ｂ）は、本発明の輪郭抽出処理を行い作成した輪郭画像（参照輪郭画像）とその元画像（参照画像）の一部、および両者の縞画像を構成する輝度値の分布（所定の空間周波数を持つ波）を比較した図である。
図１７（ａ）（ｂ）によれば、元画像（参照画像）の平均波長は０．４２ｍｍ(１７ピクセル)であったのに対し輪郭画像では平均波長０．２１ｍｍ(８．５ピクセル)と元画像に比べ半減していることが分かる。
図１７（ａ）（ｂ）は、本発明の輪郭抽出処理を行い作成した輪郭画像（参照輪郭画像）とその元画像（参照画像）の一部、および両者の縞画像を構成する輝度値の分布（所定の空間周波数を持つ波）を比較した図である。
図１７（ａ）（ｂ）によれば、元画像（参照画像）の平均波長は０．４２ｍｍ(１７ピクセル)であったのに対し輪郭画像では平均波長０．２１ｍｍ(８．５ピクセル)と元画像に比べ半減していることが分かる。
【0079】
すなわち、元画像（参照画像）の縞画像と輪郭画像（参照輪郭画像）の縞画像とを比較すると、両者の平均間隔が約１／２、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約１／２、周期も約１／２になっているので、周波数に換算した場合約２倍の周波数に相当することになり、空間分解能も約２倍になっていることが分かる。
【0080】
同様に、図１８によれば、測定画像に対して輪郭抽出処理を行った測定輪郭画像についても、元画像（測定画像）の縞画像と輪郭画像（計測輪郭画像）の縞画像とを比較すると、両者の平均間隔が約１／２、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約１／２、周期も約１／２になっているので、周波数に換算した場合約２倍の周波数に相当することになり、空間分解能も約２倍になっており、流れ場の擾乱の再現性を高めることができることが分かる。
なお、図１７、図１８では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
なお、図１７、図１８では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
【0081】
４．輪郭画像を利用した本発明の効果について
前述の「２．本発明の位相差検出型ＢＯＳ法の原理」において、測定した画像（参照画像、計測画像）に対し、Ｓ－ＢＯＳ法を適用して密度勾配を算出し、算出した密度勾配、もしくは、密度勾配の算出過程で得た位相差や変位量をもとに、流れ場の擾乱を可視化した結果画像を得ることができることはすでに説明した。
ここでは、輪郭抽出処理を行った輪郭画像のペア（参照輪郭画像、計測輪郭画像）を用いて、Ｓ－ＢＯＳ法を適用することで、より空間分解能の高い結果画像が得られることを説明する。
前述の「２．本発明の位相差検出型ＢＯＳ法の原理」において、測定した画像（参照画像、計測画像）に対し、Ｓ－ＢＯＳ法を適用して密度勾配を算出し、算出した密度勾配、もしくは、密度勾配の算出過程で得た位相差や変位量をもとに、流れ場の擾乱を可視化した結果画像を得ることができることはすでに説明した。
ここでは、輪郭抽出処理を行った輪郭画像のペア（参照輪郭画像、計測輪郭画像）を用いて、Ｓ－ＢＯＳ法を適用することで、より空間分解能の高い結果画像が得られることを説明する。
【0082】
図１９（ａ）（ｂ）は、衝撃波と混合層を含む超音速流れ場において、得られた画像に対し、Ｓ－ＢＯＳ法を適用して密度勾配を算出する過程で得た位相差をもとに結果画像を出力したものである。
図１９（ａ）は、元画像（参照画像および計測画像のペア）に対してＳ－ＢＯＳ法を適用して密度勾配を算出する過程で算出した位相差をプロットした画像である。図１９（ｂ）は、輪郭抽出処理を行った輪郭画像（参照輪郭画像および計測輪郭画像のペア）に対してＳ－ＢＯＳ法を適用して密度勾配を算出する過程で算出した位相差をプロットした画像である。
図１９（ａ）は、元画像（参照画像および計測画像のペア）に対してＳ－ＢＯＳ法を適用して密度勾配を算出する過程で算出した位相差をプロットした画像である。図１９（ｂ）は、輪郭抽出処理を行った輪郭画像（参照輪郭画像および計測輪郭画像のペア）に対してＳ－ＢＯＳ法を適用して密度勾配を算出する過程で算出した位相差をプロットした画像である。
【0083】
図１９（ａ）の結果画像では、急峻な密度勾配の変化がある衝撃波の領域（点線で囲った領域）において、現象の再現性がやや乏しく、特に、Ｘ軸の値が増加する領域（点線の右端に近い領域）において、現象が不明瞭になっていることが分かる。
これに対し、図１９（ｂ）の結果画像では、急峻な密度勾配の変化がある衝撃波の領域（点線で囲った領域）において、現象がよく再現され、Ｘ軸の値が増加する領域（点線の右端に近い領域）においても、現象が明瞭に再現できていることが分かる。
これに対し、図１９（ｂ）の結果画像では、急峻な密度勾配の変化がある衝撃波の領域（点線で囲った領域）において、現象がよく再現され、Ｘ軸の値が増加する領域（点線の右端に近い領域）においても、現象が明瞭に再現できていることが分かる。
【0084】
以上の説明では、視覚的な性能差を説明したが、グレースケールの図で表現せざるを得ないこともあって、見え方は人によって異なり、あくまで主観的な比較にならざるを得ないので、密度勾配の変化が急峻な衝撃波を含む部位について、両者の性能の違いを、密度勾配の再現率を示すグラフ（図２０）で定量的に示すこととする。
図２０（ａ）は、衝撃波と混合層を含む超音速流れ場において、Ｓ－ＢＯＳ法を適用して可視化した様子を示す図であって、密度勾配の変化が急峻な衝撃波の存在箇所を大まかに示すとともに、衝撃波を含む断面ＡＢの位置を示す図である。ＡＢのラインは、Ｙ軸に平行にＸ＝１０ｍｍの線上の密度勾配値を抜き出して示すためのラインである。
図２０（ａ）は、衝撃波と混合層を含む超音速流れ場において、Ｓ－ＢＯＳ法を適用して可視化した様子を示す図であって、密度勾配の変化が急峻な衝撃波の存在箇所を大まかに示すとともに、衝撃波を含む断面ＡＢの位置を示す図である。ＡＢのラインは、Ｙ軸に平行にＸ＝１０ｍｍの線上の密度勾配値を抜き出して示すためのラインである。
【0085】
図２０（ｂ）（ｃ）は、従来技術との性能差をよく見極められるように、密度勾配の再現率の差が出やすい密度勾配の変化が急峻な衝撃波を含む断面ＡＢについて、元画像に基づいてＳ－ＢＯＳ法を適用した場合の密度勾配と、輪郭画像に基づいてＳ－ＢＯＳ法を適用した場合の密度勾配の再現率の違いを、グラフで定量的に示した図である。
なお、２０（ｂ）（ｃ）では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
なお、２０（ｂ）（ｃ）では、背景画像として横縞を使用した場合の例を示したが、背景画像として縦縞を使用した場合も同様である。
【0086】
図２０（ｂ）（ｃ）において、Y＝０～５mmの間に見られるやや小さな密度勾配は混合層によるものであり、Y＝－３ｍｍと８ｍｍ付近に見られる大きな密度勾配は衝撃波によるものである。
ここで、図２０（ｂ）に示した矢印は、密度勾配のピーク（密度変化が急峻な個所）を示したものであり、点線の矢印が数値計算による密度勾配の値（理論値）を示し、実線の矢印が観測した元画像（参照画像と計測画像のペア）に基づいてＳ－ＢＯＳ法を適用して算出した密度勾配の値を示している。
同様に、図２０（ｃ）に示した矢印は、密度勾配のピーク（密度変化が急峻な個所）を示したものであり、点線の矢印が数値計算による密度勾配の値（理論値）を示し、実線の矢印が輪郭画像（参照輪郭画像と計測輪郭画像のペア）に基づいてＳ－ＢＯＳ法を適用して算出した密度勾配の値を示している。
数値計算による理論値は、図１４に示すような、風洞内の物理的な構成を数理モデルとして再現し、流体力学に基づくシミュレーションを行って算出した理論上の数値をいう。
ここで、図２０（ｂ）に示した矢印は、密度勾配のピーク（密度変化が急峻な個所）を示したものであり、点線の矢印が数値計算による密度勾配の値（理論値）を示し、実線の矢印が観測した元画像（参照画像と計測画像のペア）に基づいてＳ－ＢＯＳ法を適用して算出した密度勾配の値を示している。
同様に、図２０（ｃ）に示した矢印は、密度勾配のピーク（密度変化が急峻な個所）を示したものであり、点線の矢印が数値計算による密度勾配の値（理論値）を示し、実線の矢印が輪郭画像（参照輪郭画像と計測輪郭画像のペア）に基づいてＳ－ＢＯＳ法を適用して算出した密度勾配の値を示している。
数値計算による理論値は、図１４に示すような、風洞内の物理的な構成を数理モデルとして再現し、流体力学に基づくシミュレーションを行って算出した理論上の数値をいう。
【0087】
図２０（ｂ）によれば、元画像（参照画像と計測画像のペア）を利用した場合は、理論値に対して、密度変化が急峻な個所で約２８～５０％下回る程度の再現率にとどまっていることが分かる（換言するとピークの５０～７２％程度のやや低い再現率）。これは、元画像（参照画像と計測画像のペア）を利用した場合には、空間分解能の不足によって鈍ってしまい、急激な密度変化を伴う衝撃波面において、密度勾配の変化を過小評価してしまっていることを示すものである。
【0088】
これに対し、図２０（ｃ）によれば、輪郭画像（参照輪郭画像と計測輪郭画像のペア）を利用した場合には、密度変化が急峻な個所でも、約６～３５％下回るレベルまで再現率が改善されていることが分かる（換言するとピークの６５～９４％の高い再現率）。このことは、輪郭画像（参照輪郭画像と計測輪郭画像のペア）を利用した場合には空間分解能が向上し、急激な密度変化を伴う衝撃波面においても、現象をよく再現できていることを示すものである。
【0089】
なお、５０～７２％の再現率と６５～９４％の再現率との比較では、輪郭画像（参照輪郭画像と計測輪郭画像のペア）を利用した場合の方が、従来の元画像（参照画像と計測画像のペア）を利用した場合と比較して、約１．３倍性能が向上していることを示しているとみることができる。
【0090】
５．新たな精度検証方法、及びこれを用いた可視化手法の評価
前述のように、実際に撮影した計測画像に基づいて可視化手法の精度を確認することもできるが、精度を確認するための急激な密度変化を伴う衝撃波面の選択した部位によって、現象の再現率にばらつきが出て、可視化手法の精度をうまく評価できない可能性もある。
そこで、参照画像から画像処理で生成した計測画像を用いて、ばらつきが出ない精度検証の方法を開発した。
以下、本発明の新たな精度検証の手法及び、新たな精度検証方法を用いて本発明の可視化手法を用いた場合の評価について、図２１～図２３を用いて説明する。
前述のように、実際に撮影した計測画像に基づいて可視化手法の精度を確認することもできるが、精度を確認するための急激な密度変化を伴う衝撃波面の選択した部位によって、現象の再現率にばらつきが出て、可視化手法の精度をうまく評価できない可能性もある。
そこで、参照画像から画像処理で生成した計測画像を用いて、ばらつきが出ない精度検証の方法を開発した。
以下、本発明の新たな精度検証の手法及び、新たな精度検証方法を用いて本発明の可視化手法を用いた場合の評価について、図２１～図２３を用いて説明する。
【0091】
５－１．新たな精度検証方法の原理について
まず、図２１に示すように、参照画像上の幅Ａピクセル分の任意の領域を考える。この領域をＢピクセル分だけ平行移動した画像を生成し、これを計測画像とする。参照画像と生成した計測画像に対しＳ－ＢＯＳ法を適応することで変位量が求められる。そして、この変位量と、最初に与えた変位量Ｂピクセル分とを比較し計測精度を検証することにした。
ここで、１ピクセルの物理的な幅Δｈ（ｍｍ／ピクセル）を用いればＡとＢはそれぞれｍｍ単位へ変換でき、Ａが小さいほど、その小さな変位を検出するには、より高い空間分解能が必要といえる。
本発明の精度検証では、参照画像に平均波長０．４２ｍｍ（１７ピクセル）のものを用いて前述の方法で計測画像を作成した。与える変位量はＢ＝４ピクセル（約０．１ｍｍ）とした。
まず、図２１に示すように、参照画像上の幅Ａピクセル分の任意の領域を考える。この領域をＢピクセル分だけ平行移動した画像を生成し、これを計測画像とする。参照画像と生成した計測画像に対しＳ－ＢＯＳ法を適応することで変位量が求められる。そして、この変位量と、最初に与えた変位量Ｂピクセル分とを比較し計測精度を検証することにした。
ここで、１ピクセルの物理的な幅Δｈ（ｍｍ／ピクセル）を用いればＡとＢはそれぞれｍｍ単位へ変換でき、Ａが小さいほど、その小さな変位を検出するには、より高い空間分解能が必要といえる。
本発明の精度検証では、参照画像に平均波長０．４２ｍｍ（１７ピクセル）のものを用いて前述の方法で計測画像を作成した。与える変位量はＢ＝４ピクセル（約０．１ｍｍ）とした。
【0092】
この手法では実際の現象を用いないため光学系の改良に対する精度検証法として用いることはできないが、変位量Ｂを可変することで、真値がはっきりとした様々な規模の変位を再現できるため、今回提案した輪郭画像の利用などＢＯＳ法の理論を改良した際の可視化手法の精度検証方法として有用である。
また、今回与えた変位量は、ある点で不連続的に変化するという不連続面を有し、ある意味、急激な密度変化を伴う衝撃波面もしくはそれ以上に急峻な現象に類似する性質を有するので、流れ場の可視化手法の精度検証方法として有用である。
また、今回与えた変位量は、ある点で不連続的に変化するという不連続面を有し、ある意味、急激な密度変化を伴う衝撃波面もしくはそれ以上に急峻な現象に類似する性質を有するので、流れ場の可視化手法の精度検証方法として有用である。
【0093】
５－２．新たな精度検証方法を用いて本発明の可視化手法を用いた場合の評価
実験では、参照画像と、生成した計測画像を元画像とし、これらに対して輪郭抽出処理を行った輪郭画像（参照輪郭画像および計測輪郭画像）を生成し、元画像と輪郭画像それぞれにＳ－ＢＯＳ法を適用して結果の比較を行った。
図２２にＡ＝７５ピクセル（約１．９mm）、図２３にＡ＝２５ピクセル（約０．６ｍｍ）の矩形領域を設定した場合の結果を示す。
図２２、図２３において、(ａ)は幅Ａの変位量を示し、（ｂ）は（ａ）の▼印と▲印のライン上の変位量を幅Ａの近傍について抜き出した図である。図中の実線が輪郭画像を利用した場合、破線が元画像を利用した場合の結果を、点線は真値の分布をそれぞれ示す。
実験では、参照画像と、生成した計測画像を元画像とし、これらに対して輪郭抽出処理を行った輪郭画像（参照輪郭画像および計測輪郭画像）を生成し、元画像と輪郭画像それぞれにＳ－ＢＯＳ法を適用して結果の比較を行った。
図２２にＡ＝７５ピクセル（約１．９mm）、図２３にＡ＝２５ピクセル（約０．６ｍｍ）の矩形領域を設定した場合の結果を示す。
図２２、図２３において、(ａ)は幅Ａの変位量を示し、（ｂ）は（ａ）の▼印と▲印のライン上の変位量を幅Ａの近傍について抜き出した図である。図中の実線が輪郭画像を利用した場合、破線が元画像を利用した場合の結果を、点線は真値の分布をそれぞれ示す。
【0094】
図２２（ｂ）と図２３（ｂ）をみると、こうした不連続面で、輪郭画像を利用した場合の結果（Contourのライン）は真値（Trueのライン）とよく一致しているのに対し、元画像を利用した場合の結果（Originalのライン）はなだらかに分布し不連続的な変化を追えていないことがわかる。
より詳細には、図２２のように、Ａ＝７５ピクセル（約１．９mm）と背景の波長（１７ピクセル）と比較して十分に広い範囲で変位が生じている場合は、元画像を利用した場合でも、ある程度、真値に近い結果が得られているが（過小評価の程度は小）、図２３のように、Ａ＝２５ピクセル（約０．６ｍｍ）のように狭い範囲で変位が生じた場合は、真値に対し８０％程度の数値にしか到達せず、過小評価の程度は大となっていることが分かる。
より詳細には、図２２のように、Ａ＝７５ピクセル（約１．９mm）と背景の波長（１７ピクセル）と比較して十分に広い範囲で変位が生じている場合は、元画像を利用した場合でも、ある程度、真値に近い結果が得られているが（過小評価の程度は小）、図２３のように、Ａ＝２５ピクセル（約０．６ｍｍ）のように狭い範囲で変位が生じた場合は、真値に対し８０％程度の数値にしか到達せず、過小評価の程度は大となっていることが分かる。
【0095】
このように、新たな精度検証方法によっても、輪郭画像を用いるＢＯＳ法は従来のＢＯＳ法と比べてより高い空間分解能が期待できる手法といえる。
このように、本発明では、輪郭抽出処理を行って参照輪郭画像と計測輪郭画像のペアを用いてＳ－ＢＯＳ法を適用して、流れ場の擾乱の密度勾配を算出する点に最も特徴があり、適用するＳ－ＢＯＳ法は、位相反転の検出可能範囲の広狭を除けば、必ずしも、新たなＳ－ＢＯＳ法を適用する必要はなく、従来のＳ－ＢＯＳ法を適用しても構わない。
このように、本発明では、輪郭抽出処理を行って参照輪郭画像と計測輪郭画像のペアを用いてＳ－ＢＯＳ法を適用して、流れ場の擾乱の密度勾配を算出する点に最も特徴があり、適用するＳ－ＢＯＳ法は、位相反転の検出可能範囲の広狭を除けば、必ずしも、新たなＳ－ＢＯＳ法を適用する必要はなく、従来のＳ－ＢＯＳ法を適用しても構わない。
【0096】
６．密度勾配可視化システムの機能ブロック図、システム構成例
本発明の密度勾配可視化システムの機能ブロック図、システム構成の一例について、図２４～図２６を用いて説明する。
６－１．機能ブロック図
図２４は、本発明の密度勾配可視化システムにおいて画像処理などの情報処理をおこなう情報処理システム１００の機能ブロック図である。
本発明の密度勾配可視化システムの機能ブロック図、システム構成の一例について、図２４～図２６を用いて説明する。
６－１．機能ブロック図
図２４は、本発明の密度勾配可視化システムにおいて画像処理などの情報処理をおこなう情報処理システム１００の機能ブロック図である。
【0097】
入力部１１０は、外部からデータ（画像、パラメータ）を入力するための外部インターフェース、測定システム２００の各種パラメータを入力するパラメータ入力手段、キーボードで数値や文字列を入力したり画面表示されたメニューやボタンの選択を行ったりするための入力手段（図示せず）、を備えている。
【0098】
制御部１２０は、入力部１１０と連携して、撮影装置２１０によって撮影した画像を取り込む参照画像取得手段及び計測画像取得手段を備えている。
また、制御部１２０は、参照画像や計測画像に対し輪郭抽出処理を行う輪郭抽出手段、元画像や輪郭画像に基づいて密度勾配を算出する密度勾配算出手段を備えている。
また、制御部１２０は、算出した位相差や変位量ないし密度勾配に基づいて流れ場を可視化した結果画像を生成する結果画像生成手段を備えている。
また、制御部１２０は、位相反転を検出して補正する位相反転検出手段、可視化手法の精度を検証する精度検証手段などを備えている。
また、制御部１２０は、参照画像や計測画像に対し輪郭抽出処理を行う輪郭抽出手段、元画像や輪郭画像に基づいて密度勾配を算出する密度勾配算出手段を備えている。
また、制御部１２０は、算出した位相差や変位量ないし密度勾配に基づいて流れ場を可視化した結果画像を生成する結果画像生成手段を備えている。
また、制御部１２０は、位相反転を検出して補正する位相反転検出手段、可視化手法の精度を検証する精度検証手段などを備えている。
【0099】
制御部１２０の各手段は、記憶部１３０に記憶した各種の処理プログラムを、ＲＡＭ（ランダムアクセスメモリ）などの一時記憶手段にロードして、中央処理装置（ＣＰＵ）がプログラムの各ステップを実行することにより、実現している。
撮影した画像の取り込みに際しては、ネットワーク５００経由で通信部１５０の通信手段を介して取り込んでもよいし、入力部１１０の外部インターフェースを介してメモリカードなどの記憶媒体から取り込むようにしてもよい。
撮影した画像の取り込みに際しては、ネットワーク５００経由で通信部１５０の通信手段を介して取り込んでもよいし、入力部１１０の外部インターフェースを介してメモリカードなどの記憶媒体から取り込むようにしてもよい。
【0100】
記憶部１３０は、長期的な記憶を行うＨＤＤ（ハードディスク）や不揮発性メモリで構成されたＳＳＤ（ソリッドステートドライブ）のほか、ＳＲＡＭ（スタティックラム）やＤＲＡＭ（ダイナミックランダムアクセスメモリ）などの一時記憶手段などで構成される。記憶部１３０は必ずしも装置内部の記憶手段である必要はなく、装置の外部に設けられた装置、あるいはネットワーク５００を介した先のサーバーやクラウドに設けられた手段であってもよい。
【0101】
出力部１４０は、生成した結果画像などをディスプレイやプリンタに出力する結果画像出力手段、および、アイコンやメニュー画面その他の各種のユーザーインターフェースを表示するための表示手段（図示せず）を備えている。
通信部１５０は、ネットワーク５００を介して外部の装置やシステムとのデータのやり取りをするための通信手段を備えている。
通信部１５０は、ネットワーク５００を介して外部の装置やシステムとのデータのやり取りをするための通信手段を備えている。
【0102】
６－２．密度勾配可視化システムの構成例
次に、密度勾配可視化システムの構成の一例について図２５、図２６を用いて説明する。
図２５は、情報処理システム１００がスタンドアローン構成の場合の一例を示したシステム構成図である。
この場合、情報処理システム１００を構成する入力部１１０、制御部１２０、記憶部１３０、出力部１４０、通信部１５０は、端末１０の内部に構築され、いわゆるスタンドアローンの構成をとることができる。
次に、密度勾配可視化システムの構成の一例について図２５、図２６を用いて説明する。
図２５は、情報処理システム１００がスタンドアローン構成の場合の一例を示したシステム構成図である。
この場合、情報処理システム１００を構成する入力部１１０、制御部１２０、記憶部１３０、出力部１４０、通信部１５０は、端末１０の内部に構築され、いわゆるスタンドアローンの構成をとることができる。
【0103】
端末１０には、外付けで、キーボードやタブレットなどの入力部１１０、ディスプレイなどの出力部１４０を備えるようにしてもよいし、スマートフォンやタブレット端末、ノート型パソコンのように、入力部１１０や出力部１４０も含めて一体型の構成を採用することもできる。
情報処理システム１００と測定システム２００とはメモリカード又はネットワーク５００またはインターフェースケーブルを介してデータのやり取りを行うことができる。なお、ネットワーク５００には有線ＬＡＮのほか無線ＬＡＮ、その他の有線又は無線通信などを利用することができる。
情報処理システム１００と測定システム２００とはメモリカード又はネットワーク５００またはインターフェースケーブルを介してデータのやり取りを行うことができる。なお、ネットワーク５００には有線ＬＡＮのほか無線ＬＡＮ、その他の有線又は無線通信などを利用することができる。
【0104】
図２６は、情報処理システム１００がネットワーク５００を介して、クライアントサーバ、クラウド、ピアツーピア構成となっている場合の一例を示したシステム構成図である。
この場合、情報処理システム１００を構成する入力部１１０、制御部１２０、記憶部１３０、出力部１４０、通信部１５０は、端末１０や端末２０、あるいはサーバー又はクラウド３０のいずれかに存在すればよく、分散して配置されていてもよい。
この場合、情報処理システム１００を構成する入力部１１０、制御部１２０、記憶部１３０、出力部１４０、通信部１５０は、端末１０や端末２０、あるいはサーバー又はクラウド３０のいずれかに存在すればよく、分散して配置されていてもよい。
【0105】
例えば、サーバー又はクラウド３０に制御部１２０の参照画像取得手段及び計測画像取得手段、輪郭抽出手段、密度勾配算出手段、結果画像生成手段、記憶部１３０を備え、端末１０又は端末２０には制御部１２０の一部と、入力部１１０、出力部１４０、通信部１５０を備え、密度勾配可視化のためのほとんどの処理をサーバー又はクラウド３０が行い、端末１０又は端末２０では入力操作と結果のディスプレイ表示だけを行う、という構成であっても構わない（クライアントサーバ、又はクラウド構成）。
あるいは、サーバーまたはクラウド３０を備えないで、端末１０と端末２０との間で、制御部１２０の各手段、記憶部１３０の各データを分散して持ち合う構成であってもよい（ピアツーピア構成）。
あるいは、サーバーまたはクラウド３０を備えないで、端末１０と端末２０との間で、制御部１２０の各手段、記憶部１３０の各データを分散して持ち合う構成であってもよい（ピアツーピア構成）。
【0106】
なお、この場合においても、情報処理システム１００と測定システム２００とはメモリカード又はネットワーク５００またはインターフェースケーブルを介してデータのやり取りを行うことができる。なお、ネットワーク５００には有線ＬＡＮのほか無線ＬＡＮ、その他の有線又は無線通信などを利用することができる。
【0107】
【0108】
７－１．全体処理フロー
図２７は、本発明の密度勾配可視化処理の全体の一例を示すフローチャートである。
制御部１２０は、測定システム２００に対し、背景画像をセットし（ステップＳ１）、まず観測対象となる擾乱の密度分布が発生していない状態で計測領域を撮影して参照画像を取得し（ステップＳ２）、観測対象となる擾乱の密度分布が発生した状態で計測領域を撮影して計測画像を取得する（ステップＳ３）。
図２７は、本発明の密度勾配可視化処理の全体の一例を示すフローチャートである。
制御部１２０は、測定システム２００に対し、背景画像をセットし（ステップＳ１）、まず観測対象となる擾乱の密度分布が発生していない状態で計測領域を撮影して参照画像を取得し（ステップＳ２）、観測対象となる擾乱の密度分布が発生した状態で計測領域を撮影して計測画像を取得する（ステップＳ３）。
【0109】
次いで、制御部１２０は、取得した参照画像と計測画像に対して輪郭抽出処理を行い、参照輪郭画像と計測輪郭画像を生成し（ステップＳ４）、生成した参照輪郭画像と計測輪郭画像のペアに基づいて、Ｓ－ＢＯＳ法を適用して、位相差、変位量、ないし密度勾配を算出する（ステップＳ５）。
以上の処理を、背景画像として横縞および縦縞を使用した場合についてそれぞれ行う。
最終的に、算出した位相差、変位量、ないし密度勾配に基づいて計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する（ステップＳ６）。
ここで、ステップＳ４～Ｓ６は、所定の画像に対して繰り返し使用される処理なので、サブルーチンとして扱い、以下詳述する。
以上の処理を、背景画像として横縞および縦縞を使用した場合についてそれぞれ行う。
最終的に、算出した位相差、変位量、ないし密度勾配に基づいて計測領域における観測対象の密度分布の状態を可視化した結果画像を出力する（ステップＳ６）。
ここで、ステップＳ４～Ｓ６は、所定の画像に対して繰り返し使用される処理なので、サブルーチンとして扱い、以下詳述する。
【0110】
７－２．輪郭抽出処理の処理フロー
次に、制御部１２０の各手段が実行する輪郭抽出処理の処理フロー（ステップＳ４の詳細）の一例について図２８を用いて説明する。
次に、制御部１２０の各手段が実行する輪郭抽出処理の処理フロー（ステップＳ４の詳細）の一例について図２８を用いて説明する。
【0111】
輪郭抽出処理では、参照画像および計測画像に対して輪郭抽出処理を行い、参照輪郭画像および計測輪郭画像を生成するので、まず、参照画像に対する輪郭抽出処理の処理フローについて詳細に説明し、次いで、同様の処理となる、測定画像に対する輪郭抽出処理について簡単に説明する。
【0112】
輪郭抽出処理を行うための所定のフィルタを自動又は手動操作で選択し（ステップＳ４－１）、フィルタテーブルを参照してセットしたフィルタのパラメータを取得する（ステップＳ４－２）。
そして、処理対象となる画像データとして、記憶部１３０から、取得した参照画像データをワーク用のメモリ上にロードして（ステップＳ４－３）、ロードした画像データから１ピクセル分を取得し（ステップＳ４－４）、セットしたフィルタのサイズに応じて必要な周辺画素を取得して、フィルタ処理を実行し（ステップＳ４－５）、フィルタ適用後の画素値の変化率（グレースケール画像の背景画像を用いた場合において画素の輝度値を使用する場合は「輝度値の変化率」）をメモリに格納する（ステップＳ４－６）。
ここで、フィルタとして、縦方向検出用と横方向検出用の２種類あるが、背景画像として横縞をした場合には縦方向検出用を適用し、背景画像として縦縞を使用した場合には横方向検出用を適用して、ステップＳ４－４からステップＳ４－６の処理を行う。
そして、処理対象となる画像データとして、記憶部１３０から、取得した参照画像データをワーク用のメモリ上にロードして（ステップＳ４－３）、ロードした画像データから１ピクセル分を取得し（ステップＳ４－４）、セットしたフィルタのサイズに応じて必要な周辺画素を取得して、フィルタ処理を実行し（ステップＳ４－５）、フィルタ適用後の画素値の変化率（グレースケール画像の背景画像を用いた場合において画素の輝度値を使用する場合は「輝度値の変化率」）をメモリに格納する（ステップＳ４－６）。
ここで、フィルタとして、縦方向検出用と横方向検出用の２種類あるが、背景画像として横縞をした場合には縦方向検出用を適用し、背景画像として縦縞を使用した場合には横方向検出用を適用して、ステップＳ４－４からステップＳ４－６の処理を行う。
【0113】
なお、輪郭抽出処理のためのフィルタ処理としては、前述のように、一例として、線形フィルタの微分フィルタに分類されるＳＯＢＥＬフィルタを用いた積和演算（畳み込み演算）による線形のフィルタリング処理などが挙げられ、前述の式（５）～式（７）の処理を利用することができる。
この場合、例えば、プログラムによって表示されたフィルタ一覧から、３×３サイズのＳＯＢＥＬフィルタを選択した場合、図３５のフィルタテーブルのＩＤ＝Ｆ２０のフィルタ値（パラメータ）を取得して、フィルタを適用しようとする画素を含めた周辺の３×３の画素値（又は輝度値）との積和演算を行う（図１５参照）。
なお、本発明の実験では、グレースケール画像の背景画像を用いるので、画素値としては輝度値を用いるが、カラー画像を用いる場合には、ＲＧＢの画素値を用いてもよい。
この場合、例えば、プログラムによって表示されたフィルタ一覧から、３×３サイズのＳＯＢＥＬフィルタを選択した場合、図３５のフィルタテーブルのＩＤ＝Ｆ２０のフィルタ値（パラメータ）を取得して、フィルタを適用しようとする画素を含めた周辺の３×３の画素値（又は輝度値）との積和演算を行う（図１５参照）。
なお、本発明の実験では、グレースケール画像の背景画像を用いるので、画素値としては輝度値を用いるが、カラー画像を用いる場合には、ＲＧＢの画素値を用いてもよい。
【0114】
ここで、図３５に示すフィルタテーブルには、フィルタ種別、サイズ毎に所定のフィルタ値（パラメータ）が格納されており、フィルタの種類にもよるが、横方向検出用および縦方向検出用のフィルタ値を格納することができる。フィルタ値（パラメータ）としては、各フィルタの種別に応じて、一般的に知られている値を格納してもよいし、シュリーレン法によって光学的に取得した画像と、ＢＯＳ法によって生成した結果画像との差分による機械学習によって得た最適化されたフィルタ値（パラメータ）を格納して利用してもよい。
【0115】
そして、処理対象となる画像の画素について順次、輪郭抽出処理を行い（ステップＳ４－７・Ｙｅｓ）、残りの画素がなくなると（ステップＳ４－７・Ｎｏ）、フィルタ適用後の画素値の変化率（背景画像がグレースケール画像の場合は輝度値の変化率）を統合して輪郭画像（参照輪郭画像）を生成して、記憶部１３０に格納する（ステップＳ４－８）。
【0116】
次に、制御部１２０は、処理対象となる画像データとして、記憶部１３０から、取得した計測画像データをワーク用のメモリ上にロードして（ステップＳ４－３）、順次、ステップＳ４－４からＳ４－７の処理を実行し、輪郭画像（計測輪郭画像）を生成して、記憶部１３０に格納する（ステップＳ４－８）。
【0117】
７－３．密度勾配算出処理の処理フロー（Ｓ－ＢＯＳ法の適用）
次に、画像に対してＳ－ＢＯＳ法を適用して密度勾配を算出するための、制御部１２０の各手段が実行する密度勾配算出処理の処理フロー（ステップＳ５の詳細）の一例について図２９を用いて説明する。
次に、画像に対してＳ－ＢＯＳ法を適用して密度勾配を算出するための、制御部１２０の各手段が実行する密度勾配算出処理の処理フロー（ステップＳ５の詳細）の一例について図２９を用いて説明する。
【0118】
密度勾配算出処理の過程で、参照輪郭画像と計測輪郭画像との位相差や変位量なども算出でき、これらを利用して、擾乱を各種の視点で可視化した結果画像を得ることができるので以下説明する。
まず、制御部１２０は、背景２３０と観測部位３００の中心との距離ＺＤ、観測部位３００の中心と撮影装置２１０のレンズ２１の中心との距離ＺＡ、観測部位３００の幅２ΔＺＤ、レンズの焦点距離ｆ、撮影倍率ｍ、グラッドストーン・デール定数Ｇ、密度ρ０（観測部位の気体の基準密度）、背景画像を構成する所定の空間周波数を持つ波の波長λ、などの測定システムの各種パラメータをメモリ上にロードする（ステップＳ５－１）。これらのパラメータは、式（１）～式（４）の計算過程で用いるものである。
まず、制御部１２０は、背景２３０と観測部位３００の中心との距離ＺＤ、観測部位３００の中心と撮影装置２１０のレンズ２１の中心との距離ＺＡ、観測部位３００の幅２ΔＺＤ、レンズの焦点距離ｆ、撮影倍率ｍ、グラッドストーン・デール定数Ｇ、密度ρ０（観測部位の気体の基準密度）、背景画像を構成する所定の空間周波数を持つ波の波長λ、などの測定システムの各種パラメータをメモリ上にロードする（ステップＳ５－１）。これらのパラメータは、式（１）～式（４）の計算過程で用いるものである。
【0119】
次に、記憶部１３０から、処理対象となる参照輪郭画像、計測輪郭画像をメモリ上にロードし（ステップＳ５－２）、参照輪郭画像の微分画像を生成するとともに（ステップＳ５－３）、参照輪郭画像、参照輪郭画像の微分画像、及び計測輪郭画像の３種に基づいて参照輪郭画像と前記計測輪郭画像との位相差を算出する（ステップＳ５－４）。
そして、参照輪郭画像と前記計測輪郭画像との位相差、及び測定システムのパラメータから、流れ場の密度勾配∂ρ/∂rを算出し（ステップＳ５－５）、算出した各画素の密度勾配をメモリに格納する（ステップＳ５－５）。
ここで、流れ場の密度勾配∂ρ/∂rは、背景画像として横縞を使用する場合においてＹ軸方向の変位を検出しようとする場面では∂ρ/∂ｙ、背景画像として縦縞を使用する場合においてＸ軸方向の変位を検出しようとする場面では∂ρ/∂ｘが用いられる。
以上の処理を、背景画像として横縞および縦縞を使用した場合についてそれぞれ行う。
そして、参照輪郭画像と前記計測輪郭画像との位相差、及び測定システムのパラメータから、流れ場の密度勾配∂ρ/∂rを算出し（ステップＳ５－５）、算出した各画素の密度勾配をメモリに格納する（ステップＳ５－５）。
ここで、流れ場の密度勾配∂ρ/∂rは、背景画像として横縞を使用する場合においてＹ軸方向の変位を検出しようとする場面では∂ρ/∂ｙ、背景画像として縦縞を使用する場合においてＸ軸方向の変位を検出しようとする場面では∂ρ/∂ｘが用いられる。
以上の処理を、背景画像として横縞および縦縞を使用した場合についてそれぞれ行う。
【0120】
なお、上記の説明では、本発明で提案する新たなＳ－ＢＯＳ法を用いる例について説明したが、本発明の輪郭抽出処理を行った輪郭画像を利用する場合は、従来のＳ－ＢＯＳ法を適用しても、（位相反転の検出、補正を行わない限り）位相差の検出可能範囲が原則として±１／２πの範囲になること以外は、新たなＳ－ＢＯＳ法と同様の結果（位相差、変位量、密度勾配の算出結果）を得ることができるので、全体処理フロー（図２７）のステップＳ５では、様々な種類のＳ－ＢＯＳ法を適用することができる。
【0121】
また、上記の説明では、処理対象となる画像データとして、参照輪郭画像、計測輪郭画像を対象とする場合について説明したが、位相反転の検出を行う場合には、元画像である参照画像と計測画像に対しても、Ｓ－ＢＯＳ法を適用して、参照画像と計測画像の位相差、変位量、および密度勾配を算出して利用する（図３３のステップＰ１～Ｐ２を参照）。
【0122】
７－４．結果画像出力処理サブルーチン
次に、算出した位相差、変位量、又は密度勾配から、現象を可視化した結果画像を生成して出力結果画像出力処理の処理フロー（ステップＳ６の詳細）の一例について図３０を用いて説明する。
以下の説明では、一例として、密度勾配から結果画像を出力する処理について説明するが、位相差や変位量に基づいて現象を可視化した結果画像を生成、出力する場合も同様である。
なお、説明の便宜のため、背景画像としてグレースケール画像を使用した場合を例に説明することとするが、背景画像としてカラー画像を用いた場合も同様である。
次に、算出した位相差、変位量、又は密度勾配から、現象を可視化した結果画像を生成して出力結果画像出力処理の処理フロー（ステップＳ６の詳細）の一例について図３０を用いて説明する。
以下の説明では、一例として、密度勾配から結果画像を出力する処理について説明するが、位相差や変位量に基づいて現象を可視化した結果画像を生成、出力する場合も同様である。
なお、説明の便宜のため、背景画像としてグレースケール画像を使用した場合を例に説明することとするが、背景画像としてカラー画像を用いた場合も同様である。
【0123】
制御部１２０は、記憶部１３０から、画素ごとの密度勾配∂ρ/∂rをメモリ上にロードし（ステップＳ６－１）、画素ごとの密度勾配∂ρ/∂rの数値（アナログ値）に対し階調を表現するビット数に応じてデジタル値に換算する（ステップＳ６－２）。
換算した値に応じて輝度を設定し（ステップＳ６－３）、各画素を密度勾配の正負に応じて色分けし（例えば正の場合は赤色、負の場合は青色）（ステップＳ６－４）、色分けや輝度に応じて各画素をプロットして結果画像を生成（ステップＳ６－５）し、結果画像を出力する（ステップＳ６－６）。
結果画像の生成に際しては、背景画像として横縞又は縦縞の一方を使用した場合の位相差、変位量、又は密度勾配の計算結果に基づいて生成してもよいし、背景画像として横縞又は縦縞の両方を使用してそれぞれ求めた位相差、変位量、又は密度勾配の計算結果を二乗平均等で重畳した結果に基づいて生成してもよい。
換算した値に応じて輝度を設定し（ステップＳ６－３）、各画素を密度勾配の正負に応じて色分けし（例えば正の場合は赤色、負の場合は青色）（ステップＳ６－４）、色分けや輝度に応じて各画素をプロットして結果画像を生成（ステップＳ６－５）し、結果画像を出力する（ステップＳ６－６）。
結果画像の生成に際しては、背景画像として横縞又は縦縞の一方を使用した場合の位相差、変位量、又は密度勾配の計算結果に基づいて生成してもよいし、背景画像として横縞又は縦縞の両方を使用してそれぞれ求めた位相差、変位量、又は密度勾配の計算結果を二乗平均等で重畳した結果に基づいて生成してもよい。
【0124】
８．位相反転の検出と補正、および処理フロー
前述のように、Ｓ－ＢＯＳ法の種類によっては、位相差の検出可能範囲が±１／２π（従来のＳ－ＢＯＳ法）、または±π（新たなＳ－ＢＯＳ法）に限定され、それを超えると位相反転が発生し、正確な位相差、変位量、または密度勾配を算出できないという問題が生じる。
特に、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱など、観測対象となる現象が、大きな密度勾配と小さな密度勾配を同時に可視化する必要がある場合は顕著となる。
これは、部分的に、変位量が大きくなる領域と、変位量がそれほどでもない領域の両方が混在し、背景画像となる縞画像の間隔（波長）の設定との関係で、部分的に、位相反転が発生してしまうことによるものである。
前述のように、Ｓ－ＢＯＳ法の種類によっては、位相差の検出可能範囲が±１／２π（従来のＳ－ＢＯＳ法）、または±π（新たなＳ－ＢＯＳ法）に限定され、それを超えると位相反転が発生し、正確な位相差、変位量、または密度勾配を算出できないという問題が生じる。
特に、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱など、観測対象となる現象が、大きな密度勾配と小さな密度勾配を同時に可視化する必要がある場合は顕著となる。
これは、部分的に、変位量が大きくなる領域と、変位量がそれほどでもない領域の両方が混在し、背景画像となる縞画像の間隔（波長）の設定との関係で、部分的に、位相反転が発生してしまうことによるものである。
【0125】
ここで、図１８によれば、元画像（測定画像）の縞画像と、輪郭画像（計測輪郭画像）の縞画像とを比較すると、両者の平均間隔が約１／２、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約１／２、周期も約１／２になっているので、周波数に換算した場合約２倍の周波数に相当することになり、空間分解能も約２倍になっており、流れ場の擾乱の再現性を高めることができることが分かる。
【0126】
そこで、元画像に対する輪郭画像の平均間隔が約１／２、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約１／２になっていること、すなわち逆から見れば、輪郭画像に対して元画像は平均間隔が約２倍、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約２倍になっていることを利用して、輪郭画像では位相反転が発生していても、元画像では位相反転していない、という元画像に対する輪郭画像の関係を利用して、輪郭画像における位相反転の検出および補正を行うこととした。
【0127】
以下、図３１～図３３を用いて、本発明の位相反転の検出および補正方法について説明する。
図３１は、（ａ）計測画像と、（ｂ）計測画像に対して輪郭抽出処理を行った計測輪郭画像を模式的に示した図である。
図３１によれば、輪郭画像を構成する縞画像は、その間隔が計測画像の約１／２になっていることが分かる。この場合、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約１／２になっているので、πの範囲についても、（ｂ）計測輪郭画像は、（ａ）計測画像に対し約１／２になっている。
ここで、（ｂ）計測輪郭画像の一番下の部分では、上向きの矢印で示したように、変位が±πを超えているのに対し、元画像である（ａ）計測画像の一番下の上向きの矢印の部分では、変位が±πの範囲に収まっていることが分かる。
図３１は、（ａ）計測画像と、（ｂ）計測画像に対して輪郭抽出処理を行った計測輪郭画像を模式的に示した図である。
図３１によれば、輪郭画像を構成する縞画像は、その間隔が計測画像の約１／２になっていることが分かる。この場合、縞画像の輝度値の分布（所定の空間周波数を持つ波）の波長も約１／２になっているので、πの範囲についても、（ｂ）計測輪郭画像は、（ａ）計測画像に対し約１／２になっている。
ここで、（ｂ）計測輪郭画像の一番下の部分では、上向きの矢印で示したように、変位が±πを超えているのに対し、元画像である（ａ）計測画像の一番下の上向きの矢印の部分では、変位が±πの範囲に収まっていることが分かる。
【0128】
このような両者の関係を利用して、位相反転が発生した場合に、輪郭抽出処理を行った計測輪郭画像における位相差から推定される複数の変位候補について、輪郭抽出処理を行う前の元画像である計測画像における位相差から判断される変位量を基礎として、複数の変位候補から真の変位量を特定する処理の仕組みを模式的に示したものが図３２である。
なお、図中の元画像の位相差（■印）、元画像の変位量（▲印）、輪郭画像の位相差（●印）、輪郭画像の変位候補（▼印）などは、事前に、参照画像と計測画像のペアに対してＳ－ＢＯＳ法を適用して得た数値、及び参照輪郭画像と計測輪郭画像のペアに対してＳ－ＢＯＳ法を適用して得た数値を利用するものとする。
なお、図中の元画像の位相差（■印）、元画像の変位量（▲印）、輪郭画像の位相差（●印）、輪郭画像の変位候補（▼印）などは、事前に、参照画像と計測画像のペアに対してＳ－ＢＯＳ法を適用して得た数値、及び参照輪郭画像と計測輪郭画像のペアに対してＳ－ＢＯＳ法を適用して得た数値を利用するものとする。
【0129】
図３２の縦軸は位相差、横軸は変位量を示している。図３２では、便宜的に横軸も位相で表し、２πを輪郭画像（実線）の１波長分の変位量とした。したがって、輪郭画像（実線）の変位は位相差そのものを示している。
また、輪郭画像（実線）の２波長が、元画像（破線）の１波長に相当することが示されている。
また、輪郭画像（実線）の２波長が、元画像（破線）の１波長に相当することが示されている。
【0130】
この場合、高周波画像である輪郭画像（実線）の位相差がπに達すると－πに不連続に跳ぶため、位相の反転が発生する。一方、低周波画像である元画像（破線）の変位量は、位相差が±π動く間に、輪郭画像（実線）の変位で±２πの範囲を動くことになる。
【0131】
ここで、前提として、元画像（破線）から求められる変位量は、高周波の変動に追従できないためやや解像度が低く、真の変位量とは一致しない場合も多いが、この点のあいまいさがあっても、元画像（破線）の変位量からあたりをつけて、輪郭画像（実線）の変位量の候補から、真の変位量を特定することが可能となる。これを図３２の各記号を用いて説明すると次のようになる。
【0132】
まず、高周波画像である輪郭画像（実線）の位相差（●印）が定まると、低周波画像である元画像（破線）の変位の動き得る範囲で、真の変位量の候補は３点存在する（▼印）。次に、低周波画像である元画像（破線）の位相差（■印）に相当する変位量（▲印）に最も近い候補を真の変位量（★印）として特定することで、位相反転の検出と補正を行うことが可能となる。
なお、以上の説明では、背景画像として横縞を使用した場合について説明したが、背景画像として縦縞を使用した場合も同様である。
なお、以上の説明では、背景画像として横縞を使用した場合について説明したが、背景画像として縦縞を使用した場合も同様である。
【0133】
次に、位相反転検出処理の処理フローの一例について図３３を用いて説明する。
まず、元画像における位相差を算出し（ステップＰ１）、位相差から元画像における変位量を算出する（ステップＰ２）。そして、輪郭画像における位相差を算出し（ステップＰ３）、位相差から輪郭画像における複数の変位候補を出力し（ステップＰ４）、元画像における変位量に基づいて、輪郭画像における変位候補の中から真の変位量を特定する（ステップＰ５）。
以上のように、本発明によれば、輪郭画像と元画像を利用して位相反転の検出と補正が可能となる。
まず、元画像における位相差を算出し（ステップＰ１）、位相差から元画像における変位量を算出する（ステップＰ２）。そして、輪郭画像における位相差を算出し（ステップＰ３）、位相差から輪郭画像における複数の変位候補を出力し（ステップＰ４）、元画像における変位量に基づいて、輪郭画像における変位候補の中から真の変位量を特定する（ステップＰ５）。
以上のように、本発明によれば、輪郭画像と元画像を利用して位相反転の検出と補正が可能となる。
【0134】
９．精度検証の処理フロー
前述のように、図２１～図２３を用いて本発明の新たな精度検証方法、および新たな精度検証方法を利用して検証した本発明の効果を説明したが、この新たな精度検証方法の処理フローについて、図３４を用いて説明する。
前述のように、図２１～図２３を用いて本発明の新たな精度検証方法、および新たな精度検証方法を利用して検証した本発明の効果を説明したが、この新たな精度検証方法の処理フローについて、図３４を用いて説明する。
【0135】
まず、参照画像から、所定の矩形領域を観測対象として設定するための設定処理として、図２１に示すように、幅Ａ、移動距離Ｂの設定を受付け（ステップＱ１）、受け付けた設定値に基づいて、参照画像から幅Ａの矩形領域を取得し（ステップＱ２）、矩形領域を所定の距離Ｂだけ平行移動させた画像を精度検証用の計測画像として生成する（ステップＱ３）。
次に、参照画像と、生成した精度検証用の計測画像との位相差に基づいて、元画像の変位量を算出し（ステップＱ４）、参照画像および生成した精度検証用の計測画像に対し、輪郭抽出処理を行って、参照輪郭画像、精度検証用の計測輪郭画像をそれぞれ生成する（ステップＱ５）。
次に、参照画像と、生成した精度検証用の計測画像との位相差に基づいて、元画像の変位量を算出し（ステップＱ４）、参照画像および生成した精度検証用の計測画像に対し、輪郭抽出処理を行って、参照輪郭画像、精度検証用の計測輪郭画像をそれぞれ生成する（ステップＱ５）。
【0136】
そして、参照輪郭画像と、生成した精度検証用の計測輪郭画像との位相差に基づいて、輪郭画像の変位量を算出し（ステップＱ６）、所定の幅Ａの矩形領域を所定の距離Ｂだけ平行移動させたことに対応する変位量を理論値として算出する（ステップＱ７）。
最後に、変位量の理論値と、元画像の変位量と、輪郭画像の変位量とを比較して、密度勾配可視化の精度を検証する（ステップＱ８）。
最後に、変位量の理論値と、元画像の変位量と、輪郭画像の変位量とを比較して、密度勾配可視化の精度を検証する（ステップＱ８）。
【0137】
１０．小括
以上のように、本発明の密度勾配可視化方法によれば、光学系に手を入れる必要がなく、簡易な光学系のままで、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の現象を高い精度で可視化できる。
また、本発明の輪郭画像を用いる手法、すなわち、参照輪郭画像と計測輪郭画像のペアに基づいてＳ－ＢＯＳ法を適用して、位相差、変位量、密度勾配を算出する手法においては、様々なＳ－ＢＯＳ法を適用することができるので、汎用性が高いというメリットもある。
以上のように、本発明の密度勾配可視化方法によれば、光学系に手を入れる必要がなく、簡易な光学系のままで、衝撃波を含む超音速流れ場などの急激な密度変化を伴う擾乱の現象を高い精度で可視化できる。
また、本発明の輪郭画像を用いる手法、すなわち、参照輪郭画像と計測輪郭画像のペアに基づいてＳ－ＢＯＳ法を適用して、位相差、変位量、密度勾配を算出する手法においては、様々なＳ－ＢＯＳ法を適用することができるので、汎用性が高いというメリットもある。
【0138】
また、本発明で提案する新たな精度検証方法による検証結果によれば、輪郭画像は変位量の不連続的な変化にも対応できることが分かった。
また、狭い範囲に生じた変位について従来法では過小に評価されていたのに対し、参照輪郭画像と計測輪郭画像のペアを用いる本発明の手法によれば、より真値に近い結果が得られた。
特に、衝撃波と混合層を有する圧縮性流れ場の密度勾配計測ないし可視化の場面では、計測に高い空間分解能が求められる衝撃波面において輪郭画像の結果は元画像と比較して密度勾配のピーク値が約３０％向上するという結果も得られた。
また、狭い範囲に生じた変位について従来法では過小に評価されていたのに対し、参照輪郭画像と計測輪郭画像のペアを用いる本発明の手法によれば、より真値に近い結果が得られた。
特に、衝撃波と混合層を有する圧縮性流れ場の密度勾配計測ないし可視化の場面では、計測に高い空間分解能が求められる衝撃波面において輪郭画像の結果は元画像と比較して密度勾配のピーク値が約３０％向上するという結果も得られた。
【0139】
また、参照輪郭画像と計測輪郭画像のペアを用いる本発明の手法によれば、背景画像の縞画像の間隔を不用意に狭くして位相反転を招来させずに、元画像および輪郭画像における位相差や変位量を併用して、位相反転を防止しつつ、解像度の高い結果画像を提供することが可能となる。
【産業上の利用可能性】
【0140】
本発明は流体の可視化技術のみならず、非接触で変位を測定する技術や物体の形状を測定する技術にも適用可能である。例えば、材料に縞模様を貼り付けひずみ量などの変位を測定したり、物体の形状を測定して形状や位置関係を把握したうえで、産業用ロボットの停止精度を監視したりする用途に適用することが考えられる。
【符号の説明】
【0141】
１０ 端末
２０ 端末
３０ サーバー又はクラウド
１００ 情報処理システム
１１０ 機能ブロック図の入力部
１２０ 機能ブロック図の制御部
１３０ 機能ブロック図の記憶部
１４０ 機能ブロック図の出力部
１５０ 機能ブロック図の通信部
２００ 測定システム
２１０ 撮影装置
２１２ レンズ
２１４ 撮像素子
２２０ 風洞
２２２ 観測窓
２２４ 風洞センターボディ
２３０ 背景
２４０ 光源
２５０ 大気
２５２ 気流Ｍ１
２６０ 真空タンク
３００ 観測部位
３１０ 衝撃波
３２０ 混合層
３３０ 密度勾配が正の領域
５００ ネットワーク
２０ 端末
３０ サーバー又はクラウド
１００ 情報処理システム
１１０ 機能ブロック図の入力部
１２０ 機能ブロック図の制御部
１３０ 機能ブロック図の記憶部
１４０ 機能ブロック図の出力部
１５０ 機能ブロック図の通信部
２００ 測定システム
２１０ 撮影装置
２１２ レンズ
２１４ 撮像素子
２２０ 風洞
２２２ 観測窓
２２４ 風洞センターボディ
２３０ 背景
２４０ 光源
２５０ 大気
２５２ 気流Ｍ１
２６０ 真空タンク
３００ 観測部位
３１０ 衝撃波
３２０ 混合層
３３０ 密度勾配が正の領域
５００ ネットワーク
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【図26】
【図27】
【図28】
【図29】
【図30】
【図31】
【図32】
【図33】
【図34】
【図35】
