知財求人 - 知財ポータルサイト「IP Force」
▶ キヤノン株式会社の特許一覧 ▶ 国立大学法人 東京大学の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】5844884
(24)【登録日】2015年11月27日
(45)【発行日】2016年1月20日
(54)【発明の名称】画像処理装置、画像処理方法、及びプログラム
(51)【国際特許分類】
G06T 7/00 20060101AFI20151224BHJP
G06T 1/00 20060101ALI20151224BHJP
【ＦＩ】
G06T7/00 350A
G06T1/00 290A
【請求項の数】14
【全頁数】33
(21)【出願番号】特願2014-509530(P2014-509530)
(86)(22)【出願日】2012年9月26日
(65)【公表番号】特表2014-524602(P2014-524602A)
(43)【公表日】2014年9月22日
(86)【国際出願番号】JP2012075577
(87)【国際公開番号】WO2013047882
(87)【国際公開日】20130404
【審査請求日】2014年2月24日
(31)【優先権主張番号】特願2011-213381(P2011-213381)
(32)【優先日】2011年9月28日
(33)【優先権主張国】JP
(73)【特許権者】
【識別番号】000001007
【氏名又は名称】キヤノン株式会社
(73)【特許権者】
【識別番号】504137912
【氏名又は名称】国立大学法人 東京大学
(74)【代理人】
【識別番号】100076428
【弁理士】
【氏名又は名称】大塚 康徳
(74)【代理人】
【識別番号】100112508
【弁理士】
【氏名又は名称】高柳 司郎
(74)【代理人】
【識別番号】100115071
【弁理士】
【氏名又は名称】大塚 康弘
(74)【代理人】
【識別番号】100116894
【弁理士】
【氏名又は名称】木村 秀二
(74)【代理人】
【識別番号】100130409
【弁理士】
【氏名又は名称】下山 治
(74)【代理人】
【識別番号】100134175
【弁理士】
【氏名又は名称】永川 行光
(72)【発明者】
【氏名】石川 亮
(72)【発明者】
【氏名】佐藤 清秀
(72)【発明者】
【氏名】池内 克史
(72)【発明者】
【氏名】鄭 波
(72)【発明者】
【氏名】孫 永祺
【審査官】
新井 則和
(56)【参考文献】
【文献】
特開２００９−２７３５９７（ＪＰ，Ａ）
【文献】
国際公開第２００９／１３１２０９（ＷＯ，Ａ１）
【文献】
特開２００５−０８０７５７（ＪＰ，Ａ）
【文献】
国際公開第２００８／１３９８２５（ＷＯ，Ａ１）
(58)【調査した分野】（Int.Cl.，ＤＢ名）
Ｇ０６Ｔ １／００−７／６０
Ａ６１Ｂ ６／００
Ａ６１Ｂ ８／００
(57)【特許請求の範囲】
【請求項1】
画像の特徴量を算出する画像処理装置であって、
処理画像を入力する第１の入力手段と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出手段と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
処理画像を入力する第１の入力手段と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出手段と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
【請求項2】
着目次数に対応する前記寄与度は、該着目次数を持つそれぞれの単項式による前記処理画像の近似精度に関する寄与度の平方の和であることを特徴とする、請求項１に記載の画像処理装置。
【請求項3】
前記単項式による前記処理画像の近似精度に関する寄与度は、前記複数の単項式のそれぞれに対応する基底によって規定される直交空間に前記強度値を射影した際における、当該単項式に対応する座標値に対応することを特徴とする、請求項１又は２に記載の画像処理装置。
【請求項4】
画像の特徴量を算出する画像処理装置であって、
処理画像を入力する第１の入力手段と、
算出手段であって、
複数の次数のうちの何れかの次数を持つ複数の単項式で構成される近似多項式であって、前記処理画像の画素位置を入力とし且つ該画素位置における画素値を出力とする前記近似多項式について、各単項式に対応する直交ベクトルを取得する取得手段と、
前記各単項式に対応する直交ベクトルを要素とする行列と、前記処理画像中の各画素位置における画素値を要素とする行列との積演算により、各単項式に対応する特徴量を算出する第１の演算手段と、
着目次数を持つ１以上の単項式に対応する前記特徴量から、該着目次数に対応する特徴量を算出する第２の演算手段と、
を備える算出手段と、
それぞれの次数に対応する前記特徴量のセットを、前記処理画像の特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
処理画像を入力する第１の入力手段と、
算出手段であって、
複数の次数のうちの何れかの次数を持つ複数の単項式で構成される近似多項式であって、前記処理画像の画素位置を入力とし且つ該画素位置における画素値を出力とする前記近似多項式について、各単項式に対応する直交ベクトルを取得する取得手段と、
前記各単項式に対応する直交ベクトルを要素とする行列と、前記処理画像中の各画素位置における画素値を要素とする行列との積演算により、各単項式に対応する特徴量を算出する第１の演算手段と、
着目次数を持つ１以上の単項式に対応する前記特徴量から、該着目次数に対応する特徴量を算出する第２の演算手段と、
を備える算出手段と、
それぞれの次数に対応する前記特徴量のセットを、前記処理画像の特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
【請求項5】
前記処理画像は３次元画像であって、
前記取得手段が取得する各単項式に対応する前記直交ベクトルは、下記行列MをQR分解して得られる直交行列Qに含まれる行ベクトルであり、
前記第１の演算手段は、xlymznを含む単項式に対応する前記特徴量として、下式に従って値clmnを算出し、
前記第２の演算手段は、次数pに対応する前記特徴量として、下式に従って値fpを算出することを特徴とし、
前記処理画像はN個の画素を含み、該N個の画素のうちのk番目の画素の座標値はxk,yk,zkであり、k番目の画素の画素値はbkであり、列ベクトルm(xk,yk,zk)はxklykmzkn(l+m+n≦Ndegree, l,m,nは0以上の整数)を成分として有し、Ndegreeは前記近似多項式の次数である、請求項４に記載の画像処理装置。
前記取得手段が取得する各単項式に対応する前記直交ベクトルは、下記行列MをQR分解して得られる直交行列Qに含まれる行ベクトルであり、
【請求項6】
画像の特徴量を算出する画像処理装置であって、
処理画像を入力する第１の入力手段と、
算出手段であって、
複数の次数のうちの何れかの次数を持つ複数の単項式で構成される近似多項式であって、前記処理画像の画素位置を入力とし且つ該画素位置における画素値を出力とする前記近似多項式について、各単項式に対応する直交ベクトルを要素とする第１の行列の各要素の値を示す第１の情報と、第２の行列の各要素の値を示す第２の情報であって、前記第１の行列と前記第２の行列との積演算によって前記処理画像の各画素位置の座標値同士の積を要素として含む行列が得られる第２の情報と、を取得する取得手段と、
前記第１の行列と前記第２の行列と前記処理画像中の各画素位置における画素値を要素とする行列とから、前記単項式のそれぞれの係数を算出し、前記第２の行列と前記単項式のそれぞれの係数を要素とする行列との積演算により、各単項式に対応する特徴量を算出する第１の演算手段と、
着目次数を持つ１以上の単項式に対応する前記特徴量から、該着目次数に対応する特徴量を算出する第２の演算手段と、
を備える算出手段と、
それぞれの次数に対応する前記特徴量のセットを、前記処理画像の特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
処理画像を入力する第１の入力手段と、
算出手段であって、
複数の次数のうちの何れかの次数を持つ複数の単項式で構成される近似多項式であって、前記処理画像の画素位置を入力とし且つ該画素位置における画素値を出力とする前記近似多項式について、各単項式に対応する直交ベクトルを要素とする第１の行列の各要素の値を示す第１の情報と、第２の行列の各要素の値を示す第２の情報であって、前記第１の行列と前記第２の行列との積演算によって前記処理画像の各画素位置の座標値同士の積を要素として含む行列が得られる第２の情報と、を取得する取得手段と、
前記第１の行列と前記第２の行列と前記処理画像中の各画素位置における画素値を要素とする行列とから、前記単項式のそれぞれの係数を算出し、前記第２の行列と前記単項式のそれぞれの係数を要素とする行列との積演算により、各単項式に対応する特徴量を算出する第１の演算手段と、
着目次数を持つ１以上の単項式に対応する前記特徴量から、該着目次数に対応する特徴量を算出する第２の演算手段と、
を備える算出手段と、
それぞれの次数に対応する前記特徴量のセットを、前記処理画像の特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
【請求項7】
前記処理画像は３次元画像であって、
前記処理画像の各画素位置の座標値同士の積を要素として含む前記行列は下式の行列Mであり、前記第１の行列は行列MをQR分解して得られる下式の直交行列Qであり、前記第２の行列は下式の行列Rであり、
前記第１の演算手段は、下式に従ってxlymznを含む単項式の係数almnを算出し、
前記第１の演算手段はさらに、xlymznを含む単項式に対応する前記特徴量として、下式に従って値clmnを算出し、
前記第２の演算手段は、次数pに対応する特徴量として、下式に従って値fpを算出することを特徴とし、
前記処理画像はN個の画素を含み、該N個の画素のうちのk番目の画素の座標値はxk,yk,zkであり、k番目の画素の画素値はbkであり、列ベクトルm(xk,yk,zk)はxklykmzkn(l+m+n≦Ndegree, l,m,nは0以上の整数)を成分として有し、Ndegreeは前記近似多項式の次数である、請求項６に記載の画像処理装置。
前記処理画像の各画素位置の座標値同士の積を要素として含む前記行列は下式の行列Mであり、前記第１の行列は行列MをQR分解して得られる下式の直交行列Qであり、前記第２の行列は下式の行列Rであり、
【請求項8】
第１の画像と第２の画像とを入力する第２の入力手段と、
前記第１の画像の画像領域内に複数の第１の部分領域群を設定する第１の設定手段と、
前記第１の画像から前記第１の部分領域のそれぞれに含まれる画像を第１の部分画像群として抽出する第１の抽出手段と、
前記第２の画像の画像領域内に複数の第２の部分領域群を設定する第２の設定手段と、
前記第２の画像から前記第２の部分領域のそれぞれに含まれる画像を第２の部分画像群として抽出する第２の抽出手段と、
前記第１の部分画像のそれぞれについて前記算出手段が算出しかつ前記第１の出力手段が出力した特徴量を取得し、及び、前記第２の部分画像のそれぞれについて前記算出手段が算出しかつ前記第１の出力手段が出力した特徴量を取得して、前記第１の部分画像のそれぞれについて、当該第１の部分画像の特徴量に最も類似している特徴量が算出された第２の部分画像を前記第２の部分画像群から選択し当該第１の着目部分画像に対応付ける対応付け手段と、
前記対応付け手段が対応付けた前記第１の部分画像と前記第２の部分画像との組のそれぞれについて、当該第１の部分画像が抽出された第１の部分領域と、当該第２の部分画像が抽出された第２の部分領域とが近くなるように、前記第１の画像と前記第２の画像とを位置合わせする位置合わせ手段と、
前記位置合わせ手段によって位置合わせされた前記第１の画像と前記第２の画像とを出力する第２の出力手段と、
を備えることを特徴とする請求項１乃至７の何れか１項に記載の画像処理装置。
前記第１の画像の画像領域内に複数の第１の部分領域群を設定する第１の設定手段と、
前記第１の画像から前記第１の部分領域のそれぞれに含まれる画像を第１の部分画像群として抽出する第１の抽出手段と、
前記第２の画像の画像領域内に複数の第２の部分領域群を設定する第２の設定手段と、
前記第２の画像から前記第２の部分領域のそれぞれに含まれる画像を第２の部分画像群として抽出する第２の抽出手段と、
前記第１の部分画像のそれぞれについて前記算出手段が算出しかつ前記第１の出力手段が出力した特徴量を取得し、及び、前記第２の部分画像のそれぞれについて前記算出手段が算出しかつ前記第１の出力手段が出力した特徴量を取得して、前記第１の部分画像のそれぞれについて、当該第１の部分画像の特徴量に最も類似している特徴量が算出された第２の部分画像を前記第２の部分画像群から選択し当該第１の着目部分画像に対応付ける対応付け手段と、
前記対応付け手段が対応付けた前記第１の部分画像と前記第２の部分画像との組のそれぞれについて、当該第１の部分画像が抽出された第１の部分領域と、当該第２の部分画像が抽出された第２の部分領域とが近くなるように、前記第１の画像と前記第２の画像とを位置合わせする位置合わせ手段と、
前記位置合わせ手段によって位置合わせされた前記第１の画像と前記第２の画像とを出力する第２の出力手段と、
を備えることを特徴とする請求項１乃至７の何れか１項に記載の画像処理装置。
【請求項9】
前記対応付け手段は、
前記第１の部分画像群に含まれる第１の着目部分画像の前記特徴量と、前記第１の部分画像群に含まれる当該第１の着目部分画像以外の部分画像のそれぞれの前記特徴量との差異の最小値が、所定の閾値以上であるか否かを判定し、
所定の閾値以上であると判定した場合に、前記第１の着目部分画像に対して前記第２の部分画像を対応付ける
ことを特徴とする請求項８に記載の画像処理装置。
前記第１の部分画像群に含まれる第１の着目部分画像の前記特徴量と、前記第１の部分画像群に含まれる当該第１の着目部分画像以外の部分画像のそれぞれの前記特徴量との差異の最小値が、所定の閾値以上であるか否かを判定し、
所定の閾値以上であると判定した場合に、前記第１の着目部分画像に対して前記第２の部分画像を対応付ける
ことを特徴とする請求項８に記載の画像処理装置。
【請求項10】
画像の特徴量を算出する画像処理装置であって、
処理画像を入力する第１の入力手段と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数毎の誤差に対応した前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出手段と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
処理画像を入力する第１の入力手段と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数毎の誤差に対応した前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出手段と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力手段と、
を備えることを特徴とする画像処理装置。
【請求項11】
画像の特徴量を算出する画像処理装置であって、
処理画像を入力する第１の入力ユニットと、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出ユニットと、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力ユニットと、
を備えることを特徴とする画像処理装置。
処理画像を入力する第１の入力ユニットと、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出ユニットと、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力ユニットと、
を備えることを特徴とする画像処理装置。
【請求項12】
画像の特徴量を算出する画像処理装置が行う画像処理方法であって、
処理画像を入力する第１の入力工程と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出工程と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力工程と、
を含むことを特徴とする画像処理方法。
処理画像を入力する第１の入力工程と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出工程と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力工程と、
を含むことを特徴とする画像処理方法。
【請求項13】
画像の特徴量を算出する画像処理装置が行う画像処理方法であって、
処理画像を入力する第１の入力工程と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数毎の誤差に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出工程と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力工程と、
を含むことを特徴とする画像処理方法。
処理画像を入力する第１の入力工程と、
前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数毎の誤差に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出工程と、
前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力工程と、
を含むことを特徴とする画像処理方法。
【請求項14】
コンピュータを、請求項１乃至１０の何れか１項に記載の画像処理装置が有する各手段として機能させるための、コンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、画像の特徴量を算出する画像処理技術に関する。
【背景技術】
【0002】
医療の分野において、医用画像収集装置で患者を撮影して得られた医用画像を読影することにより、医師は病変部の位置や状態、経時変化などを観察する。医用画像を生成する装置としては、単純Ｘ線撮影装置、Ｘ線コンピュータ断層撮影装置（Ｘ線ＣＴ）、核磁気共鳴映像装置（ＭＲＩ）、超音波画像診断装置（ＵＳ）などがあげられる。それぞれの装置は特性が異なるため、通常医師は、撮影する部位や疾病などに適した複数の装置を選択する。例えば患者のＭＲＩを撮影し、ＭＲＩの画像を参照しながら超音波画像を撮影することにより、病変部の位置や拡がりなどの、診断に有効な情報を得ることができる。
【0003】
このように様々な時期に、様々な装置で医用画像を撮影した場合に、それらの画像を対応付けることにより、画像を用いた診断をより容易にできる。画像を対応付ける方法の１つとして、多次元画像間の対応付けに関する技術が、Warren Cheung and Ghassan Hamarneh, “n-SIFT: n-Dimensional Scale Invariant Feature Transform,” IEEE Transaction on Image Processing, Vol. 18, No. 9, September 2009（非特許文献１）に開示されている。この技術は、多次元画像の局所領域における強度勾配ヒストグラムを、局所領域における画像の特徴量として算出する。他の画像について同様に求められた特徴量と比較することにより、画像間の対応付けを行うことができる。これにより、種々の撮影機器で撮影された医用画像間の対応部位を対応づけて観察することが可能となる。
【発明の概要】
【0004】
非特許文献１に記載の技術では、高い次元数の特徴量を算出する。したがって、特徴量の算出処理及び特徴量を用いた対応付け処理に、大きな計算量を必要とする。特に画像の対応付けをリアルタイムに行う場合に、処理負荷を軽減することが求められていた。また、非特許文献１に記載の技術には、画像間の回転や画像ノイズの影響により対応付けの精度が低下するという課題があった。
【0005】
本発明は、画像の対応付けを行うために利用可能な特徴量をより高速に算出する技術を提供する。
【0006】
本発明の１つの態様によれば、画像の特徴量を算出する画像処理装置は、処理画像を入力する第１の入力手段と、前記処理画像の画素位置と該画素位置における強度値との関係を示す近似多項式であって、複数の次数のうちの何れかの次数を持つ複数の単項式で構成される前記近似多項式を用いて前記処理画像を近似した場合に、前記近似多項式のそれぞれの次数に対応した単項式による前記処理画像の近似精度に関する寄与度によって構成される数値のセットを算出する算出手段と、前記算出された数値のセットを前記処理画像の前記特徴量として出力する第１の出力手段と、を備える。
【0011】
本発明の一実施形態によれば、画像の対応付けを行うために利用可能な特徴量をより高速に算出することができる。
【0012】
添付の図面を参照して、本発明のさらなる特徴が、以下の例示的な実施形態の説明から明らかとなるだろう。
【図面の簡単な説明】
【0013】
【図1】実施例１に係る画像処理システムの機能構成を示すブロック図である。
【図2】実施例１に係る画像処理システムの装置構成を示すブロック図である。
【図3】実施例１に係る画像処理装置１０の処理手順を示すフローチャートである。
【図4】実施例１に係るウィンドウ画像を説明する図である。
【図5】実施例２に係る画像処理システムの機能構成を示すブロック図である。
【図6】実施例２に係る画像処理装置２０の処理手順を示すフローチャートである。
【図7】実施例３に係る画像処理システムの機能構成を示すブロック図である。
【図8】実施例３に係る画像処理装置３０の処理手順を示すフローチャートである。
【図9】実施例４に係る画像処理システムの機能構成を示すブロック図である。
【図10】実施例４に係る画像処理装置４０の処理手順を示すフローチャートである。
【図11A】実施例４に係るウィンドウ領域の例を示す図である。
【図11B】実施例４に係るウィンドウ領域の例を示す図である。
【図11C】実施例４に係るウィンドウ領域の例を示す図である。
【図11D】実施例４に係るウィンドウ領域の例を示す図である。
【図11E】実施例４に係るウィンドウ領域の例を示す図である。
【図11F】実施例４に係るウィンドウ領域の例を示す図である。
【図11G】実施例４に係るウィンドウ領域の例を示す図である。
【図11H】実施例４に係るウィンドウ領域の例を示す図である。
【図11I】実施例４に係るウィンドウ領域の例を示す図である。
【発明を実施するための形態】
【0014】
以下、添付図面に従って本発明に係る画像処理装置及び方法の実施例について詳説する。ただし、発明の範囲は以下の例に限定されるものではない。
【0015】
［実施例１］以下に、実施例１に係る画像処理装置１０について説明する。本実施例の第１多項式近似部１０３及び第１特徴量算出部１０４は、処理対象となる画像についての特徴量を算出することができる。算出された特徴量を用いることにより、対応算出部１１５が行う方法と同様に、画像の対応付けを行うことができる。本実施例において算出される特徴量は画像の回転にかかわらず不変であり、画像が回転している場合にもより高い精度で画像の関連付けを行うことができる。
【0016】
本実施例において画像処理装置１０は、第１多項式近似部１０３及び第１特徴量算出部１０４が算出する特徴量を用いることにより、画像の位置合わせを行う。しかしながら、本発明は画像の位置合わせを行う技術には限定されない。例えば、第１多項式近似部１０３及び第１特徴量算出部１０４以外の各部が行う処理を、画像処理装置１０とは異なる装置が実行してもよい。また、画像処理装置１０は、画像の位置合わせを行うのではなく、算出した特徴量を用いることにより複数の画像の中から類似する画像の組を抽出することもできる。
【0017】
実施例１に係る画像処理装置１０は、第１の画像及び第２の画像を取得し、位置合わせを行う。例えば第１の画像及び第２の画像は、同一被検体の略同一部位を撮影した二つの三次元画像でありうる。以下に、本実施例に係る画像処理装置が行う処理の概略を説明する。まず画像処理装置１０は、第１の画像における複数の部分領域（ウィンドウ）の特徴量を算出する。次に画像処理装置１０は、第２の画像についても同様に複数のウィンドウの特徴量を算出する。そして画像処理装置１０は、算出された特徴量に基づいて、第１の画像内のウィンドウを、第２の画像内のウィンドウに対応付ける。さらに画像処理装置１０は、この対応付けの結果に基づいて、第１の画像と第２の画像との間の位置合わせ処理を実行する。画像処理装置１０はさらに、２つの画像を比較して観察する表示画像を生成することもできる。
【0018】
次に図１を参照して、本実施例に係る画像処理装置１０の構成を説明する。画像処理装置１０は、第１画像取得部１０１、第１ウィンドウ画像取得部１０２、第１多項式近似部１０３、第１特徴量算出部１０４、及び特徴ウィンドウ選択部１０５を持つ。画像処理装置１０はさらに、第２画像取得部１１１、第２ウィンドウ画像取得部１１２、第２多項式近似部１１３、第２特徴量算出部１１４、対応算出部１１５、及び画像生成部１１６を持つ。
【0019】
第１画像取得部１０１は、第１の画像を取得する（第２の入力手段）。第１画像取得部１０１は例えば、画像記録装置１２のような記憶装置から第１の画像を取得することができる。第１画像取得部１０１はまた、画像撮影装置１１のような画像生成装置から第１の画像を取得してもよい。そして第１画像取得部１０１は、取得した第１の画像を第１ウィンドウ画像取得部１０２および画像生成部１１６へと出力する。
【0020】
第１ウィンドウ画像取得部１０２は、第１画像取得部１０１から取得した第１の画像から複数の部分領域を切り出すことにより、第１のウィンドウ画像群を生成する（第１の設定手段及び第１の抽出手段）。
【0021】
第１多項式近似部１０３は、第１ウィンドウ画像取得部１０２が生成した第１のウィンドウ画像群のそれぞれについて、ウィンドウ画像を近似する多項式を算出する。そして第１多項式近似部１０３は、算出した多項式に関する情報を第１特徴量算出部１０４へと出力する。多項式に関する情報の具体的な内容については、後に詳しく説明する。
【0022】
第１特徴量算出部１０４は、第１多項式近似部１０３から取得した多項式に関する情報に基づき、第１のウィンドウ画像群に含まれるそれぞれのウィンドウ画像（処理画像）の特徴量を算出する（算出手段）。特徴量の算出方法については、後に詳しく説明する。そして第１特徴量算出部１０４は、算出した特徴量を特徴ウィンドウ選択部１０５へと出力する。
【0023】
特徴ウィンドウ選択部１０５は、第１特徴量算出部１０４が算出した特徴量に基づき、第１のウィンドウ画像群から特徴的なウィンドウ画像を選択する。そして特徴ウィンドウ選択部１０５は、選択結果を対応算出部１１５へと出力する。
【0024】
本実施例においては、第１ウィンドウ画像取得部１０２は複数のウィンドウ画像を生成する。そして特徴ウィンドウ選択部１０５は、複数のウィンドウ画像から位置合わせに適したウィンドウ画像を選択する。しかしながら実施例はこれに限定されない。例えば、第１ウィンドウ画像取得部１０２が生成したウィンドウ画像の全てを、特徴ウィンドウとして用いてもよい。この場合第１ウィンドウ画像取得部１０２が生成するウィンドウ画像は１つ以上であればよい。画像の位置合わせのためには、第１ウィンドウ画像取得部１０２は通常、画像の次元以上の数のウィンドウ画像を生成する。例えば第１の画像が３次元画像である場合、第１ウィンドウ画像取得部１０２は３つ以上のウィンドウ画像を生成すればよい。
【0025】
第２画像取得部１１１は、第２の画像を取得する。第２画像取得部１１１は例えば、画像記録装置１２のような記憶装置から第２の画像を取得することができる。第２画像取得部１１１はまた、画像撮影装置１１のような画像生成装置から第２の画像を取得してもよい。そして第２画像取得部１１１は、取得した第２の画像を第２ウィンドウ画像取得部１１２および画像生成部１１６へと出力する。
【0026】
第１ウィンドウ画像取得部１０２は、第２画像取得部１１１から取得した第２の画像から複数の部分領域を切り出すことにより、第２のウィンドウ画像群を生成する（第２の設定手段及び第２の抽出手段）。
【0027】
第２多項式近似部１１３は、第１多項式近似部１０３と同様に、第２ウィンドウ画像取得部１１２が生成した第２のウィンドウ画像群のそれぞれについて、ウィンドウ画像を近似する多項式を算出する。そして第２多項式近似部１１３は、算出した多項式に関する情報を第２特徴量算出部１１４へと出力する。
【0028】
第２特徴量算出部１１４は、第１特徴量算出部１０４と同様に、第２多項式近似部１１３から取得した多項式に関する情報に基づき、第２のウィンドウ画像群に含まれるそれぞれのウィンドウ画像の特徴量を算出する。そして第２特徴量算出部１１４は、算出した特徴量を対応算出部１１５へと出力する。
【0029】
対応算出部１１５は特徴ウィンドウ選択部１０５および第２特徴量算出部１１４から取得した情報に基づき、第１のウィンドウ画像群と第２のウィンドウ画像群との間の対応関係を算出する（対応付け手段）。対応算出部１１５が行う処理の詳細については後述する。そして対応算出部１１５は、算出した対応関係を画像生成部１１６へ出力する。
【0030】
画像生成部１１６は、対応算出部１１５が算出した対応関係と、第１画像取得部１０１が取得した第１の画像と、第２画像取得部１１１が取得した第２の画像と、に基づいて出力画像を生成する（位置合わせ手段及び第２の出力手段）。画像生成部１１６は、生成した出力画像を表示装置（不図示）に出力してもよい。
【0031】
図２は、本実施例に係る画像処理装置１０を含む、本実施例に係る画像処理システムの構成を示す。本実施例に係る画像処理システムは、画像処理装置１０、画像記録装置１２、ローカルエリアネットワーク（ＬＡＮ）８４０、及び画像撮影装置１１を含む。
【0032】
画像撮影装置１１は画像を撮影可能な装置であって、本実施例においては超音波画像撮影装置である。本実施例における画像撮影装置１１は、超音波を送受信する不図示の超音波プローブを被検体に接触させることで、被検体の内部を超音波撮影する。本実施例で画像撮影装置１１が撮影する画像は、被検体の所定の３次元領域を撮影した３次元のＢモード超音波画像である。
【0033】
画像撮影装置１１は、撮影した画像を画像記録装置１２に格納する。画像記録装置１２は、画像を読み書きすることができる装置である。画像記録装置１２は、画像処理装置１０からの要求に従って、格納している第１の画像及び第２の画像を読み出し及び送信することができる。画像撮影装置１１及び画像記録装置１２は、図２においてはＬＡＮ８４０を介して画像処理装置１０に接続されているが、画像処理装置１０の共通バス８１８（後述）に直接又はインタフェースを介して接続されていてもよい。
【0034】
本実施例に係る画像処理装置１０は、一般的なパーソナルコンピュータ（ＰＣ）を用いて実現することができる。もっとも、画像処理装置１０は、図１に示した各部の機能を実現する専用のハードウェア要素によって実現されてもよい。例えば、第１多項式近似部１０３と第１特徴量算出部１０４とは、画像を入力する第１の入力ユニットと、特徴量の演算を行う算出ハードウェアと、特徴量を出力する第１の出力ユニットと、で実現されてもよい。
【0035】
図２は、コンピュータによって実現されている画像処理装置１０の電気的構成をさらに示す。画像処理装置１０は、中央演算処理装置（ＣＰＵ）８１１、主メモリ８１２、磁気ディスク８１３、表示メモリ８１４、モニタ８１５、マウス８１６、及びキーボード８１７を有する。ＣＰＵ８１１は、主として画像処理装置１０の各構成要素の動作を制御する。
【0036】
主メモリ８１２は、ＣＰＵ８１１が実行する制御プログラムを格納したり、ＣＰＵ８１１によるプログラム実行時の作業領域を提供したりする。磁気ディスク８１３は、オペレーティングシステム（ＯＳ）、周辺機器のデバイスドライブ、後述する位置合わせ処理等を行うためのプログラムを含む各種アプリケーションソフト等を格納する。すなわち、磁気ディスク８１３などの記憶媒体に格納されているコンピュータプログラムを主メモリ８１２に読み出し、このプログラムに従ってＣＰＵ８１１が画像処理装置１０の各部を制御することにより、本実施例に係る処理は実現される。
【0037】
表示メモリ８１４は、モニタ８１５のための表示用データを一時記憶する。モニタ８１５は、例えばＣＲＴモニタや液晶モニタ等であり、表示メモリ８１４に入力されたデータに基づいて画像を表示する。マウス８１６及びキーボード８１７は入力インタフェースであり、これらを用いてユーザはポインティング入力、文字入力、コマンド入力などを行う。これらの構成要素は、共通バス８１８により互いに通信可能に接続されている。
【0038】
次に、画像処理装置１０が行う全体の処理を、図３のフローチャートを用いて詳しく説明する。本実施例では、主メモリ８１２に格納されている各部の機能を実現するプログラムをＣＰＵ８１１が実行することにより実現される。また以下に説明する画像処理装置１０が行う各処理の結果は、主メモリ８１２に格納されることにより記録される。
【0039】
（ステップＳ１０１０）第１の画像の取得ステップＳ１０１０において第１画像取得部１０１は、第１の画像を取得する。第１の画像は、画像撮影装置１１が被検体を撮影することによって得られた画像であり、画像記録装置１２に記録されている。ここで第１の画像は、被検体の３次元領域を超音波撮影することによって得られた３次元画像である。第１の画像は、画像撮影装置１１の基準座標系における複数の３次元位置座標ベクトルX1,iと、それぞれの座標ベクトルに対応する強度値b1,iとの組の集合である。ここで、強度値は画素値と同義であり、通常は画素についての色情報である。ここで、３次元位置座標ベクトルX1,i=[x1_i y1_i z1_i]T (1≦i≦N1)とする。また、1≦i≦N1であり、N1は第１の画像を構成する画素の数である。
【0040】
第１の画像は、互いに直交する３次元格子の格子点に対応する座標と、それぞれの座標における強度値との組の集合であってもよい。また第１の画像は、３次元の位置姿勢を有する２次元の平面画像が複数集合したものであってもよい。例えば、２次元の超音波プローブに位置姿勢センサを装着して複数の２次元画像を撮影することにより、後者の３次元画像を取得することができる。
【0041】
（ステップＳ１０２０）第１ウィンドウ画像群の生成ステップＳ１０２０において、第１ウィンドウ画像取得部１０２は、第１画像取得部１０１から取得した第１の画像から複数の部分領域（ウィンドウ）を切り出すことにより、第１のウィンドウ画像群を生成する。この処理の具体的な内容を図４を使って説明する。
【0042】
図４において第１の画像の領域９００は、ステップＳ１０１０で取得した第１の画像の全体の領域である。本実施例における第１の画像は３次元画像であるが、説明の都合上、図４において第１の画像は２次元画像として表されている。
【0043】
ステップＳ１０２０では、この第１の画像の領域９００内（画像領域内）に、複数のウィンドウ領域９０１（第１の部分領域群）を設定する。以下で、複数のウィンドウ領域９０１を合わせてウィンドウ領域群と呼ぶ。図４に示される例では、ウィンドウ領域９０１のそれぞれが等間隔で配置されている。また、それぞれのウィンドウ領域９０１の大きさは等しく、それぞれのウィンドウ領域９０１は互いに接するように配置されている。しかしながら、ウィンドウ領域の配置方法は図４の例に限定されない。例えば、それぞれのウィンドウ領域９０１が互いに重なってもよいし、それぞれのウィンドウ領域９０１は互いに離れていてもよい。もっとも、ウィンドウ領域９０１の形状は任意である。
【0044】
第１ウィンドウ画像取得部１０２は、第１の画像から、それぞれのウィンドウ領域９０１内の画像を切り出す（抽出する）。切り出された画像を第１のウィンドウ画像と呼ぶ。このように、複数のウィンドウ領域９０１のそれぞれに対応する第１のウィンドウ画像を、第１ウィンドウ画像取得部１０２は生成する。以下で、複数の第１のウィンドウ画像を合わせて第１のウィンドウ画像群（第１の部分画像群）と呼ぶ。
【0045】
本実施例において、ウィンドウ領域群に含まれる各ウィンドウ領域９０１は、３次元空間において球形状を有する領域である。図４においては説明の都合上、各ウィンドウ領域９０１は円形で表されている。ここで、各ウィンドウ領域は直径がd[mm]の球であるものとする。また、j番目のウィンドウの中心位置（以降でウィンドウ中心位置と呼ぶ）を、３次元ベクトルXWC1,j=[xWC1,j yWC1,j zWC1,j]T (1≦j≦NW1)で表す。NW1は第１の画像から切り出されたウィンドウ画像の総数である。
【0046】
第１のウィンドウ画像のそれぞれは複数の画素を含む。それぞれの画素は、ウィンドウ内における画素の位置を表す３次元ベクトルと、強度値とによって表現される。ここで、j番目のウィンドウにおけるk番目の画素の位置を表す３次元ベクトルを、XW1,j,k=[xW1,j,k yW1,j,k zW1,j,k]T (1≦j≦NW1, 1≦k≦NPN1,j)とする。３次元ベクトルXW1,j,kは、各ウィンドウの中心位置XWC1,jを原点（基準点）とした位置ベクトルであるものとする。また、３次元ベクトルXW1,j,k（xk）に対応する画素の強度値を、bW1,j,k（bk）(1≦j≦NW1, 1≦k≦NPN1,j)とする。j番目のウィンドウ画像の画素数は、NPN1,j個（N個）である。
【0047】
第１ウィンドウ画像取得部１０２は、以上のように生成したそれぞれの第１のウィンドウ画像について、ウィンドウ中心位置XWC1,j(1≦j≦NW1)を記録する。また、それぞれの第１のウィンドウ画像について、画素位置XW1,j,k(1≦j≦NW1, 1≦k≦NPN1,j)と、画素の強度値bW1,j,k(1≦j≦NW1, 1≦k≦NPN1,j)との組を記録する。
【0048】
本実施例においては、３次元ベクトルXW1,j,kを、各ウィンドウの中心位置XWC1,jを原点（基準点）とした位置ベクトルとした。しかしながら、位置ベクトルの原点はウィンドウの中心に限定されない。むしろ位置ベクトルの原点は任意の位置に設定することができる。この場合、ウィンドウの中心位置の代わりに、位置ベクトルの原点の位置を、各ウィンドウ画像について記録すればよい。
【0049】
（ステップＳ１０３０）第１ウィンドウ画像群の多項式近似ステップＳ１０３０において、第１多項式近似部１０３は、ステップＳ１０２０で生成された第１のウィンドウ画像のそれぞれについて、多項式を算出する。具体的には第１多項式近似部１０３は、第１のウィンドウ画像のそれぞれを取得する（第１の入力手段）。そして第１多項式近似部は、画素位置XW1,j,kと画素の強度値bW1,j,kとの関係を近似する多項式を、第１のウィンドウ画像のそれぞれについて算出する。ここで、1≦j≦NW1、かつ1≦k≦NPN1,jとする。そして第１多項式近似部１０３は、算出した多項式を示す情報を出力する。
【0050】
以下に、ステップＳ１０３０について詳細に説明する。変数をx,y,zとするNdegree次元の多項式は、下式（１）を用いて一般的に表現できる。本実施例において画像は３次元画像であるから、この多項式を用いて画像を近似することができる。画像が２次元画像であれば変数をx,yとしてもよいし、画像が４次元画像であれば変数をx,y,z,tとしてもよい。【数1】
【0051】
式（１）においてp(x,y,z)は、座標(x,y,z)の画素について多項式によって求められる強度値（近似値）である。また、m(x,y,z)は座標値(x,y,z)に基づいて定まるベクトルを得る関数である。具体的にはm(x,y,z)は列ベクトルを示す。画像中のk番目の画素についての列ベクトルm(xk,yk,zk)は、xklykmzkn(l+m+n≦Ndegree, l,m,nは0以上の整数)を成分として有する。
【0052】
列ベクトルm(x,y,z)の各要素の値は、多項式を構成する複数の単項式について、係数を除いた（係数が乗じられていない）値である。言い換えれば列ベクトルm(x,y,z)は、ウィンドウ画像の各画素の座標値同士の積を要素として含む。さらに、almnは、多項式を構成する単項式のうち、xlymznの項の係数である。多項式を算出することは、係数almnのそれぞれを算出することに相当する。
【0053】
具体的には第１多項式近似部１０３は、第１のウィンドウ画像群のうちのj番目のウィンドウ画像についての多項式の係数ベクトルaW1,jを、式（２）の等式の二乗誤差が最小となるように算出する。【数2】
【0054】
一例として、係数ベクトルaW1,jは、以下のように求めることができる。j番目のウィンドウ画像について、1番目からNPN1,j番目までの画素についての式（２）から、式（３）を導くことができる。【数3】
【数4】
【0055】
次に第１多項式近似部１０３は式（５）のように、単項式行列MW1,jを、行列QW1,j（第１の行列）及び行列RW1,j（第２の行列）へとQR分解する。【数5】
【0056】
最後に、式（４）及び式（５）から導かれる式（６）を用いて、後退代入法により係数aW1,jを算出する。【数6】
【0057】
このようにして、係数aW1,jを持つ式（１）の近似多項式が求められる。第１多項式近似部１０３は、ステップＳ１０２０で生成した第１のウィンドウ画像群に含まれる全てのウィンドウ画像について、係数aW1,j、行列QW1,j、及び行列RW1,j(1≦j≦NW1)を以上のように算出する。そして第１多項式近似部１０３は、算出された係数aW1,j、行列QW1,j、及び行列RW1,jを、ウィンドウ画像を近似する多項式に関する情報として記録する。
【0058】
（ステップＳ１０４０）第１ウィンドウ画像群についての特徴量の算出ステップＳ１０４０において第１特徴量算出部１０４は、ステップＳ１０３０で算出された多項式に関する情報に基づいて、第１のウィンドウ画像群に含まれるそれぞれのウィンドウ画像の特徴量を算出する。具体的には第１特徴量算出部１０４は、それぞれのウィンドウ画像について、式（７）に従う積演算によりベクトルcW1,j(1≦j≦NW1)を算出する。
【数7】
【0059】
また、ベクトルcW1,jの算出方法は式（７）に限られない。例えば式（５）の関係が近似的に成り立つため、ベクトルcW1,jを式（８）に従う積演算により算出することもできる。式（７）と式（８）のいずれを用いても、略同一のベクトルcW1,j(1≦j≦NW1)を算出することができる（第１の演算手段）。【数8】
【0060】
ここで、行列QW1,jは直交行列であるから、行列QTW1,jの各行は直交基底（直交ベクトル）となる。ところで、行列QTW1,jの各行は、それぞれの単項式の係数を除いた値、すなわちxlymzn(l+m+n≦Ndegree, l,m,nは0以上の整数)のそれぞれに対応する。言い換えれば行列QTW1,jの各行は、複数の単項式のうちxlymznの項に対応する基底となる。
【0061】
ベクトルcW1,j=QTW1,jbW1,jは、行列QTW1,jによって規定される直交空間における強度値ベクトルbW1,jの座標値を示す。すなわち、xlymznの項に対応するベクトルcW1,jの成分cW1,j,lmn（値clmn）は、xlymznの項が強度値ベクトルbW1,jに与える影響を示す。すなわちcW1,j,lmnは、xlymznの項の、近似多項式によって得られる強度値への貢献度（寄与度）を示す。それゆえcW1,j,lmnは、xlymznの項の近似精度に対する貢献度を表す。この値cW1,j,lmnは、xlymznの項に対応する特徴量と考えることができる。
【0062】
さらに第１特徴量算出部１０４は、近似精度に対するp次の項の貢献度を、xlymznの項(l+m+n=p, l,m,nは0以上の整数)の近似精度に対する貢献度から求める（第２の演算手段）。着目次数についての貢献度は、着目次数を持つ１以上の単項式の貢献度から求めることができる。具体的には第１特徴量算出部１０４は、ベクトルcW1,j(1≦j≦NW1)の各要素を用いて、式（９）に従って特徴量ベクトルfW1,j(1≦j≦NW1)を算出する。本実施例においては式（９）のように、各項が近似精度に与える貢献度の平方の和を、各次数が近似精度に与える貢献度（寄与度）とする。【数9】
【0063】
特徴量ベクトルfW1,jは各次数についての貢献度を示す数値のセットとなる。この特徴量ベクトルfW1,j(1≦j≦NW1)の各要素は、ウィンドウ画像を近似する多項式の近似に対する各次数別の貢献度を意味する。具体的には、p番目の要素fW1,j,p（値fp）は、p次の項の貢献度を表す。この特徴量ベクトルfW1,jウィンドウ画像の回転に対して不変となる。そのためこの特徴量ベクトルfW1,jは、所定の点を中心に回転しているかもしれない複数の画像を位置合わせするために用いることができる。この値fW1,j,pは、次数pに対応する特徴量と考えることができる。第１特徴量算出部１０４は、ステップＳ１０３０で算出された多項式のそれぞれについて以上の処理を実行することにより、特徴量ベクトルfW1,j(1≦j≦NW1)を算出する。そして第１特徴量算出部１０４は、算出された特徴量ベクトルを記録する（第１の出力手段）。
【0064】
本実施例によれば第１多項式近似部１０３は、それぞれのウィンドウ画像について多項式近似を行った。そして上述のように、多項式の係数aW1,j、行列QW1,j、及び行列RW1,jを記録した。しかしながら式（８）のように第１特徴量算出部１０４は、行列QW1,j及び強度値ベクトルbW1,jから、ベクトルcW1,j及び特徴量ベクトルfW1,jを算出できる。したがって第１多項式近似部１０３が行列QW1,jを算出すれば、特徴量ベクトルfW1,jを求めるためには十分である。行列QW1,jは単項式行列MW1,jから導かれるため、実際には第１多項式近似部１０３が多項式の係数aW1,jを求める、すなわち多項式近似を行うことは必須ではない。
【0065】
（ステップＳ１０５０）特徴ウィンドウ画像の選択ステップＳ１０５０において特徴ウィンドウ選択部１０５は、ステップＳ１０４０で第１のウィンドウ画像のそれぞれについて算出された特徴量ベクトルに基づき、特徴ウィンドウ画像を選択する。この特徴ウィンドウ画像は、ステップＳ１１００（後述）でウィンドウ間の対応を検出する際に用いられる。
【0066】
ステップＳ１０５０の処理は以下の手順で行う。まず、それぞれのウィンドウ画像（第１の着目部分画像）について、そのウィンドウ画像と他のウィンドウ画像（第１の着目部分画像以外の部分画像）との特徴量ベクトルの差分（差異）の最小値を算出する。具体的には、j番目のウィンドウ画像について、fW1,jと、fW1,1, …,fW1,j-1, fW1,j+1, …, fW1,NW1のそれぞれとの間の差分を算出する。そして、算出した差分のうちの最小値を選択すればよい。
【0067】
特徴量ベクトルの差分の算出方法としては、任意の方法を用いることができる。たとえば、特徴量ベクトル間の差分ベクトルのL1ノルムを、差分として算出してもよい。この時、特徴量ベクトルそのものを用いてL1ノルムを求めてもよいし、予め特徴量ベクトルに正規化などの前処理を施してからL1ノルムを算出してもよい。
【0068】
そして、それぞれのウィンドウ画像について求めた差分の最小値が所定の閾値以上である場合には、そのウィンドウ画像を特徴ウィンドウ画像として選択する。特徴ウィンドウ選択部１０５は、選択された特徴ウィンドウ画像の特徴量ベクトルと、特徴ウィンドウ画像の中心位置とを記録する。以下で、k番目の特徴ウィンドウ画像についての特徴量ベクトルをfW1S,k、ウィンドウの中心位置をXWC1S,k(1≦k≦NW1S)と表記する。ただし、NW1Sは選択されたウィンドウ画像の総数である。
【0069】
ステップＳ１０５０の処理において、特徴ウィンドウ選択部１０５は、第１の画像のウィンドウ画像群のうち、他のウィンドウ画像との特徴量の差が閾値以上である特徴ウィンドウ画像を選択することができる。選択された特徴ウィンドウに対して、対応するウィンドウ画像が第２のウィンドウ画像群から後述のように選択される。このような特徴ウィンドウ画像は、第１の画像の中でも他と区別可能な画像特徴を持つことが予想される。したがって、このような特徴ウィンドウ画像に対応するウィンドウ画像を第２のウィンドウ画像群から後述のように検出することは、より高精度な位置合わせにつながることが期待される。
【0070】
本実施例では以上のように特徴量ベクトルの差分および特徴ウィンドウ画像を決定したが、他のウィンドウ画像との際が大きいウィンドウ画像を選択するための任意の方法を、ステップＳ１０５０では用いることができる。例えば特徴ウィンドウ選択部１０５は、全ての他のウィンドウ画像との特徴量の差の最小値がより大きいウィンドウ画像を、所定数選択してもよい。
【0071】
（ステップＳ１０６０）第２の画像の取得ステップＳ１０６０において第２画像取得部１１１は、第２の画像を取得する。第２の画像は、第１の画像と同様の画像である。第２の画像は、画像撮影装置１１の基準座標系における複数の３次元位置座標ベクトルX2,iと、それぞれの座標ベクトルに対応する強度値b2,iとの組の集合である。ここで、３次元位置座標ベクトルX2,i=[x2_i y2_i z2_i]T (1≦i≦N2)とする。また、1≦i≦N2であり、N2は第１の画像を構成する画素の数である。第２の画像は、ある所定の時刻前から現在までに撮影された、３次元の位置姿勢を有する複数の２次元超音波画像の集合によって構成されていてもよい。
【0072】
（ステップＳ１０７０）第２ウィンドウ画像群の生成ステップＳ１０７０において、第２ウィンドウ画像取得部１１２は、第２画像取得部１１１から取得した第２の画像から複数の部分領域（ウィンドウ、第２の部分領域群）を切り出す（抽出する）。こうして第２ウィンドウ画像取得部１１２は、第２のウィンドウ画像群（第２の部分画像群）を生成する。第２ウィンドウ画像取得部１１２の動作は、第１ウィンドウ画像取得部１０２と同様である。ここで、第２ウィンドウ画像のそれぞれの中心位置をXWC2,j=[xWC2,j yWC2,j zWC2,j]T (1≦j≦NW2)とする。ここで、NW2は第２のウィンドウ画像の総数である。
【0073】
また、j番目のウィンドウ内におけるk番目の画素の位置を表す３次元ベクトルを、XW2,j,k=[xW2,j,k yW2,j,k zW2,j,k]T (1≦j≦NW2,1≦k≦NPN2,j)とする。また、３次元ベクトルXW2,j,kに対応する画素の強度値を、bW2,j,k (1≦j≦NW2, 1≦k≦NPN2,j)とする。ここで、NPN2,jはj番目のウィンドウ画像の画素数である。
【0074】
第２ウィンドウ画像取得部１１２が取得するウィンドウ画像の数は、第１ウィンドウ画像取得部１０２が取得するウィンドウ画像の数と一致していてもよいし、異なっていてもよい。第２ウィンドウ画像取得部１１２が設定するウィンドウ領域群の位置および大きさは、第１ウィンドウ画像取得部１０２が設定するウィンドウ領域群の位置および大きさと一致していてもよいし、異なっていてもよい。
【0075】
（ステップＳ１０８０）第２ウィンドウ画像群の多項式近似ステップＳ１０８０において第２多項式近似部１１３は、ステップＳ１０７０で生成された第２のウィンドウ画像群に含まれる画像のそれぞれについて、多項式を算出する。ステップＳ１０８０における第２多項式近似部１１３の処理は、ステップＳ１０３０における第１多項式近似部１０３の処理と同様である。ステップＳ１０８０において第２多項式近似部１１３は、j番目のウィンドウ画像について、係数ベクトルaW2,j、行列QW2,j、及び行列RW2,j (1≦j≦NW2)を、ステップＳ１０７０と同様にして算出及び記録する。
【0076】
（ステップＳ１０９０）第２ウィンドウ画像群の特徴量の算出ステップＳ１０９０において第２特徴量算出部１１４は、第２のウィンドウ画像群に含まれる画像のそれぞれについて、ステップＳ１０８０で算出された多項式から特徴量を算出する。ステップＳ１０９０における特徴量の算出は、ステップＳ１０５０と同様に行うことができる。こうして第２特徴量算出部１１４は、j番目のウィンドウ画像についての特徴量ベクトルfW2,j(1≦j≦NW2)を算出して記録する。
【0077】
（ステップＳ１１００）ウィンドウ間対応の算出ステップＳ１１００において対応算出部１１５は、ステップＳ１０５０で第１のウィンドウ画像群から選択された特徴ウィンドウ画像と、第２のウィンドウ画像群に含まれるウィンドウ画像とを対応づける。具体的には対応算出部１１５は、特徴ウィンドウ画像のそれぞれについて、類似するウィンドウ画像を第２のウィンドウ画像群から選択する。例えば対応算出部１１５は、特徴ウィンドウ画像のそれぞれについて、特徴量ベクトルfW1S,k(1≦k≦NW1S)に類似する特徴量ベクトルfW2,j(1≦j≦NW2)を有する画像を第２のウィンドウ画像群から探索すればよい。
【0078】
具体的には対応算出部１１５は、特徴量ベクトルfW1S,k(1≦k≦NW1S)のそれぞれについて、所定の基準において最も類似する特徴量ベクトルを、特徴量ベクトルfW2,j(1≦j≦NW2)から探索する。例えば、２つの特徴量ベクトル間の差分ベクトルのL1ノルムがより小さい場合に、２つの特徴量ベクトルがより類似しているものとすることができる。この差分ベクトルは、特徴量ベクトルそのものの差分であってもよいし、正規化などの前処理を施された特徴量ベクトルの差分であってもよい。
【0079】
ステップＳ１１００の処理により、それぞれの特徴ウィンドウ画像に対応する画像が、第２のウィンドウ画像群から選択される。ここで、k番目の特徴ウィンドウ画像に対応する、第２のウィンドウ画像群から選択された画像の特徴量ベクトルをfW2S,k(1≦k≦NW1S)とする。また、選択された画像の中心位置をXWC2S,k(1≦k≦NW1S)とする。対応算出部１１５は、これらの特徴ウィンドウ画像に対応する画像の特徴量ベクトル及び中心位置を記録する。
【0080】
（ステップＳ１１１０）画像間の位置姿勢の算出ステップＳ１１１０において対応算出部１１５は、ステップＳ１１００で求めた対応付けに基づき、第１の画像と第２の画像との間の位置姿勢の関係を算出する。すなわち、それぞれの特徴ウィンドウ画像の位置と、第２のウィンドウ画像から選択された対応するウィンドウ画像の位置と、の関係に従って対応算出部１１５は位置姿勢関係を算出する。
【0081】
具体的には、それぞれの特徴ウィンドウ画像の中心位置XWC1S,k(1≦k≦NW1S)と、対応するウィンドウ画像の中心位置XWC2S,k(1≦k≦NW1S)と、の関係に従って対応算出部１１５は位置姿勢関係を算出できる。一例として、対応算出部１１５は、XWC1S,k(1≦k≦NW1S)とXWC2S,k(1≦k≦NW1S)との間の関係を最も良く近似する剛体変換を求めてもよい。この剛体変換を用いることにより、特徴ウィンドウ画像の中心位置と、対応するウィンドウ画像の中心位置と、が近くなるように第１の画像と第２の画像とを位置合わせすることができる。この剛体変換処理は、特異値分解を用いるなどの公知の技術によって行うことができる。
【0082】
位置姿勢関係を算出するために用いられるのは、ウィンドウ画像の中心位置には限られず、ウィンドウ画像の位置を示す任意の情報を用いることができる。こうして求められた位置姿勢関係を用いて、特徴ウィンドウ画像が切り出されたウィンドウ領域と、対応するウィンドウ画像が切り出されたウィンドウ領域と、が近くなるように第１の画像と第２の画像とを位置合わせすることができる。
【0083】
位置姿勢関係を算出する場合、対応算出部１１５は、それぞれの特徴ウィンドウ画像と対応する画像との組を平等に扱ってもよい。しかしながら位置姿勢関係を算出する際に、対応算出部１１５は、特徴ウィンドウ画像と対応する画像との組に重み付けを行ってもよい。例えば対応算出部１１５は、ステップＳ１１００で算出した特徴量ベクトル間の差分ベクトルのL1ノルムに従う重み付けを、それぞれの特徴ウィンドウ画像と対応する画像との組に適用してもよい。具体的には対応算出部１１５は、特徴量ベクトル間の差分がより小さい場合に、特徴ウィンドウ画像と対応する画像との組に大きい重みを与えてもよい。この場合、より類似度が高いウィンドウ画像の組を重視して剛体変換を求めることができる。また対応算出部１１５は、RANSACなどの公知の外れ値除去方法を用いてもよい。この場合、より安定した精度で剛体変換を求めることができる。
【0084】
また、ステップＳ１１１０で対応算出部１１５は、剛体変換の代わりに、ステップＳ１１００で求めた対応付けに基づく第１の画像と第２の画像との間の非線形な位置姿勢の関係を求めてもよい。例えば対応算出部１１５は、対応するウィンドウの中心位置XWC1S,k(1≦k≦NW1S)及びXWC2S,k(1≦k≦NW1S)を一致または略一致させる、非線形の変換関数を算出してもよい。この場合、非線形の変換関数としては放射基底関数（RBF）やB-Splineなどの公知のものを用いることができる。
【0085】
（ステップＳ１１２０）画像変換ステップＳ１１２０において画像生成部１１６は、ステップＳ１１１０で算出した第１の画像と第２の画像の位置姿勢の関係に基づき、第２の画像を座標変換することにより第３の画像を生成する。具体的には、第２の画像を第１の画像に位置合わせすることにより第３の画像が生成される。
【0086】
（ステップＳ１１３０）画像表示ステップＳ１１３０において画像生成部１１６は、第１の画像と第３の画像とを出力する。言い換えれば画像生成部１１６は、第１の画像と第２の画像とを位置合わせして出力することができる。出力された画像は、表示装置（不図示）を介してユーザに対して提示されてもよい。また出力された画像は、格納媒体（不図示）格納されてもよい。さらに出力された画像は、さらなる画像処理のために情報処理装置（不図示）に入力されてもよい。
【0087】
このとき、画像生成部１１６は第１の画像と第３の画像とをそのまま出力してもよい。また画像生成部１１６は、第１の画像及び第３の画像に対して処理を行うことにより得られる画像を出力してもよい。例えば画像生成部１１６は、第１の画像と第３の画像とのそれぞれについて、同じ切断面についての２次元の断面画像を生成してもよい。出力画像を生成する方法はこれに限られない。ユーザは、出力された画像により、第１の画像と第３の画像とを比較して観察することができる。
【0088】
（ステップＳ１１４０）終了判定ステップＳ１１４０において実時間処理部１１０は、全体の処理を終了するか否かの判定を行う。例えば処理されるべき第２の画像が残っていない場合には、実時間処理部１１０は全体の処理を終了すると判定する。終了すると判定した場合には、画像処理装置１０の処理の全体が終了する。一方、終了すると判定しなかった場合には、処理はステップＳ１０６０へと戻る。この場合、新たに取得される第２の画像を用いて、ステップＳ１０６０からステップＳ１１４０までの処理が行われる。
【0089】
以上に説明した本実施例における画像処理システムによれば、ウィンドウ画像から回転に不変な特徴量を抽出することができる。この特徴量を用いて、第１の画像及び第２の画像から抽出したウィンドウ画像を対応付けることができる。さらには、この対応付けを用いて第１の画像と第２の画像とを位置合わせすることができる。このため、第１の画像が第２の画像を回転した画像である場合でも、回転を考慮する（回転を探索的に推定する）必要なく、効率的に位置合わせを行うことができる。
【0090】
（変形例１−１）初期値がある場合に探索範囲を限定実施例１のステップＳ１１００で対応算出部１１５は、第２のウィンドウ画像群に含まれる全てのウィンドウ画像から、特徴ウィンドウ画像に類似するウィンドウ画像を検出した。しかしながら、実施例はこれに限られない。例えば第１の画像と第２の画像とが所定の誤差範囲内に位置合わせされている場合、以下のように処理を行ってもよい。すなわち、着目特徴ウィンドウ画像に対応するウィンドウ画像を検出する場合に、着目特徴ウィンドウ画像の近傍に位置するウィンドウ画像のみから対応するウィンドウ画像を検出してもよい。具体的には、着目特徴ウィンドウ画像の中心位置に対応する第２の画像の位置から所定範囲内にあるウィンドウ画像から、対応するウィンドウ画像を検出すればよい。この場合、所定範囲内にあるウィンドウ画像を対象として、着目特徴ウィンドウ画像との間で特徴量の比較を行えばよい。また同様にステップＳ１０５０において特徴ウィンドウ選択部１０５は、所定範囲内のウィンドウだけを対象として特徴量ベクトル間の比較を行うことにより、特徴ウィンドウ画像を選択してもよい。
【0091】
この方法によれば、第１の画像と第２の画像とが所定の誤差範囲内に位置合わせされている場合に、処理負荷を少なくすることができる。したがって、処理を高速化及び効率化することができる。また、互いに関連しないウィンドウ画像を対応付けする可能性が小さくなるため、より正確な位置合わせ結果が得られることが期待される。
【0092】
（変形例１−２）特徴点検出を用いてウィンドウ領域を設定実施例１において第１ウィンドウ画像取得部１０２は、第１の画像の領域９００の全体に対して、ウィンドウ領域を等間隔に設定した。しかしながら、実施例はこれに限られない。例えばステップＳ１０２０において第１ウィンドウ画像取得部１０２は、第１の画像から複数の特徴点を抽出し、それぞれの特徴点を中心とするウィンドウ領域を設定してもよい。特徴点とは例えば、強度値の変化が大きい部分、所定のパターンを有する部分、などでありうる。特徴点の抽出方法としてはHarris Corner Detectorなどの公知の方法を用いることができる。またステップＳ１０７０においても同様に、第２ウィンドウ画像取得部１１２は、第２の画像から複数の特徴点を抽出し、それぞれの特徴点を中心とするウィンドウ領域を設定してもよい。この方法によれば、画像中の特徴的な領域を含むウィンドウ画像についてのみ処理を行えばよいため、処理を高速化及び効率化できる。
【0093】
（変形例１−３）画像の多様性実施例１では３次元の超音波画像を入力画像（第１の画像及び第２の画像）として用いた。しかしながら、実施例はこれに限られない。例えば、画像は２次元又は４次元の画像であってもよい。ここで４次元の画像とは、３次元画像の時系列の集合を指す。また、画像の種類はＭＲＩ画像、ＣＴ画像、などの他の医用画像であってもよい。さらには、カメラ画像や距離画像などの、医用画像以外の画像に対して、上述の方法を適用してもよい。
【0094】
入力画像として２次元の画像を用いる場合には、ステップＳ１０２０およびステップＳ１０７０で生成するウィンドウ画像を、例えば円形のような２次元の形状とすればよい。またステップＳ１０３０及びステップＳ１０８０では、xとyとを変数とする多項式でウィンドウ画像を近似すればよい。さらに、第１の画像と第２の画像とが同一種類の画像である必要はない。例えば、第１の画像がＭＲＩ画像であり、第２の画像がＣＴ画像であってもよい。このように、入力される画像の次元及び画像の撮影方式は限定されない。
【0095】
［実施例２］ウィンドウ領域の固定実施例１では、第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについて単項式行列が生成された。また、第２のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについても単項式行列が生成された。しかし、実施例はこれに限られない。実施例２では、第１の画像に対して、複数の同じ形状のウィンドウ領域が設定される。また、第２の画像に対しても、複数の同じ形状のウィンドウ領域が設定される。第１の画像に設定されるウィンドウ領域と第２の画像に設定されるウィンドウ領域とは、同じ形状であってもよい。この場合、それぞれのウィンドウ画像について、ウィンドウ画像を構成する画素の位置（ウィンドウ中心を基準とした画素の位置）は一致する。実施例２によれば、共通の単項式行列をそれぞれのウィンドウ画像に対して用いることができるため、処理を効率化できる。
【0096】
図５は、本実施例に係る画像処理システムの構成を示す。実施例１の画像処理システムに類似する要素には同一の符号が付されており、説明は省略される。第１ウィンドウ画像取得部２０２は、第１画像取得部１０１から取得した第１の画像に対して複数のウィンドウ領域を設定する。ここで、設定される複数のウィンドウ領域は同じ形状及び同じ画素数を有する。設定されたウィンドウ領域に従って、第１ウィンドウ画像取得部２０２は、実施例１と同様に、第１のウィンドウ画像群を生成する。
【0097】
単項式行列取得部２０３は、第１ウィンドウ画像取得部２０２が生成した第１のウィンドウ画像群について、単項式行列に関する情報を取得する。単項式行列取得部２０３が行う具体的な処理については後に詳しく説明する。
【0098】
第１特徴量算出部２０４は、第１ウィンドウ画像取得部２０２が生成した第１のウィンドウ画像と、単項式行列取得部２０３が取得した単項式行列に関する情報とに基づいて、第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれの特徴量を算出する。
【0099】
第２ウィンドウ画像取得部２１２は、第２画像取得部１１１から取得した第２の画像に対して複数のウィンドウ領域を設定する。ここで、設定される複数のウィンドウ領域は同じ形状及び同じ画素数を有する。また第２ウィンドウ画像取得部２１２が設定するウィンドウ領域は、第１ウィンドウ画像取得部２０２が設定するウィンドウ領域と同じ形状及び同じ画素数を有する。設定されたウィンドウ領域に従って、第２ウィンドウ画像取得部２１２は、実施例１と同様に、第２のウィンドウ画像群を生成する。
【0100】
第２特徴量算出部２１４は、第１ウィンドウ画像取得部２１２が生成した第２のウィンドウ画像と、単項式行列取得部２１３が取得した単項式行列に関する情報とに基づいて、第２の特徴量算出部２１４は第２のウィンドウ画像群に含まれるウィンドウ画像のそれぞれの特徴量を算出する。
【0101】
次に、本実施例に係る画像処理装置２０が行う動作に関して、図６のフローチャートを参照して詳しく説明する。ステップＳ２０１０の処理は、実施例１におけるステップＳ１０１０と同様の処理であり、説明は省略する。
【0102】
（ステップＳ２０２０）共通サイズの第１のウィンドウ画像群を生成ステップＳ２０２０において、第１ウィンドウ画像取得部２０２は、第１画像取得部２０１から取得した第１の画像から複数の部分領域（ウィンドウ）を切り出すことにより、第１のウィンドウ画像群を生成する。この処理は、実施例１のステップＳ１０２０と同様であるが、それぞれのウィンドウ画像の大きさ及び形状を同一とする。ウィンドウ領域の大きさは、画像の解像度などに応じて適宜設定することができる。
【0103】
第１の画像を構成する画素の配置が、第１の画像の領域において空間的に均一でないことがある。この場合、内挿処理などの補間方法を用いて、画素の配置を均一にする。ステップＳ２０２０で生成されるウィンドウ画像のそれぞれにおいては、ウィンドウ画像を構成するそれぞれの画素は、基準点からみて同じ位置にある。したがって、それぞれのウィンドウ画像の画素数はNPNと表すことができる。また、それぞれのウィンドウ画像における、基準点（例えばウィンドウ中心位置）を原点とするk番目の画素の位置ベクトルは、XW,k(1≦k≦NPN)と表すことができる。
【0104】
以上のように設定したウィンドウ領域に従って、第１ウィンドウ画像取得部２０２は複数のウィンドウ画像を生成する。具体的には、第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについて、ウィンドウ中心位置と強度値とを記録する。j番目のウィンドウ画像の中心位置はXWC1,j(1≦j≦NW1)と表される。また、j番目のウィンドウ画像におけるk番目の画素の強度値はbW1,j,k(1≦j≦NW1, 1≦k≦NPN)と表される。上述のように、それぞれの画素の位置ベクトルXW,k(1≦k≦NPN)もまた記録される。
【0105】
（ステップＳ２０３０）共通の単項式行列からQを算出するステップＳ２０３０において単項式行列取得部２０３は、ステップＳ２０２０で生成した位置ベクトルXW1,k(1≦k≦NPN)を用いて、単項式行列を生成する。実施例１で説明したように、j番目のウィンドウ画像についての単項式行列MW1,jは、ベクトルm(xW1,j,k, yW1,j,k, zW1,j,k)(1≦k≦NPN1)によって構成される。またベクトルm(xW1,j,k, yW1,j,k, zW1,j,k)は、式（１）に示されるように、j番目のウィンドウ画像のそれぞれの画素についての位置ベクトルXW1,j,k=[xW1,j,k, yW1,j,k, zW1,j,k]で表される。本実施例においてそれぞれのウィンドウ画像のk番目の画素の位置ベクトルは、共通のベクトルXW,kを用いて表される。したがって、それぞれのウィンドウ画像についての単項式行列もまた、共通の単項式行列MWを用いて表される。
【0106】
このように単項式行列取得部２０３は、それぞれのウィンドウ画像に共通な単項式行列MWを算出する。さらに単項式行列取得部２０３は、実施例１と同様に、単項式行列MWをQR分解することにより行列QW及びRWを生成する。単項式行列取得部２０３は、以上のように生成したMW、QW、及びRWを、単項式行列に関する情報として記録する。
【0107】
（ステップＳ２０４０）第１のウィンドウ画像群の特徴量を算出ステップＳ２０４０において第１特徴量算出部２０４は、第１のウィンドウ画像群に含まれるそれぞれのウィンドウ画像について、特徴量を算出する。具体的には、各ウィンドウ画像を構成する画素の強度値と、ステップＳ２０３０で生成された単項式行列に関する情報とに基づいて、ウィンドウ画像の特徴量を算出する。
【0108】
例えば第１特徴量算出部２０４は、式（８）に従ってベクトルcW,jを計算できる。具体的には第１特徴量算出部２０４は、ステップＳ２０２０で記録されたウィンドウ画像の各画素の強度値bW1,j,k(1≦j≦NW1, 1≦k≦NPN)と、ステップＳ２０３０で生成されたQWとを用いて式（８）の計算を行う。その後第１特徴量算出部２０４は、実施例１と同様に、ベクトルcW1,jを用いてウィンドウ画像の特徴量を算出できる。以上のように第１特徴量算出部２０４は、ステップＳ２０２０で生成した全てのウィンドウ画像について、特徴量ベクトルfW1,j(1≦j≦NW1)を算出して記録する。
【0109】
ステップＳ２０５０からステップＳ２０６０の処理は、実施例１のステップＳ１０５０からステップＳ１０６０の処理と同様であるから、説明は省略する。
【0110】
（ステップＳ２０７０）共通サイズの第２のウィンドウ画像群を生成ステップＳ２０７０において、第２ウィンドウ画像取得部２１２は、第２画像取得部１１１から取得した第２の画像から複数の部分領域（ウィンドウ）を切り出すことにより、第２のウィンドウ画像群を生成する。この処理はステップＳ２０２０と同様である。ステップＳ２０７０で設定されるウィンドウ領域の大きさは、ステップＳ２０７０で設定されるウィンドウ領域の大きさに一致する。すなわち、第２のウィンドウ画像群に含まれるウィンドウ画像のそれぞれにおいて、基準点（例えばウィンドウ中心位置）を原点とするk番目の画素の位置ベクトルは、XW,k(1≦k≦NPN)と表すことができる。第２のウィンドウ画像群に含まれるウィンドウ画像の画素数NPNは、第１のウィンドウ画像群に含まれるウィンドウ画像の画素数に一致する。
【0111】
もっとも、第１の画像から生成されるウィンドウ画像の数と、第２の画像から生成されるウィンドウ画像の数とは、異なっていてもよい。ステップＳ２０２０と同様に第２ウィンドウ画像取得部２０２は、第２の画像を構成する画素X2,i(1≦i≦N2)の配置が空間的に均一でない場合、内挿処理などの補間方法を用いて、画素の配置を均一にする。
【0112】
そして第２ウィンドウ画像取得部２１２は、第２のウィンドウ画像群に含まれる第２のウィンドウ画像のそれぞれについて、ウィンドウ中心位置XWC2,j(1≦j≦NW2)、及び画素の強度値bW2,j,k(1≦j≦NW2,1≦k≦NPN2)を記録する。各ウィンドウ画像における画素の位置は、ステップＳ２０２０で記録されたXW,k(1≦k≦NPN)と同一であるから、画素位置の記録は省略できる。
【0113】
（ステップＳ２０８０）第２のウィンドウ画像群の特徴量を算出ステップＳ２０８０において第２特徴量算出部２１４は、第２のウィンドウ画像群に含まれるそれぞれのウィンドウ画像について、特徴量を算出する。具体的には、各ウィンドウ画像を構成する画素の強度値と、ステップＳ２０３０で生成された単項式行列に関する情報とに基づいて、ウィンドウ画像の特徴量を算出する。
【0114】
例えば第２特徴量算出部２１４は、式（８）に従ってベクトルcW,jを計算できる。具体的には第２特徴量算出部２１４は、ステップＳ２０７０で記録されたウィンドウ画像の各画素の強度値bW2,j,k(1≦j≦NW2, 1≦k≦NPN)と、ステップＳ２０３０で生成されたQWとを用いて式（８）の計算を行う。その後第２特徴量算出部２１４は、ステップＳ２０４０と同様に、ベクトルcW1,jを用いてウィンドウ画像の特徴量を算出できる。以上のように第２特徴量算出部２１４は、ステップＳ２０７０で生成した全てのウィンドウ画像について、特徴量ベクトルfW2,j(1≦j≦NW2)を算出して記録する。
【0115】
ステップＳ２０９０からステップＳ２１３０の処理は、実施例１におけるステップＳ１１００からステップＳ１１４０の処理と同様であり、説明は省略する。以上のように、本実施例における画像処理装置２０は処理を行う。
【0116】
以上に説明した、本実施例における画像処理システムにおいては、ウィンドウ画像の特徴量の算出において、共通の単項式行列が用いられる。本実施例によれば、単項式行列に関する情報を繰り返し生成する必要がないために、処理負荷を軽減することができる。
【0117】
（変形例２−１）実施例２においては、共通の大きさを有するウィンドウ画像が、第１の画像及び第２の画像から生成された。しかしながら、第１の画像から得られるウィンドウ画像と第２の画像から得られるウィンドウ画像とは、異なる大きさを有してもよい。また、第１の画像（及び第２の画像）からは、複数の形状のうち何れかを有するウィンドウ画像が生成されてもよい。例えば、画像の中心部から得られるウィンドウ画像と、画像の周辺部から得られるウィンドウ画像とが、異なる形状を有してもよい。
【0118】
また、ウィンドウ画像の形状は予め設定されていてもよい。単項式行列に関する情報M,Q,Rは、ウィンドウ画像の形状とウィンドウ画像における原点と、に応じて決定される。したがって、ウィンドウ画像の形状が予め設定されている場合、単項式行列に関する情報M,Q,Rもまた予め決定されていてもよい。この場合、実施例２におけるステップＳ２０３０の処理を省略することができる。例えば、画像処理装置２０が用いるウィンドウの画素位置が常に固定である場合、その画素位置によって定まる単項式行列に関する情報を磁気ディスク８１３に予め記録しておいてもよい。この場合ステップＳ２０４０及びステップＳ２０８０において、磁気ディスク８１３に記録されている単項式行列に関する情報が読み込まれる。
【0119】
［実施例３］（ウィンドウ画像の独立性に基づく対応付け）実施例１では、第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについて算出した特徴量を指標として、所定の閾値以上に類似する他のウィンドウ画像が存在しないウィンドウ画像が、特徴ウィンドウ画像として選択された。そして、特徴ウィンドウ画像のそれぞれに対して、第２のウィンドウ画像群に含まれる類似するウィンドウ画像が対応付けられた。しかしながら、実施例はこれに限られない。実施例３では、第１のウィンドウ画像群に含まれるそれぞれのウィンドウ画像の独立性が、算出された特徴量に基づいて算出される。そして、この独立性に基づいて、第１のウィンドウ画像群と第２のウィンドウ画像群とが対応付けられる。
【0121】
独立度算出部３０５は、第１特徴量算出部１０４が算出した第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについて、特徴量に基づいて独立度を算出する。対応算出部３１５は、第１特徴量算出部１０４が算出した特徴量と、第２特徴量算出部１１４が算出した特徴量と、独立度算出部３０５が算出した独立度とに基づいて、ウィンドウ画像間の対応を算出する。
【0122】
次に、本実施例に係る画像処理装置３０が行う動作に関して、図８のフローチャートを参照して詳しく説明する。ステップＳ３０１０からステップＳ３０４０の処理は、実施例１のステップＳ１０１０からステップＳ１０４０と同様であり、説明は省略する。
【0123】
（ステップＳ３０５０）独立度の算出ステップＳ３０５０において独立度算出部３０５は、第１のウィンドウ画像群に含まれるウィンドウ画像の独立度gj(1≦j≦NW1)を算出する。ここでgjは、j番目のウィンドウ画像の独立度を示す。具体的にはこの独立度は、ウィンドウ画像の特徴量fW1,j(1≦j≦NW1)に基づいて算出される。例えば独立度算出部３０５は、以下の式（１０）に従って独立度を算出できる。
【数10】
【0124】
ステップＳ３０６０からステップＳ３０９０の処理は、実施例１のステップＳ１０６０からステップＳ１０９０の処理と同様であり、説明を省略する。
【0125】
（ステップＳ３１００）独立度を考慮したウィンドウ間の対応付けステップＳ３１００において対応算出部３１５は、第１のウィンドウ画像群に関する特徴量及び独立度、並びに第２のウィンドウ画像群に関する特徴量に基づき、第１のウィンドウ画像群と第２のウィンドウ画像群とを対応付ける。具体的には以下の処理を実行する。まず対応算出部３１５は、特徴量ベクトルfW1,k(1≦k≦NW1)のそれぞれについて、所定の基準において最も類似する特徴量ベクトルを、特徴量ベクトルfW2,j(1≦j≦NW2)から探索する。例えば、２つの特徴量ベクトル間の差分ベクトルのL1ノルムがより小さい場合に、２つの特徴量ベクトルがより類似しているものとすることができる。この処理は、ステップＳ１１００と同様に行うことができる。
【0126】
さらに対応算出部３１５は、類似している特徴量の差分が、独立度と比べて大きいか否かを判定する。特徴量の差分が独立度と比べて同等か又は小さい場合には、対応算出部３１５は２つの特徴量ベクトルを対応付けない。一方で特徴量の差分が独立度と比べて大きい場合には、対応算出部３１５は２つの特徴量ベクトルを対応付ける。
【0127】
第１のウィンドウ画像群に含まれるj番目のウィンドウ画像の特徴量ベクトルが、第２のウィンドウ画像群に含まれるk番目のウィンドウ画像の特徴量ベクトルに類似している場合について、処理の一例を示す。このとき対応算出部３１５は、特徴量ベクトルfW1,jとfW2,kとの差分ベクトルのL1ノルムを算出する。L1ノルムが、独立度gj(1≦j≦NW1)以下である場合には、対応算出部３１５は第１のウィンドウ画像群に含まれるj番目のウィンドウ画像と第２のウィンドウ画像群に含まれるk番目のウィンドウ画像とを対応付ける。一方で、L1ノルムが、独立度gj(1≦j≦NW1)を超える場合には、対応算出部３１５は対応付けを行わない。ここで動作を最適化するために、独立度gjに対して所定の係数を乗じてもよい。
【0128】
以上のようにして求められた対応付けを、対応算出部３１５は記録する。具体的には対応算出部３１５は、特徴量ベクトルfW1S,k(1≦k≦NW1S)と、対応付けられた特徴量ベクトルfW2S,k(1≦k≦NW1S)と、の組を記録する。また対応算出部３１５は、ウィンドウの中心位置XWC1S,k(1≦k≦NW1S)と、対応するウィンドウの中心位置XWC2S,k(1≦k≦NW1S)と、の組を記録する。
【0129】
L1ノルムが独立度gj(1≦j≦NW1)を超える場合に対応算出部３１５は、独立度が低すぎることを示す情報とともにウィンドウ画像間の対応付けを記録してもよい。この場合、続くステップＳ３１１０において画像生成部１１６は、独立度が低すぎることを示す情報が付された対応付けを無視することができる。
【0130】
ステップＳ３１１０からステップＳ３１４０の処理は、実施例１のステップＳ１１１０からステップＳ１１４０の処理と同一であり、説明を省略する。
【0131】
本実施例においては第１のウィンドウ画像について独立度を算出した。しかしながら、第２のウィンドウ画像について独立度を算出し、対応付けを行うか否かの判断にこの独立度を用いてもよい。また、第１のウィンドウ画像と第２のウィンドウ画像との双方について独立度を算出し、双方の独立度を用いて対応付けを行うか否かの判断を行ってもよい。例えば、第１のウィンドウ画像の独立度の代わりに、第１のウィンドウ画像の独立度と第２のウィンドウ画像の独立度との和を用いてもよい。
【0132】
以上に説明した本実施例における画像処理システムによれば、第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについて独立度が算出される。そして、この独立度を参照して、対応付けを行うか否かを切り換えられる。このため、独立度の低いウィンドウ画像については、対応付けが行われる可能性が小さくなる。また、独立度の高いウィンドウ画像については、積極的に対応付けが行われる。この結果、画像内の自己類似性を考慮して、適応的にウィンドウ画像を対応付けすることができる。この結果、互いに関連しないウィンドウ画像を対応付けする可能性が小さくなるため、より正確な位置合わせ結果が得られることが期待される。
【0133】
［実施例４］ウィンドウ画像の変形実施例１では、第１の画像および第２の画像の画像領域に複数のウィンドウ中心が配置された。そして、このウィンドウ中心のそれぞれを中心とする球形の部分領域（ウィンドウ）が設定された。ウィンドウ画像群は、この設定されたウィンドウ内の画像を切り出すことによって生成された。実施例１は、第１の画像における球形の領域が、第２の画像における球形の領域に対応する場合により効果的である。すなわち、第１の画像と第２の画像との間で、被写体があまり変形していない場合に、実施例１はより有効である。しかしながら、第１の画像と第２の画像との間で、被写体が変形することがある。実施例４においては、第１の画像からに設定するウィンドウ領域の形状と、第２の画像に設定するウィンドウ領域の形状とを、異なる形状とする。本実施例によれば、第１の画像と第２の画像との間で被写体が変形している場合にも、より高い精度で画像の対応付けを行うことができる。
【0135】
第１ウィンドウ画像取得部４０２は、第１画像取得部１０１から取得した第１の画像から複数の部分領域（ウィンドウ）を切り出すことにより、第１のウィンドウ画像群を生成する。ここでウィンドウ領域の形状は必ずしも球形ではなくてもよい。例えばウィンドウ領域の形状は、球形を様々に偏平させることによって得られる楕円体の形状であってもよい。
【0136】
次に、本実施例に係る画像処理装置４０が行う動作を、図１０のフローチャートを参照して詳しく説明する。ステップＳ４０１０の処理は実施例１のステップＳ１０１０と同様であり、説明は省略する。
【0137】
（ステップＳ４０２０）様々な楕円体のウィンドウ画像を生成ステップＳ４０２０において、第１ウィンドウ画像取得部４０２は、第１画像取得部１０１から取得した第１の画像から複数の部分領域（ウィンドウ）を切り出すことにより、第１のウィンドウ画像群を生成する。ステップＳ４０２０において第１ウィンドウ画像取得部４０２が第１の画像から切り出すウィンドウ画像の形状を、図１１Ａ〜図１１Ｉを参照して説明する。
【0138】
図１１Ａは、第１ウィンドウ画像取得部４０２が第１の画像に設定するウィンドウ領域の一例を示す。図１１Ａに示されるウィンドウ領域９１０は、ウィンドウ中心９１１を中心とした球形領域である。しかしながら本実施例において第１ウィンドウ画像取得部４０２が設定するウィンドウ領域は、図１１Ａに示されるものには限られない。例えば第１ウィンドウ画像取得部４０２は、図１１Ｂ〜図１１Ｉに示される形状のウィンドウ領域を設定してもよい。図１１Ｂ及び図１１Ｃに示される形状は、球形のウィンドウをある軸方向に偏平させた形状（楕円体）である。図１１Ｄ及び図１１Ｅに示される形状は、球形のウィンドウを図１１Ｂ及び図１１Ｃとは異なる軸方向に偏平させた楕円体である。また、球が偏平される方向は第１の画像の基準座標系における座標軸方向には限定されない。例えば図１１Ｆ〜図１１Ｉに示される形状は、球径のウィンドウを任意の方向に偏平させることによって得られる。またウィンドウ形状は必ずしも楕円体である必要はない。例えば、球形の形状を変形させて生成される任意の形状でありうる。
【0139】
本ステップにおいて第１ウィンドウ画像取得部４０２は、様々な形状のウィンドウ画像を、第１の画像から切り出すことにより、第１のウィンドウ画像群を生成する。例えば第１ウィンドウ画像取得部４０２は、第１の画像の全領域に対し、等間隔にNW1center個のウィンドウの中心位置を設定する。そして第１ウィンドウ画像取得部４０２は、この位置を中心とするNangle種類の楕円体領域を設定する。それぞれの種類の楕円体領域は、球を異なる方向に偏平させることによって得られる楕円体形状でありうる。設定されたウィンドウ領域に従ってウィンドウ画像を切り出すことにより、第１ウィンドウ画像取得部４０２は、全部でNW1center×Nangle個のウィンドウ画像を生成する。
【0140】
そして第１ウィンドウ画像取得部４０２は、第１のウィンドウ画像群に含まれるウィンドウ画像のそれぞれについて、ウィンドウ中心位置と、ウィンドウに含まれる画素の位置と強度値との組と、を記録する。ここで、j番目のウィンドウ画像についてのウィンドウ中心位置はXWC1,j(1≦j≦NW1)で表される。また、j番目のウィンドウ画像におけるk番目の画素の位置はXW1,j,k(1≦j≦NW1,1≦k≦NPN1,j)で表される。さらに、j番目のウィンドウ画像についてのk番目の画素の強度値はbW1,j,k(1≦j≦NW1,1≦k≦NPN1,j)で表される。ただし、NW1は第１のウィンドウ画像群に含まれるウィンドウ画像の数であり、本実施例の場合はNW1=NW1center×Nangleとなる。
【0141】
このとき第１ウィンドウ画像取得部４０２は、画素の位置XW1,j,kを補正してもよい。例えば第１ウィンドウ画像取得部４０２は、第１のウィンドウ画像群に含まれるウィンドウ画像を、第２のウィンドウ画像群に含まれるウィンドウ画像と同じ形状となるように補正することができる。本実施例において第２のウィンドウ画像群は球形状を有する。第１ウィンドウ画像取得部４０２は、画素の位置XW1,j,kを補正することにより、第１のウィンドウ画像群を球形状に補正できる。x軸方向の長さが半分となるように偏平されたウィンドウ領域から生成されたウィンドウ画像を補正する場合について考える。ここで、ウィンドウ画像における座標系の原点は、ウィンドウ画像の中心であるものとする。この場合ウィンドウ画像取得部４０２は、ウィンドウ画像に含まれる画素の位置について、x座標を２倍にすればよい。このようにウィンドウ画像の形状を補正することにより、被写体が第１の画像と第２の画像との間で変形している場合にも、より高い精度で画像の位置合わせを行うことができる。
【0142】
ステップＳ４０３０からステップＳ４１４０までの処理は、実施例１のステップＳ１０３０からステップＳ１１４０までの処理と同様であり、説明は省略する。
【0143】
以上の本実施例における画像処理システムによれば、第１のウィンドウ画像群として、球形を様々に偏平させることによって得られる形状のウィンドウ画像を生成する。そして、これらのウィンドウ画像の特徴量が算出される。本実施例によれば、第１の画像の被写体が第２の画像においては変形している場合にも、より高精度に位置合わせを行うことができる。
【0144】
（変形例４−１）被写体の変形量に応じてウィンドウ形状を変更実施例４では、ステップＳ４０２０において、ウィンドウ領域の中心位置のそれぞれについて、所定の数のウィンドウ領域が設定された。しかしながら、実施例はこれに限られない。例えば、被写体の変形特性が予め与えられている場合には、ウィンドウ領域の形状を変形特性に応じて決定することができる。また、被写体の変形特性が推定できる場合にも、ウィンドウ領域の形状を変形特性に応じて決定することができる。
【0145】
例えば、画像中で被写体が所定の方向に圧縮又は伸張されている場合、この方向に偏平されたウィンドウ領域が設定されてもよい。圧縮又は伸張の程度に依存して、ウィンドウ領域の偏平量が決定されてもよい。さらに、画像中の部位毎に変形量が推定又は取得できる場合には、この変形量に基づいて画像中の部位毎に異なる形状のウィンドウ領域を設定してもよい。さらには、被写体に著しい変形が生じていることが推定できる部分からは、ウィンドウ画像が生成されないようにすることもできる。
【0146】
以上の方法によれば、第１の画像と第２の画像との間での被写体の変形量に応じてウィンドウ画像の形状を決定することができる。このため、生成されたウィンドウ画像を用いると、より正確に位置合わせを行うことが可能となる。
【0147】
［その他の実施例］なお、上述した本実施の形態における記述は、本発明に係る好適な画像処理装置の一例であり、本発明はこれに限定されるものではない。例えば上述の各実施例において、各構成要素は複数の構成要素に分割されていてもよい。また、複数の構成要素の機能が１つの構成要素によって実現されてもよい。上述の実施例においては、第１の画像を処理するための構成要素と、第２の画像を処理するための構成要素とは異なっている。例えば実施例１においては、第１の画像からは構成要素１０１〜１０４によって特徴量が算出される。また、第２の画像からは構成要素１１１〜１１４によって特徴量が算出される。しかしながら、構成要素１０１〜１０４が、第２の画像についての特徴量を算出してもよい。
【0148】
本発明の態様は、上述の実施形態の機能を実行するためにメモリデバイス上に記録されたプログラムを読み出し及び実行するシステム又は装置のコンピュータ（又はＣＰＵ若しくはＭＰＵのようなデバイス）によって、並びに、方法であって、方法における工程が上述の実施形態の機能を実行するために例えばメモリデバイス上に記録されたプログラムを読み出し及び実行することによってシステム又は装置のコンピュータによって実行される、方法によって、また実現されることができる。この目的のためにプログラムは、例えばネットワークを介して、又はメモリデバイスとして働く様々な形式の記録媒体（例えばコンピュータ読み取り可能な媒体）から、コンピュータへと提供される。
【0149】
本発明が例示的な実施形態を参照して説明されたが、本発明は開示された例示的な実施形態に限定されないことは理解されるであろう。以下の請求の範囲には、このような修正並びに均等な構造及び機能を包含するように、最も広い解釈が与えられる。
【0150】
本願は、日本国特許出願第２０１１−２１３３８１号の利益を主張し、その全体が本明細書に参照として組み込まれる。