hirax.net::Keywords::「多重解像度解析」のブログ

■「ウォーリーを探す出す」多重解像度解析コードを書いてみる!?　

　「ウォーリーを探し出す」Mathematica コードがとても面白く・参考になりました。そして、もしも「ウォーリーを探し出す」コードを書こうとしたならば、自分ならどんな風に書くだろう？と考えました。…そこで、試しに、書いてみることにしました。

　まずは、画像を読み込みます。

waldo=Import["hoge"]

　そして、ウォーリーを見つけるには、「やはり、赤白シャツを頼りにすべし」というわけで、赤白模様を検出するための単チャンネルを作成します。ここで、2{1,-1,-1}.#&は、R-G-Bを計算する純（無名）関数です（Stack Over Flow のコードと同処理です。カスタマイズできるように純関数にしています）。

red=ImageApply[2{1,-1,-1}.#&,waldo]

　次は、ウォーリーを探す部分です。Stack Over Flow で書かれていたコードは「横方向の赤白線の上エッジ検出」を行い、「強いエッジがある部分にウォーリーがいる」という推定を行います。
　しかし、その推定では、「ウォーリーが着ているシャツは、赤白横線が何本も入っている」とか「その赤白横線は、胴の形に添った形状をしている」といった情報は用いられていません。
　そこで、「胴の形をした赤白横線模様」を作り、さらにその「模様」を用いた多重解像度解析（ウォーリーの大きさは未知ですから）をすることによりウォーリーの位置を推定する、というコードを書いてみました。

mask=Fold[ImageAdd,#[[1]],Rest[#]]&@Table[
h=a;w=Floor[h/2];
corr = ImageApply[ 4 Abs[# - .5] &,
ImageCorrelate[red,
Image@ Table[Cos@y, {y, 1, 6 Pi, 6 Pi/h}, {h}],
NormalizedSquaredEuclideanDistance ]];
ImageApply[ #/(30 - 5) &,
Dilation[ImageApply[ If[# > .8, 1, 0] &, corr],
ConstantArray[1, 4 {h, w}]] ] , {a, 5, 30}]

　このコードでは、多重解像度解析を行った上で、ありとあらゆる解像度、すなわち、ありとあらゆる大きさのウォーリーがいる（だろう）位置を、一枚のマスク画像へと焼き付けます。

　そして、焼き付けた（推定した）「ウォーリーが立っている位置のマスク画像」を用いて、原画像からウォーリーを浮き上がらせてみます。それが、以下の（Stack Over Flowそのままの）最後のコードです。

ImageMultiply[waldo,
ImageAdd[ColorConvert[mask, "GrayLevel"], .5]]

　試しに、「ウォーリーをさがすとはどういうことか」で使われていたサンプル画像に処理をかけてみると、次のようになります。これが、「ウォーリーを探す出す」多重解像度解析コードです!?

　

---

　赤白模様の「壁」などにマッチしないようにしたい場合のために、下記のようなフィルタリングを加えると良いかもしれません（今回はコードを短くするために割愛しました）。

corr2 = ImageApply[ 4 Abs[# - .5] &,
ImageCorrelate[corr,
Image@ Table[
If[x^2 + y^2 < h, 1, 0], {y, -h, h}, {x, -h, h}],
NormalizedSquaredEuclideanDistance ] ];

■Python/OpenCVで画像多重解像度解析コードを書いてみる　

　多重解像度解析…といっても直交基底に分解するというような話ではなくて、単に各周波数帯の特性がどの程度含まれるかを眺めるといった用途なら（つまり、ガボール変換やSTFTを掛ける感じの程度の用途なら）、Python/OpenCVを使って十数行で書けるかも？と思い書いてみました。もちろん、実装は簡単第一最優先！というわけで、ガウシアンフィルタ差分で２次元のバンドパスを作成し、それを周波数軸で重ねて眺めてみるというくらいの話です。

　実際に書いてみたら、ポスト処理含めて約２０行くらいになりました。超入門的な画像処理コードですが、１次元〜２次元の多重解像度解析や周波数解析を行うことは意外に多いような気もするので、適当に貼り付けておくことにします。＊

---

＊画像処理クラスタからのコメント：
　・マルチスケールで眺めるなら、DCゲイン1同士のガウシアン差分をとり、そのL2ノルムを1に正規化しすべし。
　・周波数軸は等比的にした上で、ボリューム的表示も等比的比率で重ねたい。

import numpy as np
import cv2
from matplotlib import pyplot as plt

def DOG(img, s, r):
    img2=img.astype('uint16')
    img2=img2*128+32767
    gs = cv2.GaussianBlur(img2,(0, 0), s)
    gl = cv2.GaussianBlur(img2,(0, 0), s*r)
    return cv2.absdiff(gs, gl)

img = cv2.imread("sample2.jpg",0 )
(h,w)=img.shape
pts1 = np.float32([[0,0],[w,0],[w,h],[0,h]])
pts2 = np.float32([[0,h*1/4],[w*3/4,h*1/4],
                   [w,h*3/4],[w*1/4,h*3/4]])
M = cv2.getPerspectiveTransform(pts1,pts2)
baseImg = cv2.warpPerspective(
    img.astype('uint16'),M,(w,h))

for i in range(100,5,-10):
    pts1 = np.float32([[0,0],[w,0],[w,h],[0,h]])
    pts2 = np.float32([[0+i,h*1/4-i],
    [w*3/4+i,h*1/4-i],[w+i,h*3/4-i],[w*1/4+i,h*3/4-i]])
    M = cv2.getPerspectiveTransform(pts1,pts2)
    img2 = cv2.warpPerspective(DOG(img, i, 1.05),M,(w,h))
    baseImg = cv2.addWeighted(baseImg, 0.9, img2, 0.3, 0)

plt.figure(figsize=(6,6))
plt.imshow(np.array(baseImg)) 
plt.autoscale(False)