Rustでネコチャンを点描する

この記事はLivesense Advent Calendar 2021の15日目の記事です。

先日こんなツイートを見かけました。

https://twitter.com/artixels/status/1439909994179137537

ネコチャン！　カワイイ！

いやぁ、可愛いですね。これを見たとき、ぜひ自分でも実装してみたいと思いました。なんといってもネコチャンの柔らかい輪郭とゆるゆる動く点描の相性がとても良いので。

それで、この論文を見つけたわけですがこちらを見る限り、鍵となるアイデア自体はそこまで複雑なものではないようで、これならできるかも？　とか思っていたわけです。

もともとのアルゴリズムはこちらのリファレンス実装にある通り、GPUを利用して計算されるもので、ちょっと大変なのですが、今回はGPUを使わずにエッセンスを再現してみようと思います。そのために動画の処理は諦め、静止画の変換処理を近似的に再現します。

（一応付言しておくと、上述の論文をヒントに似ているけれど全然別のアルゴリズムを実装したというべきかもしれません。そういうものだと思って読んでください）

手法

さて、提案されているのは電荷を持った粒子のシミュレーションに基づくサンプリングです。どんな手法なのか、簡単に説明しましょう。

平面全体に電荷を持った粒子を分布させます。変換元の画像をグレースケールに変換した上で、各ピクセルにも電荷を付与します。このとき黒の強度が高いほど強く粒子を引きつけるようにピクセルの電荷を設定します。

画像の黒い部分が粒子を引きつけつつ、粒子同士は互いに反発し合うことで適度な濃淡を表現する効果を生み出すというわけです。

式で書くと

粒子$q_i$は各粒子からは、以下の力をうけ

$$F_i = q_s^2\sum_{j\neq i}^{N}\frac{1}{\left|r_i - r_j\right|^2}\hat{e}_{j,i}$$

画像のピクセルからは、各ピクセルの電荷を$q_k$として、

$$G_i = q_s^2\sum_{k=1}^{M}\frac{q_k}{\left|r_i - r_j\right|^2}\hat{e}_{k,i}$$

の力を受けます。($r_i$などは位置ベクトル、$\hat{e}_{j,i}$は方向を表す単位ベクトルです)

上記２つの式をもとに、各粒子の振る舞いを運動方程式を使ってシミュレーションすれば、点描が得られるというわけです（なお論文では粒子が分布する平面と元画像の平面の「距離」を調整することで、複数の画像を組み合わせたり、画像の特定の部分にフォーカスを当てたりする手法が提案されていますが、この記事では扱いません）

実装

さて先述の通り、リファレンス実装ではGPUを用いて全粒子×（全粒子＋全ピクセル）の作用を計算しているのですが、粒子数は（推奨が）8192、ピクセル数はちょっとした画像でも数万になるので、全く同じ計算をCPUでシミュレートしたくはありません。

それで直感的には、力は距離の二乗に反比例するので、画像のピクセルにしろ、粒子同士の作用にしろ距離の近いところだけ計算しても同じような効果が得られそうです。なので今回は、近傍の粒子とピクセルが与える力のみを考えてシミュレーションすることにしましょう。

実装にはRustおよびRust製のジェネラティブコーディング用のライブラリであるnannouを使います。

まずは画像から

まずは、画像の効果のみ実装してみます。

論文にならって、以下の計算式で各ピクセルをグレースケール、そして電荷に変換します。

let charge = BLANK_LEVEL - (0.2989 * p[0] as f32 + 0.5870 * p[1] as f32 + 0.1140 * p[2] as f32) / 255.0;

p[0]、p[1]、p[2]がそれぞれRGBに対応しています。BLANK_LEVELはパラメタの一つで、この数値を超えるピクセルは符号が反転し、粒子を弾くようになります（大体0.95くらいの数値に設定しておき、画像の白い領域に粒子が集まらないようにする効果を与えます）

ピクセルが粒子に与える力は次のとおりです

https://github.com/yskaksk/rust-stippling/blob/main/image/image_only.rs#L145-L154

let q = Q_CHARGE / d2;
let pq = model.image[[xi as usize, yi as usize]];
if d.abs() > eps {
    fx -= pq * q * dx / d;
    fy -= pq * q * dy / d;
}

fx、fyが粒子にかかる力のx成分、y成分ですね。d2はxi、yiにより指定されるピクセルと粒子の距離、Q_CHARGEは粒子がもつ電荷です。

まとめると

1. 粒子の位置から対応するピクセルの位置を計算し
2. 前後左右の予め決めた範囲のピクセルを順番に走査し
3. 粒子にかかる力を計算する

という処理で粒子にかかる力を計算しているわけです。
粒子の質量はすべて同じで単位量ということにすれば、この力はそのまま加速度とみなせます。

あとは、この加速度を使って、速度と位置をそれぞれ更新するだけです。

https://github.com/yskaksk/rust-stippling/blob/main/image/image_only.rs#L157-L176

// 速度の差分
let mut dvx = dt * fx;
let mut dvy = dt * fy;

// 速度更新
pi.vx = 0.95 * pi.vx + dvx;
pi.vy = 0.95 * pi.vy + dvy;

// 位置の差分
let mut dpx = dt * pi.vx + 0.5 * dt.pow(2) as f32 * fx;
let mut dpy = dt * pi.vy + 0.5 * dt.pow(2) as f32 * fy;

// ランダムなゆらぎを追加
pi.x += random_range(-0.005, 0.005);
pi.y += random_range(-0.005, 0.005);

// 位置を更新
pi.x = rotate(pi.x + dpx);        
pi.y = rotate(pi.y + dpy);

画像の効果だけだと、粒子が広がらないのでランダムなゆらぎを入れてみました(後で消します)
これだけでも、単純な図形ならそれなりに再現できます。

以下はランダムな初期状態から46という数字がおぼろげながら浮かび上がる様子です。

粒子同士の相互作用

続いて粒子同士の相互作用を考えましょう。
画像のピクセルは不動なので、近傍を計算するには単に対応するピクセルの位置を求めればよかったのですが、粒子の位置は常に更新されるので、近傍の計算には少し工夫が必要です。

工夫として、kd木を使います。kd木の詳しい説明は割愛しますが、粒子の位置をkd木に格納すれば、すくない計算量でk近傍が求められるので、今の目的にはぴったりです。

rustでkd木を扱うためのcrateはいくつかあるようですが、今回はこちらを使わせてもらいます。

粒子同士の相互作用は次のように計算します。

https://github.com/yskaksk/rust-stippling/blob/main/image-and-particle/image_and_particles.rs#L183-L193

for pj in model.tree.nearests(&[pi.x, pi.y], N_NEIGHBOR).iter() {
    let dx = rotate(pi.x - pj.item[0]);
    let dy = rotate(pi.y - pj.item[1]);
    let d2 = dx.powf(2.0) + dy.powf(2.0) + 0.00003;
    let d = d2.powf(0.5);
    let q = Q_CHARGE / d2;
    if d.abs() > eps {
        fx += Q_CHARGE * q * dx / d;
        fy += Q_CHARGE * q * dy / d;
    }
}

粒子同士は反発し合うので、先程と符号が逆になっていますね。model.treeは粒子の位置を格納しているkd木です。位置の更新が終わったらkd木も更新します。

let tree = KdTree::build_by_ordered_float(model.points.iter().map(|p| [p.x, p.y]).collect());
model.tree = tree;

また、先程との差分としてピクセルのもつ電荷の調整があります。リファレンス実装では（全粒子の電荷の和）=（全ピクセルの電荷の和）となるように調整されていますが、ここでは参照する近傍の粒子数とピクセル数をもとにバランスが取れるように調整します。

https://github.com/yskaksk/rust-stippling/blob/main/image-and-particle/image_and_particles.rs#L108-L113

// 参照するピクセルの電荷の和
let area_r = (PIC_R + 1).pow(2) as f32 / (w * h) as f32;
let mut bmp_charge = bmp_charge_total * area_r;

// 近傍の粒子の電荷の和
let neighbor_dot_charge = N_NEIGHBOR as f32 * Q_CHARGE;

// 調整係数
bmp_charge = neighbor_dot_charge / bmp_charge;

// グレースケールから計算した行列に調整係数をかける
let img_mul = img_mat * bmp_charge;

速度や位置の更新は先程と同じです。

結果

こちらの画像で試してみました。

結果はこの通り

うーん。かわいい。目が丸いところとかが、特に。
点描でもかわいいのがネコチャンのすごいところです。

神奈川沖浪裏

悪くないですね。

うまくいかなかった画像

いろいろパラメータを変えたけどうまくいかなかったものもあります。

こちらのネコチャンは

クリーチャーになってしまいました

胴体部分の濃淡はそんなに悪くないですが、本来目を表すべき粒子が顔に吸われてしまっています。どうも明るい狭い領域が暗い色に囲まれていると、うまく表現できないようです。

なぜそうなるかですが、私のアルゴリズムでは明るい領域（目の領域）にある粒子は囲まれる暗い領域（顔の領域）の一部のみから力を受けるためだろうと推測されます。本来のアルゴリズムならぐるりと囲まれる顔領域から均等に引き寄せられるので、ある程度は目の領域に粒子がとどまり続けるはずです。

この問題はより外側のピクセルの力をどうにかして反映させることで解消できるかもしれません。もちろんピクセルを参照する領域を広げればよいのですが、参照するピクセル数の２乗に比例して計算量が増えます。

計算量を（それほど）増やさない方法としては、例えば、画像を10x10のグリッドに切った上で、仮想的な引き寄せる粒子として画像グリッド内の加重平均をあらかじめ計算しておき、外側の領域から受ける力はその仮想粒子からの力で代替する、というようなものが考えられます。うまく機能するかどうか、いずれ試してみたいですね。

最後に

どんな画像でもうまくいくというわけではありませんし、機能するものも手動のパラメータ調整がちょっとめんどうではあるのですが、雑な近似としては良く機能しているのではないでしょうか。

今後はゆるゆる動くネコチャンの表現にも挑戦してみたいと思いました。

実装はここにまとめています。よければどうぞ。