Huモーメント不変量による類似画像検索

Googleが"画像で検索する"機能をリリースして、暫く経つ。他の人はどうか知らないけど、自分の場合、実際に使うときって、"同じ"画像を探すのにしか使わない。類似性というのは、色んな視点がありえるので、正しい答えというのもなく、そういう意味で難しい問題であるけど、同じかどうかは人間ならほぼ100%判定できる。Googleは、この機能を実現するのにDeep何とかいう手法を使って沢山のマシンで学習させたとか、どっかに書いてたけど、同じ画像を探すという目的に限定すると、もっと手軽な手法がありそうなもんではある。


同じと言っても、フォーマットが違ったり、サイズが違うと、拡大縮小やら圧縮・符号化アルゴリズムによって、データとしては違うものになってしまうので、普通のハッシュ関数に通しても仕方ない。OpenCVのcv::matchTemplateのようなものを使うという手もあるけど、サイズを合わせないといけないので、2つのデータの類似度が対称にならないし、計算量も多い。というわけで、この問題もそれほど自明ではないけど、もし、画像データから、拡大縮小(スケール変換)、並進・回転の操作によって不変な量をいくつか計算することができれば、類似度/距離の計算を簡単化・高速化できるし、データ量も大幅に減らすことができる(可能性がある)


丁度いいことに、画像のスケール不変なモーメントの多項式として、並進・回転不変な7つのHuモーメント不変量というものがよく知られていて、OpenCVにも実装されているので、簡単に使える。
Image moment
http://en.wikipedia.org/wiki/Image_moment


#上Wikipediaの記事によると、Huのoriginalのセットは、独立でない元を含むばかりか、3次までの独立な元を全て含んでないとある。OpenCVの実装は、現時点では、Huのoriginalのセットに基づいているっぽい
HuMoments
http://docs.opencv.org/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html#humoments


Huモーメントから、データ間の類似度を測る方法は、特に自然なやり方とかはないけど、cv::matchShapesとかで使われてるやり方を参考にすればいいんだと思う(ただ、ここのやり方は高次モーメントの重みがやや高くなる傾向にあって、それほどよくないかもしれない)
matchShapes
http://docs.opencv.org/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html#matchshapes


この方法では、同じ画像(拡大縮小したり、平行移動したり回転したもの)に対して、同じ不変量の値を取ることが保証されているけど、逆に、全然違う画像がたまたま近い不変量の値を取ってしまう可能性は排除されない。つまり、大規模なデータに対して、この方法で検索すると、再現率は保証できるけど、適合率はどうなるか分からない。適合率の低下を回避するには、更に多くの独立な不変量を持ってくる(Jan Flusserの論文では、4次までの独立な不変量が計算されている)とか、あるいは、別の方法(例えば、cv::matchTemplateとか)でad hocにフィルタリングしていくとか、色々やり方は考えられる。


一応、そもそも「(無限に)沢山の不変量を持ってきても原理的にすら区別不可能な"異なる"画像というのはありえるのか?」という問いがある。究極的には、定数だってスケール変換・回転・並進の不変量であるわけで、不変量を全部持ってきても、識別不能な画像が大量にあるようでは、これは使えない方法ということになってしまう。それで、まず問題となるのは、画像のモーメント(無限個ある)は、どれくらい画像を区別できるのかという話であると思う。わたしは正しい答えを知らないけど、異なる画像で全てのモーメントが等しくなることは滅多にないだろうとは思う(厳密には、モーメント問題とかいうテーマに属する問題な気がする。一次元の有限区間に対しては解かれているはず)。なので、十分に多くの不変量を持ってくれば、全然違う画像が原理的にすら区別できないという最悪の事態は避けられるのだと思う。まぁ、細かいことをいうと、不変量が、全ての"軌道"を区別可能であるとは限らないとか、色々あるし、全然違う画像が、全ての不変量で近い値を持ってしまうことがないかは別問題ではある


この方法は、3次元形状やボクセルデータにも拡張可能。動画にも使えるかどうか分からないけど、原理的には使えるはず。まぁ、理屈は純数学的なので、簡単で分かりやすい

一応、簡単なテスト

import cv2
from math import log

def hudist(f1 , f2):
   def h2m(x):
       if x==0:
          return 0.0
       elif x>0:
          return 1.0/log(x)
       else:
          return -1.0/log(-x)
   def hum(f):
       img = cv2.imread(f , cv2.CV_LOAD_IMAGE_GRAYSCALE)
       assert(img!=None),("failed to open %s" % f)
       moments = cv2.moments(img)
       return cv2.HuMoments(moments)
   hu = hum(f1)
   hu2 = hum(f2)
   d = sum([abs(h2m(x1)-h2m(x2)) for (x1,x2) in zip(hu.T[0],hu2.T[0])])
   return d

if __name__=="__main__":
   imgs = ['Lenna.png','Lena3.jpg','sendo1.jpg','Sendou2.jpg','monnalisa.jpg','apple.jpg','Red_Apple.jpg']
   for f1 in imgs:
       print f1,
       for f2 in imgs:
           v = 100*hudist(f1,f2)
           print ("%5.2f" % v),
       print ""

結果

Lenna.png Lena3.jpg sendo1.jpg Sendou2.jpg monnalisa.jpg apple.jpg Red_Apple.jpg
Lenna.png 0.00 0.04 5.93 5.87 4.19 2.43 4.83
Lena3.jpg 0.04 0.00 5.93 5.86 4.22 2.41 4.83
sendo1.jpg 5.93 5.93 0.00 0.09 7.82 4.55 1.19
Sendou2.jpg 5.87 5.86 0.09 0.00 7.86 4.47 1.11
monnalisa.jpg 4.19 4.22 7.82 7.86 0.00 3.97 7.88
apple.jpg 2.43 2.41 4.55 4.47 3.97 0.00 4.94
Red_Apple.jpg 4.83 4.83 1.19 1.11 7.88 4.94 0.00

テストに使用した画像は以下の通り。まぁ、そこそこいけそうな感じではある。








おまけ:Huモーメント不変量が代数的に独立でないことは、グレブナー基底を使って簡単に確認できる。具体的には、以下のようなRisa/Asirコードを書く

load("gr");

def relation(Exprs , Vars , NewVars){
  /* assert(length(Exprs)==length(NewVars)); */
  S = [];
  for(I = 0 ; I < length(Exprs) ; I++)S = cons(Exprs[I] - NewVars[I] , S);
  Basis = nd_gr(S , append(Vars , NewVars) , 0 , [[0, length(Vars)] , [1,length(NewVars)]]);
  Com = [];
  for(I = 0 ; I < length(Basis) ; I++){
    Check = 1;
    for(J= 0 ; J < length(Vars) ; J++){
       if( deg(Basis[I] , Vars[J]) > 0 ){
          Check = 0;
          break;
       }
    }
    if(Check == 1)Com = cons(Basis[I] , Com);
  }
  return nd_gr(Com , NewVars , 0 , 0);
}

/* Huモーメント不変量 */
I1=m20+m02;
I2=(m20-m02)^2+4*m11^2;
I3=(m30-3*m12)^2+(3*m21-m03)^2;
I4=(m30+m12)^2+(m21+m03)^2;
I5=(m30-3*m12)*(m30+m12)*((m30+m12)^2-3*(m21+m03)^2)+(3*m21-m03)*(m21+m03)*(3*(m30+m12)^2-(m21+m03)^2);
I6=(m20-m02)*((m30+m12)^2-(m21+m03)^2)+4*m11*(m30+m12)*(m21+m03);
I7=(3*m21-m03)*(m30+m12)*((m30+m12)^2-3*(m21+m03)^2)-(m30-3*m12)*(m21+m03)*(3*(m30+m12)^2-(m21+m03)^2);

relation([I1,I2,I3,I4,I5,I6,I7],[m11,m20,m02,m12,m21,m30,m03],[p1,p2,p3,p4,p5,p6,p7]);

結果は、[-p4^3*p3+p5^2+p7^2]とかなって、I5^2+I7^2=I4^3*I3であることが分かる。これは、Flusserの論文
On the independence of rotation moment invariants
http://library.utia.cas.cz/prace/20000033.pdf
に書いてある結果(Flusserは、よりad hocな手法で、これを発見している)。ついでに、スケール不変なモーメントは無限個あるけども、有限次数で打ち切れば(例えば、Huモーメント不変量は三次までしか見てない)、それはただの(合同変換群による)不変式の計算に帰着して、やっぱりグレブナー基底を使って、独立な生成元を機械的に計算することができる(はず)。


合同変換群による不変式を計算するというのは、本質的には
三角形の合同条件と不変式論
http://d.hatena.ne.jp/m-a-o/20130104#p1
でやったのと同じようなもの(これ書いた当時は、不変式を計算するアルゴリズムを知らなかったけど、Derksenのアルゴリズムを使えば計算できる)。というわけで、これはComputer Vision分野におけるグレブナー基底のちょっと面白い応用事例でもあると思う