JP3390975B2 - 画像流速推定方法、装置および画像流速推定プログラムを記録した記録媒体 - Google Patents
画像流速推定方法、装置および画像流速推定プログラムを記録した記録媒体Info
- Publication number
- JP3390975B2 JP3390975B2 JP03906798A JP3906798A JP3390975B2 JP 3390975 B2 JP3390975 B2 JP 3390975B2 JP 03906798 A JP03906798 A JP 03906798A JP 3906798 A JP3906798 A JP 3906798A JP 3390975 B2 JP3390975 B2 JP 3390975B2
- Authority
- JP
- Japan
- Prior art keywords
- image
- spatial frequency
- velocity
- energy
- estimating
- 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.)
- Expired - Lifetime
Links
Description
【０００１】
【発明の属する技術分野】本発明は、時系列画像中の物
体の移動速度を、制御量を適応的に変化させながら推定
し、気象レーダーエコー画像からの降水量の変化予測や
流体工学における流体の挙動の解析等、非剛体の動きで
生成・消滅が著しい対象の移動速度を安定に検出する方
法と装置に関する。
体の移動速度を、制御量を適応的に変化させながら推定
し、気象レーダーエコー画像からの降水量の変化予測や
流体工学における流体の挙動の解析等、非剛体の動きで
生成・消滅が著しい対象の移動速度を安定に検出する方
法と装置に関する。
【０００２】
【従来の技術】従来、画像中の剛体系の物体の移動ベク
トルを推定する場合は、物体上の照明変化がほとんどな
いモデルを採用することが多い。その検出方法には、オ
プティカルフローや相互相関法に基づいた方法が中心的
である（文献［１］：B.K.P.Horn & B.G.Schunck, "Det
ermining optical flow", Artificial Intelligence, v
ol.17, pp.185-203, 1981 ）。
トルを推定する場合は、物体上の照明変化がほとんどな
いモデルを採用することが多い。その検出方法には、オ
プティカルフローや相互相関法に基づいた方法が中心的
である（文献［１］：B.K.P.Horn & B.G.Schunck, "Det
ermining optical flow", Artificial Intelligence, v
ol.17, pp.185-203, 1981 ）。
【０００３】一方、非剛体の物体に対する適切な移動ベ
クトルはないと言える。これは、連続する画像間であっ
ても、物体の輪郭線、濃淡値等の属性が同時に変化する
ために、明確な対応づけを行えないことに起因する。即
ち、複雑な物体に対しては、統計的な類似性を追従して
いく相互相関法が適用されることが多い。例えば、気象
レーダーエコー画像を用いた場合、降水パターンは非剛
体的に変化するものの、相互相関法が適用されている。
クトルはないと言える。これは、連続する画像間であっ
ても、物体の輪郭線、濃淡値等の属性が同時に変化する
ために、明確な対応づけを行えないことに起因する。即
ち、複雑な物体に対しては、統計的な類似性を追従して
いく相互相関法が適用されることが多い。例えば、気象
レーダーエコー画像を用いた場合、降水パターンは非剛
体的に変化するものの、相互相関法が適用されている。
【０００４】オプティカルフロー法では、２枚の画像か
ら正則化と呼ばれる枠組みで対象の動き速度を推定す
る。推定式は次のように導出される。まず、２枚の画像
に基づいて誤差評価関数として定義された式から変分法
（もしくはＥｕｌｅｒ−Ｌａｇｒａｎｇｅ法）に従って
連立１次方程式を得る。この連立１次方程式を解く動作
を緩和法（ｒｅｌａｘａｔｉｏｎ）により反復誤差が小
さくなるまで繰り返す。このようにして、収束した解が
２枚の画像間における、対象の移動速度ベクトル（ｕ，
ｖ）となる。ただし、安定にかつ精度よく対象の移動速
度を推定するためには、推定式に含まれる拘束条件（制
御量α）に対する重みづけ量を試行錯誤的に調節する必
要がある。
ら正則化と呼ばれる枠組みで対象の動き速度を推定す
る。推定式は次のように導出される。まず、２枚の画像
に基づいて誤差評価関数として定義された式から変分法
（もしくはＥｕｌｅｒ−Ｌａｇｒａｎｇｅ法）に従って
連立１次方程式を得る。この連立１次方程式を解く動作
を緩和法（ｒｅｌａｘａｔｉｏｎ）により反復誤差が小
さくなるまで繰り返す。このようにして、収束した解が
２枚の画像間における、対象の移動速度ベクトル（ｕ，
ｖ）となる。ただし、安定にかつ精度よく対象の移動速
度を推定するためには、推定式に含まれる拘束条件（制
御量α）に対する重みづけ量を試行錯誤的に調節する必
要がある。
【０００５】
【発明が解決しようとする課題】本発明の目的は、試行
錯誤的に与えていた制御量αをオプティカルフローを算
出する方程式中に与え、画像処理により最適な制御量を
自動推定し、非剛体系の物体の移動ベクトルの推定精度
の向上させる、時系列画像中の流体の移動速度をオプテ
ィカルフローを用いて推定する画像流速推定方法、装置
および画像流速推定プログラムを記録した記録媒体を提
供することにある。
錯誤的に与えていた制御量αをオプティカルフローを算
出する方程式中に与え、画像処理により最適な制御量を
自動推定し、非剛体系の物体の移動ベクトルの推定精度
の向上させる、時系列画像中の流体の移動速度をオプテ
ィカルフローを用いて推定する画像流速推定方法、装置
および画像流速推定プログラムを記録した記録媒体を提
供することにある。
【０００６】
【課題を解決するための手段】本発明の画像流速推定方
法は、移動する対象を含む２次元画像を入力する画像入
力段階と、入力された画像を時系列画像として蓄積する
画像蓄積段階と、蓄積されている２次元時系列画像の空
間周波数成分を解析するとともに、画像の濃淡値の標準
偏差を同時に解析する画像解析段階と、前記画像解析段
階で得られた画像の空間周波数成分と濃淡値の標準偏差
の２つの画像特徴量を予め設計してある線形方程式に代
入して最適正則化係数αを決定する正則化係数算出段階
と、連続する２枚の画像から、各最適正則化係数αに対
する、移動する対象の速度ベクトルをオプティカルフロ
ー法を用いて推定する速度推定段階と、推定された速度
ベクトルを出力する出力段階を有する。
法は、移動する対象を含む２次元画像を入力する画像入
力段階と、入力された画像を時系列画像として蓄積する
画像蓄積段階と、蓄積されている２次元時系列画像の空
間周波数成分を解析するとともに、画像の濃淡値の標準
偏差を同時に解析する画像解析段階と、前記画像解析段
階で得られた画像の空間周波数成分と濃淡値の標準偏差
の２つの画像特徴量を予め設計してある線形方程式に代
入して最適正則化係数αを決定する正則化係数算出段階
と、連続する２枚の画像から、各最適正則化係数αに対
する、移動する対象の速度ベクトルをオプティカルフロ
ー法を用いて推定する速度推定段階と、推定された速度
ベクトルを出力する出力段階を有する。
【０００７】本発明の実施態様によれば、画像解析段階
は、画像を幾つかのサブブロックに分割して、各サブブ
ロック毎に高速フーリエ変換により、エネルギー・水平
方向の空間周波数・垂直方向の空間周波数を軸にもつ３
次元の空間周波数画像に変換し、空間周波数画像の低空
間周波数成分が画像中央に、高空間周波数成分が中央か
ら周辺に向かって分布するように、高速フーリエ変換後
の空間周波数画像の内部を配置変えし、エネルギー画像
をその中心から同心円状に抜き取り、各円ごとにエネル
ギーを求めることにより前記各円の半径に対応する空間
周波数に対するエネルギーの分布を得、各空間周波数と
そのエネルギーとの積の和をエネルギーの総和で割った
値を算出し、これを代表空間周波数とする。
は、画像を幾つかのサブブロックに分割して、各サブブ
ロック毎に高速フーリエ変換により、エネルギー・水平
方向の空間周波数・垂直方向の空間周波数を軸にもつ３
次元の空間周波数画像に変換し、空間周波数画像の低空
間周波数成分が画像中央に、高空間周波数成分が中央か
ら周辺に向かって分布するように、高速フーリエ変換後
の空間周波数画像の内部を配置変えし、エネルギー画像
をその中心から同心円状に抜き取り、各円ごとにエネル
ギーを求めることにより前記各円の半径に対応する空間
周波数に対するエネルギーの分布を得、各空間周波数と
そのエネルギーとの積の和をエネルギーの総和で割った
値を算出し、これを代表空間周波数とする。
【０００８】本発明の実施態様によれば、画像解析段階
は、画像を幾つかのサブブロックに分割して、各サブブ
ロック毎に画像の濃淡値の平均値を算出し、分散値を求
めてその平方根をとることで標準偏差を算出する。
は、画像を幾つかのサブブロックに分割して、各サブブ
ロック毎に画像の濃淡値の平均値を算出し、分散値を求
めてその平方根をとることで標準偏差を算出する。
【０００９】本発明の実施態様によれば、正則化係数算
出段階は、予めオフラインで正弦波の振幅と空間周波数
（波数）をさまざまに変化させながら一定速度で正弦波
パターンを移動させるシミュレーションパターンを予め
生成し、オプティカルフロー法における正則化係数αを
同時にさまざまに変化させながら、シミュレーションパ
ターンから速度ベクトルを推定し、ａｒｃｃｏｓｉｎｅ
関数を用いて推定された速度ベクトルと真の速度ベクト
ルとの誤差を評価して誤差が最小となった時のαを選択
することで、αと振幅と空間周波数の３者の関係を定量
づけ、さらに振幅からシミュレーションパターンの濃淡
値の標準偏差を算出する処理を行って得られる式を前記
線形方程式として用いる。
出段階は、予めオフラインで正弦波の振幅と空間周波数
（波数）をさまざまに変化させながら一定速度で正弦波
パターンを移動させるシミュレーションパターンを予め
生成し、オプティカルフロー法における正則化係数αを
同時にさまざまに変化させながら、シミュレーションパ
ターンから速度ベクトルを推定し、ａｒｃｃｏｓｉｎｅ
関数を用いて推定された速度ベクトルと真の速度ベクト
ルとの誤差を評価して誤差が最小となった時のαを選択
することで、αと振幅と空間周波数の３者の関係を定量
づけ、さらに振幅からシミュレーションパターンの濃淡
値の標準偏差を算出する処理を行って得られる式を前記
線形方程式として用いる。
【００１０】本発明の実施態様によれば、オプティカル
フロー法により定式化された速度を変数とする連立１次
方程式の解法において、加速緩和（Successive Over Re
laxation:SOR）法を適用して反復誤差が一定基準値以下
になった時点の解を推定すべき速度ベクトルとみなす。
フロー法により定式化された速度を変数とする連立１次
方程式の解法において、加速緩和（Successive Over Re
laxation:SOR）法を適用して反復誤差が一定基準値以下
になった時点の解を推定すべき速度ベクトルとみなす。
【００１１】本発明の画像流速推定装置は、移動する対
象を含む２次元画像を入力する画像入力手段と、前記画
像入力手段に入力された画像を時系列画像として蓄積す
る画像蓄積手段と、前記画像蓄積手段に蓄積されている
２次元時系列画像の空間周波数成分を解析するととも
に、画像の濃淡値の標準偏差を同時に解析する画像解析
手段と、前記画像解析手段で得られた画像の空間周波数
成分と濃淡値の標準偏差の２つの画像特徴量を予め設計
してある線形方程式に代入して最適正則化係数αを決定
する正則化係数算出手段と、連続する２枚の画像から、
各最適正則化係数αに対する、移動する対象の速度ベク
トルをオプティカルフロー法を用いて推定する速度推定
手段と、推定された速度ベクトルを出力する出力手段を
有する。
象を含む２次元画像を入力する画像入力手段と、前記画
像入力手段に入力された画像を時系列画像として蓄積す
る画像蓄積手段と、前記画像蓄積手段に蓄積されている
２次元時系列画像の空間周波数成分を解析するととも
に、画像の濃淡値の標準偏差を同時に解析する画像解析
手段と、前記画像解析手段で得られた画像の空間周波数
成分と濃淡値の標準偏差の２つの画像特徴量を予め設計
してある線形方程式に代入して最適正則化係数αを決定
する正則化係数算出手段と、連続する２枚の画像から、
各最適正則化係数αに対する、移動する対象の速度ベク
トルをオプティカルフロー法を用いて推定する速度推定
手段と、推定された速度ベクトルを出力する出力手段を
有する。
【００１２】本発明の画像流速推定プログラムを記録し
たコンピュータ読取り可能な記録媒体は、入力された、
移動する対象を含む２次元画像を時系列画像として画像
蓄積手段に蓄積する画像蓄積処理と、前記画像蓄積手段
に蓄積されている２次元時系列画像の空間周波数成分を
解析するとともに、画像の濃淡値の標準偏差を同時に解
析する画像解析処理と、前記画像解析処理で得られた画
像の空間周波数成分と濃淡値の標準偏差の２つの画像特
徴量を予め設計してある線形方程式に代入して最適正則
化係数αを決定する正則化係数算出処理と、連続する２
枚の画像から、各最適正則化係数αに対する、移動する
対象の速度ベクトルをオプティカルフロー法を用いて推
定する速度推定処理と、推定された速度ベクトルを出力
手段に出力する出力処理を有する。
たコンピュータ読取り可能な記録媒体は、入力された、
移動する対象を含む２次元画像を時系列画像として画像
蓄積手段に蓄積する画像蓄積処理と、前記画像蓄積手段
に蓄積されている２次元時系列画像の空間周波数成分を
解析するとともに、画像の濃淡値の標準偏差を同時に解
析する画像解析処理と、前記画像解析処理で得られた画
像の空間周波数成分と濃淡値の標準偏差の２つの画像特
徴量を予め設計してある線形方程式に代入して最適正則
化係数αを決定する正則化係数算出処理と、連続する２
枚の画像から、各最適正則化係数αに対する、移動する
対象の速度ベクトルをオプティカルフロー法を用いて推
定する速度推定処理と、推定された速度ベクトルを出力
手段に出力する出力処理を有する。
【００１３】
【発明の実施の形態】次に、本発明の実施の形態につい
て図面を参照して説明する。
て図面を参照して説明する。
【００１４】図１を参照すると、本発明の一実施形態の
画像流速推定装置は、移動する対象を含む２次元画像を
入力する画像入力部１００と、入力された画像を時系列
画像として蓄積する画像蓄積部１１０と、画像蓄積部１
１０に蓄積されている２次元時系列画像をいくつかのサ
ブブロックに分割し、各サブブロック内の空間周波数成
分と濃淡値の標準偏差を算出する画像解析部１２０と、
画像解析部１２０で得られた画像の空間周波数成分と濃
淡値の標準偏差の２つの画像特徴量を、予め設計してあ
る線形方程式に代入して各サブブロック毎に最適正則化
係数αを算出する正規化係数算出部１３０と、オプティ
カルフロー法を用いて連続する２つの時系列画像から画
素ごとに移動する物体から速度ベクトルを推定する速度
推定部１４０と、推定された速度ベクトルを出力する出
力部１５０で構成されている。
画像流速推定装置は、移動する対象を含む２次元画像を
入力する画像入力部１００と、入力された画像を時系列
画像として蓄積する画像蓄積部１１０と、画像蓄積部１
１０に蓄積されている２次元時系列画像をいくつかのサ
ブブロックに分割し、各サブブロック内の空間周波数成
分と濃淡値の標準偏差を算出する画像解析部１２０と、
画像解析部１２０で得られた画像の空間周波数成分と濃
淡値の標準偏差の２つの画像特徴量を、予め設計してあ
る線形方程式に代入して各サブブロック毎に最適正則化
係数αを算出する正規化係数算出部１３０と、オプティ
カルフロー法を用いて連続する２つの時系列画像から画
素ごとに移動する物体から速度ベクトルを推定する速度
推定部１４０と、推定された速度ベクトルを出力する出
力部１５０で構成されている。
【００１５】Ｈ＆Ｓ法では、ある画素値において、すな
わち画像中のある位置において対象の微小な時空間的な
ずれがあった場合は、その位置での画像の濃淡値が近似
的に等しいということを前提条件にテイラー展開を通じ
てその基本式（エネルギー誤差汎関数）が式（１）のよ
うに導かれる。
わち画像中のある位置において対象の微小な時空間的な
ずれがあった場合は、その位置での画像の濃淡値が近似
的に等しいということを前提条件にテイラー展開を通じ
てその基本式（エネルギー誤差汎関数）が式（１）のよ
うに導かれる。
【００１６】
【数１】
式（１）中、右辺第一項はペナルティー汎関数、第二項
は安定化汎関数からなる。αは正則化係数である。２次
元画像の（ｘ，ｙ）における濃淡値をＩ（ｘ，ｙ）、推
定する速度場を（ｕ，ｖ）とする。各変数の添え字は
ｘ，ｙ方向での１次微分である。
は安定化汎関数からなる。αは正則化係数である。２次
元画像の（ｘ，ｙ）における濃淡値をＩ（ｘ，ｙ）、推
定する速度場を（ｕ，ｖ）とする。各変数の添え字は
ｘ，ｙ方向での１次微分である。
【００１７】この基本式について、対象領域Ｓについて
の積分をとり、全エネルギーが最小になった時の解が推
定すべき場のオプティカルフローとなる。フローについ
て２つの成分として解が得られるように、Ｈ＆Ｓ法では
拘束条件として、速度の１次微分値が限りなくゼロであ
るという滑らか拘束が与えられている。
の積分をとり、全エネルギーが最小になった時の解が推
定すべき場のオプティカルフローとなる。フローについ
て２つの成分として解が得られるように、Ｈ＆Ｓ法では
拘束条件として、速度の１次微分値が限りなくゼロであ
るという滑らか拘束が与えられている。
【００１８】このエネルギー式を最小化する枠組みは変
分法を介して、∂Ｅ／∂χ＝０，∂Ｅ／∂ｙ＝０とおく
ことでＥＬ方程式（式（２））を得、この方程式を解を
得る。
分法を介して、∂Ｅ／∂χ＝０，∂Ｅ／∂ｙ＝０とおく
ことでＥＬ方程式（式（２））を得、この方程式を解を
得る。
【００１９】
【数２】
次に、式（２）を等間隔の格子点（画素）について差分
法に従って離散化近似する。
法に従って離散化近似する。
【００２０】
【数３】
さらに、式（２）から式（３）のように速度場（ｕ，
ｖ）に関しての連立１次方程式に陽解法形式に整理した
上で、反復的緩和法（Ｇａｕｓｓ−Ｓｅｉｄｅｌ法）に
より、一定の反復誤差基準の下で解を得ることができ
る。式（４）はそれぞれ、濃淡値Ｉ（ｘ，ｙ）の時間に
ついての前進差分と、ｘ，ｙ成分についての空間的な前
進差分である。
ｖ）に関しての連立１次方程式に陽解法形式に整理した
上で、反復的緩和法（Ｇａｕｓｓ−Ｓｅｉｄｅｌ法）に
より、一定の反復誤差基準の下で解を得ることができ
る。式（４）はそれぞれ、濃淡値Ｉ（ｘ，ｙ）の時間に
ついての前進差分と、ｘ，ｙ成分についての空間的な前
進差分である。
【００２１】
【数４】
式（３）中、ｐは反復回数である。なお、ここで用いて
いる定式化はほとんどすべて原著［１］に従っている。
△ｘ，△ｙは単位長である。
いる定式化はほとんどすべて原著［１］に従っている。
△ｘ，△ｙは単位長である。
【００２２】画像の局所的な不連続性の違いにより、最
適正規化係数αも局所的に変化させることが適当である
ので、式（１）と（３）に含まれているαをαij（式
（５））のように、各画素毎に与えられるようにしてい
る。
適正規化係数αも局所的に変化させることが適当である
ので、式（１）と（３）に含まれているαをαij（式
（５））のように、各画素毎に与えられるようにしてい
る。
【００２３】
【数５】
文献［１］では最適正規化係数αは画像全体で同じ値が
設定されているため、最適正規化係数αはスカラー量と
して使われているが、本発明ではベクトル量とする。
設定されているため、最適正規化係数αはスカラー量と
して使われているが、本発明ではベクトル量とする。
【００２４】必要に応じて、式（３）を緩和法で解く
際、収束を早めることができる。これには、加速緩和
（ＳＯＲ）法が知られており、反復時の各格子点ごとの
反復誤差を反復回数が１つ前の格子点の値に適当な重み
づけ量とともに加算すればよい。重みづけ量について
は、一般に１．５〜２．０の間の値が選択される。多く
の問題において、１．７，１．８の値が選択されてい
る。ＳＯＲ法により、収束するまでの反復回数が１／２
〜１／１００になる。
際、収束を早めることができる。これには、加速緩和
（ＳＯＲ）法が知られており、反復時の各格子点ごとの
反復誤差を反復回数が１つ前の格子点の値に適当な重み
づけ量とともに加算すればよい。重みづけ量について
は、一般に１．５〜２．０の間の値が選択される。多く
の問題において、１．７，１．８の値が選択されてい
る。ＳＯＲ法により、収束するまでの反復回数が１／２
〜１／１００になる。
【００２５】本実施形態では、空間周波数と濃淡値の標
準偏差の画像特徴量とＨ＆Ｓ法における最適正規化係数
αとの関係について解析するために、正弦波を用いて、
その振幅（ａｍｐ）・空間周波数（ｆｒｅｑ）をさまざ
まに変化させたパターンを用いたシミュレーションを行
っている。式（６）に示すように、平均濃淡値レベルを
１２７として、振幅については、１０〜１１０の間で３
０階調ずつ、空間周波数については１．０ピクセル／フ
レーム、波長（λ）については、６４，４８，３２，１
６，８，４ピクセルを選択している。
準偏差の画像特徴量とＨ＆Ｓ法における最適正規化係数
αとの関係について解析するために、正弦波を用いて、
その振幅（ａｍｐ）・空間周波数（ｆｒｅｑ）をさまざ
まに変化させたパターンを用いたシミュレーションを行
っている。式（６）に示すように、平均濃淡値レベルを
１２７として、振幅については、１０〜１１０の間で３
０階調ずつ、空間周波数については１．０ピクセル／フ
レーム、波長（λ）については、６４，４８，３２，１
６，８，４ピクセルを選択している。
【００２６】
【数６】
これらのパラメータによる組み合わせ総数は約数万であ
り、同回数だけＨ＆Ｓ法が適用される。なお、空間周波
数は速度と波長との比で定められる。また、基本波（２
π）は波長６４ピクセルとしている。この大きさは後述
するように、実画像でのサブブロックの大きさと同じく
している。したがって、例えば、波数ｋ＝１の時、波長
は６４ピクセル、ｋ＝４の時は、λ＝６４／４＝１６ピ
クセルとなる。空間周波数はｖｅｌ＝１．０ピクセル／
フレーム、ｋ＝１の時、ｆｒｅｑ＝１／６４＝０．０１
５６サークル／ピクセル、Ｋ＝４の時、ｆｒｅｑ＝１／
１６＝０．０６２５サークル／ピクセルである。ここで
行うシミュレーションでは、正弦波の移動方向は水平方
向、画面に向かって右方である。
り、同回数だけＨ＆Ｓ法が適用される。なお、空間周波
数は速度と波長との比で定められる。また、基本波（２
π）は波長６４ピクセルとしている。この大きさは後述
するように、実画像でのサブブロックの大きさと同じく
している。したがって、例えば、波数ｋ＝１の時、波長
は６４ピクセル、ｋ＝４の時は、λ＝６４／４＝１６ピ
クセルとなる。空間周波数はｖｅｌ＝１．０ピクセル／
フレーム、ｋ＝１の時、ｆｒｅｑ＝１／６４＝０．０１
５６サークル／ピクセル、Ｋ＝４の時、ｆｒｅｑ＝１／
１６＝０．０６２５サークル／ピクセルである。ここで
行うシミュレーションでは、正弦波の移動方向は水平方
向、画面に向かって右方である。
【００２７】図２に、数値解析で用いた正弦波パターン
の例を示す。１フレームの大きさが６４ｘ６４ピクセ
ル、振幅ａｍｐ＝１１０、空間周波数ｖｅｌ＝１．０、
波長λ＝４〜６４ピクセルの場合で、連続した２フレー
ムずつ組にして示す。なお、このフレーム内での振幅値
から濃淡値の標準偏差（ＳＤ）への換算式については、
図３に示す。画像における濃淡値の標準偏差（Standard
Deviation:SD) を式（７）に基づいて計算する。
の例を示す。１フレームの大きさが６４ｘ６４ピクセ
ル、振幅ａｍｐ＝１１０、空間周波数ｖｅｌ＝１．０、
波長λ＝４〜６４ピクセルの場合で、連続した２フレー
ムずつ組にして示す。なお、このフレーム内での振幅値
から濃淡値の標準偏差（ＳＤ）への換算式については、
図３に示す。画像における濃淡値の標準偏差（Standard
Deviation:SD) を式（７）に基づいて計算する。
【００２８】
【数７】
なお、実際の画像に適用するときは、画像を適当な大き
さの幾つかの小ブロックに分割して式（７）を適用す
る。式（７）中、ｓｕｂはサブブロック番号を示し、ａ
ｖｅは濃淡値の平均値を示す。ここでは、濃淡値のばら
つきについては分散を用いることと等価であるが、物理
的に直観性のある標準偏差ＳＤを適用する。式（３）を
解いて推定される速度場については、対象の動きベクト
ルが既知の場合は、式（８）を用いて定量的に評価する
ことができる。
さの幾つかの小ブロックに分割して式（７）を適用す
る。式（７）中、ｓｕｂはサブブロック番号を示し、ａ
ｖｅは濃淡値の平均値を示す。ここでは、濃淡値のばら
つきについては分散を用いることと等価であるが、物理
的に直観性のある標準偏差ＳＤを適用する。式（３）を
解いて推定される速度場については、対象の動きベクト
ルが既知の場合は、式（８）を用いて定量的に評価する
ことができる。
【００２９】
【数８】
ここでは、既知の速度ベクトルとＨ＆Ｓ法により推定さ
れる速度ベクトルとの逆余弦関数の対象領域での全格子
点による平均（全画素数ＮｘＮ ピクセル）をとってい
る。ただし、速度ベクトル
れる速度ベクトルとの逆余弦関数の対象領域での全格子
点による平均（全画素数ＮｘＮ ピクセル）をとってい
る。ただし、速度ベクトル
【００３０】
【外１】
ともに単位ベクトルである。
【００３１】式（６）に基づいて連続性の異なる時系列
パターンを生成して、既知の速度ベクトルとＨ＆Ｓ法を
適用して推定された速度ベクトルを式（８）（Ｎ＝６
４）により定量的な評価を行う。そして、最も誤差が少
なくなる時の振幅、空間周波数の組み合わせとその時の
最適正規化係数αを求めている。この結果により、推定
誤差を最小にするという規範で最適正規化係数空間を形
成できる。なお、Ｈ＆Ｓ法における反復回数は、各反復
時の速度ベクトルで最大の大きさのベクトルの反復誤差
が０．００１以下になった時点で収束したものと見なし
ている。この空間を実画像に適用する際、本実施形態で
は画像を幾つかの小ブロック（６４ｘ６４）に分割し
て、各小ブロックごとに空間周波数分布と濃淡値分布の
標準偏差（Ｎ＝６４）を同時に解析する。この２つの分
布に基づいて、最適正規化係数空間を参照することで、
各小ブロック独立に最適正規化係数αを推定する。オプ
ティカルフローはＨ＆Ｓ法に、局所的に異なるαを代入
して緩和法により得ることができる。
パターンを生成して、既知の速度ベクトルとＨ＆Ｓ法を
適用して推定された速度ベクトルを式（８）（Ｎ＝６
４）により定量的な評価を行う。そして、最も誤差が少
なくなる時の振幅、空間周波数の組み合わせとその時の
最適正規化係数αを求めている。この結果により、推定
誤差を最小にするという規範で最適正規化係数空間を形
成できる。なお、Ｈ＆Ｓ法における反復回数は、各反復
時の速度ベクトルで最大の大きさのベクトルの反復誤差
が０．００１以下になった時点で収束したものと見なし
ている。この空間を実画像に適用する際、本実施形態で
は画像を幾つかの小ブロック（６４ｘ６４）に分割し
て、各小ブロックごとに空間周波数分布と濃淡値分布の
標準偏差（Ｎ＝６４）を同時に解析する。この２つの分
布に基づいて、最適正規化係数空間を参照することで、
各小ブロック独立に最適正規化係数αを推定する。オプ
ティカルフローはＨ＆Ｓ法に、局所的に異なるαを代入
して緩和法により得ることができる。
【００３２】図４に、各最適正規化係数αに対してＨ＆
Ｓにより推定されたフロー（移動ベクトル（動きベクト
ル））の例を示す。入力した正弦波は、振幅ａｍｐ＝３
０、空間周波数ｖｅｌ＝１、波長λ＝１６の場合であ
る。推定誤差が最小となったのはα＝１６．０の時であ
り、この時のαを最適値と見なす。比較のために、非最
適値αの例を示す。αがα＝１５０．０からα＝１９
０．０まで増加すると、フローの密度がしだいに低く、
また、同時に、実際よりも速度が小さくなっていき、最
後には消えてしまうことがわかる。ここには示していな
いが、αを小さくしていくと、フローの密度は高くなる
が、実際よりも大きい速度のフローが多くなる傾向にあ
る。
Ｓにより推定されたフロー（移動ベクトル（動きベクト
ル））の例を示す。入力した正弦波は、振幅ａｍｐ＝３
０、空間周波数ｖｅｌ＝１、波長λ＝１６の場合であ
る。推定誤差が最小となったのはα＝１６．０の時であ
り、この時のαを最適値と見なす。比較のために、非最
適値αの例を示す。αがα＝１５０．０からα＝１９
０．０まで増加すると、フローの密度がしだいに低く、
また、同時に、実際よりも速度が小さくなっていき、最
後には消えてしまうことがわかる。ここには示していな
いが、αを小さくしていくと、フローの密度は高くなる
が、実際よりも大きい速度のフローが多くなる傾向にあ
る。
【００３３】図５に、式（６）におけるｖｅｌ＝１．
０、ａｍｐ＝１０，３０，１１０を入力し、式（７）で
推定誤差を評価した結果を示す。また、ａｍｐ＝１０と
０．０１〜２０．０のα（横軸）を拡大した図を示す。
ただし、予備検討の結果から、αの範囲によりαの刻み
幅を変えており、０．０１〜２０．０では０．０１毎、
２０．０〜７０．０では１．０毎、７０．０〜１５０．
０では２．０毎、１５０．０〜３００．０では６．０毎
としている。ここでは、最も不連続性が小さいものはａ
ｍｐ＝１０、最も大きいものはａｍｐ＝１１０である。
図５において、どの場合も、推定誤差（縦軸）に対して
複数の極小値が存在することなく、最小値のみが存在し
ているのがわかる。ａｍｐ＝１０からａｍｐ＝１１０に
なるに従って、推定誤差を最小とするα値はしだいに、
２０．０未満から１００．０前後まで大きくなってい
る。また、どのａｍｐにおいて、αの増加とともに推定
誤差が飽和しているのがわかる。特に、ａｍｐが小さい
時ほど、推定誤差を最小とするαの選択に関する感度は
高くなっている。各ａｍｐ中、各曲線は適用した１０種
類の空間周波数に対する結果である。
０、ａｍｐ＝１０，３０，１１０を入力し、式（７）で
推定誤差を評価した結果を示す。また、ａｍｐ＝１０と
０．０１〜２０．０のα（横軸）を拡大した図を示す。
ただし、予備検討の結果から、αの範囲によりαの刻み
幅を変えており、０．０１〜２０．０では０．０１毎、
２０．０〜７０．０では１．０毎、７０．０〜１５０．
０では２．０毎、１５０．０〜３００．０では６．０毎
としている。ここでは、最も不連続性が小さいものはａ
ｍｐ＝１０、最も大きいものはａｍｐ＝１１０である。
図５において、どの場合も、推定誤差（縦軸）に対して
複数の極小値が存在することなく、最小値のみが存在し
ているのがわかる。ａｍｐ＝１０からａｍｐ＝１１０に
なるに従って、推定誤差を最小とするα値はしだいに、
２０．０未満から１００．０前後まで大きくなってい
る。また、どのａｍｐにおいて、αの増加とともに推定
誤差が飽和しているのがわかる。特に、ａｍｐが小さい
時ほど、推定誤差を最小とするαの選択に関する感度は
高くなっている。各ａｍｐ中、各曲線は適用した１０種
類の空間周波数に対する結果である。
【００３４】図５のようにして得られた推定誤差を最小
とするαとの関係について、図６のようにまとめること
ができる。図６より、各ａｍｐにおいて、空間周波数
（ｆｒｅｑ）が０．２を境として、最適αの分布が異な
っているのがわかる。また、ａｍｐが増加すると、最適
αのオーダーも増加している。さらに、ｆｒｅｑが０．
２までは、ｆｒｅｑが高くなるに従って、最適αも増加
している。一方、ｆｒｅｑが０．２では逆に最適αは減
少している。これについては、式（６）により生成した
正弦波が、空間周波数を高くすると、分解能の低さから
十分に滑らかな濃淡値変化をもつ正弦波を生成できてお
らず、みかけ上、空間周波数が低いパターンが入力され
たことによると考えられる。
とするαとの関係について、図６のようにまとめること
ができる。図６より、各ａｍｐにおいて、空間周波数
（ｆｒｅｑ）が０．２を境として、最適αの分布が異な
っているのがわかる。また、ａｍｐが増加すると、最適
αのオーダーも増加している。さらに、ｆｒｅｑが０．
２までは、ｆｒｅｑが高くなるに従って、最適αも増加
している。一方、ｆｒｅｑが０．２では逆に最適αは減
少している。これについては、式（６）により生成した
正弦波が、空間周波数を高くすると、分解能の低さから
十分に滑らかな濃淡値変化をもつ正弦波を生成できてお
らず、みかけ上、空間周波数が低いパターンが入力され
たことによると考えられる。
【００３５】以上の結果から、パターンの不連続性が大
きい場合、即ち、振幅（ａｍｐ）が大きいか、空間周波
数が高い（０．２まで）程、推定誤差を最小とするため
には、αを大きくする必要があることが示唆された。
きい場合、即ち、振幅（ａｍｐ）が大きいか、空間周波
数が高い（０．２まで）程、推定誤差を最小とするため
には、αを大きくする必要があることが示唆された。
【００３６】図７は図５，６より、各波数毎の、濃淡値
の標準偏差と最適αの線形１次方程式を示す。実際の応
用では、画像の空間周波数を推定して、これを波数に換
算し、同時に、濃淡の標準偏差を推定することで、最適
正則化係数αを決定することができる。その結果、Ｈ＆
Ｓ法で推定精度の良い速度ベクトルを得ることができ
る。
の標準偏差と最適αの線形１次方程式を示す。実際の応
用では、画像の空間周波数を推定して、これを波数に換
算し、同時に、濃淡の標準偏差を推定することで、最適
正則化係数αを決定することができる。その結果、Ｈ＆
Ｓ法で推定精度の良い速度ベクトルを得ることができ
る。
【００３７】上述したように、正弦波パターンのａｍ
ｐ、濃淡値のＳＤをさまざまに変化させて不連続性の時
空間パターンを生成し、Ｈ＆Ｓ法を用いながら、推定誤
差が最小となる規範で最適αを導出している。正弦波パ
ターンでは速度が１ピクセル／フレームの場合を用い
た。実際の画像ではさまざまな速度をもつ物体が存在す
るために、正弦波パターンでも予めさまざまな速度の場
合についての最適α空間が設計されるべきである。しか
しながら、前に論じたように、Ｈ＆Ｓ法の差分計算では
せいぜい１ピクセル離れた格子点間での濃淡値が使わ
れ、この形式は物体の移動速度に依存せずに一定であ
る。即ち、物体の速度の違いはある格子点上での濃淡値
の差分量のみに依存するので、濃淡値の標準偏差に基づ
いた図７は他の１ピクセル／フレーム以外の対象に対し
ても有効であると考える。
ｐ、濃淡値のＳＤをさまざまに変化させて不連続性の時
空間パターンを生成し、Ｈ＆Ｓ法を用いながら、推定誤
差が最小となる規範で最適αを導出している。正弦波パ
ターンでは速度が１ピクセル／フレームの場合を用い
た。実際の画像ではさまざまな速度をもつ物体が存在す
るために、正弦波パターンでも予めさまざまな速度の場
合についての最適α空間が設計されるべきである。しか
しながら、前に論じたように、Ｈ＆Ｓ法の差分計算では
せいぜい１ピクセル離れた格子点間での濃淡値が使わ
れ、この形式は物体の移動速度に依存せずに一定であ
る。即ち、物体の速度の違いはある格子点上での濃淡値
の差分量のみに依存するので、濃淡値の標準偏差に基づ
いた図７は他の１ピクセル／フレーム以外の対象に対し
ても有効であると考える。
【００３８】画像は図８に示すように、テクスチャー解
析に広く用いられているＢｒｏｄａｔｚ（文献［２］：
P.Brodatz, "Textures: a photographic album for art
ists& designers", New York, Dover, 1966）テクスチ
ャー集の中から、さまざまな空間周波数と濃淡値をも
つ、＃９（Grass Lawn），＃７８（Oriental Straw Clo
th），＃８４（Raffia），＃９３（Fur）が選択されて
いる。４種類の画像はいずれも１２８ｘ１２８ピクセル
であり、１枚の画像（２５６ｘ２５６ピクセル）に合成
している。
析に広く用いられているＢｒｏｄａｔｚ（文献［２］：
P.Brodatz, "Textures: a photographic album for art
ists& designers", New York, Dover, 1966）テクスチ
ャー集の中から、さまざまな空間周波数と濃淡値をも
つ、＃９（Grass Lawn），＃７８（Oriental Straw Clo
th），＃８４（Raffia），＃９３（Fur）が選択されて
いる。４種類の画像はいずれも１２８ｘ１２８ピクセル
であり、１枚の画像（２５６ｘ２５６ピクセル）に合成
している。
【００３９】画像は、αを局所的に推定するために、６
４ｘ６４ピクセルのサブブロックに分割される。それぞ
れのサブブロック毎に２次元周波数のエネルギー分布を
高速フーリエ変換（fast Fourier Transform:fFT）を用
いて求める。その際、前処理として、直流成分の除去
と、Ｈａｎｎｉｎｇ窓による画像の切り取り・平滑化を
行っている。なお、サブブロックによっては画像をほと
んど含まない領域が存在するが、ＦＴは位置に不変な変
換なのでエネルギーの分布への影響は生じない。図９は
各サブブロック毎のエネルギー分布であり、各サブブロ
ックの中心部に低空間周波数成分があり、周辺にいくほ
ど、高空間周波数成分が分布する。また、明細書では考
慮しないがパターンの方向性が強い場合は、実際の方向
と直角方向にエネルギーが分布し、方向性があまり顕著
でない場合は同心円状に分布する。また、フーリエ変換
の性質から、エネルギー分布は中心に対して点対称であ
る。
４ｘ６４ピクセルのサブブロックに分割される。それぞ
れのサブブロック毎に２次元周波数のエネルギー分布を
高速フーリエ変換（fast Fourier Transform:fFT）を用
いて求める。その際、前処理として、直流成分の除去
と、Ｈａｎｎｉｎｇ窓による画像の切り取り・平滑化を
行っている。なお、サブブロックによっては画像をほと
んど含まない領域が存在するが、ＦＴは位置に不変な変
換なのでエネルギーの分布への影響は生じない。図９は
各サブブロック毎のエネルギー分布であり、各サブブロ
ックの中心部に低空間周波数成分があり、周辺にいくほ
ど、高空間周波数成分が分布する。また、明細書では考
慮しないがパターンの方向性が強い場合は、実際の方向
と直角方向にエネルギーが分布し、方向性があまり顕著
でない場合は同心円状に分布する。また、フーリエ変換
の性質から、エネルギー分布は中心に対して点対称であ
る。
【００４０】本実施形態では図９の結果において、濃淡
値と２次位置の３次元量を扱うのではなく、中心から同
心円状にエネルギーの和をとることで各周波数ごとのエ
ネルギーの分布をとらえる。
値と２次位置の３次元量を扱うのではなく、中心から同
心円状にエネルギーの和をとることで各周波数ごとのエ
ネルギーの分布をとらえる。
【００４１】エネルギー画像（フーリエ変換された画像
についてエネルギーを求めた画像で、フーリエ画像の実
部と虚部を２乗してそれぞれ和の２乗をとったもの）で
は各画素に各周波数・方向におけるエネルギーが分布し
ている。各周波数毎のエネルギーの算出式を式（９）に
示す。なお、中心に直流成分、中心から外側に向かって
低空間周波数成分〜高空間周波数成分をもつ画素が並ん
でいるものとする。中心から同一半径上に同一周波数成
分が並んでいる。
についてエネルギーを求めた画像で、フーリエ画像の実
部と虚部を２乗してそれぞれ和の２乗をとったもの）で
は各画素に各周波数・方向におけるエネルギーが分布し
ている。各周波数毎のエネルギーの算出式を式（９）に
示す。なお、中心に直流成分、中心から外側に向かって
低空間周波数成分〜高空間周波数成分をもつ画素が並ん
でいるものとする。中心から同一半径上に同一周波数成
分が並んでいる。
【００４２】
【数９】
図９のように２次元的にエネルギー分布を視覚的に比べ
ることができるが、計算機上で自動的に空間周波数を推
定する必要がある。
ることができるが、計算機上で自動的に空間周波数を推
定する必要がある。
【００４３】空間周波数領域において、周期性がある場
合はそれに対応する周波数にエネルギーのピークが立
つ。各サブブロック内で同一の帯域をもつ対象の面積が
広いほど、エネルギーは高くなる。即ち、エネルギーの
集中する帯域から代表空間周波数を推定すればよく、サ
ブブロック内の全エネルギーに対する各周波数における
エネルギーが画像特徴としての重要度を表わすことにな
る。
合はそれに対応する周波数にエネルギーのピークが立
つ。各サブブロック内で同一の帯域をもつ対象の面積が
広いほど、エネルギーは高くなる。即ち、エネルギーの
集中する帯域から代表空間周波数を推定すればよく、サ
ブブロック内の全エネルギーに対する各周波数における
エネルギーが画像特徴としての重要度を表わすことにな
る。
【００４４】このような視点から、次に述べるような手
順で図１０から代表周波数（ｗｆ）を決定する。 １）まず、低空間周波数帯域だけに直線を当てはめて、
水平軸と交わった点を低空間周波数ｆ１とする。ｆ１に
対応するエネルギーｐ１を求める。 ２）残りの帯域において、式（１０）に示すように、各
リングごとに（エネルギー画像をその中心から同心円上
に抜き取り）残りの帯域のエネルギー和を求めて（分
母）、各リング毎に各エネルギーで重み付け（分子）し
た結果との比をとる。ただし、実際に適用する際、リン
グ番号は３０までを使っている。
順で図１０から代表周波数（ｗｆ）を決定する。 １）まず、低空間周波数帯域だけに直線を当てはめて、
水平軸と交わった点を低空間周波数ｆ１とする。ｆ１に
対応するエネルギーｐ１を求める。 ２）残りの帯域において、式（１０）に示すように、各
リングごとに（エネルギー画像をその中心から同心円上
に抜き取り）残りの帯域のエネルギー和を求めて（分
母）、各リング毎に各エネルギーで重み付け（分子）し
た結果との比をとる。ただし、実際に適用する際、リン
グ番号は３０までを使っている。
【００４５】
【数１０】
３）以上により代表周波数ｆｗが推定される。
上記の方法を重み付周波数推定法と呼ぶ。
【００４６】この方法による利点としては、しきい値を
おくことなく、エネルギーの大きさ、即ち、αを左右す
る重要な画像特徴量（空間周波数）により代表周波数ｗ
ｆが自動的に可変することにある。
おくことなく、エネルギーの大きさ、即ち、αを左右す
る重要な画像特徴量（空間周波数）により代表周波数ｗ
ｆが自動的に可変することにある。
【００４７】ここでは最適正則化α空間の有効性を調べ
るために、実画像において、αに画像の全領域で人為的
に一つの値を与えた場合と、サブブロック毎に局所的に
自動的に推定したαを与えた場合のオプティカルフロー
の結果を比較する（図１１）。Ｈ＆Ｓ法に連続した２フ
レームを入力する際、平滑化のために４点の移動平均に
よるローパスをかける。Ｈ＆Ｓ法における反復回数は、
各反復時の速度ベクトルで最大の大きさのベクトルの反
復誤差が０．００１以下になった時点で収束したものと
見なしている。なお、図１１は、結果の表示は見やすく
するために、実際に推定されたベクトルの大きさの８倍
に、間隔は７ピクセル毎に表示している。
るために、実画像において、αに画像の全領域で人為的
に一つの値を与えた場合と、サブブロック毎に局所的に
自動的に推定したαを与えた場合のオプティカルフロー
の結果を比較する（図１１）。Ｈ＆Ｓ法に連続した２フ
レームを入力する際、平滑化のために４点の移動平均に
よるローパスをかける。Ｈ＆Ｓ法における反復回数は、
各反復時の速度ベクトルで最大の大きさのベクトルの反
復誤差が０．００１以下になった時点で収束したものと
見なしている。なお、図１１は、結果の表示は見やすく
するために、実際に推定されたベクトルの大きさの８倍
に、間隔は７ピクセル毎に表示している。
【００４８】表１に、各サブブロック毎の重みづけされ
た空間周波数、表２に濃淡値の標準偏差、表３に推定さ
れた最適αをまとめて示す。αは５８．９から１４０．
８まで幅広く推定されているのがわかる。なお、表１中
の波数は、後述するように、帯域補正がすべてのサブブ
ロックに１だけ加えられている。さらに、サブブロック
毎に異なるαが推定されるが、画像全体についてαの値
を４点移動平均で平滑することで収束が多少速くなるこ
とを確認している。なお、表１〜３中のｒ１〜ｒ１５は
図８中のサブブロックの番号を示している。
た空間周波数、表２に濃淡値の標準偏差、表３に推定さ
れた最適αをまとめて示す。αは５８．９から１４０．
８まで幅広く推定されているのがわかる。なお、表１中
の波数は、後述するように、帯域補正がすべてのサブブ
ロックに１だけ加えられている。さらに、サブブロック
毎に異なるαが推定されるが、画像全体についてαの値
を４点移動平均で平滑することで収束が多少速くなるこ
とを確認している。なお、表１〜３中のｒ１〜ｒ１５は
図８中のサブブロックの番号を示している。
【００４９】
【表１】
【００５０】
【表２】
【００５１】
【表３】
【００５２】図１１はα（＝０．１，１０．０，５０．
０）を一様に与えた場合、図８に示すテクスチャーに対
して推定されたフローである。αが小さい時は速度ベク
トルは大きく、また、その方向もさまざまにばらついて
いる。αを増加させるとフローはしだいに整ってくる
が、その大きさは過小になる傾向にある。
０）を一様に与えた場合、図８に示すテクスチャーに対
して推定されたフローである。αが小さい時は速度ベク
トルは大きく、また、その方向もさまざまにばらついて
いる。αを増加させるとフローはしだいに整ってくる
が、その大きさは過小になる傾向にある。
【００５３】本発明の有効性を定量的に評価するため
に、正弦波パターン実験の場合と同様に、式（７）で誤
差評価を行いながら、αを全領域で一様に、０．１から
３００．０まで、１０．０ずつ変化させてフローを画素
単位に推定誤差を求めた。
に、正弦波パターン実験の場合と同様に、式（７）で誤
差評価を行いながら、αを全領域で一様に、０．１から
３００．０まで、１０．０ずつ変化させてフローを画素
単位に推定誤差を求めた。
【００５４】図１２はαを一様に与えた場合と最適α空
間により推定されたαに基づくフローの画素ごとの平均
誤差値を比較している。図から誤差が最小となったのは
αを一様に与えた場合、推定された場合、それぞれの平
均誤差値は０．０９７４ｒａｄ／ピクセル（α＝８０．
０）、０．０６５８ｒａｄ／ピクセル（α＝５８．９〜
１４０．８）である。これから平均誤差が４７．９％向
上しているのがわかる。さらに、一様なαと最適αの場
合で、平均速度ベクトルの大きさは０．９１３ピクセル
／フレーム、速度成分は（０．００２４４５，０．９３
３６３４）と０．９５３ピクセル／フレーム、（−０．
０００３９，０．９５２８５）であり、推定されたαを
適用した時の方が実際の動き（１ピクセル／フレーム）
との誤差は、４．７％低い程度に留まった。
間により推定されたαに基づくフローの画素ごとの平均
誤差値を比較している。図から誤差が最小となったのは
αを一様に与えた場合、推定された場合、それぞれの平
均誤差値は０．０９７４ｒａｄ／ピクセル（α＝８０．
０）、０．０６５８ｒａｄ／ピクセル（α＝５８．９〜
１４０．８）である。これから平均誤差が４７．９％向
上しているのがわかる。さらに、一様なαと最適αの場
合で、平均速度ベクトルの大きさは０．９１３ピクセル
／フレーム、速度成分は（０．００２４４５，０．９３
３６３４）と０．９５３ピクセル／フレーム、（−０．
０００３９，０．９５２８５）であり、推定されたαを
適用した時の方が実際の動き（１ピクセル／フレーム）
との誤差は、４．７％低い程度に留まった。
【００５５】最適α空間を図８に示すテクスチャー画像
を並進させた場合に適用した場合、式（３）で要した反
復回数（ｐ）は３２６回であった。αを一様に与えた場
合の収束までの反復回数については、αが３００．０で
１０７１回となった。即ち、αの増加とともに収束に時
間がかかることになる。スケールスペースを用いても大
きいαから最適解を探索する枠組みではトータルの演算
コストがかかる方法と言える。
を並進させた場合に適用した場合、式（３）で要した反
復回数（ｐ）は３２６回であった。αを一様に与えた場
合の収束までの反復回数については、αが３００．０で
１０７１回となった。即ち、αの増加とともに収束に時
間がかかることになる。スケールスペースを用いても大
きいαから最適解を探索する枠組みではトータルの演算
コストがかかる方法と言える。
【００５６】図１３を参照すると、本発明の他の実施態
様の画像流速推定装置は、それぞれ図１中の画像入力部
１００、画像蓄積部１１０、出力部１５０に相当する画
像入力部２１０、画像蓄積部２２０、出力部２３０と、
画像入力部２１０から入力された画像の画像蓄積部２２
０への蓄積から推定された速度ベクトルを出力部２３０
へ出力するまでの各処理からなる画像流速推定プログラ
ムを記録した、ＦＤ、ＣＤ−ＲＯＭ、半導体メモリ等の
記録媒体２４０と、記録媒体２４０から画像流速推定プ
ログラムを読み込んで実行するデータ処理装置２５０で
構成されている。
様の画像流速推定装置は、それぞれ図１中の画像入力部
１００、画像蓄積部１１０、出力部１５０に相当する画
像入力部２１０、画像蓄積部２２０、出力部２３０と、
画像入力部２１０から入力された画像の画像蓄積部２２
０への蓄積から推定された速度ベクトルを出力部２３０
へ出力するまでの各処理からなる画像流速推定プログラ
ムを記録した、ＦＤ、ＣＤ−ＲＯＭ、半導体メモリ等の
記録媒体２４０と、記録媒体２４０から画像流速推定プ
ログラムを読み込んで実行するデータ処理装置２５０で
構成されている。
【００５７】次に、本発明の具体例について説明する。
【００５８】（流体の場合）
例１）観測される降水情報をもった流体状のパターンの
動きを気象レーダにより、本発明に基づいて検出する。
図１４に２つの例を示す。生成と消滅が激しいパターン
であるのと、局所的に動きの方向と大きさがさまざまで
あることが特徴である。また、さまざまな空間周波数と
平均濃淡値をもつ。この濃淡値の輝度は降水量の多さに
比例している。このようなパターンに対して、本発明の
方法により検出検出したフローの結果（図１５）は実際
と比べると視覚的にみて、かなり良好な精度でフローが
推定されている。最低αは画像を１６分割して、予め設
計しておいた最適α空間を参照することにより、１６個
のαを求めている。このようにして得られたフローに基
づいて、降水パターンの移動方向を予測できる。 例２）流体工学では管内の流れを解析する研究が盛んに
行われている。その中で、流体中に、流れの可視化のた
めに、粉塵や墨をいれて、その移動方向をみる実験が行
われている。粉塵や墨の動きを画像でとらえた後、本発
明により流れを検出すると、図１６のようになる。最適
αは画像を８分割して、予め設計しておいた最適α空間
を参照することにより、８個のαを求めている。例で
は、管の分岐付近に澱みが小さな渦となっているのがわ
かる。渦のある所と無い早い流れ付近では、空間周波数
と平均濃淡値が異なる。
動きを気象レーダにより、本発明に基づいて検出する。
図１４に２つの例を示す。生成と消滅が激しいパターン
であるのと、局所的に動きの方向と大きさがさまざまで
あることが特徴である。また、さまざまな空間周波数と
平均濃淡値をもつ。この濃淡値の輝度は降水量の多さに
比例している。このようなパターンに対して、本発明の
方法により検出検出したフローの結果（図１５）は実際
と比べると視覚的にみて、かなり良好な精度でフローが
推定されている。最低αは画像を１６分割して、予め設
計しておいた最適α空間を参照することにより、１６個
のαを求めている。このようにして得られたフローに基
づいて、降水パターンの移動方向を予測できる。 例２）流体工学では管内の流れを解析する研究が盛んに
行われている。その中で、流体中に、流れの可視化のた
めに、粉塵や墨をいれて、その移動方向をみる実験が行
われている。粉塵や墨の動きを画像でとらえた後、本発
明により流れを検出すると、図１６のようになる。最適
αは画像を８分割して、予め設計しておいた最適α空間
を参照することにより、８個のαを求めている。例で
は、管の分岐付近に澱みが小さな渦となっているのがわ
かる。渦のある所と無い早い流れ付近では、空間周波数
と平均濃淡値が異なる。
【００５９】（剛体の場合）
例１）車が移動している。各車の動きが本発明の方法に
より得られている（図１７）。最適αは１２個である。
より得られている（図１７）。最適αは１２個である。
【００６０】なお、図１８は、αを１０，６０，１２０
と全体で同じ値を与えた場合の結果である。車以外の領
域から多くの誤ったフローが検出されているのがわか
る。 例２）回転する台と玩具（ルービックキューブ）が共に
回転している。回転成分が本発明により得られている
（図１９）。
と全体で同じ値を与えた場合の結果である。車以外の領
域から多くの誤ったフローが検出されているのがわか
る。 例２）回転する台と玩具（ルービックキューブ）が共に
回転している。回転成分が本発明により得られている
（図１９）。
【００６１】
【発明の効果】以上説明したように、本発明によれば、
Ｈｏｒｎ ＆ Ｓｃｈｕｎｃｋの方法に基づいて対象の
移動速度を推定する際、最適な制御量を対象の空間周波
数分布の特徴に従って適応的に決定していくので、安定
した推定精度を得ることができ、また、対象の種類によ
らず、対象の空間周波数と濃淡値の標準偏差という２つ
の基本的な画像特徴量だけを用いるので、本発明の汎用
性は非常に高い。
Ｈｏｒｎ ＆ Ｓｃｈｕｎｃｋの方法に基づいて対象の
移動速度を推定する際、最適な制御量を対象の空間周波
数分布の特徴に従って適応的に決定していくので、安定
した推定精度を得ることができ、また、対象の種類によ
らず、対象の空間周波数と濃淡値の標準偏差という２つ
の基本的な画像特徴量だけを用いるので、本発明の汎用
性は非常に高い。
【図１】本発明の一実施形態の画像流速推定装置の構成
図である。
図である。
【図２】最適正則化係数α空間を設計するための正弦波
パターンを示す図である。
パターンを示す図である。
【図３】正弦波の振幅と濃淡値の標準偏差への換算式を
示すグラフである。
示すグラフである。
【図４】αを人為的に変化させたときの、推定される速
度ベクトルの違いを示す図である。
度ベクトルの違いを示す図である。
【図５】さまざまなαと、さまざまな空間周波数、振幅
に対する速度ベクトルの推定誤差を示すグラフである。
に対する速度ベクトルの推定誤差を示すグラフである。
【図６】図５で誤差が最小となった点を各空間周波数と
振幅についてまとめた結果を示すグラフである。
振幅についてまとめた結果を示すグラフである。
【図７】各波数、濃淡値の標準偏差、最適正則化係数α
との関係を示すグラフである。
との関係を示すグラフである。
【図８】並進運動させた実テクスチャー（画像が１６の
サブブロックに分割されている）を示す図である。
サブブロックに分割されている）を示す図である。
【図９】各サブブロック毎の空間周波数分布を示す図で
ある。
ある。
【図１０】代表空間周波数を決定する方法を説明する図
である。
である。
【図１１】実テクスチャー画像の動きに対する３つのα
を一様に与えたときの推定精度を示す図である。
を一様に与えたときの推定精度を示す図である。
【図１２】自動推定したαと一様にαを与えたときの、
単位画素当りの平均推定誤差を示すグラフである。
単位画素当りの平均推定誤差を示すグラフである。
【図１３】本発明の他の実施形態の画像流速推定装置の
構成図である。
構成図である。
【図１４】観測される降水情報を保った流体状のパター
ンの例を示す図である。
ンの例を示す図である。
【図１５】図１４のパターンに対して検出されたフロー
の結果を示す図である。
の結果を示す図である。
【図１６】流れの中に入れた粉塵や墨の動きを画像でと
らえた図である。
らえた図である。
【図１７】本発明の方法により得られた車の動きを示す
図である。
図である。
【図１８】図１７においてα＝１０，６０，１２０と全
体で同じ値を与えた場合のフローを示す図である。
体で同じ値を与えた場合のフローを示す図である。
【図１９】回転する台と玩具がともに回転する様子と、
本発明により回転成分を検出した結果を示す図である。
本発明により回転成分を検出した結果を示す図である。
１００，２１０ 画像入力部
１１０，２２０ 画像蓄積部
１２０ 画像解析部
１３０ 正規化係数算出部
１４０ 速度推定部
１５０，２３０ 出力部
２４０ 記録媒体
２５０ データ処理装置
フロントページの続き
(56)参考文献 特開 平11−160335（ＪＰ，Ａ)
特開 平９−297851（ＪＰ，Ａ)
(58)調査した分野(Int.Cl.7，ＤＢ名)
G06T 7/00 - 7/60
G06T 1/00
G01P 3/36
H04N 7/18
ＪＩＣＳＴファイル（ＪＯＩＳ)
Claims (11)
- 【請求項１】 時系列画像中の流体の移動速度をオプテ
ィカルフロー法を用いて推定する画像流速推定方法であ
って、 移動する対象を含む２次元画像を入力する画像入力段階
と、入力された画像を時系列画像として蓄積する画像蓄
積段階と、蓄積されている２次元時系列画像の空間周波
数成分を解析するとともに、画像の濃淡値の標準偏差を
同時に解析する画像解析段階と、前記画像解析段階で得
られた画像の空間周波数成分と濃淡値の標準偏差の２つ
の画像特徴量を予め設計してある線形方程式に代入して
最適正則化係数αを決定する正則化係数算出段階と、連
続する２枚の画像から、各最適正則化係数αに対する、
移動する対象の速度ベクトルをオプティカルフロー法を
用いて推定する速度推定段階と、推定された速度ベクト
ルを出力する出力段階を有する画像流速推定方法。 - 【請求項２】 前記画像解析段階は、画像を幾つかのサ
ブブロックに分割して、各サブブロック毎に高速フーリ
エ変換により、エネルギー・水平方向の空間周波数・垂
直方向の空間周波数を軸にもつ３次元の空間周波数画像
に変換し、空間周波数画像の低空間周波数成分が画像中
央に、高空間周波数成分が中央から周辺に向かって分布
するように、高速フーリエ変換後の空間周波数画像の内
部を配置変えし、エネルギー画像をその中心から同心円
状に抜き取り、各円ごとにエネルギーを求めることによ
り前記各円の半径に対応する空間周波数に対するエネル
ギーの分布を得、各空間周波数とそのエネルギーとの積
の和をエネルギーの総和で割った値を算出し、これを代
表空間周波数とする請求項１記載の画像流速推定方法。 - 【請求項３】 前記画像解析段階は、画像を幾つかのサ
ブブロックに分割して、各サブブロック毎に画像の濃淡
値の平均値を算出し、分散値を求めてその平方根をとる
ことで標準偏差を算出する請求項１記載の画像流速推定
方法。 - 【請求項４】 前記正則化係数算出段階は、予めオフラ
インで正弦波の振幅と空間周波数（波数）をさまざまに
変化させながら一定速度で正弦波パターンを移動させる
シミュレーションパターンを予め生成し、オプティカル
フロー法における正則化係数αを同時にさまざまに変化
させながら、シミュレーションパターンから速度ベクト
ルを推定し、ａｒｃｃｏｓｉｎｅ関数を用いて推定され
た速度ベクトルと真の速度ベクトルとの誤差を評価して
誤差が最小となった時のαを選択することで、αと振幅
と空間周波数の３者の関係を定量づけ、さらに振幅から
シミュレーションパターンの濃淡値の標準偏差を算出す
る処理を行って得られる式を前記線形方程式として用い
る請求項１記載の画像流速推定方法。 - 【請求項５】 前記速度推定段階でオプティカルフロー
法により定式化された速度を変数とする連立１次方程式
の解法において、加速緩和法を適用して反復誤差が一定
基準値以下になった時点の解を推定すべき速度ベクトル
とみなす請求項１記載の画像流速推定方法。 - 【請求項６】 時系列画像中の流体の移動速度をオプテ
ィカルフロー法を用いて推定する画像流速推定装置であ
って、 移動する対象を含む２次元画像を入力する画像入力手段
と、前記画像入力手段に入力された画像を時系列画像と
して蓄積する画像蓄積手段と、前記画像蓄積手段に蓄積
されている２次元時系列画像の空間周波数成分を解析す
るとともに、画像の濃淡値の標準偏差を同時に解析する
画像解析手段と、前記画像解析手段で得られた画像の空
間周波数成分と濃淡値の標準偏差の２つの画像特徴量を
予め設計してある線形方程式に代入して最適正則化係数
αを決定する正則化係数算出手段と、連続する２枚の画
像から、各最適正則化係数αに対する、移動する対象の
速度ベクトルをオプティカルフロー法を用いて推定する
速度推定手段と、推定された速度ベクトルを出力する出
力手段を有する画像流速推定装置。 - 【請求項７】 前記画像解析手段は、画像を幾つかのサ
ブブロックに分割して、各サブブロック毎に高速フーリ
エ変換により、エネルギー・水平方向の空間周波数・垂
直方向の空間周波数を軸にもつ３次元の空間周波数画像
に変換し、空間周波数画像の低空間周波数成分が画像中
央に、高空間周波数成分が中央から周辺に向かって分布
するように、高速フーリエ変換後の空間周波数画像の内
部を配置変えし、エネルギー画像をその中心から同心円
状に抜き取り、各円ごとにエネルギーを求めることによ
り前記各円の半径に対応する空間周波数に対するエネル
ギーの分布を得、各空間周波数とそのエネルギーとの積
の和をエネルギーの総和で割った値を算出し、これを代
表空間周波数とする請求項６記載の画像流速推定装置。 - 【請求項８】 前記画像解析手段は、画像を幾つかのサ
ブブロックに分割して、各サブブロック毎に画像の濃淡
値の平均値を算出し、分散値を求めてその平方根をとる
ことで標準偏差を算出する請求項６記載の画像流速推定
装置。 - 【請求項９】 前記正則化係数算出手段は、予めオフラ
インで正弦波の振幅と空間周波数（波数）をさまざまに
変化させながら一定速度で正弦波パターンを移動させる
シミュレーションパターンを予め生成し、オプティカル
フロー法における正則化係数αを同時にさまざまに変化
させながら、シミュレーションパターンから速度ベクト
ルを推定し、ａｒｃｃｏｓｉｎｅ関数を用いて推定され
た速度ベクトルと真の速度ベクトルとの誤差を評価して
誤差が最小となった時のαを選択することで、αと振幅
と空間周波数の３者の関係を定量づけ、さらに振幅から
シミュレーションパターンの濃淡値の標準偏差を算出す
る処理を行って得られる式を前記線形方程式として用い
る請求項６記載の画像流速推定装置。 - 【請求項１０】 前記速度推定手段でオプティカルフロ
ー法により定式化された速度を変数とする連立１次方程
式の解法において、加速緩和法を適用して反復誤差が一
定基準値以下になった時点の解を推定すべき速度ベクト
ルとみなす請求項６記載の画像流速推定装置。 - 【請求項１１】 請求項１から５のいずれか１項に記載
の画像流速推定方法をコンピュータに実行させるための
画像流速推定プログラム記録したコンピュータ読取り可
能な記録媒体。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP03906798A JP3390975B2 (ja) | 1998-02-20 | 1998-02-20 | 画像流速推定方法、装置および画像流速推定プログラムを記録した記録媒体 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP03906798A JP3390975B2 (ja) | 1998-02-20 | 1998-02-20 | 画像流速推定方法、装置および画像流速推定プログラムを記録した記録媒体 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH11238137A JPH11238137A (ja) | 1999-08-31 |
| JP3390975B2 true JP3390975B2 (ja) | 2003-03-31 |
Family
ID=12542793
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP03906798A Expired - Lifetime JP3390975B2 (ja) | 1998-02-20 | 1998-02-20 | 画像流速推定方法、装置および画像流速推定プログラムを記録した記録媒体 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3390975B2 (ja) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4266535B2 (ja) * | 2001-04-27 | 2009-05-20 | 株式会社シー・イー・デー・システム | 黒煙検知システム |
| JP4752685B2 (ja) * | 2006-08-31 | 2011-08-17 | カシオ計算機株式会社 | 撮像装置及びプログラム |
| JP5961589B2 (ja) * | 2013-07-10 | 2016-08-02 | 日本電信電話株式会社 | 映像生成装置、方法及びプログラム |
| CN108983323A (zh) * | 2018-08-08 | 2018-12-11 | 湖北河海科技发展有限公司 | 基于光流法的降水预报方法及预警平台 |
-
1998
- 1998-02-20 JP JP03906798A patent/JP3390975B2/ja not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JPH11238137A (ja) | 1999-08-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hlawatsch et al. | Hierarchical line integration | |
| Barron et al. | Tutorial: Computing 2D and 3D optical flow | |
| Jeong et al. | Adaptive determination of filter scales for edge detection | |
| JP4125786B2 (ja) | 疎配列画像の相関付け | |
| US7352386B1 (en) | Method and apparatus for recovering a three-dimensional scene from two-dimensional images | |
| KR101475382B1 (ko) | 광학적 3차원 측량의 자기 적응 윈도우 푸리에 위상추출방법 | |
| EP2420975A1 (en) | System and method for 3d wireframe reconstruction from video | |
| KR100897077B1 (ko) | 표면 모델로부터 추출된 프리미티브 특징의 조합에 의한 빔 이미지의 스캐닝 시뮬레이션을 위한 방법, 제품, 및 시스템 | |
| Zotin et al. | Techniques for medical images processing using shearlet transform and color coding | |
| Réthoré et al. | Curve and boundaries measurement using B-splines and virtual images | |
| CN115840205B (zh) | 一种基于激光雷达技术的地形面积计量方法和系统 | |
| JP2015524946A (ja) | 画像解像度が改善された超解像画像を形成するための方法及び測定装置 | |
| JP3390975B2 (ja) | 画像流速推定方法、装置および画像流速推定プログラムを記録した記録媒体 | |
| JP3953350B2 (ja) | 映像シーン予測方法、装置、プログラム、および記録媒体 | |
| Rolland-Neviere et al. | Robust diameter-based thickness estimation of 3D objects | |
| JP2015507736A (ja) | ターゲットサイズを見積もるためのシステム及び方法 | |
| Young et al. | A model-based validation framework for PIV and PTV | |
| Wang et al. | Motion estimation from noisy data with unknown distributions using multi-frame phase-preserving denoising | |
| Mageed | Fractal Dimension (D f) of Ismail’s Fourth Entropy (H IV⁽ q, a 1, a 2,.., a k⁾) with Fractal Applications to Algorithms, Haptics, and Transportation | |
| JP2020098455A (ja) | 物体識別システム、物体識別方法、並びに、画像識別プログラム | |
| JP2008008797A (ja) | カメラ挙動推定装置及びカメラ挙動推定プログラム | |
| CN114972882A (zh) | 基于多注意力机制的磨损表面损伤深度估计方法及系统 | |
| Chetverikov | Applying feature tracking to particle image velocimetry | |
| JP4538426B2 (ja) | 移動ベクトル検出装置、移動ベクトル検出方法および移動ベクトル検出プログラム | |
| Chetverikov | Particle image velocimetry by feature tracking |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080124 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090124 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090124 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100124 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110124 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110124 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120124 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130124 Year of fee payment: 10 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313531 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| EXPY | Cancellation because of completion of term |