(57)【要約】  画像処理装置は、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする。前記画像処理装置は、前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１の算出部と、前記第１の算出部によって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２の算出部と、を備える。
2017179510.tif 000002
【特許請求の範囲】
【請求項１】
カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする装置であって、
前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１の算出部と、
前記第１の算出部によって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２の算出部と、を備える画像処理装置。
【請求項２】
前記統計処理は、前記空間変化の方向を示す軸と前記時間変化の方向を示す軸とで構成される座標系において、前記画像を構成する複数の画素のそれぞれについて、前記空間変化を示す値と前記時間変化を示す値とにより座標を構成する処理と、前記座標系での前記座標の分布を近似する回帰式を求める処理と、含む請求項１に記載の画像処理装置。
【請求項３】
前記第２の算出部は、前記回帰式で示される図形が回帰平面の場合、前記空間変化の方向を示す軸と前記時間変化の方向を示す軸とで規定される平面と、前記回帰平面との交線の、前記空間変化の方向を示す軸に対する傾きを算出する請求項２に記載の画像処理装置。
【請求項４】
前記複数の画像の各画素について、前記空間変化を示す値、及び、前記画像のブレを示す値を用いて、前記カメラのブレに起因する信号変化を算出する第３の算出部をさらに備える請求項１〜３のいずれか一項に記載の画像処理装置。
【請求項５】
前記第３の算出部により算出された信号変化を用いて前記複数の画像を補正する補正部をさらに備える請求項４に記載の画像処理装置。
【請求項６】
前記複数の画像のそれぞれから、予め定められた被写体を抽出する抽出部と、
前記抽出部により抽出された被写体を、前記第３の算出部により算出された信号変化を用いて補正する補正部と、をさらに備える請求項４に記載の画像処理装置。
【請求項７】
前記第２の算出部は、前記画像を二以上の領域に分割し、二以上の前記領域のそれぞれについて、前記統計処理をすることにより、前記カメラのブレにより生じる二以上の前記領域のそれぞれのブレを示す値を算出し、
前記第２の算出部は、二以上の前記領域のそれぞれの前記ブレを示す値を用いて、前記画像に含まれる移動体の像を示す前記領域を特定し、
前記第２の算出部は、前記移動体の像を示す前記領域を除く前記領域の前記ブレを示す値を用いて、前記カメラのブレにより生じる前記画像のブレを示す値を算出する請求項１〜６のいずれか一項に記載の画像処理装置。
【請求項８】
前記空間変化を示す値は、前記画素データに対して前記空間変化の方向の方向微分をした値である請求項１〜７のいずれか一項に記載の画像処理装置。
【請求項９】
前記時間変化を示す値は、前記画素データに対して前記時間変化の方向の方向微分をした値である請求項１〜８のいずれか一項に記載の画像処理装置。
【請求項１０】
前記時間変化を示す値とは、前記画素データに含まれる交流成分である請求項１〜８のいずれか一項に記載の画像処理装置。
【請求項１１】
カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする方法であって、
前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１のステップと、
前記第１のステップによって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２のステップと、を備える画像処理方法。
【請求項１２】
カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をするプログラムであって、
前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１のステップと、
前記第１のステップによって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２のステップと、をコンピュータに実行させる画像処理プログラム。
【請求項１３】
監視対象の動画を撮影し、複数の画像が時系列に並べられた動画データを出力する赤外線カメラと、
前記赤外線カメラから入力された前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出し、
前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記赤外線カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出し、
前記複数の画像のそれぞれのブレを示す値に基づいて前記複数の画像のそれぞれを補正する画像処理部と、
前記画像処理部によって補正された前記複数の画像が時系列に並べられた動画を表示するディスプレイと、を備えるガス検知システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、カメラのブレによるノイズに対処できる画像処理技術に関する。
【背景技術】
【０００２】
ガス漏れが発生したとき、漏れたガスが漂っている空間では、わずかな温度変化が生じる。この原理を利用してガスを検知する技術として、赤外線カメラを利用したガス検知が知られている。【０００３】
赤外線カメラによって監視対象（例えば、ガス輸送管どうしが接続されている箇所）の撮影中に、赤外線カメラが揺れることにより、カメラのブレが生じることがある。これにより、赤外線カメラによって撮影された赤外画像には、カメラのブレによるノイズが含まれるので、ガスの検知精度が低下する。【０００４】
カメラのブレに対処する技術として、前フレーム画像と現フレーム画像とをパターンマッチングすることにより手ブレ量を検出し、手ブレ量を基にして、現フレーム画像から低解像度画像を形成する撮像装置が提案されている（例えば、特許文献１参照）。【０００５】
ガスの濃度を示す指標の一つとして、爆発下限界（Ｌｏｗｅｒ  Ｅｘｐｌｏｓｉｖｅ  Ｌｉｍｉｔ）がある。単位は、ｖｏｌ％である。爆発下限界は、空気と混合した可燃性ガスが着火によって爆発を起こす最低濃度である。メタンの場合、爆発下限界は、５ｖｏｌ％である。また、爆発下限界の濃度を１００とする濃度の単位として、％ＬＥＬが知られている。即ち、メタンの濃度５ｖｏｌ％は、１００％ＬＥＬと等しく、メタンの濃度１０ｖｏｌ％は、２００％ＬＥＬと等しい。【０００６】
赤外線カメラを利用したガス検知では、ガスが漂っている空間のガスの濃度を直接測定することができず、ガスの濃度厚み積を測定する。ガスの濃度厚み積とは、ガスの濃度を、ガスが漂っている空間の奥行き方向に沿って積分した値を意味する。ガスの濃度厚み積の単位は、％ＬＥＬ・ｍである。例えば、ガスの濃度が２００％ＬＥＬであり、ガスの厚みが２ｍのとき、濃度厚み積は、４００％ＬＥＬ・ｍとなる。【０００７】
赤外線カメラによって監視対象を遠隔（例えば、５０ｍ以上）から撮影し、ガス検知のアルゴリズムを用いて、比較的低い濃度（例えば、１０％ＬＥＬ・ｍ）のメタンガスを検知する場合がある。このような場合、本発明者は、カメラのブレ量が極めて小さくても（例えば、０．０１〜０．０５画素）、カメラのブレによるノイズを画像から除去しなければ、ガス検知の精度が低下することを見出した。【０００８】
カメラのブレによるノイズを画像から除去するためには、カメラのブレにより生じる画像のブレを示す値（例えば、画像全体のブレ量）を求める必要がある。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開２００６−０６７２９１号公報
【発明の概要】
【００１０】
本発明は、カメラのブレにより生じる画像のブレを示す値を算出できる画像処理装置、画像処理方法、画像処理プログラム、及び、ガス検知システムを提供することを目的とする。
【００１１】
上記目的を達成する本発明の第１の局面に係る画像処理装置は、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする装置であって、第１の算出部と、第２の算出部と、を備える。前記第１の算出部は、前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する。前記第２の算出部は、前記第１の算出部によって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する。
【００１２】
上記並びにその他の本発明の目的、特徴及び利点は、以下の詳細な記載と添付図面から明らかになるであろう。
【図面の簡単な説明】
【００１３】
【図１】現フレームと次のフレームとの間において発生したカメラのブレによる画素値の空間変化を示す値と、画素値の時間変化を示す値との関係を説明する説明図である。
【図２】次のフレームを構成する全画素のそれぞれについて、画素値のｘ軸方向の変化を示す値と画素値の時間変化を示す値との関係を示す三つのグラフである。
【図３Ａ】本実施形態に係る画像処理装置が適用されたガス検知システムの構成を示すブロック図である。
【図３Ｂ】図３Ａに示す画像処理装置のハードウェア構成を示すブロック図である。
【図４】時系列画素データを説明する説明図である。
【図５】第１実施形態に係る画像処理装置の動作を説明するフローチャートである。
【図６】時系列画素データを構成する周波数成分データを説明する説明図である。
【図７】図５のステップＳ１で用いられるバンドパスフィルターを説明する説明図である。
【図８】１番目の画素の時系列画素データから抽出された所定の周波数成分データｄを示すグラフである。
【図９】所定の周波数成分データｄをフレーム単位で説明する説明図である。
【図１０】図５のステップＳ２の処理を説明する説明図である。
【図１１】第１実施形態において、１番目のフレームのＳｘ及びＳｙを求める処理を説明する説明図である。
【図１２】第１実施形態において、２番目〜Ｋ番目のフレームのＳｘ及びＳｙを求める処理を説明する説明図である。
【図１３】第２実施形態に係る画像処理装置の動作を説明するフローチャートである。
【図１４】∂Ｇ／∂ｔを算出する処理を説明する説明図である。
【図１５】第２実施形態において、１番目のフレームのＳｘ及びＳｙを求める処理を説明する説明図である。
【図１６】第２実施形態において、２番目〜Ｋ番目のフレームのＳｘ及びＳｙを求める処理を説明する説明図である。
【図１７】第３実施形態において、Ｓｘの確定処理を説明するフローチャートである。
【図１８】Ｓｘの確定処理を説明する説明図である。
【図１９】二以上の領域に分割された一つのフレームの模式図である。
【図２０】第４実施形態において、Ｓｘを算出する処理を説明するフローチャートである。
【図２１】∂Ｇ／∂ω＝ｙ・∂Ｇ／∂ｘ＋ｘ・∂Ｇ／∂ｙが成立することを説明する説明図である。
【図２２】∂Ｇ／∂ｍ＝ｘ・∂Ｇ／∂ｘ＋ｙ・∂Ｇ／∂ｙが成立することを説明する説明図である。
【発明を実施するための形態】
【００１４】
以下、図面に基づいて本発明の実施形態を詳細に説明する。各図において、同一符号を付した構成は、同一の構成であることを示し、その構成について、既に説明している内容については、その説明を省略する。
【００１５】
本実施形態の原理について説明する。本実施形態は、カメラのブレにより生じる画像のブレ（画像全体のブレ）を示す値を、画像処理を用いて求め、この値と所定式とを用いて、カメラのブレに起因する信号変化Ｎｓ（すなわち、カメラのブレによるノイズ）を求める。図１は、現フレームと次のフレームとの間において発生したカメラのブレによる画素値Ｇの空間変化を示す値と、画素値Ｇの時間変化を示す値との関係を説明する説明図である。フレームは、動画を構成する一つの静止画像（画像）である。本実施形態は、カメラによって撮影された時系列に並ぶ複数の画像として、動画を例にして説明する。画素値Ｇは、画素の画素データの一例である。画素データは、画素が有する数値情報であり、例えば、赤外画像の場合、画素値Ｇ、又は、画素値Ｇを用いて演算された温度を意味する。
【００１６】
現フレームは、ｊ行目に位置する１番目〜５番目の画素が示され、５番目より後の画素は、省略されている。次のフレームは、ｊ行目に位置する１番目〜４番目の画素が示され、４番目より後の画素は、省略されている。ｊ行目の画素が並ぶ方向が、ｘ軸方向とする。画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）は、一つのフレームにおいて、画素値Ｇ（画素データ）の空間変化を示す値の一例である。画素値Ｇは、輝度を示しており、階調は、２５６である。現フレームのｊ行目において、１番目の画素から５番目の画素に向かって画素値Ｇが徐々に大きくなっている。これは、１番目の画素から５番目の画素に向かって、被写体が徐々に明るくなっていることを意味する。
【００１７】
被写体が変化することなく、現フレームと次のフレームとの間において、ｘ軸の正方向にカメラのブレが発生したとする。１番目〜５番目の画素に向かって、被写体が徐々に明るくなっているので、次のフレームのＮ番目の画素の画素値Ｇは、現フレームのＮ番目の画素の画素値Ｇ以上になると考えられる。例えば、次のフレームの１番目の画素の画素値Ｇは、現フレームの１番目の画素の画素値Ｇ以上になると考えられる。
【００１８】
次のフレームのＮ番目の画素の画素値Ｇが、現フレームのＮ番目の画素の画素値Ｇと比べて、どのくらい大きくなるか考える。カメラのブレにより、現フレームと次のフレームとにおいて、ブレ量が＋０．１画素のブレが発生したとする。ブレ量は、カメラのブレにより生じるフレーム（画像）のブレを示す値の一例である。
【００１９】
ブレ量が１画素以内であれば、画素値Ｇの変化量とブレ量とは、正比例すると考えることができる。従って、ブレ量が、＋０．１画素の場合、画素値Ｇは、以下のように変化すると考えられる。現フレームの１番目の画素の画素値Ｇは、０である。現フレームにおいて、１番目の画素と２番目の画素とにおいて、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）は、１０である。従って、次のフレームにおいて、１番目の画素の画素値Ｇは、１（＝０＋０．１×１０）になる。
【００２０】
現フレームの２番目の画素の画素値Ｇは、１０である。現フレームにおいて、２番目の画素と３番目の画素とにおいて、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）は、０である。従って、次のフレームにおいて、２番目の画素の画素値Ｇは、１０（＝１０＋０．１×０）になる。
【００２１】
現フレームの３番目の画素の画素値Ｇは、１０である。現フレームにおいて、３番目の画素と４番目の画素とにおいて、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）は、２０である。従って、次のフレームにおいて、３番目の画素の画素値Ｇは、１２（＝１０＋０．１×２０）になる。
【００２２】
現フレームの４番目の画素の画素値Ｇは、３０である。現フレームにおいて、４番目の画素と５番目の画素とにおいて、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）は、１０である。従って、次のフレームにおいて、４番目の画素の画素値Ｇは、３１（＝３０＋０．１×１０）になる。
【００２３】
連続する２つのフレームにおいて、画素値Ｇ（画素データ）の時間変化を示す値について説明する。現フレームと次のフレームとにおいて（連続する２つのフレームにおいて）、１番目の画素の画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）は、１（＝１−０）であり、２番目の画素の画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）は、０（＝１０−１０）であり、３番目の画素の画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）は、２（＝１２−１０）であり、４番目の画素の画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）は、１（＝３１−３０）である。
【００２４】
以上より、カメラのブレによる画素値Ｇ（画素データ）のｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）と、時間変化を示す値（∂Ｇ／∂ｔ）との関係は、以下の式１で表されると推測できる。
【００２５】
∂Ｇ／∂ｔ＝Ｓｘ・∂Ｇ／∂ｘ・・・（式１）
【００２６】
Ｓｘは、カメラのブレにより生じるフレーム（画像）のｘ軸方向のブレを示す値である。図１の場合、Ｓｘは、連続する２つのフレーム（画像）間のｘ軸方向のブレ量である。
【００２７】
カメラのブレによる画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）は、カメラのブレに起因する信号変化Ｎｓであり、従って、カメラのブレに起因する信号変化Ｎｓは、以下の式２で表されると推測できる。
【００２８】
Ｎｓ＝∂Ｇ／∂ｔ＝Ｓｘ・∂Ｇ／∂ｘ・・・（式２）
【００２９】
図２は、次のフレームを構成する全画素のそれぞれについて、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）と画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）との関係を示す三つのグラフである。グラフ１、グラフ２及びグラフ３の横軸は、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）を示し、縦軸は、画素値Ｇの時間変化を示す値（∂Ｇ／∂ｔ）を示す。直線ＳＬで表される式は、式２である。×は、次のフレームを構成する各画素に対応する座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｔ）を示す。例えば、図１を参照して、次のフレームの１番目の画素に対応する座標は、（９，１）となり、２番目の画素に対応する座標は、（２，０）となり、３番目の画素に対応する座標は、（１９，２）となる。
【００３０】
図２を参照して、グラフ１は、被写体が変化しておらず、ｘ軸方向にカメラのブレが発生した場合を示している。座標には、被写体変化に起因する信号変化が含まれておらず、カメラのブレに起因する信号変化Ｎｓが含まれている。このため、座標は、直線ＳＬ上に分布している。被写体変化に起因する信号変化とは、例えば、ガス漏れの監視対象から漏れたガスの像を示す。
【００３１】
グラフ２は、被写体が変化し、かつ、ｘ軸方向にカメラのブレが発生した場合を示している。座標には、カメラのブレに起因する信号変化Ｎｓ、及び、被写体変化に起因する信号変化が含まれている。このため、直線ＳＬより上に分布する座標、及び、直線ＳＬより下に分布する座標が存在している。
【００３２】
被写体の変化は、一般に、ｘ軸方向のカメラのブレと相関しないので、カメラのブレに起因する信号変化Ｎｓと被写体変化に起因する信号変化とは、分離することができる。そこで、本実施形態は、フレームからカメラのブレに起因する信号変化Ｎｓを除去する補正をする。これにより、カメラのブレに起因する信号変化Ｎｓをフレームから除去し、被写体変化に起因する信号変化をフレームに残すことができる。グラフ３は、グラフ２に対して、カメラのブレに起因する信号変化Ｎｓを除去する補正をした場合を示す。
【００３３】
Ｓｘ（カメラのブレにより生じるフレームのｘ軸方向のブレを示す値）を求める方法としては、平均値を利用する方法と統計処理を利用する方法がある。平均値を利用する方法は、グラフ２に示す座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｔ）において、∂Ｇ／∂ｘの平均値、及び、∂Ｇ／∂ｔの平均値を求め、これらの平均値と式２とからＳｘを求める。但し、∂Ｇ／∂ｘ＝０付近では、Ｓｘの値が不安定となるので、∂Ｇ／∂ｘ＝０付近の座標は除外する。
【００３４】
統計処理を利用する方法は、グラフ２に示す座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｔ）を統計処理（例えば、回帰直線を求める処理）することで、Ｓｘを求める。グラフ２に示す座標を用いて求めた回帰直線は、直線ＳＬとなるので、グラフ２に示す座標を用いて回帰直線を求めることにより、Ｓｘが求まる。
【００３５】
統計処理は、空間変化の方向を示す軸と時間変化の方向を示す軸とで構成される座標系において、フレーム（画像）を構成する複数の画素のそれぞれについて、空間変化を示す値と時間変化を示す値とにより座標を構成する処理と、座標系での座標の分布を近似する回帰式を求める処理と、含む。回帰式で示される図形は、回帰直線、回帰曲線、回帰平面や回帰曲面がある。図２は、回帰直線の場合を示している。
【００３６】
以上において、説明を簡単にするために、カメラのブレがｘ軸方向に発生した場合で説明した。カメラのブレがｘ軸方向及びｙ軸方向に発生した場合、カメラのブレに起因する信号変化Ｎｓは、式１及び式２の場合と同様に考えると、以下の式３で表されると推測できる。
【００３７】
Ｎｓ＝∂Ｇ／∂ｔ＝Ｓｘ・∂Ｇ／∂ｘ＋Ｓｙ・∂Ｇ／∂ｙ・・・（式３）
【００３８】
∂Ｇ／∂ｙは、画素値Ｇ（画素データ）のｙ軸方向の変化を示す値（画素値の空間変化を示す値）である。Ｓｙは、カメラのブレにより生じるフレーム（画像）のｙ軸方向のブレを示す値である。このような値として、例えば、連続する２つのフレーム間のｙ軸方向のブレ量がある。
【００３９】
以上が本実施形態の原理の説明である。
【００４０】
図３Ａを用いて、本実施形態に係る画像処理装置３が適用されたガス検知システム１について説明する。図３Ａは、このシステム１の構成を示すブロック図である。ガス検知システム１は、赤外線カメラ２と画像処理装置３とを備える。
【００４１】
赤外線カメラ２は、ガス漏れの監視対象及び背景の赤外画像の動画（監視対象の動画）を撮影し、動画を示す動画データＤ１を生成する。動画に限らず、赤外線カメラ２によって、ガス漏れの監視対象及び背景の赤外画像を複数の時刻で撮影してもよい。赤外線カメラ２は、光学系４、フィルター５、二次元イメージセンサー６及び信号処理部７を備える。
【００４２】
光学系４は、被写体（監視対象及び背景）の赤外画像を二次元イメージセンサー６上で結像させる。フィルター５は、光学系４と二次元イメージセンサー６との間に配置され、光学系４を通過した光のうち、特定波長の赤外線のみを通過させる。赤外の波長帯のうち、フィルター５を通過させる波長帯は、検知するガスの種類に依存する。例えばメタンの場合、３．２〜３．４μｍの波長帯を通過させるフィルター５が用いられる。二次元イメージセンサー６は、例えば、冷却型インジウムアンチモン（ＩｎＳｂ）イメージセンサーであり、フィルター５を通過した赤外線を受光する。信号処理部７は、二次元イメージセンサー６から出力されたアナログ信号を、デジタル信号に変換し、公知の画像処理をする。このデジタル信号が、動画データＤ１となる。
【００４３】
画像処理装置３は、パーソナルコンピュータ、スマートフォン、タブレット端末等であり、機能ブロックとして、画像データ入力部８、画像処理部９、表示制御部１０、ディスプレイ１１及び入力部１２を備える。
【００４４】
画像データ入力部８は、赤外線カメラ２の通信部（不図示）と通信する通信インターフェイスである。画像データ入力部８には、赤外線カメラ２の通信部から送られてきた動画データＤ１が入力される。画像データ入力部８は、動画データＤ１を画像処理部９へ送る。
【００４５】
動画データＤ１は、フレームが時系列に複数並べられた構造を有する。複数のフレーム（複数の赤外画像）において、同じ位置にある画素の画素データを時系列に並べたデータを、時系列画素データとする。時系列画素データを具体的に説明する。図４は、時系列画素データを説明する説明図である。赤外画像の動画のフレーム数をＫとする。一つのフレームがＭ個の画素、すなわち、１番目の画素、２番目の画素、・・・、Ｍ−１番目の画素、Ｍ番目の画素で構成されている。画素データは、上述したように、画素が有する数値情報であり、例えば、画素の画素値Ｇである。赤外画像の場合、画素値Ｇ、又は、画素値Ｇを用いて演算された温度である。
【００４６】
複数（Ｋ個）のフレームの同じ位置にある画素とは、同じ順番の画素を意味する。例えば、１番目の画素で説明すると、１番目のフレームに含まれる１番目の画素の画素データ、２番目のフレームに含まれる１番目の画素の画素データ、・・・、Ｋ−１番目のフレームに含まれる１番目の画素の画素データ、Ｋ番目のフレームに含まれる１番目の画素の画素データを、時系列に並べたデータが、１番目の画素の時系列画素データとなる。また、Ｍ番目の画素で説明すると、１番目のフレームに含まれるＭ番目の画素の画素データ、２番目のフレームに含まれるＭ番目の画素の画素データ、・・・、Ｋ−１番目のフレームに含まれるＭ番目の画素の画素データ、Ｋ番目のフレームに含まれるＭ番目の画素の画素データを、時系列に並べたデータが、Ｍ番目の画素の時系列画素データとなる。時系列画素データの数は、一つのフレームを構成する画素の数と同じである。
【００４７】
図３Ａの説明に戻る。画像処理部９は、第１の算出部９１、第２の算出部９２、第３の算出部９３及び補正部９４を備える。これらについては、後で説明する。
【００４８】
表示制御部１０は、動画データＤ１で示される動画、及び、画像処理部９で所定の処理がされた動画（別の表現を用いれば、画像処理部９によって補正された複数の画像が時系列に並べられた画像）を、ディスプレイ１１に表示させる。
【００４９】
入力部１２は、ガス検知に関連する各種入力がされる。
【００５０】
図３Ｂは、図３Ａに示す画像処理装置３のハードウェア構成を示すブロック図である。画像処理装置３は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）３ａ、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）３ｂ、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）３ｃ、ＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）３ｄ、液晶ディスプレイ３ｅ、通信インターフェイス３ｆ、キーボード等３ｇ、及び、これらを接続するバス３ｈを備える。液晶ディスプレイ３ｅは、ディスプレイ１１を実現するハードウェアである。液晶ディスプレイ３ｅの替わりに、有機ＥＬディスプレイ（Ｏｒｇａｎｉｃ  Ｌｉｇｈｔ  Ｅｍｉｔｔｉｎｇ  Ｄｉｏｄｅ  ｄｉｓｐｌａｙ）、プラズマディスプレイ等でもよい。通信インターフェイス３ｆは、画像データ入力部８を実現するハードウェアである。キーボード等３ｇは、入力部１２を実現するハードウェアである。キーボードの替わりに、タッチパネルでもよい。
【００５１】
ＨＤＤ３ｄには、画像処理部９及び表示制御部１０について、これらの機能ブロックをそれぞれ実現するためのプログラム、及び、各種データ（例えば、動画データＤ１）が格納されている。画像処理部９を実現するプログラムは、動画データＤ１（画像データ）を取得し、動画データＤ１に上記所定の処理をする画像処理プログラムである。表示制御部１０を実現するプログラムは、例えば、動画データＤ１で示される動画をディスプレイ１１に表示させたり、画像処理部９によって上記所定の処理がされた動画をディスプレイ１１に表示させたりする表示制御プログラムである。これらのプログラムは、ＨＤＤ３ｄに予め記憶されているが、これに限定されない。例えば、これらのプログラムを記録している記録媒体（例えば、磁気ディスク、光学ディスクのような外部記録媒体）が用意されており、この記録媒体に記憶されているプログラムがＨＤＤ３ｄに記憶されてもよい。また、これらのプログラムは、画像処理装置３とネットワーク接続されたサーバに格納されており、ネットワークを介して、これらのプログラムがＨＤＤ３ｄに送られ、ＨＤＤ３ｄに記憶されてもよい。これらのプログラムは、ＨＤＤ３ｄの替わりにＲＯＭ３ｃに記憶してもよい。画像処理装置３は、ＨＤＤ３ｄの替わりに、フラッシュメモリを備え、これらのプログラムはフラッシュメモリに記憶してもよい。
【００５２】
ＣＰＵ３ａは、これらのプログラムを、ＨＤＤ３ｄから読み出してＲＡＭ３ｂに展開させ、展開されたプログラムを実行することによって、画像処理部９及び表示制御部１０が実現される。但し、画像処理部９の機能及び表示制御部１０の機能について、各機能の一部又は全部は、ＣＰＵ３ａによる処理に替えて、又は、これと共に、ＤＳＰ（Ｄｉｇｉｔａｌ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｏｒ）による処理によって実現されてもよい。又、同様に、各機能の一部又は全部は、ソフトウェアによる処理に替えて、又は、これと共に、専用のハードウェア回路による処理によって実現されてもよい。
【００５３】
なお、画像処理装置３は、次に説明するように、第１実施形態から第４実施形態がある。これらの実施形態は、それぞれ、複数の要素によって構成される。従って、ＨＤＤ３ｄには、これらの要素を実現するためのプログラムが格納されている。例えば、画像処理装置３の第１実施形態は、要素として、第１の算出部９１、第２の算出部９２、第３の算出部９３及び補正部９４を含む。ＨＤＤ３ｄには、これらの要素のそれぞれを実現するためのプログラムが格納されている。これらのプログラムは、第１の算出プログラム、第２の算出プログラム、第３の算出プログラム、補正プログラムと表現される。
【００５４】
これらのプログラムは、要素の定義を用いて表現される。第１の算出部９１及び第１の算出プログラムを例にして説明する。第１の算出部９１は、赤外線カメラ２によって撮影された時系列に並ぶ複数の画像（動画）について、画素データの空間変化を示す値、及び、画素データの時間変化を示す値を算出する。第１の算出プログラムは、赤外線カメラ２によって撮影された時系列に並ぶ複数の画像（動画）について、画素データの空間変化を示す値、及び、画素データの時間変化を示す値を算出するプログラムである。
【００５５】
ＣＰＵ３ａによって実行されるこれらのプログラム（第１の算出プログラム、第２の算出プログラム、第３の算出プログラム、補正プログラム）のフローチャートが、後で説明する図５である。
【００５６】
上述したように、本実施形態に係る画像処理装置３には、第１実施形態〜第４実施形態がある。いずれの実施形態も、カメラのブレがｘ軸方向（第１の変化方向）及びｙ軸方向（第２の変化方向）に発生した場合を例にして説明する。この場合、カメラのブレに起因する信号変化Ｎｓは、上記式３で表される。第１実施形態〜第４実施形態で用いられる主な用語について説明する。
【００５７】
画素値Ｇ（画素データ）の空間変化を示す値として、画素値Ｇのｘ軸方向の変化を示す値（∂Ｇ／∂ｘ）や、画素値Ｇ（画素データ）のｙ軸方向の変化を示す値（∂Ｇ／∂ｙ）がある。∂Ｇ／∂ｘは、画素値Ｇのｘ軸方向の変化を方向微分することにより求められる。言い換えれば、∂Ｇ／∂ｘは、ｘ軸方向に隣り合う画素の画素値Ｇの差分である。∂Ｇ／∂ｙは、画素値Ｇのｙ軸方向の変化を方向微分することにより求められる。言い換えれば、∂Ｇ／∂ｙは、ｙ軸方向に隣り合う画素の画素値Ｇの差分である。∂Ｇ／∂ｘ、∂Ｇ／∂ｙのいずれも、画素度Ｇ（画素データ）に対して空間変化の方向の方向微分をした値である。
【００５８】
画素値Ｇ（画素データ）の時間変化を示す値として、連続する２つのフレーム（画像）において、同じ順番の画素の画素値Ｇ（画素データ）を差分した値（∂Ｇ／∂ｔ）や、図８及び図９に示す所定の周波数成分データｄがある。前者は、第２及び第３実施形態で用いられ、言い換えれば、画素値Ｇに対して、時間変化の方向の方向微分をした値である。後者は、第１実施形態で用いられ、言い換えれば、１番目〜Ｍ番目の画素の画素データに含まれる交流成分である。
【００５９】
カメラのブレにより生じる画像のブレを示す値として、カメラのブレにより生じるフレーム（画像）のｘ軸方向のブレを示す値（Ｓｘ）や、カメラのブレにより生じるフレーム（画像）のｙ軸方向のブレを示す値（Ｓｙ）がある。
【００６０】
Ｓｘとして、カメラのブレにより生じる、連続する２つのフレーム（画像）間のｘ軸方向のブレ量や、カメラのブレにより生じる、各フレーム（画像）のｘ軸方向のブレの周波数成分がある。前者は、第２及び第３実施形態で用いられ、後者は、第１実施形態で用いられる。
【００６１】
Ｓｙとして、カメラのブレにより生じる、連続する２つのフレーム（画像）間のｙ軸方向のブレ量や、カメラのブレにより生じる、各フレーム（画像）のｙ軸方向のブレの周波数成分がある。前者は、第２及び第３実施形態で用いられ、後者は、第１実施形態で用いられる。
【００６２】
第１実施形態から順番に説明する。図５は、第１実施形態に係る画像処理装置３の動作を説明するフローチャートであり、画像処理装置３のＣＰＵ３ａが画像処理プログラムに従って実行するものである。図３Ａ及び図５を参照して、第１の算出部９１は、バンドパス処理部の機能を有する。第１の算出部９１は、図４に示す１番目〜Ｋ番目のフレームにより構成されるオリジナルの動画に対して、バンドパス処理をする（ステップＳ１）。
【００６３】
これについて説明する。図６は、時系列画素データを構成する周波数成分データを説明する説明図である。時系列画素データは、以下の（１）〜（５）の周波数成分データにより構成される。
【００６４】
（１）背景の温度変化を示す信号Ｂ（言い換えれば、背景の温度変化を示す周波数成分データ）
（２）カメラのブレに起因する信号変化Ｎｓ（言い換えれば、カメラのブレによるノイズを示す周波数成分データ）
（３）漏れたガスによる温度変化を示す信号ＮＬ（言い換えれば、漏れたガスによる温度変化を示す周波数成分データ）
（４）その他のノイズＮｃ（言い換えれば、その他のノイズを示す周波数成分データ）
（５）高周波ノイズＮｈ（言い換えれば、高周波ノイズを示す周波数成分データ）
【００６５】
（３）によって、漏れたガスの像（ガス像）が形成される。（２）、（３）、（４）は、（１）より周波数が高く、（５）より周波数が低い。（１）は、（３）と比べて、かなり大きく、かつ、空間変化の方向の方向微分に対して値の偏りが発生するので、除去する必要がある。（２）及び（５）は、（１）に比べて小さいが、無視できない大きさなので、除去する必要がある。（４）は、（３）と比べてかなり小さいので、無視する。
【００６６】
第１の算出部９１は、バンドパスフィルターを用いて、図４に示すＭ個の時系列画素データを処理することにより、Ｍ個の時系列画素データのそれぞれについて、背景の温度変化を示す信号Ｂ、及び、高周波ノイズＮｈを除去する。言い換えれば、第１の算出部９１は、バンドパスフィルターを用いて、Ｍ個の時系列画素データを処理することにより、Ｍ個の時系列画素データのそれぞれから、所定の周波数成分データｄを抽出する。所定の周波数成分データｄは、カメラのブレに起因する信号変化Ｎｓ、漏れたガスによる温度変化を示す信号ＮＬ、及び、その他のノイズＮｃにより構成される。所定の周波数成分データｄは、オリジナルの動画から、背景の温度変化を示す信号Ｂ、及び、高周波ノイズＮｈが除去された動画である。従って、所定の周波数成分データｄは、ガス像の動画を含むが、背景（監視対象）の動画を含まない。このように、第１の算出部９１は、抽出部（不図示）として機能し、抽出部は、フレームから予め定められた被写体（ここでは、ガス像）を抽出することを被写体抽出とし、１番目〜Ｋ番目のフレームのそれぞれに対して、被写体抽出をする。別の表現を用いれば、抽出部は、複数の画像（複数のフレーム）のそれぞれから、予め定められた被写体を抽出する。
【００６７】
漏れたガスによる温度変化を示す信号ＮＬの周波数を、０．３〜３Ｈｚとする。例えば、０．３〜３Ｈｚを抽出できる重み付け係数を用いて、図４に示すＫ個のフレーム数より少ない所定数のフレームを単位とする加重移動平均を算出する。これにより、時系列画素データから、所定の周波数成分データｄが抽出される。
【００６８】
図７は、ステップＳ１で用いられるバンドパスフィルターを説明する説明図である。横軸は、フレームを示し、縦軸は、重み付け係数を示している。所定数のフレームは、例えば、９９フレームである。内訳は、ターゲットとなるフレーム、これより前の連続する４９フレーム、これより後の連続する４９フレームである。
【００６９】
ステップＳ１の処理により、Ｍ個の所定の周波数成分データｄが生成される。内訳は、１番目の画素の時系列画素データから抽出された所定の周波数成分データｄ、２番目の画素の時系列画素データから抽出された所定の周波数成分データｄ、・・・、Ｍ番目の画素の時系列画素データから抽出された所定の周波数成分データｄである。
【００７０】
図８は、１番目の画素の時系列画素データから抽出された所定の周波数成分データｄを示すグラフである。グラフの縦軸は、画素値Ｇを示し、横軸は、フレームの順番を示す。所定の周波数成分データｄは、上述したように、カメラのブレに起因する信号変化Ｎｓ、漏れたガスによる温度変化を示す信号ＮＬ、及び、その他のノイズＮｃにより構成される。
【００７１】
図９は、所定の周波数成分データｄをフレーム単位で説明する説明図である。各フレームにおいて、１番目の画素の所定の周波数成分データｄの値をｄ（１）、２番目の画素の所定の周波数成分データｄの値をｄ（２）、・・・、Ｍ番目の画素の所定の周波数成分データｄの値をｄ（Ｍ）とする。各フレームの１番目の画素のｄ（１）を時系列に並べたデータが、１番目の画素の時系列画素データから抽出された所定の周波数成分データｄとなり、各フレームの２番目の画素のｄ（２）を時系列に並べたデータが、２番目の画素の時系列画素データから抽出された所定の周波数成分データｄとなり、・・・、各フレームのＭ番目の画素のｄ（Ｍ）を時系列に並べたデータが、Ｍ番目の画素の時系列画素データから抽出された所定の周波数成分データｄとなる。
【００７２】
ｄ（１）〜ｄ（Ｍ）は、それぞれ、１番目〜Ｍ番目の画素の画素データに含まれる交流成分（言い換えれば、周波数成分）である。ｄ（１）〜ｄ（Ｍ）は、連続する２フレームにおいて、画素値Ｇ（画素データ）の時間変化を示す値の例である。このように、第１の算出部９１は、画素値Ｇ（画素データ）の時間変化を示す値を算出する。
【００７３】
図３Ａ及び図５を参照して、第１の算出部９１は、１番目〜Ｋ番目のフレームのそれぞれに対して、∂Ｇ／∂ｘ及び∂Ｇ／∂ｙを算出する処理をする（ステップＳ２）。図１０は、この処理を説明する説明図である。∂Ｇ／∂ｘは、一つのフレームにおいて、画素値Ｇのｘ軸方向の変化を示す値である。∂Ｇ／∂ｙは、一つのフレームにおいて、画素値Ｇのｙ軸方向の変化を示す値である。いずれも、一つのフレームにおいて、画素値Ｇ（画素データ）の空間変化を示す値である。
【００７４】
ステップＳ１の処理及びステップＳ２の処理の順番は、いずれが先でもよいし、同時に処理をしてもよい。以上説明したように、第１の算出部９１は、複数の画像（複数のフレーム）について、画素データの空間変化を示す値（ステップＳ２）、及び、画素データの時間変化を示す値を算出する（ステップＳ１）。
【００７５】
図３Ａ及び図５を参照して、第２の算出部９２は、ステップＳ１の結果及びステップＳ２の結果を用いて、１番目〜Ｋ番目のフレームのそれぞれについて、Ｓｘ及びＳｙを算出する（ステップＳ３）。これについて説明する。
【００７６】
図１１は、１番目のフレームのＳｘ及びＳｙを求める処理を説明する説明図である。図１１に示す二つのグラフは、∂Ｇ／∂ｘを示す軸、∂Ｇ／∂ｙを示す軸、及び、ｄを示す軸を有する。これらの軸は、互いに直交している。∂Ｇ／∂ｘ（画素値Ｇのｘ軸方向の変化を示す値）は、図１０に示す１番目のフレームの∂Ｇ／∂ｘである。∂Ｇ／∂ｙ（画素値Ｇのｙ軸方向の変化を示す値）は、図１０に示す１番目のフレームの∂Ｇ／∂ｙである。ｄ（画素値Ｇの時間変化を示す値）は、図９に示す所定の周波数成分データｄのうち、１番目のフレームの周波数成分データｄ（１）〜ｄ（Ｍ）である。∂Ｇ／∂ｘを示す軸は、第１の変化方向を示す軸であり、∂Ｇ／∂ｙを示す軸は、第２の変化方向を示す軸であり、ｄを示す軸は、時間方向の変化を示す軸である。
【００７７】
第２の算出部９２は、図１０に示す１番目のフレームを構成する１番目〜Ｍ番目の画素の∂Ｇ／∂ｘ及び∂Ｇ／∂ｙ、並びに、図９に示す１番目のフレームを構成する１番目〜Ｍ番目の画素のｄ（１）〜ｄ（Ｍ）を用いて、１番目のフレームを構成する１番目〜Ｍ番目の画素のそれぞれについて、座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｙ，ｄ（１））〜座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｙ，ｄ（Ｍ））を生成し、これらの座標に対する回帰面を算出する。回帰面の算出は、統計処理の一例である。
【００７８】
第２の算出部９２は、∂Ｇ／∂ｘを示す軸とｄを示す軸とで規定される平面と、１番目のフレームの回帰面との交線を求め、∂Ｇ／∂ｘを示す軸に対するその交線の傾きを、１番目のフレームのＳｘとする（Ｓｘ（１））。第２の算出部９２は、∂Ｇ／∂ｙを示す軸とｄを示す軸とで規定される平面と、１番目のフレームの回帰面との交線を求め、∂Ｇ／∂ｙを示す軸に対するその交線の傾きを、１番目のフレームのＳｙとする（Ｓｙ（１））。すなわち、第２の算出部９２は、第１の変化方向を示す軸と時間変化の方向を示す軸とで規定される平面と、回帰平面との交線を第１の交線とし、第２の変化方向を示す軸と時間変化の方向を示す軸とで規定される平面と、回帰平面との交線を第２の交線とし、第１の変化方向を示す軸に対する第１の交線の傾きを、第１の変化方向に沿ったフレーム（画像）のブレを示す値とし、第２の変化方向を示す軸に対する第２の交線の傾きを、第２の変化方向に沿ったフレーム（画像）のブレを示す値とする。これを上位概念で表現すると、第２の算出部９２は、回帰式で示される図形が回帰平面の場合、空間変化の方向を示す軸と時間変化の方向を示す軸とで規定される平面と、回帰平面との交線の、空間変化の方向を示す軸に対する傾きを算出する。
【００７９】
図１２は、２番目のフレームのＳｘ（Ｓｘ（２））及び２番目のフレームのＳｙ（Ｓｙ（２））〜Ｋ番目のフレームのＳｘ（Ｓｘ（Ｋ））及びＫ番目のフレームのＳｙ（Ｓｙ（Ｋ））を求める処理を説明する説明図である。第２の算出部９２は、２番目〜Ｋ番目のフレームのＳｘ及びＳｙを、１番目のフレームと同様にして求める。以上説明したように、第２の算出部９２は、第１の算出部９１によって算出された空間変化を示す値、及び、時間変化を示す値を統計処理して、カメラのブレにより生じる複数の画像（複数のフレーム）のそれぞれのブレを示す値を算出する。
【００８０】
図３Ａ及び図５を参照して、第３の算出部９３は、ステップＳ２の結果及びステップＳ３の結果を用いて、カメラのブレに起因する信号変化Ｎｓを算出する（ステップＳ４）。別の表現を用いれば、第３の算出部９３は、複数の画像（複数のフレーム）の各画素について、空間変化を示す値、及び、画像のブレを示す値を用いて、カメラのブレに起因する信号変化Ｎｓを算出する。これについて説明する。第３の算出部９３は、１番目〜Ｍ番目の画素のそれぞれに関して、カメラのブレに起因する信号変化Ｎｓを算出する。これを説明する表を、以下に示す。Ｎｓは、上記式３を用いて算出する。
【００８１】
【表１】
2017179510.tif 000003
【００８２】
１番目の画素に関するＮｓから説明する。１番目の画素に関するＮｓは、１番目のフレームに含まれる１番目の画素に関するＮｓ、２番目のフレームに含まれる１番目の画素に関するＮｓ、・・・、Ｋ番目のフレームに含まれる１番目の画素に関するＮｓを時系列に並べたデータである。
【００８３】
１番目のフレームに含まれる１番目の画素に関するＮｓは、図１０に示す１番目のフレームに含まれる１番目の画素の∂Ｇ／∂ｘ及び∂Ｇ／∂ｙ、並びに、図１１に示す１番目のフレームのＳｘ（Ｓｘ（１））及びＳｙ（Ｓｙ（１））を、式３に代入することにより求められる。２番目のフレームに含まれる１番目の画素に関するＮｓは、図１０に示す２番目のフレームに含まれる１番目の画素の∂Ｇ／∂ｘ及び∂Ｇ／∂ｙ、並びに、図１２に示す２番目のフレームのＳｘ（Ｓｘ（２））及びＳｙ（Ｓｙ（２））を、式３に代入することにより求められる。Ｋ番目のフレームに含まれる１番目の画素に関するＮｓは、図１０に示すＫ番目のフレームに含まれる１番目の画素の∂Ｇ／∂ｘ及び∂Ｇ／∂ｙ、並びに、図１２に示すＫ番目のフレームのＳｘ（Ｓｘ（Ｋ））及びＳｙ（Ｓｙ（Ｋ））を、式３に代入することにより求められる。
【００８４】
２番目の画素に関するＮｓは、１番目のフレームに含まれる２番目の画素に関するＮｓ、２番目のフレームに含まれる２番目の画素に関するＮｓ、・・・、Ｋ番目のフレームに含まれる２番目の画素に関するＮｓを時系列に並べたデータである。Ｍ番目の画素に関するＮｓは、１番目のフレームに含まれるＭ番目の画素に関するＮｓ、２番目のフレームに含まれるＭ番目の画素に関するＮｓ、・・・、Ｋ番目のフレームに含まれるＭ番目の画素に関するＮｓを時系列に並べたデータである。
【００８５】
図３Ａ及び図５を参照して、補正部９４は、所定の周波数成分データｄとカメラのブレに起因する信号変化Ｎｓとの差分を算出する（ステップＳ５）。これについて説明する。図８及び表１を参照して、補正部９４は、１番目の画素の時系列画素データから抽出された所定の周波数成分データｄと、１番目の画素に関するＮｓとの差分、２番目の画素の時系列画素データから抽出された所定の周波数成分データｄと、２番目の画素に関するＮｓとの差分、・・・、Ｍ番目の画素の時系列画素データから抽出された所定の周波数成分データｄと、Ｍ番目の画素に関するＮｓとの差分を算出する。
【００８６】
図６を参照して、ステップＳ５の差分により得られた周波数成分データは、漏れたガスによる温度変化を示す信号ＮＬ、及び、その他のノイズＮｃにより構成される。その他のノイズＮｃは、漏れたガスによる温度変化を示す信号ＮＬと比べてかなり小さいので、ステップＳ５の差分により得られた周波数成分データは、漏れたガスによる温度変化を示す信号ＮＬと見なすことができる。従って、ステップＳ５は、ガス像の動画に対して、カメラのブレを補正する処理となる。
【００８７】
以上説明したように、補正部９４は、カメラのブレに起因する信号変化Ｎｓを用いて被写体（ここでは、ガス像）を補正することを第２の補正とし、動画を構成する１番目〜Ｋ番目のフレームのそれぞれから抽出された被写体に対して、第２の補正をする。別の表現を用いれば、補正部９４は、上記抽出部により抽出された被写体を、第３の算出部９３により算出された信号変化Ｎｓを用いて補正する。
【００８８】
図３Ａ及び図５を参照して、表示制御部１０は、ステップＳ５の差分により得られた周波数成分データで示される動画（すなわち、カメラのブレが補正されたガス像の動画）を、ディスプレイ１１に表示させる（ステップＳ６）。
【００８９】
第１実施形態の主な効果を説明する。図１を参照して、（ａ）カメラのブレが発生すれば、画素値Ｇの空間方向の変化と画素値Ｇの時間方向の変化とが発生する。（ｂ）カメラのブレ量が１画素以内であれば、画素値Ｇの変化量とブレ量とは、正比例すると考えられる。本発明者は、これらを前提にして、赤外線カメラ２によって撮影された動画を構成する各フレーム（画像）について、フレームを構成する複数の画素のそれぞれの空間変化を示す値、及び、時間変化を示す値に対して、所定の統計処理（回帰平面の算出）をすれば、各フレームのＳｘ及びＳｙが求まることを見出した（図１１、図１２）。
【００９０】
図５を参照して、第１実施形態は、上記（ａ）及び（ｂ）を基にして、各フレームのＳｘ及びＳｙを算出するので（ステップＳ３）、カメラのブレが極めて小さくても（例えば、０．０１〜０．０５画素）、各フレームのＳｘ及びＳｙを算出することができる。そして、第１実施形態は、Ｓｘ及びＳｙを用いて、カメラのブレに起因する信号変化Ｎｓを算出し（ステップＳ４）、信号変化Ｎｓを用いて、ガス像を補正する（ステップＳ５）。従って、第１実施形態によれば、カメラのブレが極めて小さくても、それによるノイズを除去できるので、ガスの検知精度を向上させることができる。
【００９１】
以上の効果は、第２実施形態〜第４実施形態でも言えることである。
【００９２】
第２実施形態について説明する。第１実施形態は、カメラのブレが補正されたガス像の動画を生成している。これに対して、第２実施形態は、オリジナルの動画に対して、カメラのブレが補正された動画（すなわち、監視対象の像及びガス像を含む動画）を生成する。
【００９３】
図１３は、第２実施形態に係る画像処理装置３の動作を説明するフローチャートであり、画像処理装置３のＣＰＵ３ａが画像処理プログラムに従って実行するものである。図３Ａ及び図１３を参照して、第１の算出部９１は、図４に示す１番目〜Ｋ番目のフレームにより構成されるオリジナルの動画に対して、∂Ｇ／∂ｘ及び∂Ｇ／∂ｙを算出する処理、並びに、∂Ｇ／∂ｔを算出する処理をする（ステップＳ１１）。∂Ｇ／∂ｘ及び∂Ｇ／∂ｙを算出する処理は、第１実施形態のステップＳ２（図５）と同じなので、説明を省略する。
【００９４】
∂Ｇ／∂ｔを算出する処理について説明する。図１４は、これを説明する説明図である。∂Ｇ／∂ｔは、画素値Ｇ（画素データ）に対して時間変化の方向の方向微分をした値である。第１実施形態では、時間変化を示す値として、図８及び図９で説明した所定の周波数成分データｄを用いたが、第２実施形態では、時間変化を示す値として、∂Ｇ／∂ｔを用いる。∂Ｇ／∂ｔは、連続する２つのフレームにおいて、同じ順番の画素の画素値Ｇの差分である。
【００９５】
連続する２つのフレームは、２番目のフレームの場合、１番目と２番目のフレームであり、３番目のフレームの場合、２番目と３番目のフレームであり、・・・、Ｋ番目のフレームの場合、Ｋ−１番目とＫ番目のフレームである。１番目のフレームの場合、連続する２つのフレームがないので、１番目のフレームを構成する１番目〜Ｍ番目の画素の∂Ｇ／∂ｔには、所定値が割り当てられる。
【００９６】
例えば、２番目のフレームにおいて、１番目の画素の∂Ｇ／∂ｔは、１番目のフレームに含まれる１番目の画素の画素値Ｇと、２番目のフレームに含まれる１番目の画素の画素値Ｇとの差分であり、２番目の画素の∂Ｇ／∂ｔは、１番目のフレームに含まれる２番目の画素の画素値Ｇと、２番目のフレームに含まれる２番目の画素の画素値Ｇとの差分であり、・・・、Ｍ番目の画素の∂Ｇ／∂ｔは、１番目のフレームに含まれるＭ番目の画素の画素値Ｇと、２番目のフレームに含まれるＭ番目の画素の画素値Ｇとの差分である。
【００９７】
図３Ａ及び図１３を参照して、第２の算出部９２は、ステップＳ１１の結果を用いて、１番目〜Ｋ番目のフレームのそれぞれについて、Ｓｘ及びＳｙを算出する（ステップＳ１２）。これについて説明する。図１５は、１番目のフレームのＳｘ及びＳｙを求める処理を説明する説明図である。図１５に示す二つのグラフは、∂Ｇ／∂ｘを示す軸、∂Ｇ／∂ｙを示す軸、及び、∂Ｇ／∂ｔを示す軸を有する。これらの軸は、互いに直交している。∂Ｇ／∂ｘ（画素値Ｇのｘ軸方向の変化を示す値）は、図１０に示す１番目のフレームの∂Ｇ／∂ｘである。∂Ｇ／∂ｙ（画素値Ｇのｙ軸方向の変化を示す値）は、図１０に示す１番目のフレームの∂Ｇ／∂ｙである。∂Ｇ／∂ｔ（画素値Ｇの時間変化を示す値）は、１番目のフレームの∂Ｇ／∂ｔである。
【００９８】
第２の算出部９２は、１番目のフレームを構成する１番目〜Ｍ番目の画素の∂Ｇ／∂ｘ、∂Ｇ／∂ｙ、及び、∂Ｇ／∂ｔを用いて、１番目のフレームを構成する１番目〜Ｍ番目の画素のそれぞれについて、座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｙ，∂Ｇ／∂ｔ）〜座標（∂Ｇ／∂ｘ，∂Ｇ／∂ｙ，∂Ｇ／∂ｔ）を生成し、これらの座標に対する回帰面を算出する。回帰面の算出は、統計処理の一例である。
【００９９】
第２の算出部９２は、∂Ｇ／∂ｘを示す軸と∂Ｇ／∂ｔを示す軸とで規定される平面と、１番目のフレームの回帰面との交線を求め、∂Ｇ／∂ｘを示す軸に対するその交線の傾きを、１番目のフレームのＳｘとする（Ｓｘ１）。第２の算出部９２は、∂Ｇ／∂ｙを示す軸と∂Ｇ／∂ｔを示す軸とで規定される平面と、１番目のフレームの回帰面との交線を求め、∂Ｇ／∂ｙを示す軸に対するその交線の傾きを、１番目のフレームのＳｙとする（Ｓｙ１）。
【０１００】
図１６は、２番目のフレームのＳｘ（Ｓｘ（２））及び２番目のフレームのＳｙ（Ｓｙ（２））〜Ｋ番目のフレームのＳｘ（Ｓｘ（Ｋ））及びＫ番目のフレームのＳｙ（Ｓｙ（Ｋ））を求める処理を説明する説明図である。第２の算出部９２は、２番目〜Ｋ番目のフレームのＳｘ及びＳｙを、１番目のフレームと同様にして求める。
【０１０１】
図３Ａ及び図１３を参照して、第３の算出部９３は、ステップＳ１１の結果及びステップＳ１２の結果を用いて、カメラのブレに起因する信号変化Ｎｓを算出する（ステップＳ１３）。これは、表１を用いて説明した第１実施形態のステップＳ４（図５）と同じなので、説明を省略する。
【０１０２】
図３Ａ及び図１３を参照して、第３の算出部９３は、ステップＳ１３の結果を用いて、カメラのブレに起因する信号変化Ｎｓを積算する処理をする（ステップＳ１４）。これを説明する表を、以下に示す。
【０１０３】
【表２】
2017179510.tif 000004
【０１０４】
表１及び表２を参照して、１番目のフレームの場合、１番目の画素に関するＮｓ積算値は、１番目のフレームに含まれる１番目の画素に関するＮｓであり、２番目の画素に関するＮｓ積算値は、１番目のフレームに含まれる２番目の画素に関するＮｓであり、・・・、Ｍ番目の画素に関するＮｓ積算値は、１番目のフレームに含まれるＭ番目の画素に関するＮｓである。
【０１０５】
２番目のフレームの場合、１番目の画素に関するＮｓ積算値は、１番目のフレームに含まれる１番目の画素に関するＮｓ＋２番目のフレームに含まれる１番目の画素に関するＮｓであり、２番目の画素に関するＮｓ積算値は、１番目のフレームに含まれる２番目の画素に関するＮｓ＋２番目のフレームに含まれる２番目の画素に関するＮｓであり、・・・、Ｍ番目の画素に関するＮｓ積算値は、１番目のフレームに含まれるＭ番目の画素に関するＮｓ＋２番目のフレームに含まれるＭ番目の画素に関するＮｓである。
【０１０６】
Ｋ番目のフレームの場合、１番目の画素に関するＮｓ積算値は、１番目のフレームに含まれる１番目の画素に関するＮｓ＋２番目のフレームに含まれる１番目の画素に関するＮｓ＋・・・＋Ｋ番目のフレームに含まれる１番目の画素に関するＮｓであり、２番目の画素に関するＮｓ積算値は、１番目のフレームに含まれる２番目の画素に関するＮｓ＋２番目のフレームに含まれる２番目の画素に関するＮｓ＋・・・＋Ｋ番目のフレームに含まれる２番目の画素に関するＮｓであり、・・・、Ｍ番目の画素に関するＮｓ積算値は、１番目のフレームに含まれるＭ番目の画素に関するＮｓ＋２番目のフレームに含まれるＭ番目の画素に関するＮｓ＋・・・＋Ｋ番目のフレームに含まれるＭ番目の画素に関するＮｓである。
【０１０７】
図３Ａ及び図１３を参照して、補正部９４は、オリジナルの動画とステップＳ１４の結果との差分を算出する（ステップＳ１５）。これについて詳しく説明する。補正部９４は、オリジナル動画の１番目のフレームを構成する１番目〜Ｍ番目の画素のそれぞれの画素値Ｇと、表２に示す１番目のフレームを構成する１番目〜Ｍ番目の画素に関するそれぞれのＮｓ積算値との差分を算出する。補正部９４は、オリジナル動画の２番目のフレームを構成する１番目〜Ｍ番目の画素のそれぞれの画素値Ｇと、表２に示す２番目のフレームを構成する１番目〜Ｍ番目の画素に関するそれぞれのＮｓ積算値との差分を算出する。補正部９４は、オリジナル動画の３番目〜Ｋ番目のフレームについても同様の処理をする。ステップＳ１５の処理により、オリジナルの動画に対して、カメラのブレを補正することができる。
【０１０８】
以上説明したように、補正部９４は、カメラのブレに起因する信号変化Ｎｓを用いてフレームを補正することを第１の補正とし、動画を構成する１番目〜Ｋ番目のフレームのそれぞれに対して、第１の補正をする。別の表現を用いれば、補正部９４は、第３の算出部９３により算出された信号変化Ｎｓを用いて複数の画像（複数のフレーム）を補正する。
【０１０９】
図３Ａ及び図１３を参照して、表示制御部１０は、ステップＳ１５の差分により得られた動画（すなわち、カメラのブレが補正された動画）を、ディスプレイ１１に表示させる（ステップＳ１６）。
【０１１０】
ガス像だけの動画を生成する場合、画像処理部は、ステップＳ１５の差分により得られた動画に対して、バンドパス処理をする（ステップＳ１７）。この処理は、図５のステップＳ１と同様なので、説明を省略する。表示制御部１０は、ステップＳ１７の処理がされた動画（すなわち、カメラのブレが補正されたガス像の動画）を、ディスプレイ１１に表示させる（ステップＳ１８）。
【０１１１】
第３実施形態について説明する。第３実施形態は、ブレ量が大きい場合に対処することができる。図１を参照して、例えば、現フレームの２番目の画素と３番目の画素とは、画素値Ｇが同じである。このことからに、ブレ量が１画素を超えると、画素値Ｇの変化量とブレ量とは、正比例しなくなり、Ｓｘの誤差が大きくなる（Ｓｙについても同様に誤差が大きくなる）。第３実施形態は、連続する２つのフレーム間のブレ量を１画素以下にしてＮｓを算出する。
【０１１２】
第３実施形態が、第２実施形態と異なるのは、図１３に示すステップＳ１２の替わりに、Ｓｘ及びＳｙの確定処理をすることである。図１７は、第３実施形態において、Ｓｘの確定処理を説明するフローチャートである。図１８は、Ｓｘの確定処理を説明する説明図である。図３Ａ、図１７及び図１８を参照して、第２の算出部９２は、図１３のステップＳ１２で算出した１番目のフレームのＳｘ（Ｓｘ（１））〜Ｋ番目のフレームのＳｘ（Ｓｘ（Ｋ））を、１番目のフレームのＳｘ（Ｓｘ（１））〜Ｋ番目のフレームのＳｘ（Ｓｘ（Ｋ））の初期値とする（ステップＳ２１）。１番目のフレームのＳｘ（Ｓｘ（１））の初期値は、０画素とする。
【０１１３】
第２の算出部９２は、Ｓｘ積算値を算出する（ステップＳ２２）。具体的に説明する。２番目のフレームのＳｘ（Ｓｘ（２））の初期値が、例えば、０．３画素とする。１番目のフレームのＳｘ積算値は、０画素であり、２番目のフレームのＳｘ積算値は、０．３画素（＝０＋０．３）である。第２の算出部９２は、１番目のフレームのＳｘ積算値、及び、２番目のフレームのＳｘ積算値をそれぞれ、整数化する（ステップＳ２３）。これは、Ｓｘ積算値を四捨五入する処理と同様である。１番目のフレームのＳｘ積算値は、０画素なので、整数化すれば、０画素となる。２番目のフレームのＳｘ積算値は、０．３画素なので、整数化すれば、０画素となる。
【０１１４】
第２の算出部９２は、連続する２つのフレームにおいて、整数化された積算値が、いずれも０画素か否かを判定する（ステップＳ２４）。ここでは、連続する２つのフレームとは、１番目と２番目のフレームである。これらのフレームにおいて、整数化された積算値が、いずれも０画素なので（ステップＳ２４でＹｅｓ）、第２の算出部９２は、１番目のフレームのＳｘ（Ｓｘ（１））の初期値（０画素）を、１番目のフレームのＳｘ（Ｓｘ（１））と確定し、２番目のフレームのＳｘ（Ｓｘ（２））の初期値（０．３画素）を、２番目のフレームのＳｘ（Ｓｘ（２））と確定する（ステップＳ２５）。
【０１１５】
第２の算出部９２は、最後のフレームのＳｘを確定したか否かを判断する（ステップＳ２６）。最後のフレームのＳｘが確定していないので（ステップＳ２６でＮｏ）、第２の算出部９２は、次のフレームのＳｘ積算値を算出する（ステップＳ２２）。ここでは、３番目のフレームのＳｘ積算値が算出される。３番目のフレームのＳｘ（Ｓｘ（３））の初期値が、例えば、０．１画素とする。３番目のフレームのＳｘ積算値は、０．４画素（＝０．３＋０．１）となる。
【０１１６】
第２の算出部９２は、３番目のフレームのＳｘ積算値を整数化する（ステップＳ２３）。３番目のフレームのＳｘ積算値は、０．４画素なので、整数化すれば、０画素となる。
【０１１７】
第２の算出部９２は、連続する２つのフレームにおいて、整数化された積算値が、いずれも０画素か否かを判定する（ステップＳ２４）。ここでは、連続する２つのフレームとは、２番目と３番目のフレームである。これらのフレームにおいて、整数化された積算値が、いずれも０画素なので（ステップＳ２４でＹｅｓ）、第２の算出部９２は、３番目のフレームのＳｘ（Ｓｘ（３））の初期値（０．１画素）を、３番目のフレームのＳｘ（Ｓｘ（３））と確定する（ステップＳ２５）。
【０１１８】
第２の算出部９２は、最後のフレームのＳｘを確定したか否かを判断する（ステップＳ２６）。最後のフレームのＳｘが確定していないので（ステップＳ２６でＮｏ）、第２の算出部９２は、次のフレームのＳｘ積算値を算出する（ステップＳ２２）。ここでは、４番目のフレームのＳｘ積算値が算出される。４番目のフレームのＳｘ（Ｓｘ（４））の初期値が、例えば、０．７画素とする。４番目のフレームのＳｘ積算値は、１．１画素（＝０．４＋０．７）となる。
【０１１９】
第２の算出部９２は、４番目のフレームのＳｘ積算値を整数化する（ステップＳ２３）。４番目のフレームのＳｘ積算値は、１．１画素なので、整数化すれば、１画素となる。第２の算出部９２は、連続する２つのフレームにおいて、整数化された積算値が、いずれも０画素か否かを判定する（ステップＳ２４）。ここでは、連続する２つのフレームとは、３番目と４番目のフレームである。４番目のフレームについて、整数化されたＳｘ積算値が、１画素である（ステップＳ２４でＮｏ）。よって、第２の算出部９２は、３番目のフレームと、ｘ軸方向に−１画素ずらした４番目のフレームとを用いて、４番目のフレームのＳｘ（Ｓｘ（４））を算出する（ステップＳ２７）。この例では、Ｓｘ（４）は、−０．３画素となる。第２の算出部９２は、４番目のフレームのＳｘ（Ｓｘ（４））を−０．３画素と確定する（ステップＳ２８）。
【０１２０】
第２の算出部９２は、最後のフレームのＳｘを確定したか否かを判断する（ステップＳ２６）。最後のフレームのＳｘが確定していないので（ステップＳ２６でＮｏ）、第２の算出部９２は、次のフレームのＳｘ積算値を算出する（ステップＳ２２）。ここでは、５番目のフレームのＳｘ積算値が算出される。５番目のフレームのＳｘ（Ｓｘ（５））の初期値が、例えば、−０．８画素とする。５番目のフレームのＳｘ積算値は、０．３画素（＝１．１−０．８）となる。
【０１２１】
第２の算出部９２は、５番目のフレームのＳｘ積算値を整数化する（ステップＳ２３）。５番目のフレームのＳｘ積算値は、０．３画素なので、整数化すれば、０画素となる。第２の算出部９２は、連続する２つのフレームにおいて、整数化された積算値が、いずれも０画素か否かを判定する（ステップＳ２４）。ここでは、連続する２つのフレームとは、４番目と５番目のフレームである。４番目のフレームについて、整数化されたＳｘ積算値が、１画素である（ステップＳ２４でＮｏ）。よって、第２の算出部９２は、５番目のフレームと、ｘ軸方向に−１画素ずらした４番目のフレームとを用いて、５番目のフレームのＳｘ（Ｓｘ（５））を算出する（ステップＳ２７）。この例では、Ｓｘ（５）は、０．２画素となる。第２の算出部９２は、５番目のフレームのＳｘ（Ｓｘ（５））を０．２画素と確定する（ステップＳ２８）。
【０１２２】
第２の算出部９２は、６番目以降のフレームに対して、上記説明と同様の処理をする。第２の算出部９２は、Ｋ番目のフレームのＳｘを確定したと判断したとき（ステップＳ２６でＹｅｓ）、Ｓｘの確定処理が終了する。
【０１２３】
第２の算出部９２は、Ｓｙについても、Ｓｘと同様にして確定処理をする。画像処理部９は、確定処理がされたＳｘ及びＳｙを用いて、図１３に示すステップＳ１３以降の処理をする。
【０１２４】
第４実施形態について説明する。図３Ａに示す赤外線カメラ２による動画の撮影中に、カメラのブレの発生に加えて、被写体（例えば、自動車、人、木の枝等）の動きが発生することがある。第１実施形態〜第３実施形態は、フレームを構成する複数の画素のそれぞれの空間変化を示す値及び時間変化を示す値を統計処理して、フレームのＳｘ及びＳｙを算出する。この場合に、動く被写体の像を示す領域を構成する画素を含めて、フレームのＳｘ及びＳｙが算出されれば、これらの値の精度が低下する。そこで、第４実施形態は、被写体の動きが発生しても、フレームのＳｘ及びＳｙの精度が低下しないようにする。
【０１２５】
動く被写体は、動画において、移動体の像として表れる。二以上の領域に分割されたフレームにおいて、移動体の像を示す領域のブレを示す値は、これ以外の領域のブレを示す値と異なる。そこで、第４実施形態は、二以上の領域に分割されたフレームにおいて、各領域のブレを示す値を用いて、フレームに含まれる移動体の像を示す領域を特定し、移動体の像を示す領域を除く二以上の領域のブレを示す値を用いて、カメラのブレにより生じる画像のブレを示す値を算出する。
【０１２６】
図１９は、二以上の領域Ｒに分割された一つのフレームの模式図である。二以上の領域Ｒのそれぞれに含まれる矢印は、その領域Ｒのブレ量の大きさ及び方向を示す。二以上の領域Ｒの中に、ブレ量の大きさ及び方向が他の領域Ｒと異なる領域Ｒが存在する場合がある。例えば、図１９では、領域Ｒａ、領域Ｒｂ、領域Ｒｃ、領域Ｒｄ、領域Ｒｅ、領域Ｒｆのブレ量の大きさ及び方向が、他の領域Ｒのブレ量の大きさ及び方向と異なる。第４実施形態は、領域Ｒａ〜Ｒｆを移動体の像を示す領域とし、一つのフレームのブレ量の算出に、領域Ｒａ〜Ｒｆのような領域Ｒを除外する。
【０１２７】
図２０は、第４実施形態において、Ｓｘ，Ｓｙを算出する処理を説明するフローチャートである。第２の算出部９２は、ｉ＝１に設定する（ステップＳ４１）。これは、１番目のフレームに対する処理を意味する。第１実施形態〜第３実施形態では、一つのフレームの全体でＳｘ，Ｓｙを算出したが（例えば、図１３のステップＳ１２）、第４実施形態の第２の算出部９２は、１番目のフレームを、二以上の領域Ｒに分割し、二以上の領域Ｒのそれぞれについて、Ｓｘ，Ｓｙを算出する（ステップＳ４２）。各領域ＲのＳｘ，Ｓｙを算出する方法は、一つのフレームのＳｘ，Ｓｙを算出する方法と同じである（例えば、図１３のステップＳ１２）。すなわち、第２の算出部９２は、フレームのＳｘ，Ｓｙの算出方法と同様の方法を用いて、二以上の領域ＲのそれぞれのＳｘ，Ｓｙを算出する。
【０１２８】
第２の算出部９２は、ステップＳ４２で算出したベクトル（Ｓｘ，Ｓｙ）の平均値を算出する（ステップＳ４３）。
【０１２９】
第２の算出部９２は、ステップＳ４２で算出された複数の領域Ｒのそれぞれのベクトル（Ｓｘ，Ｓｙ）の中で、ステップＳ４３で算出されたベクトル（Ｓｘ，Ｓｙ）の平均値と最も離れているベクトル（Ｓｘ，Ｓｙ）を探索する（ステップＳ４４）。
【０１３０】
第２の算出部９２は、ステップＳ４４で探索されたベクトル（Ｓｘ，Ｓｙ）とステップＳ４３で算出された平均値との差の絶対値が、予め定められた値以上か否かを判断する（ステップＳ４５）。
【０１３１】
第２の算出部９２が、ステップＳ４４で探索されたベクトル（Ｓｘ，Ｓｙ）とステップＳ４３で算出された平均値との差が、予め定められた値以上と判断したとき（ステップＳ４５でＹｅｓ）、ステップＳ４４で探索されたベクトル（Ｓｘ，Ｓｙ）を有する領域Ｒを、移動体の像を示す領域として特定し、ステップＳ４４で探索されたベクトル（Ｓｘ，Ｓｙ）を、平均値の算出に用いない決定をする（ステップＳ４６）。すなわち、第２の算出部９２は、二以上の領域Ｒのそれぞれのベクトル（Ｓｘ，Ｓｙ）を用いて、フレームに含まれる移動体の像を示す領域Ｒを特定する。そして、第２の算出部９２は、ステップＳ４３の処理に戻る。ここでは、第２の算出部９２は、ステップＳ４２で算出されたベクトル（Ｓｘ，Ｓｙ）の中から、ステップＳ４６で決定されたベクトル（Ｓｘ，Ｓｙ）を除外して、ベクトル（Ｓｘ，Ｓｙ）の平均値を算出する（ステップＳ４３）。そして、第２の算出部９２は、ステップＳ４４以降の処理をする。
【０１３２】
第２の算出部９２が、ステップＳ４４で探索されたベクトル（Ｓｘ，Ｓｙ）とステップＳ４３で算出された平均値との差が、予め定められた値より小さいと判断したとき（ステップＳ４５でＮｏ）、第２の算出部９２は、ｉ＝ｋか否かを判断する（ステップＳ４７）。ここでは、ｉ＝１なので（ステップＳ４７でＮｏ）、第２の算出部９２は、ステップＳ４３で算出された直近の平均値を、ｉ番目のフレームのベクトル（Ｓｘ，Ｓｙ）と確定する（ステップＳ４８）。ここでは、１番目のフレームのベクトル（Ｓｘ，Ｓｙ）、すなわち、ベクトル（Ｓｘ（１），Ｓｙ（１））と確定される。すなわち、第２の算出部９２は、移動体の像を示す領域Ｒを除く二以上の領域Ｒのベクトル（Ｓｘ，Ｓｙ）を用いて（ベクトル（Ｓｘ，Ｓｙ）の平均値）、フレームのベクトル（Ｓｘ，Ｓｙ）を算出する。
【０１３３】
第２の算出部９２は、ｉ＝ｉ＋１に設定する（ステップＳ４９）。ここでは、ｉ＝２に設定される。そして、第２の算出部９２は、ステップＳ４２の処理に戻る。ここでは、第２の算出部９２は、２番目のフレームを複数の領域Ｒに分割し、複数の領域Ｒのそれぞれについて、ベクトル（Ｓｘ，Ｓｙ）を算出する（ステップＳ４２）。
【０１３４】
第２の算出部９２は、ｉ＝ｋと判断したとき（ステップＳ４７でＹｅｓ）、第２の算出部９２は、ステップＳ４３で算出された直近の平均値を、ｉ番目のフレームのベクトル（Ｓｘ，Ｓｙ）と確定する。ここでは、Ｋ番目のフレームのベクトル（Ｓｘ，Ｓｙ）、すなわち、ベクトル（Ｓｘ（Ｋ），Ｓｙ（Ｋ））と確定される。
【０１３５】
以上が第４実施形態によるフレームのベクトル（Ｓｘ，Ｓｙ）の算出処理である。これは、第１実施形態〜第３実施形態のいずれのベクトル（Ｓｘ，Ｓｙ）の算出処理に置き換えることができる。
【０１３６】
本実施形態は、図１１及び図１２、並びに、図１５及び図１６で説明したように、回帰面を求め、回帰面を用いてＳｘ及びＳｙを求めている。これは、以下の処理をしていることを意味する。
【０１３７】
図６を参照して、時系列画素データから、背景の温度変化を示す信号Ｂ及び高周波ノイズＮｈを引いた信号をｇとする。ｇは、以下の式で表される。
【０１３８】
ｇ＝Ｎｓ＋ＮＬ＋Ｎｃ
【０１３９】
上述したように、Ｎｓは、カメラのブレに起因する信号変化であり、ＮＬは、漏れたガスによる温度変化を示す信号であり、Ｎｃは、その他のノイズである。
【０１４０】
上記ｇの式は、以下の式に変形できる。
【０１４１】
ＮＬ＋Ｎｃ＝ｇ−Ｎｓ
【０１４２】
上述したように、カメラのブレがｘ軸方向及びｙ軸方向に発生した場合、カメラのブレに起因する信号変化Ｎｓは、式３で表されると推測できる。カメラのブレに起因する信号変化Ｎｓは、信号ｇ、漏れたガスによる温度変化ＮＬ、及び、その他のノイズＮｃのそれぞれに対して、独立であり、相関関係を有さない。一つのフレーム（一つの画像）では、Ｓｘ及びＳｙがそれぞれ一定と仮定する。様々な値のＳｘ及びＳｙの中で、信号値（ＮＬ＋Ｎｃ＝ｇ−Ｎｓ）が最小となるＳｘ及びＳｙの組み合わせが、Ｓｘ及びＳｙの正解値と思われる。これは、以下の式４で表される信号値の二乗ノルムが最小となるＳｘ及びＳｙの組み合わせを求めればよいことを意味する。
【０１４３】
【数１】
2017179510.tif 000005
【０１４４】
式４を計算すると、式５となる。
【０１４５】
【数２】
2017179510.tif 000006
【０１４６】
信号値の二乗ノルムが最小となる点では、式５を、Ｓｘ及びＳｙで偏微分した値が、ともに０になると思われる。式５を、Ｓｘで偏微分する式は、式６となり、式５を、Ｓｙで偏微分する式は、式７となる。
【０１４７】
【数３】
2017179510.tif 000007
【０１４８】
式６を計算（すなわち、式６をＳｘで偏微分）すると、式８となり、式７を計算（すなわち、式７をＳｙで偏微分）すると、式９となる。
【０１４９】
【数４】
2017179510.tif 000008
【０１５０】
式８及び式９から、Ｓｘは式１０となり、Ｓｙは式１１となる。
【０１５１】
【数５】
2017179510.tif 000009
【０１５２】
第１実施形態〜第４実施形態では、式３を用いて、カメラのブレに起因する信号変化Ｎｓを算出している。これは、カメラのブレがｘ軸方向及びｙ軸方向に発生した場合のＮｓである。カメラのブレが回転ブレ等の場合のＮｓについて説明する。
【０１５３】
カメラのブレが回転ブレの場合、カメラのブレに起因する信号変化Ｎｓは、式１及び式２の場合と同様に考えると、以下の式１２で表されると推測できる。
【０１５４】
Ｎｓ＝∂Ｇ／∂ｔ＝Ｓω・∂Ｇ／∂ω＝Ｓω（ｙ・∂Ｇ／∂ｘ＋ｘ・∂Ｇ／∂ｙ）・・・（式１２）
【０１５５】
Ｓωは、カメラのブレが回転ブレの場合に生じるフレームの回転ブレを示す値である。このような値として、例えば、回転軸方向のフレームのブレ量がある。∂Ｇ／∂ωは、画素値Ｇの回転方向の変化を示す値を示す。式４から、∂Ｇ／∂ｘ及び∂Ｇ／∂ｙが求まれば、∂Ｇ／∂ωを求めることなく、Ｎｓが求まる。
【０１５６】
式１２によれば、∂Ｇ／∂ω＝ｙ・∂Ｇ／∂ｘ＋ｘ・∂Ｇ／∂ｙである。これが成立することを説明する。図２１は、これらを説明する説明図である。カメラが原点を中心として、角度Δωの回転をしたとき、座標（ｘ１，ｙ１）に位置する画素のｘ座標値の変化量Δｘ、及び、ｙ座標値の変化量Δｙは、以下の通りである。
【０１５７】
Δｘ≒ｙ１・Δω
Δｙ≒ｘ１・Δω
【０１５８】
上記回転による座標（ｘ１，ｙ１）に位置する画素の画素値Ｇの変化量ΔＧは、以下の通りである。
【０１５９】
ΔＧ≒∂Ｇ／∂ｘ・Δｘ＋∂Ｇ／∂ｙ・Δｙ
≒∂Ｇ／∂ｘ・ｙ１・Δω＋∂Ｇ／∂ｙ・ｘ１・Δω
【０１６０】
ΔＧをΔωで割った式は、以下の通りである。
【０１６１】
ΔＧ／Δω≒ｙ１・∂Ｇ／∂ｘ・＋ｘ１・∂Ｇ／∂ｙ
【０１６２】
この式は、Δωが十分に小さい場合に成立する。従って、∂Ｇ／∂ω＝ｙ・∂Ｇ／∂ｘ＋ｘ・∂Ｇ／∂ｙが成立することになる。
【０１６３】
拡大縮小ブレ（フォーカス時等でレンズが駆動されるときに生じる像の大きさの変動、又は、被写体までの距離が変化するときに生じる像の大きさの変動）の場合、カメラのブレに起因する信号変化Ｎｓは、式１及び式２の場合と同様に考えると、以下の式１３で表されると推測できる。
【０１６４】
Ｎｓ＝∂Ｇ／∂ｔ＝Ｓｍ・∂Ｇ／∂ｍ＝Ｓｍ（ｘ・∂Ｇ／∂ｘ＋ｙ・∂Ｇ／∂ｙ）・・・（式１３）
【０１６５】
Ｓｍは、拡大縮小ブレの場合に生じるフレームの拡大縮小ブレを示す値である。∂Ｇ／∂ｍは、拡大縮小ブレによる画素値Ｇの変化を示す値を示す。式１３から、∂Ｇ／∂ｘ、及び、∂Ｇ／∂ｙが求まれば、∂Ｇ／∂ｍを求めることなく、Ｎｓが求まる。
【０１６６】
式１３によれば、∂Ｇ／∂ｍ＝ｘ・∂Ｇ／∂ｘ＋ｙ・∂Ｇ／∂ｙである。これが成立することを説明する。図２２は、これらを説明する説明図である。カメラのズームが原点を中心として、（１＋Δｍ）の拡大をしたとき、座標（ｘ１，ｙ１）に位置する画素のｘ座標値の変化量Δｘ、及び、ｙ座標値の変化量Δｙは、以下の通りである。
【０１６７】
Δｘ＝ｘ１・Δｍ
Δｙ＝ｙ１・Δｍ
【０１６８】
上記拡大による座標（ｘ１，ｙ１）に位置する画素の画素値Ｇの変化量ΔＧは、以下の通りである。
【０１６９】
ΔＧ≒∂Ｇ／∂ｘ・Δｘ＋∂Ｇ／∂ｙ・Δｙ
≒∂Ｇ／∂ｘ・ｘ１・Δｍ＋∂Ｇ／∂ｙ・ｙ１・Δｍ
【０１７０】
ΔＧをΔｍで割った式は、以下の通りである。
【０１７１】
ΔＧ／Δｍ≒ｘ１・∂Ｇ／∂ｘ・＋ｙ１・∂Ｇ／∂ｙ
【０１７２】
この式は、Δｍが十分に小さい場合に成立する。従って、∂Ｇ／∂ｍ＝ｘ・∂Ｇ／∂ｘ＋ｙ・∂Ｇ／∂ｙが成立することになる。
【０１７３】
図３Ａを参照して、第１実施形態〜第４実施形態では、説明を簡単にするために、二次元イメージセンサー６が受光した光量に対して、二次元イメージセンサー６の出力がリニアであることを前提として説明した。リニアでない場合、以下の２つの方法のいずれかで対処することができる。一つ目の対処方法は、信号処理部７が、二次元イメージセンサー６の出力をリニアに補正した後、動画データＤ１を生成する。二つ目の対処方法は、第２の算出部９２が、統計処理として回帰曲面を算出し、この回帰曲面を用いて、Ｓｘ及びＳｙを算出する。
【０１７４】
図３Ａを参照して、第１実施形態〜第４実施形態では、画像処理部９が、動画に対して所定の画像処理をし、表示制御部１０が、その画像処理がされた動画をディスプレイ１１に表示させている。本発明は、この構成に限定されず、画像処理部９を備えるが、表示制御部１０及びディスプレイ１１を備えない構成でもよいし、画像処理部９及び表示制御部１０を備えるが、ディスプレイ１１を備えない構成でもよい。
【０１７５】
第１実施形態〜第４実施形態では、カメラとして赤外線カメラ２を例に説明したが、本発明は、可視光カメラについても適用することができる。
【０１７６】
（実施形態の纏め）
上記目的を達成する本実施形態の第１の局面に係る画像処理装置は、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする装置であって、前記画像を構成する複数の画素のそれぞれについて、前記画素の画素データの空間変化を示す値、及び、前記画素の画素データの時間変化を示す値を算出することを第１の算出とし、複数の前記画像のそれぞれに対して、前記第１の算出をする第１の算出部と、前記画像を構成する複数の前記画素のそれぞれについて算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記画像のブレを示す値を算出することを第２の算出とし、複数の前記画像のそれぞれに対して、前記第２の算出をする第２の算出部と、を備える。別の表現を用いると、本実施形態の第１の局面に係る画像処理装置は、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする装置であって、前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１の算出部と、前記第１の算出部によって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２の算出部と、を備える。
【０１７７】
カメラのブレが発生すれば、画素データの空間方向の変化と画素データの時間方向の変化とが発生する。カメラのブレ量が１画素以内であれば、画素データの変化量とブレ量とは、正比例すると考えられる。本発明者は、これらを前提にして、カメラによって撮影された時系列に並ぶ複数の画像について、画像を構成する複数の画素のそれぞれの空間変化を示す値、及び、時間変化を示す値に対して、所定の統計処理をすれば、カメラのブレにより生じる画像のブレを示す値が求まることを見出した。そこで、本実施形態の第１の局面に係る画像処理装置において、第１の算出部は、画像を構成する複数の画素のそれぞれについて、画素の画素データの空間変化を示す値、及び、画素の画素データの時間変化を示す値を算出し、第２の算出部は、画像を構成する複数の画素のそれぞれについて、第１の算出部によって算出された空間変化を示す値、及び、時間変化を示す値に対して、所定の統計処理をする。従って、本実施形態の第１の局面に係る画像処理装置によれば、カメラのブレにより生じる画像のブレを示す値を求めることができる。上記別の表現の場合についても、同様のことが言える。
【０１７８】
画素データとは、画素が有する数値情報であり、例えば、画素の画素値（輝度値）である。赤外画像の場合、画素値、又は、画素値を用いて算出された温度を意味する。
【０１７９】
空間変化を示す値として、画素の画素データに対して空間変化の方向の方向微分をした値がある。
【０１８０】
時間変化を示す値として、画素の画素データに対して時間変化の方向の方向微分をした値がある。また、時間変化を示す値として、複数の画素のそれぞれの画素データに含まれる交流成分がある。
【０１８１】
上記構成において、前記統計処理は、前記空間変化の方向を示す軸と前記時間変化の方向を示す軸とで構成される座標系において、前記画像を構成する複数の画素のそれぞれについて、前記空間変化を示す値と前記時間変化を示す値とにより座標を構成する処理と、前記座標系での前記座標の分布を近似する回帰式を求める処理と、含む。
【０１８２】
この構成は、統計処理の一例である。回帰式で示される図形は、回帰直線、回帰曲線、回帰平面や回帰曲面がある。
【０１８３】
上記構成において、前記第２の算出部は、前記回帰式で示される図形が回帰平面の場合、前記空間変化の方向を示す軸を、第１の変化方向を示す軸、及び、前記第１の変化方向と直交する第２の変化方向を示す軸とし、前記第１の変化方向を示す軸と前記時間変化の方向を示す軸とで規定される平面と、前記回帰平面との交線を第１の交線とし、前記第２の変化方向を示す軸と前記時間変化の方向を示す軸とで規定される平面と、前記回帰平面との交線を第２の交線とし、前記第１の変化方向を示す軸に対する前記第１の交線の傾きを、前記第１の変化方向に沿った前記画像のブレを示す値とし、前記第２の変化方向を示す軸に対する前記第２の交線の傾きを、前記第２の変化方向に沿った前記画像のブレを示す値とする。別の表現を用いると、前記第２の算出部は、前記回帰式で示される図形が回帰平面の場合、前記空間変化の方向を示す軸と前記時間変化の方向を示す軸とで規定される平面と、前記回帰平面との交線の、前記空間変化の方向を示す軸に対する傾きを算出する。
【０１８４】
回帰式で示される図形が回帰平面の場合、カメラのブレにより生じる画像のブレを示す値は、第１の変化方向に沿った画像のブレを示す値と、第２の変化方向に沿った画像のブレを示す値とに分けられる。第１の変化方向に沿った画像のブレを示す値は、第１の変化方向を示す軸に対する第１の交線の傾きである。第２の変化方向に沿った画像のブレを示す値は、第２の変化方向を示す軸に対する第２の交線の傾きである。上記別の表現の場合についても、同様のことが言える。
【０１８５】
上記構成において、前記画像を構成する複数の前記画素のそれぞれについて、前記空間変化を示す値、及び、前記画像のブレを示す値を用いて、前記カメラのブレに起因する信号変化を算出することを第３の算出とし、複数の前記画像のそれぞれに対して、前記第３の算出をする第３の算出部をさらに備える。別の表現を用いると、前記第３の算出部は、前記複数の画像の各画素について、前記空間変化を示す値、及び、前記画像のブレを示す値を用いて、前記カメラのブレに起因する信号変化を算出する。
【０１８６】
この構成によれば、カメラのブレに起因する信号変化、すなわち、カメラのブレによるノイズを算出することができる。上記別の表現の場合についても、同様のことが言える。
【０１８７】
上記構成において、前記信号変化を用いて前記画像を補正することを第１の補正とし、複数の前記画像のそれぞれに対して、前記第１の補正をする補正部をさらに備える。別の表現を用いると、前記補正部は、前記第３の算出部により算出された信号変化を用いて前記複数の画像を補正する。
【０１８８】
この構成によれば、時系列に並ぶ複数の画像（例えば、ガス漏れの監視対象の画像）に対して、カメラのブレを補正することができる。上記別の表現の場合についても、同様のことが言える。
【０１８９】
上記構成において、前記画像に含まれる予め定められた被写体を抽出することを被写体抽出とし、複数の前記画像のそれぞれに対して、前記被写体抽出をする抽出部と、前記信号変化を用いて前記被写体を補正することを第２の補正とし、複数の前記画像のそれぞれから抽出された前記被写体に対して、前記第２の補正をする補正部と、をさらに備える。別の表現を用いると、抽出部は、前記複数の画像のそれぞれから、予め定められた被写体を抽出する。補正部は、前記抽出部により抽出された被写体を、前記第３の算出部により算出された信号変化を用いて補正する。
【０１９０】
この構成によれば、時系列に並ぶ複数の画像に含まれる被写体（例えば、ガス像）を抽出し、その被写体に対して、カメラのブレを補正することができる。上記別の表現の場合についても、同様のことが言える。
【０１９１】
上記構成において、前記第２の算出部は、前記画像を二以上の領域に分割し、二以上の前記領域のそれぞれについて、前記第２の算出（前記統計処理）をすることにより、前記カメラのブレにより生じる二以上の前記領域のそれぞれのブレを示す値を算出し、前記第２の算出部は、二以上の前記領域のそれぞれの前記ブレを示す値を用いて、前記画像に含まれる移動体の像を示す前記領域を特定し、前記第２の算出部は、前記移動体の像を示す前記領域を除く前記領域の前記ブレを示す値を用いて、前記カメラのブレにより生じる前記画像のブレを示す値を算出する。
【０１９２】
カメラによる画像の撮影中に、カメラのブレの発生に加えて、被写体（例えば、自動車、人、木の枝等）の動きが発生することがある。本実施形態の第１の局面に係る画像処理装置は、画像を構成する複数の画素のそれぞれの空間変化を示す値及び時間変化を示す値を統計処理して、カメラのブレにより生じる画像のブレを示す値を算出する。この場合に、動く被写体の像を示す領域を構成する画素を含めて、カメラのブレにより生じる画像のブレを示す値が算出されれば、この値の精度が低下する。
【０１９３】
動く被写体は、時系列に並べられた複数の画像において、移動体の像として表れる。二以上の領域に分割された画像において、移動体の像を示す領域のブレを示す値は、これ以外の領域のブレを示す値と異なる。そこで、この構成は、二以上の領域に分割された画像において、各領域のブレを示す値を用いて、画像に含まれる移動体の像を示す領域を特定し、移動体の像を示す領域を除く領域のブレを示す値を用いて、カメラのブレにより生じる画像のブレを示す値を算出する。
【０１９４】
本実施形態の第２の局面に係る画像処理方法は、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする方法であって、前記画像を構成する複数の画素のそれぞれについて、前記画素の画素データの空間変化を示す値、及び、前記画素の画素データの時間変化を示す値を算出することを第１の算出とし、複数の前記画像のそれぞれに対して、前記第１の算出をする第１のステップと、前記第１のステップによって、前記画像を構成する複数の前記画素のそれぞれについて算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記画像のブレを示す値を算出することを第２の算出とし、複数の前記画像のそれぞれに対して、前記第２の算出をする第２のステップと、を備える。別の表現を用いると、本実施形態の第２の局面に係る画像処理方法は、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をする方法であって、前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１のステップと、前記第１のステップによって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２のステップと、を備える。
【０１９５】
本実施形態の第２の局面に係る画像処理方法は、本実施形態の第１の局面に係る画像処理装置を方法の観点から規定しており、本実施形態の第１の局面に係る画像処理装置と同様の作用効果を有する。上記別の表現の場合についても、同様のことが言える。
【０１９６】
本実施形態の第３の局面に係る画像処理プログラムは、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をするプログラムであって、前記画像を構成する複数の画素のそれぞれについて、前記画素の画素データの空間変化を示す値、及び、前記画素の画素データの時間変化を示す値を算出することを第１の算出とし、複数の前記画像のそれぞれに対して、前記第１の算出をする第１のステップと、前記第１のステップによって、前記画像を構成する複数の前記画素のそれぞれについて算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記画像のブレを示す値を算出することを第２の算出とし、複数の前記画像のそれぞれに対して、前記第２の算出をする第２のステップと、をコンピュータに実行させる。別の表現を用いると、本実施形態の第３の局面に係る画像処理プログラムは、カメラによって撮影された時系列に並ぶ複数の画像に対して画像処理をするプログラムであって、前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出する第１のステップと、前記第１のステップによって算出された前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出する第２のステップと、をコンピュータに実行させる。
【０１９７】
本実施形態の第３の局面に係る画像処理プログラムは、本実施形態の第１の局面に係る画像処理装置をプログラムの観点から規定しており、本実施形態の第１の局面に係る画像処理装置と同様の作用効果を有する。上記別の表現の場合についても、同様のことが言える。
【０１９８】
本実施形態の第４の局面に係るガス検知システムは、監視対象の動画を撮影し、複数の画像が時系列に並べられた動画データを出力する赤外線カメラと、前記赤外線カメラから入力された前記複数の画像について、画素データの空間変化を示す値、及び、前記画素データの時間変化を示す値を算出し、前記空間変化を示す値、及び、前記時間変化を示す値を統計処理して、前記赤外線カメラのブレにより生じる前記複数の画像のそれぞれのブレを示す値を算出し、前記複数の画像のそれぞれのブレを示す値に基づいて前記複数の画像のそれぞれを補正する画像処理部と、前記画像処理部によって補正された前記複数の画像が時系列に並べられた動画を表示するディスプレイと、を備える。
【０１９９】
本実施形態の第４の局面に係るガス検知システムは、本実施形態の第１の局面に係る画像処理装置をガス検知システムの観点から規定しており、本実施形態の第１の局面に係る画像処理装置と同様の作用効果を有する。
【０２００】
この出願は、２０１６年４月１１日に出願された日本国特許出願特願２０１６−０７８５７８を基礎とするものであり、その内容は、本願に含まれるものである。
【０２０１】
本発明を表現するために、上述において図面を参照しながら実施形態を通して本発明を適切且つ十分に説明したが、当業者であれば上述の実施形態を変更および／または改良することは容易に為し得ることであると認識すべきである。したがって、当業者が実施する変更形態または改良形態が、請求の範囲に記載された請求項の権利範囲を離脱するレベルのものでない限り、当該変更形態または当該改良形態は、当該請求項の権利範囲に包括されると解釈される。
【産業上の利用可能性】
【０２０２】
本発明によれば、画像処理装置、画像処理方法、画像処理プログラム、及び、ガス検知システムを提供することができる。
【図１】
2017179510.tif 000010
【図２】
2017179510.tif 000011
【図３Ａ】
2017179510.tif 000012
【図３Ｂ】
2017179510.tif 000013
【図４】
2017179510.tif 000014
【図５】
2017179510.tif 000015
【図６】
2017179510.tif 000016
【図７】
2017179510.tif 000017
【図８】
2017179510.tif 000018
【図９】
2017179510.tif 000019
【図１０】
2017179510.tif 000020
【図１１】
2017179510.tif 000021
【図１２】
2017179510.tif 000022
【図１３】
2017179510.tif 000023
【図１４】
2017179510.tif 000024
【図１５】
2017179510.tif 000025
【図１６】
2017179510.tif 000026
【図１７】
2017179510.tif 000027
【図１８】
2017179510.tif 000028
【図１９】
2017179510.tif 000029
【図２０】
2017179510.tif 000030
【図２１】
2017179510.tif 000031
【図２２】
2017179510.tif 000032
【国際調査報告】
2017179510.tif 000033
2017179510.tif 000034
2017179510.tif 000035
2017179510.tif 000036