(57)【特許請求の範囲】
【請求項１】
フラクタル次元を測定して画像処理を行う画像処理装置において、
多値画像から一定領域の画像を抽出して抽出画像を生成する抽出画像生成手段と、
前記抽出画像に解像度を設定し、前記解像度を階層的に変えて、複数の第１の解像度変換画像及び複数の第２の解像度変換画像を生成する解像度変換画像生成手段と、
前記第１の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最大輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とする最大輝度値設定手段と、
前記第２の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最小輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とする最小輝度値設定手段と、
前記第１の解像度変換画像と前記第２の解像度変換画像の同一座標の画素の輝度値の差分の総和である輝度差分値を計算する輝度差分値計算手段と、
前記解像度及び前記輝度差分値の座標軸上に、前記輝度差分値をプロットして生成したグラフの傾きを求め、前記傾きを前記フラクタル次元として前記多値画像全体の各画素にわたって測定するフラクタル次元測定手段と、
前記フラクタル次元の画像を表示するフラクタル次元画像表示手段と、
を有することを特徴とする画像処理装置。
【請求項２】
フーリエ変換を行って、前記多値画像をフーリエ変換画像に変換する画像変換手段をさらに有することを特徴とする請求項１記載の画像処理装置。
【請求項３】
複数の連続画像の処理を行う場合は、複数の２次元多値画像をＹ軸の方向に切って、複数のＸ−Ｚ軸からなる２次元のフラクタル次元画像を生成するフラクタル次元画像生成手段をさらに有することを特徴とする請求項１記載の画像処理装置。
【請求項４】
前記最大輝度値設定手段は、下位の階層の前記第１の解像度変換画像の画素（ｉ，ｊ）に対し、ｉ／２及びｊ／２の数値を越えない最大の整数〔ｉ／２〕及び〔ｊ／２〕を、上位の階層の前記第１の解像度変換画像の画素（〔ｉ／２〕，〔ｊ／２〕）に対応させて、前記最大輝度値を設定することを特徴とする請求項１記載の画像処理装置。
【請求項５】
前記最小輝度値設定手段は、下位の階層の前記第２の解像度変換画像の画素（ｉ，ｊ）に対し、ｉ／２及びｊ／２の数値を越えない最大の整数〔ｉ／２〕及び〔ｊ／２〕を、上位の階層の前記第２の解像度変換画像の画素（〔ｉ／２〕，〔ｊ／２〕）に対応させて、前記最小輝度値を設定することを特徴とする請求項１記載の画像処理装置。
【請求項６】
前記フラクタル次元画像表示手段は、フラクタル次元画像または前記フラクタル次元の統計的分布を表示することを特徴とする請求項１記載の画像処理装置。
【請求項７】
複数の連続画像の処理を行う場合は、２枚目からの画像は前画像からの差分を記憶し、前記差分のある箇所のみ前記フラクタル次元を測定することを特徴とする請求項１記載の画像処理装置。
【請求項８】
フラクタル次元を測定してＸ線医療画像の処理を行うＣＴスキャナにおいて、
被検体の前記Ｘ線医療画像である多値画像を収集する多値画像収集手段と、
前記多値画像から一定領域の画像を抽出して抽出画像を生成する抽出画像生成手段と、前記抽出画像に解像度を設定し、前記解像度を階層的に変えて、複数の第１の解像度変換画像及び複数の第２の解像度変換画像を生成する解像度変換画像生成手段と、前記第１の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最大輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とする最大輝度値設定手段と、前記第２の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最小輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とする最小輝度値設定手段と、前記第１の解像度変換画像と前記第２の解像度変換画像の同一座標の画素の輝度値の差分の総和である輝度差分値を計算する輝度差分値計算手段と、前記解像度及び前記輝度差分値の座標軸上に、前記輝度差分値をプロットして生成したグラフの傾きを求め、前記傾きを前記フラクタル次元として前記多値画像全体の各画素にわたって測定するフラクタル次元測定手段と、前記フラクタル次元の画像を表示するフラクタル次元画像表示手段と、から構成される画像処理手段と、
を有することを特徴とするＣＴスキャナ。
【請求項９】
フラクタル次元を高速測定するフラクタル次元測定方法において、
多値画像から一定領域の画像を抽出して抽出画像を生成し、
前記抽出画像に解像度を設定し、前記解像度を階層的に変えて、複数の第１の解像度変換画像及び複数の第２の解像度変換画像を生成し、
前記第１の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最大輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とし、
前記第２の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最小輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とし、
前記第１の解像度変換画像と前記第２の解像度変換画像の同一座標の画素の輝度値の差分の総和である輝度差分値を計算し、
前記解像度及び前記輝度差分値の座標軸上に、前記輝度差分値をプロットして生成したグラフの傾きを求め、前記傾きを前記フラクタル次元として前記多値画像全体の各画素にわたって測定することを特徴とするフラクタル次元測定方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は画像処理装置、ＣＴ（Computerized Tomography)スキャナ及びフラクタル次元計算方法に関し、特にフラクタル次元を測定して画像処理を行う画像処理装置、フラクタル次元を測定してＸ線医療画像の処理を行うＣＴスキャナ及びフラクタル次元を高速測定するフラクタル次元測定方法に関する。
【０００２】
【従来の技術】
従来、医療分野では広くＣＴスキャナが用いられており、ＣＴスキャナから得られた２次元画像のシークエンスから医師が解剖学的知識にもとづいて、腫瘍の存在の検出を行っている。
【０００３】
しかしながら、一人の医師が読むＸ線画像の時間と数は限られているため、機械による診断支援により、癌などの早期発見率を上げることが望まれている。そのためには、より高速に画像処理を行う技術が必要である。
【０００４】
また、近年では、２次元画像から腫瘍を検出する方法として、腫瘍部のフラクタル次元の変化を捉えて検出する方法が考えられている。フラクタル次元による腫瘍検出の例としては、肝臓癌、乳癌、腎臓癌、肺癌などがある（V. Swarnakar et al., Fractal Based Characterization of Structural Changes in Biomedical Images, SPIE, Vol.2709, pp 444-455、及びJ. Levy Vehel, Using fractal and morphological criteria for automatic classification of diseases,SPIE Vol.1199, Visual Communications and Image Processing IV,1989,pp903-912）。
【０００５】
このフラクタル次元を測定することで、腫瘍部の形を量的に表現し、その評価を客観的にしかも定量的に行うことで、より精度の高い腫瘍検出を行うことができる。
【０００６】
フラクタル次元の測定では、ボックスカウンティング法が最もよく用いられている。図２４は２次元の画像にボックスカウンティング法を適用した場合の図である。
【０００７】
（Ａ）は２次元画像を示す。（Ｂ）は２次元画像の境界線を示す。（Ｃ）は２次元画像の境界線及び境界線を覆うボックスを示す。（Ｄ）は（Ｃ）のボックスを細分化したものである。
【０００８】
このように、２値で表現された２次元画像を正方形のボックスで分割し、その正方形の一辺の長さを変化させ、正方形と画像が重なり合う正方形の数から、フラクタル次元を求める。
【０００９】
図２５は３次元の画像にボックスカウンティング法を適用した場合の図である。多値で表現された３次元画像に対し、輝度の値を２次元平面の曲面として、その曲面をスケールを変えながら立方体のボックスで覆い、その立方体の数からフラクタル次元を求める。
【００１０】
一方、フラクタル次元は、テクスチャの認識、識別及び分類等に役立つために、医療分野以外でも幅広く応用されつつある。
例えば、特開平３−２６９７８２号公報では、文書の文字部分の認識にフラクタル次元を利用している。また、特開平８−１９７１００号公報では、汚泥を凝集させるための添加物の量を制御するために、フラクタル次元を利用して汚泥の画像を解析している。
【００１１】
【発明が解決しようとする課題】
しかし、上記のような従来のボックスカウンティング法を用いた２次元及び３次元画像のフラクタル次元の測定では、画像の一端から走査してボックスと画像が重なり合うボックスの数を検出するといった処理を、ボックスの一辺を変化させる度に行っているので、長い測定時間を要するといった問題があった。
【００１２】
また、上記のような従来技術は、特定の分野にフラクタル次元を応用するだけで、フラクタル次元を高速に測定することに関しては、何ら考慮されていない。フラクタル次元を利用した技術を実用化する際には、フラクタル次元を効率よく測定し、測定処理時間を短縮する必要がある。
【００１３】
特に集団検診などで医療画像を処理する場合には、Ｘ線画像を効率よく処理して、リアルタイムあるいはそれに近い短い時間でフラクタル次元を測定しなければならない。
【００１４】
本発明はこのような点に鑑みてなされたものであり、フラクタル次元を高速に測定して、画像処理を行う画像処理装置を提供することを目的とする。
また、本発明の他の目的は、フラクタル次元を高速に測定して、Ｘ線医療画像の処理を行うＣＴスキャナを提供することである。
【００１５】
さらに、本発明の他の目的は、フラクタル次元を高速に測定して、処理速度の向上を図るフラクタル次元測定方法を提供することである。
【００１６】
【課題を解決するための手段】
本発明では上記課題を解決するために、フラクタル次元を測定して画像処理を行う画像処理装置において、多値画像に対し、設定したしきい値以上の輝度値を高輝度値、前記しきい値以下の輝度値を低輝度値として、２値画像に変換する画像変換手段と、前記２値画像から一定領域の画像を抽出して抽出画像を生成する抽出画像生成手段と、前記抽出画像に解像度を設定し、前記解像度を階層的に変えて複数の解像度変換画像を生成する解像度変換画像生成手段と、下位の階層の前記解像度変換画像を構成する隣接画素のうち、一つでも前記高輝度値があれば、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値を前記高輝度値とし、前記隣接画素がすべて低輝度値の場合は対応する画素の輝度値を前記低輝度値とする輝度値設定手段と、前記高輝度値を持つ画素数の総和である高輝度画素数を前記解像度変換画像毎にそれぞれ計算する高輝度画素数計算手段と、前記解像度及び前記高輝度画素数の座標軸上に、前記高輝度画素数をプロットして生成したグラフの傾きを求め、前記傾きを前記フラクタル次元として前記２値画像全体の各画素にわたって測定するフラクタル次元測定手段と、前記フラクタル次元の画像を表示するフラクタル次元画像表示手段と、を有することを特徴とする画像処理装置が提供される。
【００１７】
ここで、画像変換手段は、多値画像に対し、設定したしきい値以上の輝度値を高輝度値、しきい値以下の輝度値を低輝度値として、２値画像に変換する。抽出画像生成手段は、２値画像から一定領域の画像を抽出して抽出画像を生成する。解像度変換画像生成手段は、抽出画像に解像度を設定し、解像度を階層的に変えて複数の解像度変換画像を生成する。輝度値設定手段は、下位の階層の解像度変換画像を構成する隣接画素のうち、一つでも高輝度値があれば、隣接画素に対応する上位の階層の解像度変換画像の画素の輝度値を高輝度値とし、隣接画素がすべて低輝度値の場合は対応する画素の輝度値を低輝度値とする。高輝度画素数計算手段は、高輝度値を持つ画素数の総和である高輝度画素数を解像度変換画像毎にそれぞれ計算する。フラクタル次元測定手段は、解像度及び高輝度画素数の座標軸上に、高輝度画素数をプロットして生成したグラフの傾きを求め、傾きをフラクタル次元として２値画像全体の各画素にわたって測定する。フラクタル次元画像表示手段は、フラクタル次元の画像を表示する。
【００１８】
また、フラクタル次元を測定して画像処理を行う画像処理装置において、多値画像から一定領域の画像を抽出して抽出画像を生成する抽出画像生成手段と、前記抽出画像に解像度を設定し、前記解像度を階層的に変えて、複数の第１の解像度変換画像及び複数の第２の解像度変換画像を生成する解像度変換画像生成手段と、前記第１の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最大輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とする最大輝度値設定手段と、前記第２の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最小輝度値を、前記隣接画素に対応する上位の階層の前記解像度変換画像の画素の輝度値とする最小輝度値設定手段と、前記第１の解像度変換画像と前記第２の解像度変換画像の同一座標の画素の輝度値の差分の総和である輝度差分値を計算する輝度差分値計算手段と、前記解像度及び前記輝度差分値の座標軸上に、前記輝度差分値をプロットして生成したグラフの傾きを求め、前記傾きを前記フラクタル次元として前記多値画像全体の各画素にわたって測定するフラクタル次元測定手段と、前記フラクタル次元の画像を表示するフラクタル次元画像表示手段と、を有することを特徴とする画像処理装置が提供される。
【００１９】
ここで、抽出画像生成手段は、多値画像から一定領域の画像を抽出して抽出画像を生成する。解像度変換画像生成手段は、抽出画像に解像度を設定し、解像度を階層的に変えて、複数の第１の解像度変換画像及び複数の第２の解像度変換画像を生成する。最大輝度値設定手段は、第１の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最大輝度値を、隣接画素に対応する上位の階層の解像度変換画像の画素の輝度値とする。最小輝度値設定手段は、第２の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最小輝度値を、隣接画素に対応する上位の階層の解像度変換画像の画素の輝度値とする。輝度差分値計算手段は、第１の解像度変換画像と第２の解像度変換画像の同一座標の画素の輝度値の差分の総和である輝度差分値を計算する。フラクタル次元測定手段は、解像度及び輝度差分値の座標軸上に、輝度差分値をプロットして生成したグラフの傾きを求め、傾きをフラクタル次元として多値画像全体の各画素にわたって測定する。フラクタル次元画像表示手段は、フラクタル次元の画像を表示する。
【００２０】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して説明する。図１は画像処理装置の原理図である。画像処理装置１０は、フラクタル次元を測定して画像処理を行う。
【００２１】
画像変換手段１は、多値画像Ｉに対し、輝度の分布を調べて、設定したしきい値以上の輝度値（画素値）を高輝度値（輝度値＝１）、しきい値以下の輝度値を低輝度値（輝度値＝０）として、２値画像Ｉａに変換する。
【００２２】
抽出画像生成手段２は、２値画像Ｉａから一定領域の画像を抽出して抽出画像Ｉｂを生成する。
解像度変換画像生成手段３は、抽出画像Ｉｂに解像度（ボックス）を設定し、解像度を階層的に変えて複数の解像度変換画像を生成する。なお、解像度を変えるとは、ボックスの大きさを変えることである。
【００２３】
例えば、解像度変換画像Ａの解像度は１６×１６ビット、解像度変換画像Ｂの解像度は８×８ビット、解像度変換画像Ｃの解像度は４×４ビットというように、階層毎に解像度を下位から上位に向かって粗くして、解像度変換画像を生成する。
【００２４】
また、解像度の変換を、この例では倍ずつに低解像度（ボックスの大きさを２倍ずつ）としたが、３倍、４倍、あるいは解像度毎に異なった倍率で次の低解像度を定義してもよい。また、画像内の領域毎に異なった倍率を与えるようにしてもよい。これらのことは、以降の実施の形態すべてに当てはまる。このように解像度を倍毎に変えることが計算量削減につながることは後述する。
【００２５】
また、抽出画像生成手段２は、実際には最初に設定した解像度で一定領域の画像を抽出するので、抽出画像Ｉｂと解像度変換画像Ａは同一（１６×１６ビット）の画像である。
【００２６】
輝度値設定手段４は、下位の階層の解像度変換画像を構成する隣接画素のうち、一つでも高輝度値があれば、隣接画素に対応する上位の階層の解像度変換画像の画素の輝度値を高輝度値とし、隣接画素がすべて低輝度値の場合は対応する画素の輝度値を低輝度値とする。
【００２７】
すなわち、下位の階層の解像度変換画像Ａを構成する隣接画素（隣接４画素とする）Ｐａ１には、高輝度値が１つあるので、隣接４画素Ｐａ１に対応する上位の階層の解像度変換画像Ｂの画素Ｐｂ１の輝度値は高輝度値となる。
【００２８】
また、隣接４画素Ｐａ２はすべて低輝度値であるから、対応する画素Ｐｂ２の輝度値は低輝度値となる。なお、各階層における画素の対応づけ等の説明は後述する。
【００２９】
高輝度画素数計算手段５は、高輝度値を持つ画素数の総和である高輝度画素数を解像度変換画像毎にそれぞれ計算する。
図では、高輝度画素数をＳ（Ａ）、Ｓ（Ｂ）、Ｓ（Ｃ）としている。ここで、解像度変換画像Ａの高輝度画素数が例えば５個（１の数が５個）あれば、Ｓ（Ａ）＝５である。
【００３０】
フラクタル次元測定手段６は、解像度をＸ軸、高輝度画素数をＹ軸にとり（実際にはそれぞれ対数をとって）、この座標軸上に、高輝度画素数（ここでは、Ｓ（Ａ）、Ｓ（Ｂ）、Ｓ（Ｃ））をプロットして生成したグラフの傾きａを例えば、最小２乗法等で求める。
【００３１】
そして、この傾きａをフラクタル次元として、２値画像Ｉａ全体の各画素にわたって測定する。この場合、抽出画像Ｉｂを２値画像Ｉａ上を移動させることによって、すべての画素のフラクタル次元を測定する。
【００３２】
フラクタル次元画像表示手段７は、測定されたフラクタル次元の画像を表示する。また、後述のフラクタル次元の統計的分布を表示してもよい。
次にフラクタル次元の定義及びフラクタル次元の算出について説明する。
【００３３】
平面が一辺ｒの正方形（ボックス）に分割された画面内に画像を取り込む。そして、画像の少なくとも一部を含む正方形の数をＮ（ｒ）とし、Ｎ（ｒ）の一辺ｒの大きさを種々の大きさに変えると、
【００３４】
【数１】
Ｎ（ｒ）・ｒD ＝Ｃ  （Ｃは定数）                          …（１）
となる時、Ｄをフラクタル次元と呼ぶ。
【００３５】
また、式（１）の両辺の対数をとると、
【００３６】
【数２】
ｌｏｇＮ（ｒ）＝−Ｄｌｏｇｒ＋ｌｏｇＣ                    …（２）
となる。
【００３７】
図２は線形回帰直線を示す図である。式（２）が線形回帰直線である。Ｘ軸にｌｏｇｒ、Ｙ軸にｌｏｇＮ（ｒ）をとる。図のようにｌｏｇｒとｌｏｇＮ（ｒ）のグラフは直線関係である。
【００３８】
したがって、ｒとＮ（ｒ）との関係をプロットして、線形回帰直線を最小２乗法により算出すれば、その傾きがフラクタル次元Ｄとして求まる。本発明では、このフラクタル次元を高速に算出して画像処理を行うものである。
次に画像処理装置１０のフラクタル次元の測定を図３〜図５を用いてさらに詳しく説明する。図３、図４はフラクタル次元測定の処理手順を示すフローチャートである。
〔Ｓ１〕画像変換手段１は、多値画像Ｉを２値画像Ｉａに変換する。以降、ステップＳ２からステップＳ１２は、解像度を変換しながら行う繰り返し処理を順次説明するものである。
〔Ｓ２〕抽出画像生成手段２は、２値画像Ｉａから抽出画像Ｉｂを生成し、解像度変換画像生成手段３は、抽出画像Ｉｂに１６×１６ビットの解像度を設定する。すなわち、解像度変換画像Ａは、１６×１６ビット解像度の画像である。
〔Ｓ３〕高輝度画素数計算手段５は、１６×１６ビットの解像度変換画像Ａ内の輝度値１の総和を求め、それをＳ（Ａ）とする。
【００３９】
ここでは、この解像度変換画像Ａの中心の隣接する４画素（７，７）、（７，８）、（８，７）、（８，８）のフラクタル次元を求めるものとする。
なお、１６×１６ビットのうち、フラクタル次元はどのビットでも対応させることができるが、一番意味のある対応づけがこの中心部の隣接する４ビットに対応させることである（この例は、エプソン社のEndeavor ET-6500(CPU Pentium 133MHz)上でインプリメントされた。また、人間の目で見て、その視覚に訴えるには、１画素よりも隣接４画素で同じフラクタル次元を表わす方が効果的であることが知られている。）。
〔Ｓ４〕解像度変換画像生成手段３は、１６×１６ビットの解像度変換画像Ａに対して、８×８ビットの解像度変換画像Ｂを用意する。そして、解像度変換画像Ａの座標（ｉ，ｊ）の画素に解像度変換画像Ｂの座標（〔ｉ／２〕，〔ｊ／２〕）の画素が対応づけられる。
【００４０】
ここで〔ｎ〕は、記号（ガウス記号）内の数値を越えない最大の整数とする。したがって、解像度変換画像Ａの隣接する４画素が解像度変換画像Ｂの１画素に対応することになる。
〔Ｓ５〕輝度値設定手段４は、この隣接４画素のうち一つでも輝度値が１ならば、解像度変換画像Ｂに対応する画素の輝度値は１にする。また、隣接４画素すべての輝度値が０ならば、解像度変換画像Ｂに対応する画素の輝度値は０にする。〔Ｓ６〕高輝度画素数計算手段５は、８×８ビットの解像度変換画像Ｂ内の輝度値１の総和を求め、Ｓ（Ｂ）とする。
〔Ｓ７〕解像度変換画像生成手段３は、８×８ビットの解像度変換画像Ｂに対して、４×４ビットの解像度変換画像Ｃを用意する。解像度変換画像Ｂの（ｉ，ｊ）の画素に解像度変換画像Ｃの画素（〔ｉ／２〕，〔ｊ／２〕）の画素が対応づけられる。したがって、解像度変換画像Ｂの隣接する４画素が解像度変換画像Ｃの１画素に対応する。
〔Ｓ８〕輝度値設定手段４は、この隣接する４画素のうち一つでも輝度値が１ならば、解像度変換画像Ｃに対応する画素の輝度値は１にする。また、隣接４画素すべての輝度値が０ならば、解像度変換画像Ｃに対応する画素の輝度値は０にする。
〔Ｓ９〕高輝度画素数計算手段５は、４×４ビットの解像度変換画像Ｃ内の輝度値１の総和を求め、Ｓ（Ｃ）とする。
〔Ｓ１０〕解像度変換画像生成手段３は、４×４ビットの解像度変換画像Ｃに対して、２×２ビットの解像度変換画像Ｄを用意する。また、対応する画素を上記同様に定義する。
〔Ｓ１１〕輝度値設定手段４は、この隣接する４画素のうち一つでも輝度値が１ならば、解像度変換画像Ｃに対応する画素の輝度値は１にする。また、隣接４画素すべての輝度値が０ならば、解像度変換画像Ｃに対応する画素の輝度値は０にする。
〔Ｓ１２〕高輝度画素数計算手段５は、２×２ビットの解像度変換画像Ｄ内の輝度値１の総和を求め、Ｓ（Ｄ）とする。
【００４１】
ここで、Ｓ（Ａ）、Ｓ（Ｂ）、Ｓ（Ｃ）及びＳ（Ｄ）は、各対応する解像度で１を持つボックスの総和である。それぞれの解像度の各ボックスは、もとの１６×１６ビットの解像度変換画像Ａに対し、輝度値が１であるものが少なくとも１個はあり、その時に限って１である。
【００４２】
ボックスの大きさはそれぞれ、１×１ビット、２×２ビット、４×４ビット、８×８ビットであるので、そのボックスの１辺の対数（底２）をとれば０、１、２、３である。
〔Ｓ１３〕フラクタル次元測定手段６は、Ｘ軸にボックスの１辺の大きさの対数、Ｙ軸に各辺に対応するボックスの数Ｓ（Ａ）、Ｓ（Ｂ）、Ｓ（Ｃ）及びＳ（Ｄ）の対数をとる。そして、各点を最も近似する直線を最小２乗法により求め、その傾きからフラクタル次元を求める。
【００４３】
なお、この場合のフラクタル次元は、解像度変換画像Ａの中心隣接４画素のフラクタル次元と見なし、画像処理装置１０の記憶領域の同じ座標上の４画素に記憶する。
〔Ｓ１４〕１６×１６ビットの抽出画像Ｉｂ（解像度変換画像Ａ）を２値画像Ｉａに対してスライドさせながら（この例では、２ビットずつずらす。すなわち、２ビットずつずらした画像を抽出するのと同じ）画像全体の各画素におけるフラクタル次元を測定する。
【００４４】
そして、すべての各画素のフラクタル次元が求まった時点で、フラクタル次元画像は、フラクタル次元画像表示手段７により表示される。
図５は輝度値が設定された解像度変換画像の生成過程の流れを示す図である。解像度変換画像Ａから解像度変換画像Ｃまでを示している。
【００４５】
解像度変換画像Ａの隣接画素Ｐａ１は、４画素で構成され、輝度値１の画素と輝度値０の３つの画素からなる。隣接画素Ｐａ２〜Ｐａ４は、輝度値０の４画素からなる。
【００４６】
解像度変換画像Ｂの画素Ｐｂ１は、隣接画素Ｐａ１に対応し、上述したように輝度値が１となる。解像度変換画像Ｂの画素Ｐｂ２〜Ｐｂ４は、それぞれ隣接画素Ｐａ２〜Ｐａ４に対応し、輝度値が０となる。
【００４７】
解像度変換画像Ｃの画素Ｐｃは、画素Ｐｂ１〜〜Ｐｂ４の４つの画素からなる隣接４画素Ｐｂに対応し、輝度値が１となる。
なお、上記の説明では、最初の設定として解像度変換画像Ａを１６×１６ビットとしたが、別のビットのサイズでもよい。このサイズは画像処理装置１０のアプリケーションによって決まる。
【００４８】
この例では、フラクタル次元を画像のローカルの値として求める必要上、またディスプレイ上での画像の認識度を考慮して１６×１６ビットを採用した。
一般にｎ×ｎビットの画像全体のフラクタル次元を求めるには、ｎ／２×ｎ／２、ｎ／２2 ×ｎ／２2 、ｎ／２3 ×ｎ／２3 、…、２×２のサイズの補助のデータ領域があればよい。
【００４９】
また、抽出画像Ｉｂをスライドさせる際には、中心部の２×２ビットのフラクタル次元を求めているので、次の抽出画像は２ビットＸ軸またはＹ軸方向にずれた場所にとり、以降繰り返していく。
【００５０】
なお、同じ１６×１６の画像でも次の抽出画像を４ビットずらしたり、８ビットずつずらしたり、あるいはオーバーラップなしにして行ってもよい。
次にフラクタル次元測定の詳細アルゴリズムについて、図６、図７を用いて説明する。図６は抽出画像Ｉｂの生成手順を示すフローチャートである。
【００５１】
フラクタル次元測定時に、抽出画像生成手段２が２値画像Ｉａ上に抽出画像Ｉｂを数ビットスライドさせて、新たな抽出画像を生成する際の処理手順を示している。なお、以降では抽出画像をウインドとも呼ぶ。
【００５２】
ここでは、処理される入力はサイズがｎ×ｎの２値画像である。出力はサイズがｎ×ｎのフラクタル次元の２次元データＤ（ｉ，ｊ）である。簡単のため、ｎは２のべき乗の数であるとする。また、対数の底は２とする。
〔Ｓ２０〕アルゴリズムの始めに初期値として、Ｘ軸方向とＹ軸方向へのフラクタル次元を求めている画素の座標値の各ループでの増加量(increment＿x,increment ＿y)を与える。ここでは２とする。
【００５３】
もとの画像のうち、ある任意の画素を中心付近に含む２次元の解像度変換画像ＳL （ｉ，ｊ）の大きさwindow＿sizeを与える。ここでは１６(=24) とする。
すなわち、１６×１６サイズのウインドの単位でフラクタル次元を求める。この例の１６という数は、一般のｎに置き換えてもよい。その際の必要なメモリ領域の数（ウインドを含めて）をmax ＿level と呼ぶ。ここでは４の値を初期値として持つ。
【００５４】
なお、異なった階層（以降、レベル）が異なった解像度、あるいは異なったサイズのボックスに対応していると考えてよい。
ウインドを２次元の画像上でスライドさせながら、フラクタル次元を求めるが、そのウインドの中心座標をｉ，ｊで表し、初期値として、Ｘ座標、Ｙ座標とも７をとる。なお、画像の座標は、Ｘ軸、Ｙ軸とも０から始まるものとする。
〔Ｓ２１〕ｉとｊがその最大値 i＿max とj ＿max かどうかを判断する。そうであれば、フラクタル次元を必要とする画素すべてのフラクタル次元が求まったことになるので終了し、そうでなければステップＳ２２へ行く。
〔Ｓ２２〕Ｘ座標またはＹ座標をincrement ＿x またはincrement ＿y だけ増やす。なお、ここではＹ座標優先で行う。すなわち、一定のｉに対して、ｊを増やしてその最大値までにする。そこでｊをリセットして、かつｉをincrement ＿x だけ増やす。すなわち、２次元データを行毎に上から下へスキャンする。
〔Ｓ２３〕ある特定のｉ、ｊに対して、すなわち（ｉ，ｊ）座標を中心にもつ１６×１６のウインド領域からフラクタル次元を求める。なお、フラクタル次元ｄを求める詳細な処理は図７で説明する。
〔Ｓ２４〕フラクタル次元ｄを出力データＤ（ｉ，ｊ）上のウインドの中心部に対応する隣接４座標に記録する。すなわち、Ｄ（ｉ，ｊ）、Ｄ（ｉ，ｊ＋１）、Ｄ（ｉ＋１，ｊ）、Ｄ（ｉ＋１，ｊ＋１）に今計算されたｄが代入される。そして、ステップＳ２１に戻る。
【００５５】
次に解像度変換画像の生成、輝度値の設定及び高輝度画素数の計算を行って、フラクタル次元を測定する際の処理手順について説明する。図７はフラクタル次元測定の処理手順を示すフローチャートである。
〔Ｓ３０〕フラクタル次元の計算に必要な入力は、ウィンドウで囲まれた輝度値、ウインドウの元の画像での位置を表す（ｉ，ｊ）、及びレベルの最小値Lmin( 普通は１) と最大値Lmaxである。
〔Ｓ３１〕始めに現在のレベルＬを１(=Lmin) に初期化する。
〔Ｓ３２〕レベルが最大Lmaxに達しているかをテストする。達していればステップＳ３３へ、そうでなければステップＳ３５へ行く。
〔Ｓ３３〕ｄ1 ，ｄ2 ，…ｄLmaxから最小２乗法によってフラクタル次元ｄを求める。これらのサフィックスはレベルを表し、ｄ1 ，ｄ2 ，…ｄLmaxは各長さのボックスに対し、その中に含まれるレベル１のボックスのうち、少なくとも一つの値が１であるボックスの数である。そして、これらの点の集まりを最もよく近似する直線を最小２乗法で求める。この直線の傾きからフラクタル次元が求まる。
〔Ｓ３４〕フラクタル次元が求まったなら、元のルーチンに戻る。
〔Ｓ３５〕レベル番号を１増やす（ＬをＬ＋１）。さらに、各レベルで使われるポインタｉ、ｊをリセットする。
〔Ｓ３６〕未処理の隣接４画素があるかどうかを調べる。なければこのレベルの画素の計算が終わったことになり、ｄL にはこのレベルでの１の総和が入っているはずである。そして、ステップＳ３２に戻る。未処理の隣接４画素があればステップＳ３７へ行く。
〔Ｓ３７〕レベルの数Ｌが１増えたのに注意して、これまでに計算したボックスカウンティング法での値ＳL-1 （ｉ，ｊ）、ＳL-1 （ｉ，ｊ＋１）、ＳL-1 （ｉ＋１，ｊ）、ＳL-1 （ｉ＋１，ｊ＋１）のうち一つでも１があるかどうかを判断する。あればステップＳ３８へ、なければステップＳ４０へ行く。
〔Ｓ３８〕次のレベルでのボックスカウンティングのためデータ領域ＳL （〔ｉ／２〕，〔ｊ／２〕）の場所に１を記憶する。
〔Ｓ３９〕このデータ領域ＳL （ｉ，ｊ）の１の画素の総和を表すｄL の値を１増やす。そして、ステップＳ３６へ戻る。
〔Ｓ４０〕ＳL （〔ｉ／２〕，〔ｊ／２〕）の場所に０を記憶する。そして、ステップＳ３６へ戻る。
【００５６】
その後、解像度を変えることによって得られるデータから生成される直線を最小２乗法で求め、フラクタル次元を測定する。
次にフラクタル次元測定のアルゴリズムの正しさ(the validity of algorithm）について説明する。データが有限であり、かつレベル毎に計算に利用されるデータが４分の１に減っているので、有限の大きさのデータに対して必ず有限ステップで停止する。すなわち、上記のステップＳ３２、Ｓ３５、Ｓ３６から構成されるループは有限回で停止する。
【００５７】
一方、各レベルの各ボックスに含まれる元の画像の画素（レベル１のボックス）が一つでも１なら１であり、各ボックスに含まれる元の画像の画素全てが０の時０になるようにした。
【００５８】
この主張はステップＳ３７、３８、４０に対し、もし、ＳL-1 （ｉ，ｊ）、ＳL-1 （ｉ，ｊ＋１）、ＳL-1 （ｉ＋１，ｊ）、ＳL-1 （ｉ＋１，ｊ＋１）がすべて０であるならば、ＳL （〔ｉ／２〕，〔ｊ／２〕）＝０であり、また前者の一つでも１があれば、ＳL （〔ｉ／２〕，〔ｊ／２〕）＝１であるので、サフィックスＬについての数学的帰納法により簡単に証明される。このこととボックスカウンティング方法の正しさにより、このアルゴリズムの正当性は証明される。
【００５９】
次にこのアルゴリズムに必要なメモリーの大きさ（space complexity）について説明する。必要なメモリは計算に有効な２次元のメモリがレベル毎に４分の１になること、もとの入力ｎ2 データをのぞいて、
【００６０】
【数３】
（ｎ×２×ｎ／２＋ｎ×２2 ×ｎ／２2 ＋ｎ×２3 ×ｎ／２3 ＋…＋２×２）
＝ｎ2 （１／４＋１／４2 ＋１／４3 ＋…＋（ｎ×ｎ）／（ｎ／２）2 ）
＜ｎ2 （１／４＋１／４2 ＋１／４3 ＋…）
＝（１／３）×ｎ2                                             …（３）
の余分のスペースが必要となる。ここで２次元データの各要素には０または１の値をとる。なお、このスペースはｎ2 ／２2 で充分である（データ領域最小化の説明部分で後述）。
【００６１】
次にアルゴリズムの時間の複雑さ（time complexity ）について説明する。各レベルに対し、次のレベルの画素の値は２×２の小画像から計算される。
すなわち、隣接４画素の一つでも１があれば、次のレベルの対応する画素は１となり、隣接４画素のすべてが０の場合０となる。
【００６２】
これは３回の論理和をとるだけで充分である。これはコンピュータ内部では実質比較と同じオペレーションである（４画素の和が正かどうかで判断するとすれば、３回の和と１回の比較ですることもできる）。
【００６３】
そのような２×２の小画像は全部でｎ2 ／２2 ある。このことは各レベルでも成り立ち、かつ各レベルでの画像の大きさは前レベルの１／２2 になるので、すべてのレベルでの小画像の総和は、
【００６４】
【数４】
ｎ2 （１／２2 ＋１／２2 ×１／２2 ＋１／２2 ×１／２2 ×１／２2 ＋…＋１／（ｎ／２）2 ）＜ｎ2 （１／４＋１／４2 ＋１／４3 ＋…）
＝（１／３）×ｎ2                                             …（４）
である。したがって、全部で多くて３×（１／３）×ｎ2 ＝ｎ2 の比較が必要である。
【００６５】
また、加算については各レベルでの画素を足しあわせる必要があるので、最初の画素の加算ｎ2 を考慮すると全体で（４／３）ｎ2 （＝ｎ2 ＋（１／３）ｎ2 ）の加算が必要になる。
【００６６】
次にフラクタル次元の測定速度について説明する。従来の方法では、ボックスの各サイズにおいて、最悪の場合ｎ2 の比較が必要である。したがって、Ｌを異なったボックスの数（異なったスケールの数、あるいは本発明のアルゴリズムでのレベルの値に相当) とすれば、全体の比較の数はＬ×ｎ2 ある。同様に、ボックスの各サイズにおいて、加算は最悪の場合ｎ2 であり、全体の加算の数はＬ×ｎ2 である。
【００６７】
一方、本発明の計算アルゴリズムでは、可変させる異なるボックスの種類の数にかかわらずｎ2 の比較、また加算は（４／３）ｎ2 ですむ。
したがって、本アルゴリズムの速さＴ1 と従来のアルゴリズムの速さＴ2 の比＝Ｔ1/Ｔ2 ＝１／Ｌ（比較の場合）、あるいはＴ1/Ｔ2 ＝（４／３）・１／Ｌ（加算の場合）となる。数Ｌは最低３から４であり、実用上では１０から４０であるので、従来のアルゴリズムに比べて著しく優越性があることが分かる。
【００６８】
次にフラクタル次元画像表示手段７について説明する。図８は人間の肺のＸ線画像を示す図である。ＣＴスキャナから得られた肺のＸ線画像である。５１２×５１２ピクセルのモノクロ画像で各画素は８ビットのグレースケールのデータを保持する。
【００６９】
図９は図８のＸ線画像を処理したフラクタル次元画像を示す図である。各画素のフラクタル次元を求めたものである。両図ともコントラストとエッジを強調処理して見やすくしてある。
【００７０】
Ｘ線画像１００には、細かな暗影部ｍと画像領域１１０に大きな暗影部Ｍが映し出されており、大きな暗影部Ｍは癌腫瘍であろうと推測できる。
フラクタル次元画像１００ｆには、画像領域１１０に画像Ｍｆが島状になって現れている。したがって、Ｘ線画像１００の推測に踏まえて、さらに確実に画像領域１１０の部分に異常が発生していると判断することができる。
【００７１】
なお、もとのＸ線画像１００の細かく黒く現れた箇所（暗影部ｍ）が、必ずしもフラクタル次元が高く島状となって現れないことに注意する。
図１０、図１１はフラクタル次元の統計的分布を示す図である。人間の脳のＸ線画像を解析し、Ｘ軸にとったフラクタル次元に対する統計値をプロットした図である。さらに、図１１は図１０の３ヵ月後に撮った図である。
【００７２】
図１１では３ヵ月前にはなかった画像Ｍａが現れており、撮影した箇所に異常が発生していると判断できる。
このように、フラクタル次元画像表示手段７は、フラクタル次元の画像として、フラクタル次元画像またはフラクタル次元の統計的分布を表示するので、画像解析者にとって有効な情報を提供することができる。
【００７３】
次に多値画像（グレースケール画像）のフラクタル次元を測定する画像処理装置について説明する。図１２は多値画像のフラクタル次元を測定する画像処理装置の原理図である。画像処理装置１０ａは、多値画像のフラクタル次元を測定して画像処理を行う。
【００７４】
抽出画像生成手段１２は、多値画像Ｉから一定領域の画像を抽出して抽出画像Ｉｃを生成する。
解像度変換画像生成手段１３は、抽出画像Ｉｃに解像度を設定し、解像度を階層的に変えて、複数の第１の解像度変換画像及び複数の第２の解像度変換画像を生成する。
【００７５】
ここで、第１の解像度変換画像及び第２の解像度変換画像は、同一解像度の画像である。解像度変換画像生成手段１３は、階層毎に解像度を下位から上位に向かって粗くして、２系統の解像度変換画像を生成する。
【００７６】
また、抽出画像生成手段１２は、実際には最初に設定した解像度で一定領域の画像を抽出するので、抽出画像Ｉｃと、第１の解像度変換画像及び第２の解像度変換画像の最初の解像度変換画像とは同一である。
【００７７】
最大輝度値設定手段１４ａは、第１の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最大輝度値を、隣接画素に対応する上位の階層の解像度変換画像の画素の輝度値とする。
【００７８】
最小輝度値設定手段１４ｂは、第２の解像度変換画像に対し、下位の階層の解像度変換画像を構成する隣接画素の中の最小輝度値を、隣接画素に対応する上位の階層の解像度変換画像の画素の輝度値とする。最大輝度値及び最小輝度値の設定の詳細については後述する。
【００７９】
輝度差分値計算手段１５は、第１の解像度変換画像と第２の解像度変換画像の同一座標の画素の輝度値の差分の総和である輝度差分値を計算する。
フラクタル次元測定手段１６は、解像度及び輝度差分値の座標軸上に、輝度差分値をプロットして生成したグラフの傾きを最小２乗法により求め、この傾きをフラクタル次元として多値画像全体の各画素にわたって測定する。この場合、抽出画像Ｉｃを多値画像Ｉ上を移動させることによって、すべての画素のフラクタル次元を測定する。
【００８０】
フラクタル次元画像表示手段１７は、上述したようなフラクタル次元画像またはその統計的分布を表示する。
次に画像処理装置１０ａのフラクタル次元の測定を図１３〜図１５を用いてさらに詳しく説明する。図１３はフラクタル次元測定の処理手順を示すフローチャートである。
【００８１】
なお、以降の説明で、セルとはｎ×ｎの２次元の小画像のことであり、ボックスとは、このセル上に定義されるｎ×ｎ×ｎの３次元の小画像を意味する。
ここでの処理は、解像度変換画像Ａ１の中を２×２ビットの重なりのない小画像（２次元のボックス）で分割する。隣接する４画素のうち、左上の画素の座標を（ｉ，ｊ）とする。
【００８２】
そして、それぞれの小画像の４画素に対し、輝度の最大値を求める。同様に、解像度変換画像Ａ２の中を２×２ビットの小画像で分割する。４画素のうち、左上の画素の座標を（ｉ，ｊ）とする。そして、それぞれの小画像の４画素に対し、輝度の最小値を求める。
〔Ｓ５０〕抽出画像生成手段１２は、多値画像Ｉから抽出画像Ｉｃを生成し、解像度変換画像生成手段１３は、抽出画像Ｉｃに１６×１６ビットの解像度を設定して、第１の解像度変換画像Ａ１と第２の解像度変換画像Ａ２を生成する。
【００８３】
すなわち、第１の解像度変換画像Ａ１と第２の解像度変換画像Ａ２は、１６×１６ビット解像度の画像である。
ここでは、抽出画像Ｉｃの中心の隣接する４画素、（７，７）、（７，８）、（８，７）、（８，８）のフラクタル次元を求めるものとする。以降、ステップＳ５１からステップＳ５９までは、画像の解像度の変換を順次繰り返しながら行うことを説明する。
〔Ｓ５１〕解像度変換画像生成手段１３は、１６×１６ビットの解像度変換画像Ａ１、Ａ２に対して、８×８ビットの解像度変換画像Ｂ１、Ｂ２を用意する。そして、解像度変換画像Ａ１、Ａ２の画素（ｉ，ｊ）に解像度変換画像Ｂ１、Ｂ２の座標（〔ｉ／２〕，〔ｊ／２〕）の画素がそれぞれ対応づけられる。〔ｎ〕は上述したので説明は省略する。
〔Ｓ５２〕最大輝度値設定手段１４ａは、解像度変換画像Ａ１の隣接４画素内の最大輝度値を解像度変換画像Ｂ１の対応する画素に設定する。最小輝度値設定手段１４ｂは、解像度変換画像Ａ２の隣接４画素内の最小輝度値を解像度変換画像Ｂ２の対応する画素に設定する。
〔Ｓ５３〕輝度差分値計算手段１５は、解像度変換画像Ｂ１、Ｂ２の同一アドレスの差分の総和を求め、それをＳ（Ｂ）とする。
〔Ｓ５４〕解像度変換画像生成手段１３は、８×８ビットの解像度変換画像Ｂ１、Ｂ２に対して、４×４の解像度変換画像Ｃ１、Ｃ２を用意する。画素の対応づけは、上記と同様である。
〔Ｓ５５〕最大輝度値設定手段１４ａは、解像度変換画像Ｂ１の隣接４画素内の最大輝度値を解像度変換画像Ｃ１の対応する画素に設定する。最小輝度値設定手段１４ｂは、解像度変換画像Ｂ２の隣接４画素内の最小輝度値を解像度変換画像Ｃ２の対応する画素に設定する。
〔Ｓ５６〕輝度差分値計算手段１５は、解像度変換画像Ｃ１、Ｃ２の同一アドレスの差分の総和を求め、それをＳ（Ｃ）とする。
〔Ｓ５７〕解像度変換画像生成手段１３は、４×４ビットの解像度変換画像Ｃ１、Ｃ２に対して、２×２の解像度変換画像Ｄ１、Ｄ２を用意する。画素の対応づけは、上記と同様である。
〔Ｓ５８〕最大輝度値設定手段１４ａは、解像度変換画像Ｃ１の隣接４画素内の最大輝度値を解像度変換画像Ｄ１の対応する画素に設定する。最小輝度値設定手段１４ｂは、解像度変換画像Ｃ２の隣接４画素内の最小輝度値を解像度変換画像Ｄ２の対応する画素に設定する。
〔Ｓ５９〕輝度差分値計算手段１５は、解像度変換画像Ｄ１、Ｄ２の同一アドレスの差分の総和を求め、それをＳ（Ｄ）とする。
【００８４】
ここで、Ｓ（Ｂ）、Ｓ（Ｃ）及びＳ（Ｄ）は、各対応する解像度での各セルにおいて、輝度の最大値と最小値の差を各セルの一辺の長さで割った数である。
したがって、この数はボックスがセルの辺の長さだけの高さが合ったとした場合の立方体で輝度の最大と最大を含むように覆った場合の立方体の数である。
【００８５】
なお、立方体の高さをボックスと同じ長さにできるのは、輝度の最大値が画像の一辺の長さと同じ場合である。そうでない場合は補正を行う。
ボックスの大きさはそれぞれ、２×２ビット、４×４ビット、８×８ビットであるので、そのボックスの１辺の対数をとれば、１、２、３である。
〔Ｓ６０〕フラクタル次元測定手段１６は、Ｘ軸のボックスの１辺の大きさの対数、Ｙ軸に異なった辺の大きさに対応するボックスの数の対数をとる。
【００８６】
そして、各点をよく近似する直線を最小２乗法により求め、その傾きからフラクタル次元を求める。
なお、この場合のフラクタル次元は、解像度変換画像Ａの中心隣接４画素のフラクタル次元と見なし、画像処理装置１０ａの記憶領域の同じ座標の４画素に記憶される。
〔Ｓ６１〕１６×１６ビットの抽出画像Ｉｃを多値画像Ｉに対してスライドさせながら（この例では２ビットずつずらした画像を抽出する）画像全体の各画素におけるフラクタル次元を求める。
【００８７】
そして、すべての各画素のフラクタル次元が求まった時点で、フラクタル次元画像は、フラクタル次元画像表示手段１７により表示される。
図１５は最大輝度値及び最小輝度値が設定された解像度変換画像の生成過程の流れを示す図である。解像度変換画像Ａ１、Ａ２〜解像度変換画像Ｄ１、Ｄ２までを示している。
【００８８】
図では例えば、解像度変換画像Ａ１の隣接４画素の中の最大輝度値２４６が解像度変換画像Ｂ１の対応する画素に設定される。以降、解像度変換画像Ｄ１まで同様である。
【００８９】
また、解像度変換画像Ａ２の隣接４画素の中の最小輝度値８８が解像度変換画像Ｂ２の対応する画素に設定される。以降、解像度変換画像Ｄ２まで同様である。
【００９０】
なお、上記の説明では、最初の設定として解像度変換画像Ａ１、Ａ２を１６×１６ビットとしたが、別のビットのサイズでもよい。このサイズは画像処理装置１０ａのアプリケーションによって決まる。
【００９１】
この例では、フラクタル次元を画像のローカルの値として求める必要上、またディスプレイ上での画像の認識度を考慮して１６×１６ビットを採用した。
一般にｎ×ｎビットの画像全体のフラクタル次元を求めるには、ｎ／２×ｎ／２、ｎ／２2 ×ｎ／２2 、ｎ／２3 ×ｎ／２3 、…、２×２のサイズの補助のデータ領域があればよい。
【００９２】
また、この例ではフラクタル次元を求めるために、８×８，４×４，２×２の画像からＳ（Ｂ），Ｓ（Ｃ），Ｓ（Ｄ）を求めている。Ｓ（Ａ）を利用するために、もとの画像から抽出する小画像の大きさを３２×３２にしてもよいが、実験の結果、Ｘ線写真の場合、３つのＳ（＊）値でフラクタル次元を求めるのに充分であることが分かった。
【００９３】
さらに、抽出画像Ｉｃをスライドさせる際には、中心部の２×２ビットのフラクタル次元を求めているので、次の抽出画像は２ビットＸ軸またはＹ軸方向にずれた場所にとり、以降繰り返していく。
【００９４】
ただし、同じ１６×１６の画像でも次の抽出画像を４ビットずらしたり、８ビットずつずらしたり、あるいはオーバーラップなしにしてもよい。
以上説明したように、画像処理装置１０ａは、多値画像Ｉから第１と第２の解像度変換画像を構成し、それぞれの解像度を階層的に変化させて解像度変換画像を生成し、解像度変換画像に最大、最小輝度値を設定して、輝度値の差の総和を計算した後、フラクタル次元を測定する構成とした。これにより、解像度を変化させる度に画像をスキャンさせる必要がないので、フラクタル次元を高速に測定することが可能になる。
【００９５】
次にフラクタル次元測定の詳細アルゴリズムについて説明する。抽出画像Ｉｃの生成手順は図６と同様であるので、フラクタル次元を計算する部分のみを説明する。図１６はフラクタル次元の測定処理手順を示すフローチャートである。
〔Ｓ７０〕フラクタル次元計算ルーチンは、インプットとしてウインドウの元の画像上の所在を示す座標ｉ、ｊ（これによって、ウインドウ小画像を表す）、解像度のレベルの最大値を表すLmax、レベルの最小を表すLminがある。
【００９６】
Lminは普通１である。このルーチンのアウトプットはウインドウ小画像のフラクタル次元ｄである。
〔Ｓ７１〕LminをＬに入れる。また、ウインドウで囲まれた画素値Ｗ（ｉ，ｊ）をＳmax,L （ｉ，ｊ）、Ｓmin,L （ｉ，ｊ）にそれぞれにコピーする。すなわち、第一レベル（最下位のレベル）では各座標値の最大値と最小値は同じであるとする。これはアルゴリズムを簡単にするために行われる。
〔Ｓ７２〕ＬがLmaxに達したかどうかテストする。達した場合はステップＳ７３へ、そうでなければステップＳ７５へ行く。
〔Ｓ７３〕ｄ2 ，ｄ3 ，…ｄLmaxから（それぞれの対数を求めた後）最小２乗法によってフラクタル次元ｄを求める。
〔Ｓ７４〕フラクタル次元が求まったなら、もとのルーチンに戻る。なお、この例の場合は、図７のステップＳ３３のように、フラクタル次元の計算にｄ1 を利用しないことに注意する。
〔Ｓ７５〕Ｌを１だけ増やし、小画像のポインタｉ，ｊ及び、imax,jmax(現在の小画像の１辺の画素数（下記の説明でａとする）をリセットする。
〔Ｓ７６〕未処理の隣接４画素があるかどうかを調べる。なければステップＳ７７へ行き、あればステップＳ７８へ行く。
〔Ｓ７７〕現在のレベルの小画像の処理が終わったので、ｄL を画像の一辺の長さと輝度の比率（Ｇ／Ｍ）で補正したボックスの一辺の長さでもって割る。
【００９７】
【数５】
ｄL ←ｄL ／（２L ・Ｇ／Ｍ）                                …（５）
この時点で、ｄL ＝１／（２L ・Ｇ／Ｍ）max i,j=0...,(a-1)｛Ｓmax,L （ｉ，ｊ）―Ｓmin,L （ｉ，ｊ）｝を計算したことと同じになる。
【００９８】
ここでａ＝Ｍ／２L である。すなわち、ｉ、ｊは各レベルでのすべての座標値をとる。また、Ｇは可能な輝度の最大値、Ｍは元の画像の１辺の長さ（画素の数）である。そしてステップＳ７２へ戻る。
〔Ｓ７８〕Ｌが１だけ増えたので、Ｌ―１のレベルの隣接する４画素Ｓmax,L-1 （ｉ，ｊ）、Ｓmax,L-1 （ｉ，ｊ＋１）、Ｓmax,L-1 （ｉ＋１，ｊ）、Ｓmax,L-1 （ｉ＋１，ｊ＋１）から最大値（smax）を選ぶ。
【００９９】
同様にＬ―１のレベルの隣接する４画素Ｓmin,L-1 （ｉ，ｊ）、Ｓmin,L-1 （ｉ，ｊ＋１）、Ｓmin,L-1 （ｉ＋１，ｊ）、Ｓmin,L-1 （ｉ＋１，ｊ＋１），から最小値（smin）を選ぶ。
〔Ｓ７９〕ここで次のレベルの輝度のmax とmin の差分をこれまでの差分の部分和に加算する。さらに、Ｌのレベルに対し、smaxをＳmax,L （〔ｉ／２〕，〔ｊ／２〕）に、sminをＳmin,L （〔ｉ／２〕，〔ｊ／２〕）に記憶させる。その後、ステップＳ７６に戻る。
【０１００】
次にアルゴリズムの正当性(the validity of algorithm) について説明する。データが有限でありかつレベル毎に計算に利用されるデータが４分の１に減っているので、有限の大きさのデータに対して必ず有限ステップで停止する。
【０１０１】
今、ある座標（ｘ，ｙ）における輝度値ｆ（ｘ，ｙ）からなる擬似曲面を考える。するとある解像度で一つのセル内の最大値と最小値を求めて、それをセルの１辺の大きさで割るとその擬似曲面を上下から囲む３次元のボックスの数が求められる（正確にはセルの１辺×Ｇ／Ｍで補正したもの、ここで、Ｇは輝度の最大値、Ｍは元の画像の１辺の長さ）。
【０１０２】
したがって、すべてのセルの３次元ボックスの数の総和は、現在のセルに対応する解像度で見た時の擬似曲面の近似値を与える。なぜなら、ボックスは曲面を下と上から覆いつくし、また、画素同士の間隔が無限に小さくなれば、画素値の集まりが近似的になだらかな曲面を表し、無限に輝度からなる曲面に近づくからである。以上の説明及びボックスカウンティング法の正しさから、このアルゴリズムの正しさが証明される。
【０１０３】
次にアルゴリズムに必要なメモリの大きさ（space complexity）について説明する。必要なメモリは計算に有効な２次元のメモリがレベル毎に４分の１になること、及び各レベルで最大値と最小値を保持するために２倍のメモリが必要なことから、もとの入力ｎ2 データをのぞいて
【０１０４】
【数６】
２×（ｎ／２×ｎ／２＋ｎ／２2 ×ｎ／２2 ＋ｎ／２3 ＋…＋２×２）
＝２×ｎ2 （１／４＋１／４2 ＋１／４3 ＋…＋（ｎ×ｎ）／（ｎ／２）2 ）
＜２×ｎ2 （１／４＋１／４2 ＋１／４3 ＋…）＝（２／３）×ｎ2   …（６）
の余分のスペースが必要となる。なお、このスペースは、ｎ2 ／２（＝（ｎ2 ×２）／２2 ）で充分である（データ領域最小化の説明部で後述する）。
【０１０５】
次にアルゴリズムの時間の複雑さ（time complexity ）について説明する。各レベルで最大値のみを保持する２次元データと最小値のみを保持する２次元データがある。
【０１０６】
あるレベルの隣接する４画素から最大値を求めるには３回の比較ですむ。また、そのような４画素からなる小画像は全部でｎ2 ／２2 個ある。
次のレベルでの解像度変換画像１枚につき隣接４画素の数は常に４分の１になっていることを考慮にすると、各レベルでの隣接４画素の数の総和は、最大値と最小値を保持するデータ領域が必要であることを考慮して、
【０１０７】
【数７】
２×ｎ2 （１／２2 ＋１／４2 ＋１／８2 ＋…１／（ｎ／２）2 ）
＜２×ｎ2 （１／４＋１／４2 ＋１／４3 ＋…）
＝（２／３）×ｎ2                                             …（７）
したがって、全部で多くて３×（２／３）×ｎ2 ＝２×ｎ2 の比較が必要である。乗算はまったく必要としない。加算と減算、及び解像度が異なるために対する補正（すなわち、ボックスの大きさが変わるため）のための割り算はステップＳ７９で生じて、それぞれ補助画像データの異なったインデックスをもつ画素の総和回数だけ行う。すなわち、（１／３）×ｎ2 である。
【０１０８】
次にフラクタル次元の測定速度について説明する。本発明のアルゴリズムは、Sarkar and Chaudhuriらによるこれまで最速と知られているアルゴリズムの時間の複雑さより、はるかによいことがわかる。
【０１０９】
この論文によれば、画像の大きさをｎ×ｎとして、スケールダウンの比率をｒとすれば（本件のボックスの一辺のサイズは１／ｒに対応する) 、各解像度、すなわち、レベル（あるいは比率ｒ）に対して、２×ｎ2 の比較、４／ｒ2 の加算、２／ｒ2 の減算、３／ｒ2 の乗算、３／ｒ2 の除算を必要とする。
【０１１０】
全体でＬ異なる解像度で計算したとすると、全体の計算量は、単純に同じ解像度で計算が行われたとして、それぞれＬ×２×ｎ2 、Ｌ×４／ｒ2 、Ｌ×２／ｒ2 、Ｌ×３／ｒ2 、Ｌ×３／ｒ2 である（なお、最小２乗法を使うため、Ｌ＞２である）。
【０１１１】
２つのアルゴリズムを比較するために、条件を同じにしてSarkar and Chaudhuriらのアルゴリズムでも、異なる解像度を２、４、８... 画素間ずつ取ったとすると、簡単な計算の後、最後の４つは、それぞれ（４／３）×ｎ2 、（２／３）×ｎ2 、ｎ2 、ｎ2 となる。
【０１１２】
本発明のアルゴリズムでは２×ｎ2 の比較を必要とするが、加算、減算、乗算、除算はそれぞれ、（１／３）×ｎ2 、（１／３）×ｎ2 、０、（１／３）×ｎ2 である。これはすべての解像度に合わせての時間のコストであり、また、乗算は本発明で使われていないのは注目に値する。
【０１１３】
同じ解像度でフラクタル次元の計算がなされたと仮定して、本発明のアルゴリズムと従来のSarkar and Chaudhuriのアルゴリズムの比をとると、比較、加算、減算、乗算、除算の時間に対して、それぞれ、１／Ｌ、１／４、１／２、０、１／３となる。よって、全体としてこの新しいアルゴリズムは実用的な見地からも優れたアルゴリズムである。
【０１１４】
以上説明したように、画像処理装置１０ａは、解像度を階層的に変化させて第１及び第２の解像度変換画像を生成し、第１の解像度変換画像に最大輝度値の設定、第２の解像度変換画像に最小輝度値を設定し、輝度差分値を計算した後、最小２乗法によりフラクタル次元を測定する構成とした。
【０１１５】
これにより、解像度を変化させる度に画像をスキャンさせる必要がないので、フラクタル次元を高速に測定することが可能になる。
次にフラクタル次元画像生成手段について説明する。
【０１１６】
３次元データは、次のようにして作成する。ＣＴスキャンやＭＲＩから得られた人体の画像は、通常連続する２次元画像である場合が多い。
これをＧ０（ｘ，ｙ）、Ｇ１（ｘ，ｙ）、....Ｇｎ−１（ｘ，ｙ）、ｘ，ｙ＝０、...(画像の１辺の長さをｎで表わす）。これらを３次元のデータに宣言されたメモリ上の｛（ｘ，ｙ，０）, （ｘ，ｙ，１）,.....（ｘ，ｙ，ｎ−１）｝の座標値に入れると自動的に３次元画像が構成される。
【０１１７】
また、この３次元画像から、Ｚ軸方向についてスライスした画像は、固定したx'に対して、Ｙ−Ｚ軸からなる２次元平面｛(x',y,z)｝を考えればよい。これも、例えば、コンピュータのハードディスクを利用すれば容易に構成できる。
【０１１８】
データ構造としてＸ、Ｙ、Ｚ軸からなる直方体を考える。図１７はデータ構造の３次元座標を示す図である。各点の座標は（ｘ、ｙ、ｚ）で表現できる。各座標は２値、０または１の値を持つ。それをｆ（ｘ、ｙ、ｚ）で表わす。
【０１１９】
座標は直方体の奥下の角の一つを（０、０、０）原点とし、手前にＸ座標値が増える方向にＸ軸があり、右手にＹ座標値が増える方向にＹ軸があり、また上にＺ座標値が増える方法にＺ軸があるものとする。
【０１２０】
図１８はＸ−Ｙ平面の画像からＸ−Ｚ平面の画像を構成する様子を示す図である。複数のＸ−Ｙ軸からなる２次元画像（１からｎまで）を３次元のデータ領域からなる記憶部に読み込む。各２次元平面はＸ軸とＹ軸がそれぞれ同じ長さの辺からなるものとする。すなわち、ｋ番目の画像の座標（ｘ，ｙ）は３次元記憶部の座標（ｘ，ｙ，ｋ）に記憶される。
【０１２１】
したがって、この記憶部からＹ座標を固定してあるＹ＝ｍについて、（ｘ，ｍ，ｚ）のすべての値を読み出すと、ｍ番目の（ｘ，ｚ）軸からなる画像が構成される。なお、このような変換は通常のコンピュータのハードディスク上に３次元データを定義することによって容易に構成できる。
【０１２２】
３次元２値画像のフラクタル次元の測定については、まず、直方体（ボックス）は１６×１６×１６の大きさとして説明する。この直方体の構造をＡと呼ぶ。
最初に直方体Ａのなかの要素（画素に対応する）の値１の総和Ｓ（Ａ）を求める。次に８×８×８の大きさの直方体Ｂを考える。Ａにおいてf(i,j,k)、f(i+1,j,k)、f(i,j+1,k)、f(i+1,j+1,k)、 f(i,j,k+1) 、f(i+1,j,k+1)、f(i,j+1,k+1)、f(i+1,j+1,k+1)の８つの値の内一つでも１があれば、Ｂにおけるf(〔i/2 〕，〔j/2 〕, 〔k/2 〕）の値は１となる。また、８つの値の全てが０の時、 f（〔i/2 〕, 〔j/2 〕, 〔k/2 〕) の値は０となる。
【０１２３】
次に直方体Ｂのなかの要素の値１の総和Ｓ（Ｂ）を求める。同様に、４×４×４の直方体Ｃを考え、Ｃの各要素の値は、Ｂの各要素の値がＡの要素から決められたように、Ｂの隣接する８つの要素から定義される。
【０１２４】
次に直方体Ｃのなかの要素の値１の総和Ｓ（Ｃ）を求める。さらに、２×２×２の直方体Ｄを考え、Ｄの各要素の値は、これまで記述したように、Ｃの隣接する８つの要素から定義される。
【０１２５】
次に直方体Ｄのなかの要素の値１の総和Ｓ（Ｄ）を求める。最後に、Ｘ軸にボックスの１辺の大きさの対数、Ｙ軸に異なる長さの各辺に対応するボックスの数、Ｓ（Ａ），Ｓ（Ｂ），Ｓ（Ｃ），Ｓ（Ｄ）の対数をとれば、各点を最も近似する直線を最小２乗法により求め、その傾きからフラクタル次元を求める。
【０１２６】
一般的な形でのアルゴリズムはここまでの説明から容易に作成できる。また、アルゴリズムの解析はｎ×ｎ×ｎの入力データでもって行う。アルゴリズムの正当性は上述したものと同様に証明されるので省略する。
【０１２７】
次にアルゴリズムに必要なメモリの大きさ（space complexity）について説明する。あるレベルの３次元的に隣接する８要素の値が、次のレベルの１個の要素の値を決める。
【０１２８】
したがって、必要なメモリは計算に有効な３次元のメモリスペースがレベル毎に８分の１になること、もとの入力データを保持する領域ｎ3 をのぞいて、ほとんど上述したものと同様である。
【０１２９】
【数８】
（ｎ／２×ｎ／２×ｎ／２＋ｎ／２2 ×ｎ／２2 ×ｎ／２2 ＋ｎ／２3 ×ｎ／２3 ×ｎ／２3 ＋…＋２×２×２）
＝ｎ3 （１／８＋１／８2 ＋１／８3 ＋…＋（（ｎ×ｎ×ｎ）／（ｎ／２）3 ）＜ｎ3 （１／８＋１／８2 ＋１／８3 ＋…）
＝（１／７）×ｎ3                                             …（８）
に比例する余分のスペースが必要となる。ここで２次元データの各要素には０または１の値をとる。なお、このスペースはｎ3 ／８で充分である（データ領域最小化の説明部で後述する）。
【０１３０】
次にアルゴリズムの時間の複雑さ（time complexity ）について説明する。各レベルにおいて、次のレベルの要素の値は前のレベルの２×２×２の小画像から計算される。
【０１３１】
すなわち、８要素の一つでも１があれば、次のレベルの対応する要素は１となり、８要素の全てが０の場合０となる。これは多くとも７回の論理和オペレーションをするだけで充分である。これはコンピュータ内部では実質比較と同じオペレーションである。
【０１３２】
そのような小画像は全部でｎ3 ／２3 ある。したがって、各レベルにおいて、３次元画像の各辺の長さは半分になるので、全てのレベルを通して、
【０１３３】
【数９】
ｎ3 （１／２3 ＋１／４3 ＋１／８3 ＋…１／（ｎ／２）3 ）
＜ｎ3 （１／２3 ＋１／４3 ＋１／８3 ＋…）
＝（１／７）×ｎ3                                             …（９）
の小画像に比較処理をしたことになる。よって、全部で多くて７×( １／７) ×ｎ3 ＝ｎ3 の比較が必要である。
【０１３４】
また、加算については各レベルでの要素を足しあわせる必要があるので、最初の要素の加算ｎ2 を考慮すると全体で（８／７）ｎ2 （＝ｎ2 ＋（１／７）ｎ2 ）の加算が必要になる。
【０１３５】
アルゴリズムの測定速度についてはこの場合、３Ｄの対象物の形状そのもののフラクタル次元を、２次元のデータに落とすことなく効率良く測れる。その計算速度は２Ｄの場合に比べて遜色がない。
【０１３６】
次に画像処理装置１０ａで３次元の多値画像を扱う場合について説明する。すなわち、自然のカラー３Ｄオブジェクトを全て扱うものである。
始めに、元の多値画像から１６×１６×１６ビットの画像を抽出する。この画像のコピーを二つ作る。これらの解像度変換画像をＡ１とＡ２と呼ぶ。この例では１６×１６×１６ビットの中心部の座標のフラクタル次元を求める。
【０１３７】
ここで、１６は偶数なので、この小画像の中心部は隣接する８画素からなるものとする。すなわち、中心部とは、この例の場合、（７、７、７）、（８、７、７）、（７、８、７）、（８、８、７）…（８、８、８）からなる座標で示された隣接８要素のことである。
【０１３８】
この画像Ａ１の中を２×２×２ビットの小画像で分割する。隣接する８画素の内、左下角の画素の座標を（ｉ，ｊ，ｋ）とする。それぞれの小画像の８画素に対し、輝度の最大値を求める。同様に、この画像Ａ２の中を２×２×２ビットの小画像で分割する。４画素の内、左下角の画素の座標を（ｉ，ｊ，ｋ）とする。それぞれの小画像の８画素に対し、輝度の最小値を求める。
【０１３９】
１６×１６×１６ビットの画像Ａ１、Ａ２に対して、それぞれ８×８×８ビットの画像Ｂ１、Ｂ２を用意する。画像Ａ１またはＡ２の画素（ｉ，ｊ，ｋ）に対して、Ｂ１またはＢ２の画素（〔ｉ／２〕，〔ｊ／２〕，〔ｋ／２〕）が対応する。したがって、画像Ａ１（Ａ２）の隣接する８画素が画像Ｂ１（Ｂ２）の１画素に対応する。
【０１４０】
画像Ｂ１の座標（〔ｉ／２〕，〔ｊ／２〕，〔ｋ／２〕）にＡ１の隣接する８画素の最大値を記憶させる。同様にして、画像Ｂ２の座標（〔ｉ／２〕，〔ｊ／２〕，〔ｋ／２〕）にＡ２の隣接する８画素の最小値を記憶させる。
【０１４１】
そして、画像Ｂ１とＢ２の同一アドレスの差分の総和を求め、それをＳ（Ｂ）とする。同様に８×８×８ビットの画像Ｂ１、Ｂ２に対して、それぞれ４×４×４の画像Ｃ１、Ｃ２を用意する。
【０１４２】
画像Ｂ１（Ｂ２）の（ｉ，ｊ，ｋ）の画素がＣ１（Ｃ２）の座標（〔ｉ／２〕，〔ｊ／２〕，〔ｋ／２〕）に対応する。したがって、画像Ｂ１（Ｂ２）の８画素が画像Ｃ１（Ｃ２）の１画素に対応する。
【０１４３】
このＢ１の８画素の内、最大値を画像Ｃ１の対応する画素が受け取る。同様に、Ｂ２の８画素のうち、最小値を画像Ｃ２の対応する画素が受け取る。
画像Ｃ１とＣ２の同一アドレスの差分の総和を求め、それをＳ（Ｃ）とする。同様に４×４×４ビットの画像Ｃ１、Ｃ２に対して、２×２×２の画像Ｄ１、Ｄ２を用意して、その画素値を上記同様に定義する。
【０１４４】
画像Ｄ１とＤ２の同一アドレスの差分の総和を求め、それをＳ（Ｄ）とする。Ｓ（Ｂ），Ｓ（Ｃ），Ｓ（Ｄ）は各対応する解像度での、すなわち、大きさがそれぞれ、２×２×２ビット、４×４×４ビット、８×８×８ビットのボックスで測った時の、各ボックスでの輝度の最大値と最小値の差分の総和である。
【０１４５】
なお、そのボックスの１辺の対数（２を底数とする）をとれば、１、２、３である。
したがって、Ｘ軸のボックスの１辺の大きさの対数、Ｙ軸に長さの異なった辺に対応する輝度差分の総和の対数をとる。すると、各点を最も近似する直線を最小２乗法により求めると、その傾きからフラクタル次元が求まる。このフラクタル次元は画像Ａの中心８画素のフラクタル次元と見なす。
【０１４６】
そして、１６×１６×１６ビットの３Ｄウインドウを元の画像に対してスライドさせながら（例えばこの例では２ビットずつ）画像全体の各画素におけるフラクタル次元を求める。
【０１４７】
上記の１６×１６×１６ビットの小ウインドは１６×１６×１６ビットである必要はなく、別のサイズでもよい。サイズはアプリケーションによる。この例ではフラクタル次元を画像のローカルの値として求める必要から１６×１６×１６ビットになった。
【０１４８】
一般にｎ×ｎ×ｎビットの画像全体のフラクタル次元を求めるには、多くともｎ／２×ｎ／２×ｎ／２、ｎ／２2 ×ｎ／２2 ×ｎ／２2 、ｎ／２3 ×ｎ／２3 ×ｎ／２3 …２×２×２のサイズのデータ領域の２倍の領域があればよい。
【０１４９】
また、この例ではフラクタル次元を求めるために、８×８×８、４×４×４、２×２×２の画像からＳ（Ｂ）、Ｓ（Ｃ）、Ｓ（Ｄ）を求めている。
Ｓ（Ａ）を利用するために、もとの画像から抽出する小画像の大きさを３２×３２×３２にすることもできる。
【０１５０】
コンピュータによる実験によれば、Ｘ線写真の場合、３つのＳ（＊）値でフラクタル次元を求めれば充分である。別法としてＢ，Ｃ，ＤにつづいてＥを定義して、Ｓ（Ｃ），Ｓ（Ｄ），Ｓ（Ｅ）からフラクタル次元を計算するようにしてもよい。
【０１５１】
また、同じ１６×１６×１６の小ウインドウでも次のウインドを４ビットずらしたり、８ビットずつずらしたり、オーバーラップなしにしてもよい。
次にアルゴリズムに必要なメモリの大きさ（space complexity）について説明する。なお、アルゴリズムの正当性については上記と同様なので省略する。
【０１５２】
必要なメモリは計算に有効な３次元のメモリスペースがレベル毎に８分の１になること、及び各レベルで最大値と最小値を保持するために２倍のメモリが必要なことから、もとの入力ｎ3 データをのぞいて
【０１５３】
【数１０】
２×（ｎ／２×ｎ／２×ｎ／２＋ｎ／２2 ×ｎ／２2 ×ｎ／２2 ＋×ｎ／２3 ×ｎ／２3 ×ｎ／２3 ＋…＋２×２×２）
＝２×ｎ3 （１／８＋１／８2 ＋１／８3 ＋…＋（ｎ×ｎ×ｎ）／（ｎ／２）3 ＜２×ｎ３（１／８＋１／８2 ＋１／８3 ＋…）
＝（２／７）×ｎ3                                           …（１０）
の余分のスペースが必要となる。このスペースはｎ3 ／４（＝（ｎ2 ×２）／８）で充分である（データ領域最小化の説明部で後述する）。
【０１５４】
次にアルゴリズムの時間の複雑さ（time complexity ）について説明する。各レベルの各要素は、一つ手前のレベルの隣接する８要素の最大値のみを保持する３次元データと最小値のみを保持する３次元データがある。
【０１５５】
あるレベルの隣接する８要素から最大値を求めるには多くて７回の比較ですむ。また、そのような８要素からなる小画像は全部でｎ3 ／２3 ある。同様にして、３次元の画像の一辺は各レベル毎に半分になるので、
【０１５６】
【数１１】
２×ｎ3 （１／２3 ＋１／４3 ＋１／８3 ＋…＋１／（ｎ／２）3 ）
＜２×ｎ2 （１／２3 ＋１／４3 ＋１／８3 ＋…）
＝（２／７）×ｎ3                                             …（１１）
の小画像に対して比較計算を行うことになる。従って、全部で多くて７×（２／７）×ｎ3 ＝ｎ3 ＝２×ｎ3 の比較で処理が終わる。
【０１５７】
乗算は全く必要としない。加算と減算、及び、ボックスの大きさによる補正のための除算が補助面積の異なったインデックスの画素の総和だけ、すなわち、（１／７）×ｎ3 行われる。
【０１５８】
フラクタル次元測定速度については、従来は画像の次元が上がる度にボックスの大きさを増やせば増やすほど、多大の計算を必要としたが、本発明では、３次元データを２次元データに分解することなく、多値３次元画像のフラクタル次元の高速処理ができるようになった。その計算量は異なったサイズのボックスの数に依存しない。
【０１５９】
次にフラクタル次元を求めるためデータ領域を最小化する方法について説明する。あるレベルＬの２次元データの内隣接する４要素が次のレベルＬ＋１の一つの要素の値を決める。すべての隣接する４要素はオーバーラップがない背反の要素なので、次のレベルの要素は独立に（あるいは並列に）計算することができる。
【０１６０】
図１９はメモリスペースを節約する構造を示す図である。円で示された隣接４要素が矢印の先の要素の値を決める。
ＳL （ｉ，ｊ）、ＳL （ｉ，ｊ＋１）、ＳL （ｉ＋１，ｊ），ＳL （ｉ＋１，ｊ＋１）がＳL+1 （〔ｉ／２〕, 〔ｊ／２〕）の要素を決める。したがって、ｉ, ｊが行優先で増加するならば、次のレベルの要素を別の領域に入れずに、一旦、上記隣接する４要素の値を０にリセットしてからＳL+1 （〔ｉ／２〕, 〔ｊ／２〕）に入れればよい。行と列番号の最大値は各レベルが進む毎に半分になるので、各レベルでデータのオーバラップを生じることはない。
【０１６１】
上述した式（３）、（６）、（８）、（１０）のメモリスペースも同様にデータ領域を最小化することができる。
次にフラクタル次元の測定で複数の連続画像を取り扱う場合の効率的なアルゴリズムについて図２０〜図２２を用いて説明する。
【０１６２】
図２０はフラクタル次元の測定で複数の連続画像を取り扱う場合の概念を示す図である。
画像の量を減らすために、最初の画像を除いて、２枚目からの画像は前画像からの差分のみを記憶するようにする。このアルゴリズムで２枚目以降のフラクタル次元計算に差分のあった所のみ計算し直すことによって計算されることを示す。なお、ここで説明してあるアルゴリズムのプロトタイプはAkia Mircobook 55(MMX pentium 233MHz) 上でその動作を確認した。
【０１６３】
１枚目と２枚目の画像の間に差分がある場合には２枚目の画像はその差分がある座標に１の値を持たせる。そうでない場合は０にする。２枚目の元の画像は、１枚目の各要素をＨ１（ｉ，ｊ）で２枚目をＨ２（ｉ，ｊ）で表わし、差分の画像をＨＤ２（ｉ，ｊ）で表現したとすれば、次の関係を満たす。
【０１６４】
【数１２】
Ｈ２（ｉ，ｊ）＝Ｈ１（ｉ，ｊ）○ＨＤ２（ｉ，ｊ）          …（１２）
ここで○は排他的論理和(Exclusive or)( すなわち、もしx ≠y であれば、ｘ○ｙ=1、もしx=y であれば、ｘ○ｙ=0) である。
【０１６５】
なお、ここでの計算はブール代数（Boolean Algebra ）であるので、容易に
【０１６６】
【数１３】
ＨＤ２（ｉ，ｊ）＝Ｈ１（ｉ，ｊ）○Ｈ２（ｉ，ｊ）          …（１３）
であることがわかる。すなわち差分を求める方法と差分を１枚目に重ねて２枚目を得る方法は同じである。
【０１６７】
３枚目の画像は１枚目に２枚目の変更を重ねて得られる画像（２枚目の元の画像）) と３枚目の画像の差を記録させる。４枚目以降も同様にする。
図２１、図２２は複数の連続画像を取り扱う場合のフラクタル次元測定の処理手順を示すフローチャートである。
〔Ｓ９０〕始めに、１枚目の画像の処理を図６、図７で説明したアルゴリズムを使って行う。ただし、中間補助データＳ0 （ｉ，ｊ）、ｄ0 、ＳL （ｉ，ｊ）、ｄ1 、Ｓ2 （ｉ，ｊ）、ｄ2 、Ｓ3 （ｉ，ｊ）、ｄ3 …はすべてそのままにしておく。なお、Ｓ0 （ｉ，ｊ）には前の画像のインプットが入っている。
〔Ｓ９１〕次の差分画像があるかないかを調べる。なければ終了し、そうでなければ、ステップＳ９２へ行く。
〔Ｓ９２〕Ｌに１を設定する。
〔Ｓ９３〕次の画像、すなわち、最初の差分画像を読み込む。これをＳ-L（ｉ，ｊ）で表わす。
〔Ｓ９４〕未処理の隣接４要素の選択において、もし隣接する４要素Ｓ-L（ｉ，ｊ）の内一つでも０でないものがあるならば、すなわち、一つでも１があるならば、今問題となっている隣接４要素について、排他的論理和Ｓ0 （ｉ，ｊ）＝Ｓ0 （ｉ，ｊ）○Ｓ-L（ｉ，ｊ）を行って、Ｓ0 （ｉ，ｊ）の補正を行う。
〔Ｓ９５〕次に、現在のレベルＬがレベルの最大値Ｌmax より大きいかを調べ、大きければ次に差分画像をとる為に、ステップＳ９１に戻る。そうでなければ、ステップＳ９６へ行く。
〔Ｓ９６〕初期設定の後の繰り返し処理について最初に説明する。補正後のＳ0 （ｉ，ｊ）の隣接４要素の論理和をＡと置く。そのＡとＳ1 （〔ｉ／２〕、〔ｊ／２〕）の値が違えば、Ｓ1 （〔ｉ／２〕、〔ｊ／２〕）をＡで置き換えると同時に、Ｓ-2（〔ｉ／２〕、〔ｊ／２〕）＝１とする。
【０１６８】
そして、今度はＳ-2（ｉ，ｊ）＝１であるようなインデックス（ｉ，ｊ）についてＳ1 （ｉ，ｊ）の隣接４要素の論理和を求める。Ｓ2(〔ｉ／２〕、〔ｊ／２〕）の値との比較を行う。 両者に違いがあれば、Ｓ-3（〔ｉ／２〕、〔ｊ／２〕）＝１とする。
【０１６９】
したがって、Ｓ0 （ｉ，ｊ）の補正の後の一般の繰り返しループは次のようになる。一般にＳ-K（ｉ，ｊ）はＳK-1 （ｉ，ｊ）のどの要素を更新するかを決めるパイロットの役目をする。別の言葉で言えば、Ｓ-K（ｉ，ｊ）＝１である要素のインデックスについてそのインデックスを含む隣接４要素ＳK-1 （ｉ，ｊ）の更新が行われる。
【０１７０】
補正が行われて隣接４要素Ｓ-L（ｉ，ｊ）について、少なくとも１の要素があるか調べる。
そのような要素がなければ、ステップＳ９７へ行き、そうでなければステップＳ９８へ行く。
〔Ｓ９７〕レベルを１増やし、ステップＳ９５に戻る。
〔Ｓ９８〕この（ｉ，ｊ）を含む隣接４要素ＳL-1 （ｉ，ｊ）の論理和をＡと記す。隣接４要素は２次元のデータを２×２ずつのオーバーラップのない小画像に分割することによって定義される。
〔Ｓ９９〕そして、一つ上のレベルの対応要素の値ＳL （〔ｉ／２〕，〔ｊ／２〕）と比べ、変化があれば（Ａ≠ＳL （〔ｉ／２〕，〔ｊ／２〕）ステップＳ１００へ、なければステップＳ１０４へ行く。
〔Ｓ１００〕ＳL （〔ｉ／２〕，〔ｊ／２〕）にＡを入れると同時にＳ-(L+1)（〔ｉ／２〕，〔ｊ／２〕）= １とする。
〔Ｓ１０１〕Ａが１ならステップＳ１０２へ、そうでなければステップＳ１０３へ行く。
〔Ｓ１０２〕ＳL （〔ｉ／２〕，〔ｊ／２〕）が０から１にリセットされるのであれば（リセットさえる前に４要素のうちあるｉ，ｊが存在して、Ａ＝ＳL-1 （ｉ，ｊ）＝１≠０＝ＳL （〔ｉ／２〕，〔ｊ／２〕）が成り立っている）、ｄL を１増やす。
〔Ｓ１０３〕ＳL （〔ｉ／２〕，〔ｊ／２〕）が１から０にリセットされるのであれば（リセットさえる前に４要素の全てのインデックスｉ，ｊについて、Ａ＝ＳL-1 （ｉ，ｊ）＝０≠１＝ＳL （〔ｉ／２〕，〔ｊ／２〕）が成り立っている）、ｄｌを１減らす。
〔Ｓ１０４〕Ｓ-L（〔ｉ／２〕，〔ｊ／２〕）＝０とする。したがって  変化のない場合これ以降のレベルでの補正は不必要となる。そして、ステップＳ９６へ戻る。
【０１７１】
なお、全体のアルゴリズムは以下の通りである。（１）１枚目の画像のフラクタル次元を求める。次元を求めるための中間のデータを保持しておく。（２）２枚目（ｎ枚目) の差分画像を読みこむ。読み込まれた差分画像はＳ-L（ｉ，ｊ）に保持される。（３）ｎ−１枚目（１枚目）の中間データをｎ枚目（２枚目）の差分画像でもって、上記に説明したアルゴリズムを利用して更新し、ｎ枚目（２枚目）の画像のフラクタル次元を求める。
【０１７２】
このアルゴリズムにおいて、差分の部分だけ再計算が行われるようにするためには、サフィックスがマイナスであるＳ-L（ｉ，ｊ）の要素値が１であるものをリンクで繋げればよい。但し、単純にリンクで繋げた場合、対応する隣接４要素ＳL （ｉ，ｊ）の値を同時に調べることが出来ないのでインプリメンテーションが複雑になる。
【０１７３】
これを回避する方法として、Ｓ-L（ｉ，ｊ）に対応してサイズが行列とも２分の１ずつの２次元データＳ-L' （ｉ，ｊ）を用意する。これはＳ-(L+1)（ｉ，ｊ）のデータ構造と同じであるが、区別する為にＳ-L' と表記する。このＳ-L' （ｉ，ｊ）は、Ｓ-L（ｉ，ｊ）の隣接４要素の内一つでも１があれば１、そうでなければ０とする。そして、Ｓ-L' （ｉ，ｊ）＝１なる要素をリンクで結ぶ。よって、このリンクを辿りながら、隣接４要素ＳL-1 （２＊ｉ，２＊ｊ）、ＳL-1 （２＊ｉ＋１，２＊ｊ）、ＳL-1 （２＊ｉ，２＊ｊ＋１）、ＳL-1 （２＊ｉ＋１，２＊ｊ＋１）の要素の１の存在をこれまで説明したように検査すればよい。なお、このようなリンクづけは差分画像を得られる時に同時に作ればよい。
【０１７４】
更に、２枚目以降の画像のフラクタル次元を求めるのに差分があった所だけを、中間データを更新していけば充分である。よって、２枚目以降のフラクタル次元の計算を１枚ずつ初めからやり直して計算する時に比べて、差分があったところだけ（Ｓ-L（ｉ，ｊ）＝１となるようなｉ，ｊについてのみ）計算し直せばよいので大幅に計算時間を減らすことができる。また、２枚目からの画像は一つ手前の画像からの差分だけを記憶すればよく、連続する画像全体を保持する為の記憶スペースを大幅に節約することが出来る。
【０１７５】
差分画像の差分のある部分の大きさは、元の画像に比べて約１／１０から１／１００（＝ｋ）の大きさである。従って、例えば、１枚のオリジナルの画像のフラクタル処理時間をＣ１とするとｎ枚のオリジナル画像を処理した時間のコストはｎ×Ｃ１であるが、１枚目だけオリジナルを２枚目以降を差分画像を利用した場合の時間のコストはＣ１＋（ｎ−１）×ｋ×Ｃ１である。
【０１７６】
従って、この場合のアルゴリズムの時間のコストと上述した２次元２値画像でのアルゴリズムの時間のコストは（Ｃ１＋（ｎ−１）×ｋ×Ｃ１）／（ｎ×Ｃ１）である。この値はｎが大きくなればなるほどｋに近づき、ｋは通常の動画等では約１／１０から１／１００であるので、その効果は大となる。同様なことがメモリスペースについて言える。
【０１７７】
なお、ここでは連続する２値画像について説明したが、連続する多値画像の場合も０，１のかわりに、ｍａｘとｍｉｎを扱うため、計算が複雑になるが全体のアルゴリズムの流れは同様にして処理することができる。
【０１７８】
特にＭＲＩやＣＴのように人体の輪切り画像の場合、連続画像である場合が多く、隣接する画像の相関(coherence) が非常に高いのでデータの省スペース、ひいては計算時間の削減を大幅に達成できる。このことは、遠隔医療診断において多量の画像を送らねばならないことを考えると非常に重要な効果である。
【０１７９】
次にフーリエ変換画像のフラクタル次元測定について説明する。解析される元の画像は２値画像でも多値画像Ｉ（ｉ，ｊ）でもよい。この画像はフラクタルブラウン関数であると仮定すると、そのフーリエ変換された画像のパワースぺクトラム（fourier power spectrum）はこの関数の密度関数は１／ｆ２Ｈ＋１に従うことが知られている。ここで、ｆは帯域を表わす。このＨを求めることにより、フラクタル次元Ｄ＝２−Ｈが求められる。
【０１８０】
次にＣＴスキャナについて説明する。図２３はＣＴスキャナの構成を示す図である。ＣＴスキャナ２０は、フラクタル次元を測定して、Ｘ線医療画像の処理を行う。
【０１８１】
ＣＴスキャナ２０は、ガントリ２００と寝台２０１を有する。多値画像収集手段２１は、被検体をＸ線撮影した医療画像である多値画像を収集する。画像処理手段２２は、画像処理装置１０、１０ａを適用したものであり、多値画像からフラクタル次元を測定し表示する。詳細は上述したので説明は省略する。
【０１８２】
以上説明したように、ＣＴスキャナ２０は、フラクタル次元を高速に測定する画像処理手段２２を設ける構成とした。これにより、Ｘ線画像を効率よく高速に処理することが可能になる。
【０１８３】
【発明の効果】
以上説明したように、本発明の画像処理装置は、多値画像から２値画像に変換し、解像度を階層的に変化させて解像度変換画像を生成し、解像度変換画像に輝度値を設定して、高輝度画素数を計算した後、フラクタル次元を測定する構成とした。これにより、解像度を変化させる度に画像全体をスキャンする必要がなくなるので、効率よくボックスと画像が重なり合うボックス数の検出を行うことができ、フラクタル次元を高速に測定することが可能になる。
【０１８４】
また、本発明の画像処理装置は、多値画像に対し、解像度を階層的に変化させて第１及び第２の解像度変換画像を生成し、第１の解像度変換画像に最大輝度値の設定、第２の解像度変換画像に最小輝度値を設定し、輝度差分値を計算した後、フラクタル次元を測定する構成とした。これにより、解像度を変化させる度に画像全体をスキャンする必要がなくなるので、効率よくボックスと画像が重なり合うボックス数の検出を行うことができ、フラクタル次元を高速に測定することが可能になる。
【図面の簡単な説明】
【図１】  画像処理装置の原理図である。
【図２】  線形回帰直線を示す図である。
【図３】  フラクタル次元測定の処理手順を示すフローチャートである。
【図４】  フラクタル次元測定の処理手順を示すフローチャートである。
【図５】  輝度値が設定された解像度変換画像の生成過程の流れを示す図である。
【図６】  抽出画像の生成手順を示すフローチャートである。
【図７】  フラクタル次元測定の処理手順を示すフローチャートである。
【図８】  人間の肺のＸ線画像を示す図である。
【図９】  図８のＸ線画像を処理したフラクタル次元画像を示す図である。
【図１０】  フラクタル次元の統計的分布を示す図である。
【図１１】  フラクタル次元の統計的分布を示す図である。
【図１２】  多値画像のフラクタル次元を測定する画像処理装置の原理図である。
【図１３】  フラクタル次元測定の処理手順を示すフローチャートである。
【図１４】  フラクタル次元測定の処理手順を示すフローチャートである。
【図１５】  最大輝度値及び最小輝度値が設定された解像度変換画像の生成過程の流れを示す図である。
【図１６】  フラクタル次元の測定処理手順を示すフローチャートである。
【図１７】  データ構造の３次元座標を示す図である。
【図１８】  Ｘ−Ｙ平面の画像からＸ−Ｚ平面の画像を構成する様子を示す図である。
【図１９】  メモリスペースを節約する構造を示す図である
【図２０】  フラクタル次元の測定で複数の連続画像を取り扱う場合の概念を示す図である。
【図２１】  複数の連続画像を取り扱う場合のフラクタル次元測定の処理手順を示すフローチャートである。
【図２２】  複数の連続画像を取り扱う場合のフラクタル次元測定の処理手順を示すフローチャートである。
【図２３】  ＣＴスキャナの構成を示す図である。
【図２４】  ２次元の画像にボックスカウンティング法を適用した場合の図である。（Ａ）は２次元画像を示す。（Ｂ）は２次元画像の境界線を示す。（Ｃ）は２次元画像の境界線及び境界線を覆うボックスを示す。（Ｄ）は（Ｃ）のボックスを細分化したものである。
【図２５】  ３次元の画像にボックスカウンティング法を適用した場合の図である。
【符号の説明】
１  画像変換手段
２  抽出画像生成手段
３  解像度変換画像生成手段
４  輝度値設定手段
５  高輝度画素数計算手段
６  フラクタル次元測定手段
７  フラクタル次元画像表示手段
１０  画像処理装置
Ｉ  多値画像
Ｉａ  ２値画像
Ｉｂ  抽出画像
Ａ、Ｂ、Ｃ  解像度変換画像
Ｓ（Ａ）、Ｓ（Ｂ）、Ｓ（Ｃ）  高輝度画素数
【図１】
0003837921.tif 000002
【図２】
0003837921.tif 000003
【図３】
0003837921.tif 000004
【図４】
0003837921.tif 000005
【図５】
0003837921.tif 000006
【図６】
0003837921.tif 000007
【図７】
0003837921.tif 000008
【図８】
0003837921.tif 000009
【図９】
0003837921.tif 000010
【図１０】
0003837921.tif 000011
【図１１】
0003837921.tif 000012
【図１２】
0003837921.tif 000013
【図１３】
0003837921.tif 000014
【図１４】
0003837921.tif 000015
【図１５】
0003837921.tif 000016
【図１６】
0003837921.tif 000017
【図１７】
0003837921.tif 000018
【図１８】
0003837921.tif 000019
【図１９】
0003837921.tif 000020
【図２０】
0003837921.tif 000021
【図２１】
0003837921.tif 000022
【図２２】
0003837921.tif 000023
【図２３】
0003837921.tif 000024
【図２４】
0003837921.tif 000025
【図２５】
0003837921.tif 000026