JP2004139219A

JP2004139219A - Image processing method and image processor

Info

Publication number: JP2004139219A
Application number: JP2002301470A
Authority: JP
Inventors: Koji Kakiuchi; 垣内　宏司
Original assignee: Seiko Instruments Inc
Current assignee: Seiko Instruments Inc
Priority date: 2002-10-16
Filing date: 2002-10-16
Publication date: 2004-05-13

Abstract

<P>PROBLEM TO BE SOLVED: To provide a composition processing method for combining continuously photographed digital images, and for compositing those digital images into one image by solving the problem that the luminance of the whole region of a screen is uncorrectable, or that an image is destroyed in luminance unevenness correcting method. <P>SOLUTION: This image processing method for continuously combining the plurality of digital images photographed while the photographic places or photographic directions are continuously moved so that a portion of the photographed digital images can be overlapped, and for compositing those digital images into one panorama image is characterized to express luminance unevenness correction coefficients as the linear combination of the known functions of image face coordinates, and to calculate linear combination constants by a statistical means. <P>COPYRIGHT: (C)2004,JPO

Description

【０００１】
【発明の属する技術分野】
本発明は、ビデオカメラによって取り込んだ連続画像や複数の静止画像をつなぎ合わせて一つのパノラマ画像として画像を合成処理する画像処理方法、及びこの画像処理方法を用いた画像処理装置に関する。
【０００２】
【従来の技術】
デジタルビデオの普及とパーソナルコンピュータの高性能化に伴って、かつては不可能であったデジタル画像データの大量保存や大量処理が可能になってきた。これに伴い、デジタル画像データを処理する新しい方法が行われるようになってきた。その一つとして、移動しながら撮影したビデオ画像や、同じ領域を場所や角度を変えて撮影した複数枚の静止画をつなぎ合わせて一つの合成画像（パノラマ画像）を生成する技術がある。これをモザイキングと呼ぶこともある。
【０００３】
この技術はトンネルや壁のひび割れなどの検査を行う場合に有効に使われている。画像処理装置においてこのモザイキング処理を行うときには、まず各画像が相互に重なる重複領域を求め、次にこの重複領域を「のりしろ」として画像を貼り付けるという処理を行っている。
【０００４】
このような画像合成を行う場合、画像毎の照明状態の変化などによって、合成画像の画素の輝度が画像のつなぎ目に相当する部分で急激に変化し、そのために合成画像内に濃淡および色調の境目ができることがある。
【０００５】
図１１は従来のモザイキング処理における撮影状態を説明する図である。図１１における領域１０１が一つの画像として撮影され、同様に領域１０２が一つ画像として撮影され、さらに領域１０３が一つ画像として撮影されたとする。境界線１１１から境界線１１２が領域１０１の画像の範囲となり、境界線１１３から境界線１１４が領域１０２の画像の範囲となり、境界線１１５から境界線１１６が領域１０３の画像の範囲となる。図１２は照度むらがあるときの図１１における各画像の輝度分布状態を説明するための図である。
【０００６】
図１２では、領域１０１を撮影した画像の輝度分布１２１に示したが、このような分布になる理由は、人工照明されている場合には、中央が明るくて端が暗くなるためである。同じ理由によって領域１０２、領域１０３を撮影した画像における輝度分布１２２、１２３においても同じような輝度分布が発生する。このとき重なり領域１０４を「のりしろ」として二つの領域１０１と１０２の画像を貼り付ける。ここで、領域１０１の画像の中で図１２の１１１から１１３の部分を切り取り、この切り取った画像を領域１０２の画像の左に貼り付ける。次に領域１０２の画像から図１２の１１３から１１５の領域を切り取って、領域１０３の画像の左に貼り付ける。
【０００７】
図１３はこのような従来の方法により図１２の画像をモザイキングしたときの結果を示す図である。ここで説明したような画面の貼り付けを行うと、図１３の境界線１１３と境界線１１５で輝度の不連続的な「飛び」１３１と１３２がここで発生し、パノラマ画像に不自然な線として見える。
【０００８】
この問題を解決するための一つの方法として、異なる画像間において重なっている領域を何らかの方法で検出し、そこでの輝度の違いを用いて照度むら効果を補正しようとする方法がある。即ち、重なっている領域は本来同じ場所であるから輝度は同じはずであるが、照度むら効果により画像に輝度が異なっている。そこでその領域での輝度を比較して同じになるように補正すればよいというのが基本的な考え方である。
【０００９】
このような考え方に基づく方法として、ヒストグラム法・線形濃度変換法・平均濃度差補正法などがある。ヒストグラム法と線形濃度変換法では画像全体の輝度補正はできるが、場所に応じた細かい輝度補正ができないという欠点がある。
【００１０】
人間の目は不連続的な変化には非常に敏感であるため、これらの方法では目立たない程度にまで輝度の不連続さをなくしてしまうことが難しい。一方、平均濃度差補正法では場所に応じた細かい補正が可能となるものの、重なっていない領域を含めた画像全体の補正はできないという欠点がある（例えば、非特許文献１参照）。
【００１１】
また別の種類の考え方として、画像の輝度データに平滑化フィルターを掛けて照度の画像内分布をカットするという方法もある（例えば、非特許文献２参照）。しかしこの方法では照度むらを検出しているのではないため、照度むら以外に元々の画像が持っていた情報までカットしている可能性がある。実際このような処理をモザイキング画像に対して行うことにより画像がぼけてしまったり、全体に白っぽくなったり黒っぽくなったりといった画像の劣化現象が発生していた。画像合成時の照度むら補正に関する従来技術としては、照明装置が固定されている場合に、照明からの位置により輝度が変化することを補正するための技術がある（例えば、特許文献１参照）。しかしこの技術は本発明が想定しているような照明装置が撮影装置ともに移動している状態には使えない。
【００１２】
また画像合成時の照度むら補正に関する別の従来技術としてスキャナなどで画像を取り込むときに、照明装置の揺らぎなどによって輝度が変化することを補正しようとするものがある（例えば、特許文献３参照。）。しかしそこでの補正手段は、基本的に画像全体の輝度平均の違いを補正しようとするものであり、本発明が考えているような場所ごとの細かい補正を考えたものではない。
【００１３】
また以上のような処理を行う場合に、画像列の間の画像移動量（オプティカルフローと通常呼ぶ）を画像パターンから計算するための方法として、マッチング法とグラディエント法が使われてきた（例えば、非特許文献３参照）。しかしマッチング法では画像データの量子化単位（ピクセル）までの計算精度しか出ないという問題があり、グラディエント法では大きな移動量は安定に計算できないという問題も付帯的問題としてあった。
【００１４】
【特許文献１】
特許第３２３３６０１号公報
【００１５】
【特許文献２】
特開平１１−１６４１３３号公報
【００１６】
【非特許文献１】
高木幹雄、下田陽久監修「画像解析ハンドブック」東京大学出版会　１９９１年１月、ｐ．４６３〜４６５　および　ｐ．４７８〜４７９
【００１７】
【非特許文献２】
Ｐｌａｎｅｔｒｏｎ，Ｉｎｃホームページ、Ｉｍａｇｅ−Ｐｒｏ　Ｐｌｕｓスタートアップマニュアル、［ｏｎｌｉｎｅ］、２００２年５月１５日、［平成１４年９月５日検索］、インターネット　＜ＵＲＬ：ｈｔｔｐ：／／ｗｗｗ．ｐｌａｎｅｔｒｏｎ．ｃｏ．ｊｐ／ｄｌ＿ｄｏｃｓ．ｈｔｍ＞第五章＜ＵＲＬ：ｈｔｔｐ：／／ｗｗｗ．ｐｌａｎｅｔｒｏｎ．ｃｏ．ｊｐ／ｐｄｆ／Ｖ４−Ｓｔ０５．ｐｄｆ＞第５−６頁
【００１８】
【非特許文献３】
三池秀敏　他　著、「パソコンによる動画像処理」森北出版株式会社　１９９３年７月、ｐ．１３５〜１７８　および　ｐ．２４１〜２４４
【００１９】
【発明が解決しようとする課題】
本発明は、上記のようなモザイキング時に発生する輝度の不連続な変化を人間の目には目立たない程度にまで精密に補正できて、しかも元の画像を劣化させることのない画像処理方法及びこの画像処理方法を用いた画像処理装置を提供することを目的としている。
【００２０】
【課題を解決するための手段】
本発明に係る画像処理方法は、例えば、画面の中央は明るく、画面の端は暗いというように、人工照明による照度分布は画像内の位置によってほぼ決まっていることに注目してなされたものであり、平均的な輝度変化を求めるのではなく、この画面内照度分布を画像面内座標値の関数として求めることを基本とするものである。
【００２１】
そのために画像面内座標値のべき乗関数の線形結合関数（多項式関数）を補正関数として構成し、その線形結合係数を最小二乗法により計算することによって、平均的な方法では除去できなかった場所ごとの細かい補正が比較的簡単な計算により可能となり、接合面における輝度の不連続性をかなり低減できる。また補正係数が座標値の関数として得られるので補正係数を画面全体にわたって計算することもできるので、重なった領域だけしか補正できないという欠点もない。さらに画像の照度分布を計算した上で照度むらの補正をしているので、照度むら以外の要因をカットしてしまうこともないため、画像の劣化も防ぐことができる。
【００２２】
また本発明に関わる画像処理方法は、撮影された画像を直接使うのではなく、その微分画像を使うものである。このように微分画像を使うと画像パターンのエッジが強調される。つまり、画像の中で特徴的なパターンに注目してオプティカルフローを算出することによってオプティカルフローの推定を安定に行うことができる。
【００２３】
また本発明に関わる画像処理方法では、従来から行われてきた二つのオプティカルフローの算出方法（マッチング法とグラディエント法）を組み合わせて、粗い算出と細かい算出の二段階で行うことにより、オプティカルフローの算出精度と安定性を高めている。
【００２４】
また本発明に関わる画像処理方法は、一つの貼り付け画像が持つ左右の境界線（図１１での１１３と１１５）での照度むら補正関数をそれぞれ計算し、画像の位置に応じてそれら二つの照度むら補正係数を連続的に変化・補間して補正するものであり、これにより更にスムーズな画像を得ることができる。
【００２５】
【発明の実施の形態】
以下では本発明の典型的な実施形態について説明する。
【００２６】
図１は、本発明に係る画像処理装置を含む合成画像作成装置１１の電気的構成を示すブロック図である。合成画像作成装置１１は、カメラで撮影した画像データの取り込み装置１２、画像処理装置１６、表示装置１７を含む。画像処理装置１６は、メモリ１３、外部記憶装置１４、中央演算処理回路１５を含む。画像データ取り込み装置１２でカメラから取り込んだ画像データは一度外部記憶装置１４に保存される。これらの画像データを入力画像データ列と以下呼ぶ。この入力画像データ列は、中央演算処理回路１５により順次取り出されながら処理を受け、再び外部記憶装置１４に保存される。ここで中央演算処理回路１５が処理を行うときにメモリ１３に一時的にデータを保存する処理をする。
【００２７】
図２は、本発明に係る画像処理装置１６の画像処理動作を説明するための機能的ブロック図である。機能的ブロック図では単一のブロックが中央演算処理回路１５の動作プログラムのサブルーチンを表す。中央演算処理回路１５が動作プログラムを実行した場合、まず中央演算処理回路１５は画像移動量算出手段２１として動作し、入力画像データ列各画像の相互の対応関係を調べて画像移動量を算出する。次いで中央演算処理回路９は照度むら補正手段として動作し、画像移動量算出手段２１が算出した画像移動量に基づいて照度むらの補正係数を算出する。次いで中央演算処理回路９は合成手段２３として動作し、画像移動量を参照して照度むら補正後の画像を合成する。この画像が出力画像として外部記憶装置１４に保存される。
【００２８】
図３は、本発明のモザイキング処理における撮影状態を説明する図である。図３における領域３１が一つの画像として撮影され、同様に領域３２が一つ画像として撮影され、さらに領域３３が一つ画像として撮影される。境界線４１から境界線４２が領域３１の画像の範囲となり、境界線４３から境界線４４が領域３２の画像の範囲となり、境界線４５から境界線４６が領域３３の画像の範囲となる。図４は、図３における各画像の輝度分布状態を説明するための図である。図４では、輝度分布５１は領域３１を撮影した画像の輝度分布であり、輝度分布５２は領域３２を撮影した画像の輝度分布であり、輝度分布５３は領域３３を撮影した画像の輝度分布である。
【００２９】
輝度分布グラフ５１に対応する輝度分布を画像面上の座標値の関数として表してＣ１（ｘ、ｙ）、とする。同様に輝度分布５２の輝度分布をＣ２（ｘ、ｙ）、輝度分布５３の輝度分布をＣ３（ｘ、ｙ）とする。ここでｘ、ｙは画像面での座標値である。以下ではモノクロ画像の場合を説明する。カラーの場合にはＲ（赤），Ｇ（緑），Ｂ（青）の色成分ごとに別々に照度むら補正係数を計算する方法や、Ｒ、Ｇ、Ｂ系をＸ、Ｙ、Ｚ系に変換した上で、輝度成分に対応するＹについて照度むら補正係数を求めてその照度むら補正係数をＲ、Ｇ，Ｂのすべてに適用すれば同様に適用可能である。
【００３０】
図５は、本発明の実施例における照度むら補正処理全体をまとめた図である。まず、ステップ６１では元画像の微分画像を計算する。次にステップ６２でオプティカルフローを粗く計算するためにマッチング法でのオプティカルフロー計算をする。次にステップ６３でオプティカルフローをさらに精密に計算するために、マッチング法の結果を出発点としてグラディエント法を用いることによりオプティカルフローを計算する。
【００３１】
次にステップ６４ではステップ６２とステップ６３の結果を組み合わせて最終的なオプティカルフローを計算する。次にステップ６５では照度むら補正係数を画像面の上の座標値ｘ、ｙの多項式（二次式）の形に展開して、多項式の係数を未知数として設定する。次にステップ６６ではステップ６５で未知数として設定した多項式係数の値を計算する。このとき重なり領域（図４の３４）では照度むら補正後の輝度が二つの輝度分布５１と輝度分布５２で一致するように決める。
【００３２】
以上のようにして二つの画像間での照度むら補正係数を計算する。この処理は図４における輝度分布５１と輝度分布５２の間で行われるとともに、輝度分布５２と輝度分布５３の間でも行われる。輝度分布５２と輝度分布５３の間での照度むら補正係数の計算が終了するとステップ６７に進み、ここでは、輝度分布５１と輝度分布５２の間での照度むら補正係数と、輝度分布５２と輝度分布５３の間での照度むら補正係数を、境界線４３からの距離に応じて重み付けを変えながら補間することによって、各ピクセルでの照度むら補正係数を計算する。
【００３３】
次にステップ６８ではステップ６７で計算した各ピクセルでの照度むら補正係数で補正しながら領域３２の画像を貼り付けていく。以下では図５の各ステップを詳細に説明していく。
【００３４】
ステップ６１：Ｃ１（ｘ、ｙ）とＣ２（ｘ、ｙ）の微分画像を計算して画像の輝度変化程度を計算する。微分画像を計算するためには画像データに微分フィルターを作用させる。微分フィルターとしては、Ｓｏｂｅｌ、Ｐｒｅｗｉｔｔ、Ｒｏｂｅｒｔｓ、Ｋｉｒｓｈ、Ｒｏｂｉｎｓｏｎなどの微分フィルターを使うことができる。本実施例では以下の微分フィルターを画像データに適用する。
【００３５】
【数１】

【００３６】
上のフィルターの場合には前後左右に隣り合ったピクセルの輝度データとの間で差を取ることに対応する。画像面座標（ｘ、ｙ）における微分画像輝度値をｄＣ１（ｘ、ｙ）とすると、上のフィルターは次の式で表すことができる。
【００３７】
【数２】

【００３８】
となる。ただし画面の端の部分（ｘ±１　あるいは　ｙ±１　がない場所）では若干変更する。同じ計算をＣ２についても行い、ｄＣ２（ｘ、ｙ）を計算する。
【００３９】
ステップ６２：次にマッチング法によるオプティカルフローの推定値を計算するために、画像全体を一定の大きさの区画（以下ブロックと呼ぶ）に分けて、その中からマッチング対象となるブロックを選び出す。例えば、１６ピクセルｘ１６ピクセルの区画を考えてこれをブロックとする。図６は、領域３１での微分画像ｄＣ１（ｘ、ｙ）をブロックに分割し、マッチング対象ブロック７０を選び出した状態を示す図である。マッチング対象ブロックを選び出す方法は次のようにする。
【００４０】
これらの各ブロックについて微分画像値の大きさを評価する。これはなるべく輝度変化の著しい特徴的なパターンを持つ領域をマッチングの対象とすることによりマッチング精度を高めるためである。このために各ブロックについて、ブロック内の各ピクセルにおける微分画像輝度値ｄＣ１（ｘ、ｙ）の２乗和を計算する。この計算結果をｄＣ１＊２と表すと次のような式で表すことができる（＾はべき乗を表す）。
【００４１】
【数３】

【００４２】
上式ではブロック内のピクセルについての和をΣ（ブロック内）と表している。ここで、大きな輝度変化を相対的に強調して小さな輝度変化を相対的に無視するために、微分画像値の４乗和を計算する方法を使うこともできる。
【００４３】
このようにして各ブロックについてｄＣ１＊２を計算し、それらの中で最も大きいｄＣ１＊２値を持つブロックを選び、このブロックをパターンマッチングの対象パターンとする。図６では斜線を引いたブロック７０がマッチング対象パターンとして選ばれている。
【００４４】
次に、マッチング対象パターン７０と領域３２の画像ｄＣ２（ｘ、ｙ）を重ねて相関係数を計算する。図７は、領域３２の微分画像ｄＣ２（ｘ、ｙ）にマッチング対象パターン７０を重ねて相関係数を計算するときの状態を示す図である。このときブロックの位置を元の位置から７３＝（ｄｘｗ、ｄｙｗ）だけ移動させた位置でｄＣ２（ｘ、ｙ）と重ねて相関を求める。このとき相関値は以下の式で計算する。
【００４５】
【数４】

【００４６】
ここで移動量（ｄｘｗ、ｄｙｗ）のｄｘｗ、ｄｙｗを全画面に渡って１ずつ変化させながら相関値を求める。そして相関値が最小となる移動量（ｄｘｗ、ｄｙｗ）を求める。これがマッチング法でのオプティカルフロー推定値となる。このようにして求めたオプティカルフロー推定値を以下ではｄｘ、ｄｙとする。
【００４７】
ステップ６３：次にステップ６２で得たオプティカルフロー推定値の精度をさらに高めるために、サブピクセル単位の微小なずれに対してどのように画像が変化するかを予測する方法を用いてオプティカルフローを計算するグラディエント法を適用する。
【００４８】
この方法では輝度の空間的な変化率から予測される輝度変化が、実際の輝度変化に対応するように空間変化量を計算し、それをオプティカルフローの推定値とする。即ち、或る点での輝度Ｉの空間微分が（∂Ｉ／∂ｘ、∂Ｉ／∂ｙ）であるとき、Δｘ、Δｙだけそこから離れた点における輝度変化量ΔＩは次式で与えられると推定する。
【００４９】
【数５】

【００５０】
このようにして得た輝度の変化量ΔＩが画像間での輝度変化と一致するようにΔｘ、Δｙを決めるというのがこの方法である。
本実施例の場合にはステップ６２で得たｄＣ１（ｘ＋ｄｘ、ｙ＋ｄｙ）とｄＣ２（ｘ、ｙ）との間での輝度変化を比較する。すると以下の式からΔｘとΔｙが推定される。
【００５１】
【数６】

【００５２】
ここで画像全体が一定量だけずれているとするならば、すべてのピクセルに対してΔｘとΔｙ　は同じ量になる。したがって二つの未知数ΔｘとΔｙ　に対してピクセルの数だけの方程式が存在することになる。これらの過剰な方程式を取り込むために最小二乗法により未知数ΔｘとΔｙ　を決める。即ち、以下の費用関数を定義して、
【００５３】
【数７】

【００５４】
これが最小になるようにΔｘとΔｙを決める。この方法によれば多数のデータを統計的に処理していることになるため、本来の単位であるピクセル以下の精度で運動量を決定することができる。ステップ６２とステップ６３を総合することによって画像Ｃ１（ｘ、ｙ）と画像Ｃ２（ｘ、ｙ）の間のオプティカルフローが（ｄｘ＋Δｘ、ｄｙ＋Δｙ）と計算される。このとき以下の式が成り立つ。
【００５５】
【数８】

【００５６】
（ｄｘ＋Δｘ、ｄｙ＋Δｙ）を（−Ｄｘ、−Ｄｙ）とすると、上の式は以下のようになる。
【００５７】
【数９】

【００５８】
これをＣ１（ｘ、ｙ）を中心に書き直して以下のように書くこともできる。
【００５９】
【数１０】

【００６０】
即ち、Ｃ１（ｘ、ｙ）において（ｘ、ｙ）であった点は、Ｃ２（ｘ、ｙ）では（ｘ＋Ｄｘ，ｙ＋Ｄｙ）に移動している。以下ではこの書き方を使う。
ステップ６４：次に照度むら補正のステップに進む。このために照度むら補正関数の形を仮定し、いくつかのパラメータを用いた関数形によって表現する。照度
むらの効果は
【課題を解決するための手段】でも述べたように画像上の座標で決まっていると考えられる。そこで画像面内の照度分布Ｌ（ｘ，ｙ）がスムーズな関数であると仮定する。このとき照度分布を均一にするための補正関数Ｈ（ｘ，ｙ）もスムーズな関数になる。ここでＨ（ｘ，ｙ）とＬ（ｘ，ｙ）の関係は次の式で与えられる。
【００６１】
【数１１】

【００６２】
ここで照度分布の関数を物理的に決める。図８は、本実施例における照明の位置と撮影対象の位置との関係を示す図である。８０は照明の位置を示し、８１が撮影対象点の位置である。８１の座標を（ｘ、ｙ）とする。照明と撮影平面との距離８２をｈ、照明の真正面に当たる場所８３と撮影対象点との距離８４をｄとする。このとき撮影対象点８１と照明８０との距離ｒは次式で与えられる。
【００６３】
【数１２】

【００６４】
一般に照度は光源からの距離の二乗に反比例するとされている。したがって、撮影対象点８１での照度Ｌ（ｘ，ｙ）は、Ｌ０を或る定数として次式で与えられる。
【００６５】
【数１３】

【００６６】
これから補正関数Ｈ（ｘ，ｙ）は次式のようになる。
【００６７】
【数１４】

【００６８】
これを一般化して、補正関数はｘ、ｙの二次関数を第一近似として採用することができる。一般にｘ、ｙ座標に関するより高次の多項式関数を補正関数に用いることにより、補正の精度を順次上げて行けば、必要な精度に合わせて補正をすることができて、実用上有効である。以下の説明では二次関数を仮定する。このとき、補正関数を以下のように６個のパラメータ（ａ、ｂ、ｃ、ｄ、ｅ、ｆ）を用いて表現することができる。
【００６９】
【数１５】

【００７０】
ステップ６５：次に画像Ｃ１（ｘ、ｙ）とＣ２（ｘ、ｙ）で対応する点（物理的に同じ点）は同じ輝度になるように照度むら補正係数を決める。Ｃ１（ｘ、ｙ）とＣ２（ｘ、ｙ）において対応する点はオプティカルフローを用いて知ることができる。即ち、Ｃ１（ｘ、ｙ）において（ｘ、ｙ）であった点は、Ｃ２（ｘ、ｙ）では（ｘ＋Ｄｘ，ｙ＋Ｄｙ）に移動している。そこで、照度むら補正係数をＨ１（ｘ，ｙ）と書いたとき、Ｃ１（ｘ、ｙ）に照度むら補正を施すと
【００７１】
【数１６】

【００７２】
になる。一方Ｃ２（ｘ＋Ｄｘ，ｙ＋Ｄｙ）に照度むら補正を施すと、
【００７３】
【数１７】

【００７４】
になる。これら二つの点は本来同じ点であり同じ輝度であるので次式が成り立つ。
【００７５】
【数１８】

【００７６】
ここでＨ１（ｘ、ｙ）６個のパラメータ（ａ、ｂ、ｃ、ｄ、ｅ、ｆ）で表されているから未知数が６個である。これに対して上の方程式はピクセルの数だけある。そこでステップ６３と同じように，最小二乗法によって次の費用関数を最小化するように６個のパラメータ（ａ、ｂ、ｃ、ｄ、ｅ、ｆ）を求める。
【００７７】
【数１９】

【００７８】
ここで画像Ｃ１（ｘ、ｙ）と画像Ｃ２（ｘ、ｙ）が重なっている領域（図３の３４）内のピクセルに関する和をΣ（重なり領域）と表現した。
しかしこれだけでは（ａ、ｂ、ｃ、ｄ、ｅ、ｆ）が求められない。なぜならば、もし式（ａ）を満足するＨ（ｘ、ｙ）が求まったとすると、それを定数倍したｋ＊Ｈ１（ｘ、ｙ）もやはり式（ａ）を満足するから、結局式（ａ）だけでは一意に（ａ、ｂ、ｃ、ｄ、ｅ、ｆ）を決めることができないからである。このためにもう一つ条件式を追加する。それは上で述べた不定係数ｋを決めるためのもので、照度むら補正後の輝度の画面平均が照度むら補正前と変わらないという条件である。これを数式で表現すると次式のようになる。
【００７９】
【数２０】

【００８０】
以上をまとめると、式（Ｂ）の制限条件の下で式（Ａ）を満足するように（ａ、ｂ、ｃ、ｄ、ｅ、ｆ）を定めればよい。制限条件付の最大値・最小値問題は、ラグランジュ（Ｌａｇｒａｎｇｅ）の未定係数法により解くことができることが解析学で知られている。
ステップ６６：次にステップ６５で決まった照度むら補正関数Ｈ１（ｘ、ｙ）を用いてＣ１（ｘ、ｙ）とＣ２（ｘ、ｙ）の輝度補正を行う。輝度補正後のＣ１（ｘ、ｙ）、Ｃ２（ｘ、ｙ）をｃＣ２（ｘ、ｙ）、ｃＣ２（ｘ、ｙ）と書くと次式のようになる。
【００８１】
【数２１】

【００８２】
図９は、照度むら補正を行った状態を示す図である。領域３１での画像と領域３２での画像の間の照度むら補正係数９０を領域３１での画像の輝度分布５１に適用すると、照度むらが補正されて輝度分布９４が得られる。領域３１での画像と領域３２での画像の間の照度むら補正係数９１を領域３２での画像の輝度分布５２に適用すると、照度むらが補正されて輝度分布９５が得られる。両者は境界４３において一致していて、不連続的な輝度の変化がない。
【００８３】
以上のプロセスにより図３で領域３１での画像と領域３２での画像の間のモザイキングと照度むら補正ができる。そして同じことを図３の領域３２と領域３３との間でも行うことにより、領域３２での画像と領域３３での画像の間の照度むら補正係数Ｈ２（ｘ、ｙ）を計算することができる。この状態を図９に示した。領域３２での画像と領域３３での画像の間の照度むら補正係数９２を領域３２での画像の輝度分布５２に適用すると、照度むらが補正されて輝度分布９６が得られる。領域３１での画像と領域３２での画像の間の照度むら補正係数９３を領域３３での画像の輝度分布５３に適用すると、照度むらが補正されて輝度分布９７が得られる。両者は境界４５において一致していて、不連続的な輝度の変化がない。
【００８４】
ところが、一般には領域３１と領域３２の画像間で計算した照度むら補正係数と領域３２と領域３３の画像間で計算した照度むら補正係数は異なる。このために、領域３１と領域３２の画像間で計算した照度むら補正係数を領域３２での画像に使うと、境界線４５での不連続的な輝度変化９９が発生する。一方、領域３２と領域３３の画像間で計算した照度むら補正係数を領域３２での画像に使うと、境界線４３での不連続的な輝度変化９８が発生する。
【００８５】
この問題を解決するために、照度むら補正係数の値を領域３２での画像の中で連続的に変化させる方法を使う。図１０は、領域３２の画像面内での照度むら補正係数を補間した状態を示す図である。領域３２での画像において、境界４３では領域３１と領域３２の画像間で計算した照度むら補正係数を使い、境界４５では領域３２と領域３３の画像間で計算した照度むら補正係数を使う。その間は１００のように連続的に照度むら補正係数が変わるようにする。
【００８６】
この方法を、数式を用いて表現すると次のようになる。領域３１と領域３２の画像間で計算した照度むら補正係数をＨ１（ｘ、ｙ）とし、領域３２と領域３３の画像間で計算した照度むら補正係数をＨ２（ｘ、ｙ）とする。境界４３の、領域３２での画像上のｘ座標をｘ１，境界４５の領域３２での画像上のｘ座標をｘ２とするとき、座標値がｘの場所での補正係数として以下の係数を使う。
【００８７】
【数２２】

【００８８】
この方法により境界４３においても境界４５においても不連続的な輝度変化をなくすことができる。この手順の場合には、領域３２と領域３３の画像間で計算した照度むら補正係数が求まった段階で、領域３２での画像を照度むら補正しながらパノラマ画像に貼り付けることになる。
【００８９】
【発明の効果】
本発明によれば、画像面内の照度分布により発生するモザイキング画像内の不連続的な輝度変化を除去し、自然なパノラマ画像を得ることができる。しかも照度むら効果のみを抽出して本来の画像輝度を復元しているため、周波数フィルタリングのように本来の画像を壊してしまうことがない。また、画像面内の照度分布を考えているので場所毎の細かい補正が可能となり、不連続的な輝度変化の除去能力が高い。したがって極めて有効な方法である。
【図面の簡単な説明】
【図１】本発明の画像処理装置を含む合成画像作成装置の構成を示すブロック図である。
【図２】本発明に係る画像処理装置の画像処理動作を説明するための機能的ブロック図である。
【図３】本発明のモザイキング処理における撮影状態を説明する図である。
【図４】図３における各画像の輝度分布状態を説明するための図である。
【図５】本発明に係る画像処理方法を説明するためのフローチャートである。
【図６】領域３１での微分画像ｄＣ１（ｘ、ｙ）をブロックに分割し、マッチング対象パターン７０を選び出した状態を示す図である。
【図７】領域３２の微分画像ｄＣ２（ｘ、ｙ）にマッチング対象パターン７０を重ねて相関係数を計算するときの状態を示す図である。
【図８】本実施例における照明の位置と撮影対象の位置との関係を示す図である。
【図９】照度むら補正を行った状態を示す図である。
【図１０】領域３２の画像面内での照度むら補正係数を補間した状態を示す図である。
【図１１】従来のモザイキング処理における撮影状態を説明する図である。
【図１２】照度むらがあるときの図１１における各画像の輝度分布状態を説明するための図である。
【図１３】従来の方法により図１２の画像をモザイキングしたときの結果を示す図である。
【符号の説明】
１１　合成画像作成装置
１２　画像データ取り込み装置
１３　メモリ
１４　外部記憶装置
１５　中央演算処理回路
１６　画像処理装置
１７　表示装置
５１，５２，５３　画像列における各画像の輝度分布
９０，９１，９２，９３　画像列における各画像の照度むら補正係数の計算値
９４，９５，９６，９７　照度むら補正係数で画像列における画像の照度むらを補正した結果の画像の輝度分布[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to an image processing method for combining continuous images or a plurality of still images captured by a video camera to synthesize an image as one panoramic image, and an image processing apparatus using the image processing method.
[0002]
[Prior art]
With the spread of digital video and the high performance of personal computers, it has become possible to store and process a large amount of digital image data, which was once impossible. Accordingly, new methods for processing digital image data have been introduced. As one of the techniques, there is a technique of generating one composite image (panoramic image) by connecting a video image captured while moving or a plurality of still images captured in the same area at different locations and angles. This is sometimes called mosaiking.
[0003]
This technology is used effectively when inspecting tunnels and walls for cracks. When performing the mosaicing process in the image processing apparatus, first, an overlapping area where the images overlap with each other is obtained, and then the image is pasted using the overlapping area as a “paste”.
[0004]
When such image synthesis is performed, the luminance of the pixels of the synthesized image rapidly changes at a portion corresponding to the seam of the image due to a change in the lighting state of each image, and therefore, the boundary between the shading and the tone in the synthesized image. Can be done.
[0005]
FIG. 11 is a view for explaining a photographing state in the conventional mosaicing process. It is assumed that the region 101 in FIG. 11 is photographed as one image, the region 102 is photographed as one image, and the region 103 is photographed as one image. The boundary line 111 to the boundary line 112 are the range of the image of the region 101, the boundary line 113 to the boundary line 114 are the range of the image of the region 102, and the boundary line 115 to the boundary line 116 are the range of the image of the region 103. FIG. 12 is a diagram for explaining the luminance distribution state of each image in FIG. 11 when there is uneven illuminance.
[0006]
In FIG. 12, the luminance distribution 121 of the image obtained by photographing the area 101 is shown. The reason why such a distribution is obtained is that when artificial illumination is performed, the center is bright and the edges are dark. For the same reason, similar luminance distributions occur in the luminance distributions 122 and 123 in the images obtained by photographing the

regions

102 and 103. At this time, the images of the two

regions

101 and 102 are pasted together with the overlap region 104 as the “paste”. Here, portions 111 to 113 in FIG. 12 are cut out of the image of the area 101, and the cut-out image is pasted to the left of the image of the area 102. Next, regions 113 to 115 in FIG. 12 are cut out from the image of the region 102 and pasted to the left of the image of the region 103.
[0007]
FIG. 13 is a diagram showing a result when the image of FIG. 12 is mosaiked by such a conventional method. When the screen is pasted as described here, discontinuous “jumps” 131 and 132 in luminance occur at the boundary line 113 and the boundary line 115 in FIG. 13 and an unnatural line appears in the panoramic image. Looks like.
[0008]
As one method for solving this problem, there is a method of detecting an overlapping area between different images by some method, and correcting the illuminance unevenness effect using a difference in luminance there. That is, since the overlapping areas are originally the same place, the luminance should be the same, but the luminance differs in the image due to the uneven illuminance effect. Therefore, the basic idea is to compare the luminances in the areas and correct them so that they are the same.
[0009]
Methods based on such a concept include a histogram method, a linear density conversion method, and an average density difference correction method. The histogram method and the linear density conversion method can correct the luminance of the entire image, but have the disadvantage that fine luminance correction cannot be performed according to the location.
[0010]
Because the human eye is very sensitive to discontinuous changes, it is difficult to eliminate the discontinuity in brightness to a lesser degree with these methods. On the other hand, although the average density difference correction method enables fine correction according to the location, it has a drawback that it is not possible to correct the entire image including non-overlapping regions (for example, see Non-Patent Document 1).
[0011]
As another type of idea, there is a method of applying a smoothing filter to luminance data of an image to cut the distribution of illuminance in the image (for example, see Non-Patent Document 2). However, since this method does not detect uneven illuminance, there is a possibility that information other than uneven illuminance may have been cut off from the original image. Actually, when such processing is performed on the mosaicing image, the image is degraded, and the entire image becomes whitish or dark. As a conventional technique for correcting uneven illuminance at the time of image synthesis, there is a technique for correcting a change in luminance depending on a position from illumination when a lighting device is fixed (for example, see Patent Document 1). However, this technique cannot be used in a state where the illumination device is moving together with the image capturing device as assumed by the present invention.
[0012]
Further, as another conventional technique relating to correction of uneven illuminance at the time of image synthesis, there is a technique which attempts to correct a change in luminance due to fluctuation of an illumination device when an image is captured by a scanner or the like (for example, see Patent Document 3). ). However, the correction means there is basically intended to correct the difference in the average luminance of the entire image, and is not a fine correction for each location as considered by the present invention.
[0013]
In performing the above-described processing, a matching method and a gradient method have been used as a method for calculating an image movement amount (commonly referred to as an optical flow) between image sequences from an image pattern (for example, Non-Patent Document 3). However, the matching method has a problem that only the calculation accuracy up to the quantization unit (pixel) of the image data is obtained, and the gradient method cannot calculate a large moving amount stably.
[0014]
[Patent Document 1]
Japanese Patent No. 3233601
[0015]
[Patent Document 2]
JP-A-11-164133
[0016]
[Non-patent document 1]
Mikio Takagi and Hirohisa Shimoda, “Image Analysis Handbook,” University of Tokyo Press, January 1991, p. 463-465 and p. 478-479
[0017]
[Non-patent document 2]
Planetron, Inc homepage, Image-Pro Plus Startup Manual, [online], May 15, 2002, [searched September 5, 2002], Internet <URL: http: // www. plantron. co. jp / dl_docs. htm > Chapter 5 <URL: http: // www. plantron. co. jp / pdf / V4-St05. pdf > Pages 5-6
[0018]
[Non-Patent Document 3]
Miike Hidetoshi et al., "Moving image processing by personal computer", Morikita Publishing Co., Ltd., July 1993, p. 135-178 and p. 241-244
[0019]
[Problems to be solved by the invention]
The present invention provides an image processing method that can accurately correct discontinuous changes in luminance that occur during mosaicking as described above to such an extent that they are inconspicuous to human eyes, and that does not degrade the original image. It is an object of the present invention to provide an image processing device using an image processing method.
[0020]
[Means for Solving the Problems]
The image processing method according to the present invention focuses on the fact that the illuminance distribution due to artificial illumination is substantially determined by the position in the image, for example, the center of the screen is bright and the edges of the screen are dark. In other words, instead of calculating an average luminance change, the illuminance distribution in the screen is obtained as a function of the coordinate value in the image plane.
[0021]
For this purpose, a linear combination function (polynomial function) of the power function of the coordinate values in the image plane is configured as a correction function, and the linear combination coefficient is calculated by the least-squares method. Correction can be made by relatively simple calculation, and the discontinuity of luminance at the joint surface can be considerably reduced. Further, since the correction coefficient is obtained as a function of the coordinate value, the correction coefficient can be calculated over the entire screen, so that there is no drawback that only the overlapping area can be corrected. Further, since the illuminance unevenness is corrected after calculating the illuminance distribution of the image, factors other than the illuminance unevenness are not cut off, so that the deterioration of the image can be prevented.
[0022]
The image processing method according to the present invention does not directly use a captured image but uses a differential image thereof. The use of the differential image emphasizes the edges of the image pattern. That is, the optical flow can be stably estimated by calculating the optical flow while paying attention to the characteristic pattern in the image.
[0023]
Also, in the image processing method according to the present invention, the two optical flow calculation methods (matching method and gradient method) that have been conventionally performed are combined and performed in two stages of coarse calculation and fine calculation, so that the optical flow Calculation accuracy and stability are improved.
[0024]
Further, the image processing method according to the present invention calculates the illuminance unevenness correction functions at the left and right border lines (113 and 115 in FIG. 11) of one pasted image, and calculates the two illuminance unevenness functions according to the position of the image. The correction is performed by continuously changing and interpolating the illuminance unevenness correction coefficient, whereby a smoother image can be obtained.
[0025]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, a typical embodiment of the present invention will be described.
[0026]
FIG. 1 is a block diagram showing an electrical configuration of a composite image creating device 11 including an image processing device according to the present invention. The composite image creation device 11 includes a device 12 for capturing image data captured by a camera, an image processing device 16, and a display device 17. The image processing device 16 includes a memory 13, an external storage device 14, and a central processing circuit 15. The image data captured from the camera by the image data capturing device 12 is once stored in the external storage device 14. These image data are hereinafter referred to as an input image data sequence. The input image data sequence is processed while being sequentially taken out by the central processing circuit 15 and stored again in the external storage device 14. Here, when the central processing circuit 15 performs the processing, the processing for temporarily storing data in the memory 13 is performed.
[0027]
FIG. 2 is a functional block diagram for explaining an image processing operation of the image processing device 16 according to the present invention. In the functional block diagram, a single block represents a subroutine of the operation program of the central processing unit 15. When the central processing circuit 15 executes the operation program, first, the central processing circuit 15 operates as the image moving amount calculating means 21 to calculate the image moving amount by examining the mutual correspondence of each image of the input image data sequence. . Next, the central processing circuit 9 operates as an illuminance unevenness correcting unit, and calculates a correction coefficient for the illuminance unevenness based on the image moving amount calculated by the image moving amount calculating unit 21. Next, the central processing circuit 9 operates as the synthesizing means 23, and synthesizes the image after the illuminance unevenness correction with reference to the image movement amount. This image is stored in the external storage device 14 as an output image.
[0028]
FIG. 3 is a diagram for explaining a photographing state in the mosaicking process of the present invention. The region 31 in FIG. 3 is photographed as one image, the region 32 is photographed as one image, and the region 33 is photographed as one image. The boundary line 41 to the boundary line 42 become the image range of the region 31, the boundary line 43 to the boundary line 44 become the image range of the region 32, and the boundary line 45 to the boundary line 46 become the image range of the region 33. FIG. 4 is a diagram for explaining a luminance distribution state of each image in FIG. In FIG. 4, a luminance distribution 51 is a luminance distribution of an image of the area 31, a luminance distribution 52 is a luminance distribution of an image of the area 32, and a luminance distribution 53 is a luminance distribution of an image of the area 33. is there.
[0029]
The brightness distribution corresponding to the brightness distribution graph 51 is expressed as a function of the coordinate value on the image plane, and is represented by C1 (x, y). Similarly, let the luminance distribution of the luminance distribution 52 be C2 (x, y) and the luminance distribution of the luminance distribution 53 be C3 (x, y). Here, x and y are coordinate values on the image plane. Hereinafter, a case of a monochrome image will be described. In the case of color, a method of separately calculating the illuminance unevenness correction coefficient for each of the R (red), G (green), and B (blue) color components, or changing the R, G, B system to the X, Y, Z system After conversion, an uneven illuminance correction coefficient is obtained for Y corresponding to the luminance component, and the uneven illuminance correction coefficient is applied to all of R, G, and B, and the same applies.
[0030]
FIG. 5 is a diagram summarizing the entire illuminance unevenness correction process according to the embodiment of the present invention. First, in step 61, a differential image of the original image is calculated. Next, in step 62, an optical flow is calculated by a matching method in order to roughly calculate the optical flow. Next, in step 63, in order to calculate the optical flow more precisely, the optical flow is calculated by using the gradient method with the result of the matching method as a starting point.
[0031]
Next, in step 64, the final optical flow is calculated by combining the results of

steps

62 and 63. Next, in step 65, the uneven illuminance correction coefficient is developed into a polynomial (quadratic) of coordinate values x and y on the image plane, and the coefficient of the polynomial is set as an unknown. Next, in step 66, the value of the polynomial coefficient set as an unknown in step 65 is calculated. At this time, in the overlapping area (34 in FIG. 4), the luminance after the illuminance non-uniformity correction is determined so as to match between the two luminance distributions 51 and 52.
[0032]
The illuminance unevenness correction coefficient between the two images is calculated as described above. This processing is performed between the luminance distribution 51 and the luminance distribution 52 in FIG. When the calculation of the illuminance unevenness correction coefficient between the luminance distribution 52 and the luminance distribution 53 is completed, the process proceeds to step 67, where the illuminance unevenness correction coefficient between the luminance distribution 51 and the luminance distribution 52, the luminance distribution 52 and the luminance The illuminance unevenness correction coefficient at each pixel is calculated by interpolating the uneven illuminance correction coefficient between the distributions 53 while changing the weight according to the distance from the boundary line 43.
[0033]
Next, in step 68, the image of the area 32 is pasted while correcting with the illuminance unevenness correction coefficient for each pixel calculated in step 67. Hereinafter, each step of FIG. 5 will be described in detail.
[0034]
Step 61: A differential image of C1 (x, y) and C2 (x, y) is calculated to calculate a degree of luminance change of the image. To calculate a differential image, a differential filter is applied to the image data. As a differential filter, a differential filter such as Sobel, Prewitt, Roberts, Kirsh, and Robinson can be used. In the present embodiment, the following differential filter is applied to image data.
[0035]
(Equation 1)

[0036]
In the case of the above filter, it corresponds to taking the difference between the luminance data of the pixels adjacent to the front, rear, left and right. Assuming that the differential image luminance value at the image plane coordinates (x, y) is dC1 (x, y), the above filter can be expressed by the following equation.
[0037]
(Equation 2)

[0038]
It becomes. However, slight changes are made at the edge of the screen (where there is no x ± 1 or y ± 1). The same calculation is performed for C2, and dC2 (x, y) is calculated.
[0039]
Step 62: Next, in order to calculate an estimated value of the optical flow by the matching method, the whole image is divided into blocks (hereinafter, referred to as blocks) of a fixed size, and a block to be matched is selected from the divided blocks. For example, a block of 16 pixels × 16 pixels is considered and is set as a block. FIG. 6 is a diagram showing a state where the differential image dC1 (x, y) in the area 31 is divided into blocks, and a matching target block 70 is selected. The method of selecting a matching target block is as follows.
[0040]
The magnitude of the differential image value is evaluated for each of these blocks. This is to improve the matching accuracy by targeting a region having a characteristic pattern with a remarkable luminance change as much as possible. For this, for each block, the sum of squares of the differential image luminance value dC1 (x, y) at each pixel in the block is calculated. If this calculation result is expressed as dC1 * 2, it can be expressed by the following equation (＾ represents a power).
[0041]
[Equation 3]

[0042]
In the above equation, the sum of the pixels in the block is represented as Σ (in the block). Here, in order to relatively emphasize large luminance changes and relatively ignore small luminance changes, a method of calculating the sum of squares of differential image values may be used.
[0043]
In this way, dC1 * 2 is calculated for each block, and a block having the largest dC1 * 2 value is selected from these, and this block is set as a target pattern for pattern matching. In FIG. 6, a hatched block 70 is selected as a matching target pattern.
[0044]
Next, a correlation coefficient is calculated by superimposing the matching target pattern 70 and the image dC2 (x, y) of the region 32. FIG. 7 is a diagram illustrating a state in which the matching coefficient is calculated by superimposing the matching target pattern 70 on the differential image dC2 (x, y) of the region 32. At this time, a correlation is obtained by superimposing dC2 (x, y) at a position where the position of the block has been shifted by 73 = (dxw, dyw) from the original position. At this time, the correlation value is calculated by the following equation.
[0045]
(Equation 4)

[0046]
Here, the correlation value is obtained while changing dxw and dyw of the movement amount (dxw, dyw) one by one over the entire screen. Then, the movement amount (dxw, dyw) at which the correlation value becomes minimum is obtained. This is an optical flow estimation value by the matching method. The optical flow estimation values obtained in this manner are hereinafter referred to as dx and dy.
[0047]
Step 63: Next, in order to further improve the accuracy of the optical flow estimation value obtained in step 62, the optical flow is estimated using a method of predicting how an image changes with respect to a slight shift in subpixel units. Apply the gradient method to calculate.
[0048]
In this method, a spatial change amount is calculated so that a luminance change predicted from a spatial change rate of luminance corresponds to an actual luminance change, and the calculated amount is used as an estimated value of an optical flow. That is, when the spatial derivative of the luminance I at a certain point is (∂I / ∂x, ∂I / ∂y), the luminance change amount ΔI at a point away from it by Δx and Δy is given by the following equation. It is estimated.
[0049]
(Equation 5)

[0050]
In this method, Δx and Δy are determined so that the amount of change in luminance ΔI obtained in this way matches the change in luminance between images.
In the case of this embodiment, the luminance change between dC1 (x + dx, y + dy) obtained in step 62 and dC2 (x, y) is compared. Then, Δx and Δy are estimated from the following equations.
[0051]
(Equation 6)

[0052]
Here, if the entire image is shifted by a certain amount, Δx and Δy are the same for all pixels. Therefore, for the two unknowns Δx and Δy, there are as many equations as the number of pixels. In order to take in these excess equations, the unknowns Δx and Δy are determined by the least squares method. That is, by defining the following cost function,
[0053]
(Equation 7)

[0054]
Δx and Δy are determined so as to minimize this. According to this method, a large amount of data is statistically processed, so that the momentum can be determined with an accuracy of a pixel or less, which is the original unit. By combining step 62 and step 63, the optical flow between the image C1 (x, y) and the image C2 (x, y) is calculated as (dx + Δx, dy + Δy). At this time, the following equation is established.
[0055]
(Equation 8)

[0056]
If (dx + Δx, dy + Δy) is (−Dx, −Dy), the above equation is as follows.
[0057]
(Equation 9)

[0058]
This can be rewritten with C1 (x, y) as the center and written as follows.
[0059]
(Equation 10)

[0060]
That is, the point that was (x, y) in C1 (x, y) has moved to (x + Dx, y + Dy) in C2 (x, y). We use this notation below.
Step 64: Next, proceed to the step of correcting uneven illuminance. For this purpose, the form of the illuminance unevenness correction function is assumed and expressed by a function form using some parameters. Illumination
The effect of unevenness
As described above, it is considered that the coordinates are determined by the coordinates on the image. Therefore, it is assumed that the illuminance distribution L (x, y) in the image plane is a smooth function. At this time, the correction function H (x, y) for making the illuminance distribution uniform also becomes a smooth function. Here, the relationship between H (x, y) and L (x, y) is given by the following equation.
[0061]
[Equation 11]

[0062]
Here, the function of the illuminance distribution is physically determined. FIG. 8 is a diagram illustrating a relationship between the position of the illumination and the position of the imaging target in the present embodiment. Numeral 80 indicates the position of the illumination, and 81 indicates the position of the photographing target point. The coordinates of 81 are (x, y). The distance 82 between the illumination and the imaging plane is h, and the distance 84 between the location 83 directly in front of the illumination and the imaging target point is d. At this time, the distance r between the imaging target point 81 and the illumination 80 is given by the following equation.
[0063]
(Equation 12)

[0064]
Generally, the illuminance is inversely proportional to the square of the distance from the light source. Therefore, the illuminance L (x, y) at the imaging target point 81 is given by the following equation, using L0 as a certain constant.
[0065]
(Equation 13)

[0066]
From this, the correction function H (x, y) is as follows.
[0067]
[Equation 14]

[0068]
Generalizing this, a quadratic function of x and y can be adopted as a first approximation as the correction function. In general, if a higher-order polynomial function relating to the x and y coordinates is used as a correction function, and the correction accuracy is sequentially increased, the correction can be performed to the required accuracy, which is practically effective. In the following description, a quadratic function is assumed. At this time, the correction function can be expressed using the following six parameters (a, b, c, d, e, f).
[0069]
[Equation 15]

[0070]
Step 65: Next, an illuminance unevenness correction coefficient is determined so that the corresponding points (physically the same points) in the images C1 (x, y) and C2 (x, y) have the same luminance. The corresponding points in C1 (x, y) and C2 (x, y) can be known using an optical flow. That is, the point that was (x, y) in C1 (x, y) has moved to (x + Dx, y + Dy) in C2 (x, y). Therefore, when the illuminance unevenness correction coefficient is written as H1 (x, y), the illuminance unevenness correction is performed on C1 (x, y).
[0071]
(Equation 16)

[0072]
become. On the other hand, when illuminance unevenness correction is performed on C2 (x + Dx, y + Dy),
[0073]
[Equation 17]

[0074]
become. Since these two points are originally the same point and have the same luminance, the following equation holds.
[0075]
(Equation 18)

[0076]
Here, H1 (x, y) is represented by six parameters (a, b, c, d, e, f), so that there are six unknowns. In contrast, the above equation has only the number of pixels. Thus, as in step 63, six parameters (a, b, c, d, e, f) are obtained by the least squares method so as to minimize the following cost function.
[0077]
[Equation 19]

[0078]
Here, the sum of the pixels in the area where the image C1 (x, y) and the image C2 (x, y) overlap (34 in FIG. 3) is expressed as Σ (overlap area).
However, (a, b, c, d, e, f) cannot be obtained by this alone. This is because, if H (x, y) satisfying the equation (a) is obtained, k * H1 (x, y) obtained by multiplying H (x, y) by a constant also satisfies the equation (a). ) Alone cannot uniquely determine (a, b, c, d, e, f). For this purpose, another conditional expression is added. This is for determining the indefinite coefficient k described above, and is a condition that the screen average of the luminance after the correction of the uneven illuminance is not different from that before the correction of the uneven illuminance. This can be expressed by the following equation.
[0079]
(Equation 20)

[0080]
To summarize the above, (a, b, c, d, e, f) may be determined so as to satisfy Expression (A) under the restriction condition of Expression (B). It is known from analytics that the maximum / minimum problem with a constraint can be solved by the Lagrange's undetermined coefficient method.
Step 66: Next, the luminance correction of C1 (x, y) and C2 (x, y) is performed using the illuminance unevenness correction function H1 (x, y) determined in step 65. When C1 (x, y) and C2 (x, y) after luminance correction are written as cC2 (x, y) and cC2 (x, y), the following equation is obtained.
[0081]
(Equation 21)

[0082]
FIG. 9 is a diagram showing a state in which uneven illuminance has been corrected. When the uneven illuminance correction coefficient 90 between the image in the area 31 and the image in the area 32 is applied to the luminance distribution 51 of the image in the area 31, the uneven illuminance is corrected and a luminance distribution 94 is obtained. When the uneven illuminance correction coefficient 91 between the image in the area 31 and the image in the area 32 is applied to the luminance distribution 52 of the image in the area 32, the uneven illuminance is corrected, and the luminance distribution 95 is obtained. Both coincide at the boundary 43 and there is no discontinuous change in luminance.
[0083]
By the above process, the mosaicing between the image in the region 31 and the image in the region 32 in FIG. By performing the same operation between the region 32 and the region 33 in FIG. 3, it is possible to calculate the illuminance unevenness correction coefficient H2 (x, y) between the image in the region 32 and the image in the region 33. . This state is shown in FIG. When the uneven illuminance correction coefficient 92 between the image in the area 32 and the image in the area 33 is applied to the luminance distribution 52 of the image in the area 32, the uneven illuminance is corrected, and the luminance distribution 96 is obtained. When the uneven illuminance correction coefficient 93 between the image in the area 31 and the image in the area 32 is applied to the luminance distribution 53 of the image in the area 33, the uneven illuminance is corrected, and the luminance distribution 97 is obtained. Both coincide at the boundary 45, and there is no discontinuous change in luminance.
[0084]
However, generally, the illuminance unevenness correction coefficient calculated between the images of the area 31 and the area 32 is different from the illuminance unevenness correction coefficient calculated between the images of the area 32 and the area 33. For this reason, if the illuminance unevenness correction coefficient calculated between the images in the region 31 and the region 32 is used for the image in the region 32, a discontinuous luminance change 99 at the boundary 45 occurs. On the other hand, when the illuminance unevenness correction coefficient calculated between the images of the region 32 and the region 33 is used for the image of the region 32, a discontinuous luminance change 98 at the boundary 43 occurs.
[0085]
In order to solve this problem, a method of continuously changing the value of the illuminance unevenness correction coefficient in the image in the area 32 is used. FIG. 10 is a diagram illustrating a state in which the illuminance unevenness correction coefficient in the image plane of the region 32 is interpolated. In the image of the region 32, the boundary 43 uses the uneven illuminance correction coefficient calculated between the images of the region 31 and the region 32, and the boundary 45 uses the uneven illuminance correction coefficient calculated between the images of the region 32 and the region 33. In the meantime, the illuminance unevenness correction coefficient is continuously changed like 100.
[0086]
This method is expressed as follows using mathematical expressions. The illuminance unevenness correction coefficient calculated between the images of the area 31 and the area 32 is H1 (x, y), and the illuminance unevenness correction coefficient calculated between the images of the area 32 and the area 33 is H2 (x, y). When the x-coordinate of the boundary 43 on the image in the region 32 is x1 and the x-coordinate of the boundary 45 on the image in the region 32 is x2, the following coefficient is used as the correction coefficient at the position of x. .
[0087]
(Equation 22)

[0088]
With this method, discontinuous luminance change can be eliminated at both the boundary 43 and the boundary 45. In the case of this procedure, when the illuminance unevenness correction coefficient calculated between the images of the area 32 and the area 33 is obtained, the image in the area 32 is pasted onto the panoramic image while correcting the illuminance unevenness.
[0089]
【The invention's effect】
According to the present invention, a natural panoramic image can be obtained by removing a discontinuous change in luminance in a mosaicing image caused by an illuminance distribution in an image plane. Moreover, since the original image luminance is restored by extracting only the uneven illuminance effect, the original image is not broken unlike frequency filtering. In addition, since the illuminance distribution in the image plane is considered, fine correction can be performed for each location, and the ability to remove discontinuous luminance changes is high. Therefore, this is an extremely effective method.
[Brief description of the drawings]
FIG. 1 is a block diagram illustrating a configuration of a composite image creating device including an image processing device of the present invention.
FIG. 2 is a functional block diagram for explaining an image processing operation of the image processing apparatus according to the present invention.
FIG. 3 is a diagram illustrating a photographing state in the mosaicing process of the present invention.
FIG. 4 is a diagram for explaining a luminance distribution state of each image in FIG. 3;
FIG. 5 is a flowchart illustrating an image processing method according to the present invention.
FIG. 6 is a diagram showing a state in which a differential image dC1 (x, y) in an area 31 is divided into blocks, and a matching target pattern 70 is selected.
FIG. 7 is a diagram showing a state when a matching coefficient is calculated by superimposing a matching target pattern 70 on a differential image dC2 (x, y) in an area 32.
FIG. 8 is a diagram illustrating a relationship between a position of an illumination and a position of a shooting target in the present embodiment.
FIG. 9 is a diagram showing a state in which uneven illuminance has been corrected.
FIG. 10 is a diagram showing a state where an uneven illuminance correction coefficient in an image plane of an area 32 is interpolated.
FIG. 11 is a diagram illustrating a photographing state in a conventional mosaicing process.
12 is a diagram for explaining a luminance distribution state of each image in FIG. 11 when there is uneven illuminance.
FIG. 13 is a diagram showing a result when the image of FIG. 12 is mosaiked by a conventional method.
[Explanation of symbols]
11 Composite image creation device
12 Image data capture device
13 memory
14 External storage device
15 Central processing circuit
16 Image processing device
17 Display device
51, 52, 53 Luminance distribution of each image in image sequence
90, 91, 92, 93 Calculated value of correction coefficient for illuminance unevenness of each image in image sequence
94, 95, 96, 97 Image luminance distribution as a result of correcting the illuminance unevenness of the image in the image sequence by the illuminance unevenness correction coefficient

Claims

An image processing method for continuously combining a plurality of digital images taken while moving a shooting place or a shooting direction so that a part of the taken digital images overlaps to synthesize a single panoramic image. hand,
An image moving amount calculating step of calculating a moving amount of the first digital image and the second digital image, each of which partially overlaps among the plurality of digital images taken while moving, on the image; Image points predicted to correspond to the same point based on the movement amount obtained in the movement amount calculation step are obtained for the first image and the second image, respectively, and the second image is obtained by using the brightness of those image points. An illuminance unevenness correction step of calculating an illuminance unevenness correction coefficient for correcting an illuminance difference between the first image and the second image;
When calculating the illuminance unevenness correction coefficient in the illuminance unevenness correction step, the illuminance unevenness correction coefficient is expressed by a linear combination of known functions using coordinate values on an image plane as variables, and a linear combination constant in the linear combination is statistically calculated. An image processing method characterized by a dynamic method.

2. The image processing method according to claim 1, wherein in the illuminance unevenness correction step, a polynomial function is used as the known function.

2. The image processing method according to claim 1, wherein in the illuminance unevenness correction step, a least square method is used as the statistical method for obtaining the linear combination constant.

2. The image processing method according to claim 1, wherein in the moving amount calculating step, the moving amount is estimated using a differential image.

The moving amount calculating step includes a first calculating step by a matching method for estimating a moving amount in a pixel unit which is a quantization unit of a screen, and a moving amount in a sub-pixel unit having higher accuracy than the quantization unit by a statistical method. 2. The image processing method according to claim 1, further comprising: a second calculating step based on a gradient method for estimating the gradient.

In the illuminance unevenness correction step, the illuminance unevenness correction coefficient in the screen is continuously changed based on a plurality of illuminance unevenness correction coefficients calculated by comparing two or more images with one image. 2. The image processing method according to claim 1, wherein:

An image processing apparatus that continuously joins a plurality of digital images taken while moving the shooting place or shooting direction so that a part of the shot digital images overlaps and synthesizes a single panoramic image. hand,
A central processing circuit for performing image processing for continuously connecting the plurality of digital images by the image processing method according to claim 1, and storing the digital image processed by the central processing circuit. An image processing apparatus comprising: