MeanShiftFilteringをProcessingでアニメーションしたかも

月に1回は更新したかった。GANばっかも飽きるので別のこと書く
openCVにcv2.pyrMeanShiftFilteringっていう画像変換がある。何というか入力画像の色合いが均一になってのっぺりした感じのものが出力として返ってくる。

f:id:Owatank:20181029130442j:plain f:id:Owatank:20181029130356j:plain
左が入力画像で右がcv2.pyrMeanShiftFilteringを実行した出力結果

もうどこのサイトだったか忘れたけど、この手法を3次元の色空間にプロットして直感的に理解しやすい画像を前に見た。分散の気持ちになれて感動した。
せっかくなので自分もProcessingでMeanShiftFilteringを実装してアニメーションというか変化をプロットをしてようと思った。ちょっと失敗した

コード

github.com

デザインはこんな感じで作ってみる
f:id:Owatank:20181029131925p:plain

入出力画像を表示するのは簡単なんだけど、RGB値を3次元空間にプロットするための軸はどうしようか最初に悩んだ
P3Dというのがあるらしく、これを使えばキャンバス内で3D表現できるらしい。
このモードで面が塗られていない線だけの長方形を用意すれば3次元空間もどきが作れた。rotateY()使ってマウスでぐりぐり立方体を動かせるようにもできた。

f:id:Owatank:20181029133142p:plain

入力画像を (256, 256) サイズにリサイズして、ピクセル情報をこの立体内にプロットするとこんな感じになった
f:id:Owatank:20181029133547p:plain
雑だけどメインの部分じゃないからいいや

次は難しいMeanShiftFilteringの部分の実装。次のリンクとかを参考にした。

http://visitlab.jp/pdf/MeanShiftSegmentation.pdf
http://seiya-kumada.blogspot.com/2013/05/mean-shift-filtering.html

うーん。入力画像の各ピクセル $\boldsymbol{p_i} = (x_i, y_i, r_i, g_i, b_i)$ について、次の条件式を満たすようなピクセル $(x^{\prime}, y^{\prime}, r^{\prime}, g^{\prime}, b^{\prime})$ の集合をまずは求めていくのか

$|x^{\prime} - x_i| \leq h_s \ \ |y^{\prime} - y_i| \leq h_s , \ \ \sqrt{ {(r^{\prime} - r_i)}^2 + {(g^{\prime} - g_i)}^2 + {(b^{\prime} - b_i)}^2 \leq h_r}$

$h_s$ と $h_r$ はそれぞれ画像(の位置)の空間、色空間のカーネルサイズ。公式だと窓の半径とか書かれてた。これを指定しないと、あるピクセルの近傍 $h_s$ 内に存在してかつそのピクセルのRGB値の近傍 $h_r$ 内にも存在しているピクセルを集めることができないからかな。

よく見ると画像の位置の方の条件式と、色の条件式は式の形が違う。後者はユークリッド距離に見える
これって前に本で読んだこの部分を意味してるのかな
f:id:Owatank:20181029144948j:plain f:id:Owatank:20181029144921j:plain

前者はL1距離かなと思ったけど、それならXとYの和の式になるはずだから違うっぽい。絶対値で判断してるだけか
色の条件式で球のイラストが載ってたのは多分こういうこと・・なのかなあ・・・

次にこの条件式に当てはまったピクセルの集合に対して位置情報と色情報の重心を求める。ある場所に要素が密集してたらそっちに重心が偏るはず。これが出力画像がのっぺりする理由かな

位置情報と色情報の重心をベクトルで $(x^\ast , y^\ast , r^\ast, g^\ast, b^\ast )$ 表現したら参考文献によると $\boldsymbol{p_i} = (x_i, y_i, r_i, g_i, b_i)$ とこの重心が以下の条件のいずれかを満たせばいいらしい

$x_i = x^\ast \cap y_i = y^\ast$

$n \ge N$

$|x_i - x^\ast | + |y_i - y^\ast | + {(r_i - r^\ast )}^2 + {(g_i - g^\ast )}^2 + {(b_i - b^\ast)}^2 \leq \varepsilon$

3つ目の $\varepsilon$ は閾値を表す。思ったのが $|x_i - x^\ast | + |y_i - y^\ast |$ と ${(r_i - r^\ast )}^2 + {(g_i - g^\ast )}^2 + {(b_i - b^\ast)}^2$ の二つは単位が違くね？距離ってわけじゃ無いのかな。数学は難しい

2つ目の $n$ と $N$ は何かというと上の条件式のどれにも満たなかった場合、 $\boldsymbol{p_i} = (x_i, y_i, r_i, g_i, b_i)$ を $\boldsymbol{p_i} = (x^\ast, y^\ast, r^\ast, g^\ast, b^\ast)$ にして新たにまた重心を求めるステップから繰り返す。 $n$ と $N$ は反復回数と反復の上限数を表す。何か勾配となるようなものを手がかりに最適な値へと近づいていく感じか

どれかの条件式に当てはまったなら、出力ピクセルを $\boldsymbol{p'_i} = (x_i, y_i, r^\ast , g^\ast , b^\ast)$ として出力用の画像配列に代入する。これを入力の各ピクセルに対して行うとのっぺりした出力が得られる

次のように最初は実装した

for(int iter=0; iter < 5; iter++){
    center_p = get_center(filter_img, x_dash, y_dash, r_dash, g_dash, b_dash);
    condition = abs(x_dash - center_p.x) + abs(y_dash - center_p.y) 
    + pow(r_dash - center_p.r, 2) + pow(g_dash - center_p.g, 2) + pow(b_dash - center_p.b, 2);
    if(condition <= eps || (int(x_dash) == int(center_p.x) && int(y_dash) == int(center_p.y))){
      r = center_p.get_r();
      g = center_p.get_g();
      b = center_p.get_b();
      break;
    }
    x_dash = center_p.get_x();
    y_dash = center_p.get_y();
    r_dash = center_p.get_r();
    g_dash = center_p.get_g();
    b_dash = center_p.get_b();
}
output_color = color(r, g, b);

get_center()はこんな風にした

Pixel_vec get_center(PImage data, int x_dash, int y_dash, float r_dash, float g_dash, float b_dash){
  
  color c;
  float r = 0, g = 0, b = 0;
  float rgb_distance = 0;  
  
  float len = 0;  
  float sum_x = 0, sum_y = 0;
  float sum_r = 0, sum_g = 0, sum_b = 0;
  float center_x = 0, center_y = 0;
  float center_r = 0, center_g = 0, center_b = 0;
  
  Pixel_vec center_p;
  for(int x = 0; x < data.width; x++){
    for(int y = 0; y < data.height; y++){
      
        c = data.pixels[x*data.width+y];
        r = red(c);
        g = green(c);
        b = blue(c);

        rgb_distance = sqrt(pow(r-r_dash, 2) + pow(g-g_dash, 2) + pow(b-b_dash, 2));
        if(abs(x-x_dash) <= hs && abs(y-y_dash) <= hs && rgb_distance <= hr){
        len += 1;
        sum_x += x;
        sum_y += y;
        sum_r += r;
        sum_g += g;
        sum_b += b;
        
        center_x = sum_x / len;
        center_y = sum_y / len;
        center_r = sum_r / len;
        center_g = sum_g / len;
        center_b = sum_b / len;
        }
    }
  }
  center_p = new Pixel_vec(int(center_x), int(center_y), center_r, center_g, center_b);
  return center_p;
}

画像空間のカーネルサイズ hs = 16、色空間のカーネルサイズ hr = 16にして実行してみた
卑劣なので反復上限を超えても最適なピクセルが得られなかった場合はその位置の入力画像のピクセル情報を出力として代入するズルをした

f:id:Owatank:20181029130442j:plain f:id:Owatank:20181029164954p:plain

f:id:Owatank:20181029165055g:plain

画像だとなんか赤みが増したくらいにしか見えないけど、出力画像の各ピクセルのRGB値のプロットを見るとどことなく色の分散が小さくなっているのがわかる。でも冒頭に貼ったopenCVのcv2.pyrMeanShiftFiltering結果とは全然似てない。openCVの方は色味はあまり変化せずにのっぺりしてるし。

原因としては重心を求めるためのピクセルを集めるための条件式とか、重心を求める式の部分が違うのかな。微分の気持ちになれなくてここら辺の理解が足りなかった。またチャレンジしよう

そういえば、どっかで計算幾何学にいる人は数式で証明して終わりってだけじゃなくて、それをルールとしてコンピュータに与えられるところまでやる責任があるって本で読んだ気がする。まだまだわからんことのが多いけどかっこいいなあ