プログラミングの備忘録

プログラムをつくる過程を残すもの

processingの備忘録 -マーブリング-

こんにちは。
今回は、processing で「マーブリング」をやってみます。


以前にこんなツイートを見ました。

流体シミュレータ的なものに興味があったので、こんなに単純化したモデルもあったのか、と感心していました。


そしてつい最近、我らが Shiffman 先生が Coding Challenge で似たものを取りあげてくださいました。


興味があってやり方もわかったとなれば、やるしかありません。


目次


マーブリングとは

マーブリングは、水面に絵の具を垂らしたり少しかき混ぜたりすることにより、いわゆる「マーブル模様」をつくり出す手法です。
日本では「墨流し」が有名ですね。
紙や布などに写しとって作品とされることが多いようです。

フルイドアートやラテアートなどもこの一種といえるのではないでしょうか。


実装

絵の具の滴下と押しのけ処理

まずは、クリックで円が描かれるようにしてみます。


ArrayList<Drop> drops = new ArrayList<Drop>();

void setup(){
  size(500, 500);
  noStroke();
}

void draw(){
  background(255);
  
  for(Drop d: drops){
    d.display();
  }
}

void mousePressed(){
  Drop drop = new Drop(mouseX, mouseY);
  drops.add(drop);
}

class Drop{
  PVector center; 
  float r;
  color col;
  int n;
  PVector[] edges;
  
  Drop(float _x, float _y){
    center = new PVector(_x, _y);
    r = 50;
    col = int(random(50, 200));
    
    n = 100;
    edges = new PVector[n];
    for(int i = 0; i < n; i ++){
      float theta = 2*PI*i/n;
      edges[i] = new PVector(center.x + r * cos(theta), 
                             center.y + r * sin(theta));
    }
  }
  
  void display(){
    fill(col);    
    beginShape();
      for(int i = 0; i < n; i ++){
        vertex(edges[i].x, edges[i].y);
      }
    endShape();
  }
}

垂らした絵の具を「drop」とよぶことにし、クラスをつくりました。
また、今後計算にあたって円の外周の点が必要なので「ある点から一定距離にある点の集まり edges」として円を表現しています。


これだけでは単に円が重なっていくだけのプログラムですが、ここに以下のような関係を導入することで絵の具が押しのけられるような挙動を表現できます。

$$ P' = C + (P - C) \times \sqrt{ 1 + \frac{ r^{2} }{ |P-C|^{2} } } $$

ここで、P'P は絵の具の外周の点の位置 (P は移動前、P' は移動後)、C は新たに垂らしている絵の具の中心の位置、r は新たに垂らしている絵の具の半径を指します。
ある色の絵の具の範囲を表すにはその外周の点がわかれば良いので、外周の点 P が新たに垂らされた絵の具によってどのくらい動くか、を計算すれば良いということになります。

$C$ に $P - C$ を足すと $P$ になりますが、既存の絵の具は新たに垂らされている絵の具によって押し広げられるので、その分を $\displaystyle \sqrt{ 1 + \frac{ r^{2} }{ |P-C|^{2} } }$ として $P - C$ にかけることにより $C$ を中心として若干遠くに動かしています。
既存の絵の具の移動量は、新たな絵の具が遠い (= $|P - C|$ が大きい) ほど小さく、新たな絵の具が多い (= $r$ が大きい) ほど大きいと考えられるので $\displaystyle \frac{ r }{ |P-C| }$ のという形となっていると思われます。


というわけで、この数式を実装します。

クラスに以下のメソッドを追加し、

void isPushedBy(Drop drop){
  for(PVector edge: edges){
    PVector c = drop.center;
    PVector p = edge;
    float r = drop.r;
      
    PVector s = PVector.sub(p, c);
    float m = s.mag();
    float root = sqrt(1 + (r*r)/(m*m));
      
    edge.set(PVector.add(c, s.mult(root)));
  }
}

呼び出しは mousePressed() 内で行うことにします。

void mousePressed(){
  Drop drop = new Drop(mouseX, mouseY);
  drops.add(drop);
  
  for(int i = 0; i < drops.size()-1; i ++){
    Drop other = drops.get(i);
    other.isPushedBy(drop);
  }
}

新たな絵の具 drop により既存の絵の具 other が押しのけられるので、isPushedBy() という名前にしてみました。
isPushedBy() 関数では先程挙げた数式から、edge の情報を更新しています。

実行するとこんな感じ。

押しのけられる挙動が見えるようになりました。


滴下処理をいじる

より実際のマーブリングらしく、長押しで絵の具の面積が大きくなっていくようにしてみます。

mousePressed() では押した瞬間のみの実行になってしまいますが、draw() 内で mousePressed とif文を組み合わせればマウスボタンが押されている間ずっと実行させることができるようになります。
そこで、mousePressed() には新たな絵の具の追加の処理を残したまま、mousePressed を使って押しのける処理を行うことにしました。


ArrayList<Drop> drops = new ArrayList<Drop>();

void setup(){
  size(500, 500);
  noStroke();
}

void draw(){
  background(255);
  
  for(Drop d: drops){
    d.display();
  }
  
  if(mousePressed == true){
    Drop drop = drops.get(drops.size()-1);
    
    for(int i = 0; i < drops.size(); i ++){
      Drop other = drops.get(i);
      other.isPushedBy(drop);
    }
  }
}

void mousePressed(){
  Drop drop = new Drop(mouseX, mouseY);
  drops.add(drop);
}

class Drop{
  PVector center; 
  float r;
  color col;
  int n;
  PVector[] edges;
  
  Drop(float _x, float _y){
    center = new PVector(_x, _y);
    r = 10;
    col = int(random(0, 255));
    
    n = 100;
    edges = new PVector[n];
    for(int i = 0; i < n; i ++){
      float theta = 2*PI*i/n;
      edges[i] = new PVector(center.x + r * cos(theta), 
                             center.y + r * sin(theta));
    }
  }
  
  void display(){
    fill(col);    
    beginShape();
      for(int i = 0; i < n; i ++){
        vertex(edges[i].x, edges[i].y);
      }
    endShape();
  }
  
  void isPushedBy(Drop drop){
    for(PVector edge: edges){
      PVector c = drop.center;
      PVector p = edge;
      float r = drop.r;
      
      PVector s = PVector.sub(p, c);
      float m = s.mag();
      float root = sqrt(1 + (r*r)/(m*m));
      
      edge.set(PVector.add(c, s.mult(root)));
    }
  }
}

処理的には、新たに追加された絵の具による押しのけをクリックしている間ずっと続ける、という形なので押しのけ処理を適用する絵の具は drops.size()-1 までにする必要があります。


実行してみると、


乱し処理

ここまででなんとなくできてきましたが、マーブリングでは乱す操作も重要な要素です。

直線

まずは直線状に乱します。なぞるような感じでしょうか。


軸に沿って

y軸に沿った方向になぞる場合、移動前の $P = (x, y)$ はなぞることにより $P' = (x, y + z u ^ {| x - x_{L} |})$ に移動します。
ここで、$u$ は双曲線に近似するための値だというようなことが書かれており (The Mathematics of Marbling)、$u = \left( \frac{1}{2} \right) ^ {\frac{1}{c}}$ で表されるようです。
また、$z$ は歪みの大きさ、$c$ は歪みの鋭さ、$x_{\mathrm{L}}$ はなぞる位置のx座標を表します。

これを適用すると、$z$ を係数として $x _{\mathrm{L}}$ に近い点は大きく、遠い点は小さく移動する、という感じになります。


コード上では、Dropクラスに tineLine() という関数を追加して処理を行います。

void tineLine(float xl, float z, float c){
  for(PVector edge: edges){
    float u = 1 / pow(2, 1/c);
    
    edge.x = edge.x;
    edge.y = edge.y + z*pow(u, abs(edge.x-xl));
  }
}

変数として $x_{\mathrm{L}}$、$z$、$c$ がありますので、引数として受け取るようにしました。
そして計算式に従い、円周上の点である egdes の座標を変更していきます。

呼び出しは draw() 内で、

if(keyPressed == true){
  for(Drop d: drops){
    d.tineLine(mouseX, 5, 10);
  }
}

として、キーを押したときに実行されるようにしてみました。

zc を変えることで、若干ですが見た目も変わると思います。


x軸に沿った場合は、$x$ と $y$ の関係が入れ替わるだけです。
つまり、$P = (x, y)$ が $P' = (x + z u ^ {| y - y_{L} |}, y)$ に移動します。


一般化

これまでの計算方法では、軸に沿った向きにしか歪みを加えられませんでした。
実際には棒か何かで任意の方向になぞっていく感じになると思いますので、上の式を一般化します。

といっても、やること自体は軸に沿った場合のものを任意の角度に回転させただけです。

つまり、

$$ P' = P + z u^{|(P - B) \cdot N|} M $$

$P$ に対して、移動方向に $z u^{|(P - B) \cdot N|}$ だけ動かす分のベクトルを足す、という計算になります。
ここで、大文字はベクトルだと考え、$N$ と ${M}$ はそれぞれの方向の単位ベクトルとします。


ちなみに、距離が $|(P - B) \cdot N|$ と表されていることについては、

と見たとき、赤い点線の長さは $|P - B| \mathrm{cos} ~ \theta$ と表されることから、$N$ は単位ベクトルなので $|P - B| |N| \mathrm{cos} ~ \theta$、内積の定義から $(P - B) \cdot N$ と変形できることから説明できます。


ではプログラムに組み込んでいきます。
先程つくった tineLine() を一般化したものにつくりかえれば良いでしょう。

void tineLine(PVector start, PVector direction, float z, float c){
  for(PVector edge: edges){
    PVector l = direction;
    PVector p = edge;
    PVector b = start;
    PVector n = l.copy().rotate(PI/2).normalize();
    PVector m = l.copy().normalize();
        
    float u = 1 / pow(2, 1/c);
    float d = abs(PVector.dot(PVector.sub(p, b), n));
        
    edge.set(PVector.add(p, m.mult(z*pow(u, d))));
  }
}

ベクトルの演算に関して、PVector.sub(p, b) のように呼び出した場合は問題ありません。
しかし、何か特定のベクトルのメソッドとして、p.sub(b) のように呼び出した場合は p の情報を変更しながら計算が行われてしまうので、予期せぬ挙動になってしまいます。
この問題は、p.copy().sub(b) のように copy() を介することで回避できます。

呼び出しは draw() 内で、キーを押したときに。

if(keyPressed == true){
  for(Drop d: drops){
    PVector s = new PVector(mouseX, mouseY);
    PVector v = new PVector(0, 1);
    d.tineLine(s, v, 5, 10);
  }
}

マウスの位置 s を始点として、方向 v に歪みを加えます。
(0, 1) なら、y軸方向に、(1, 1) なら45°の方向に、など指定できますし、うまくやればマウスでなぞった部分に適用するようにもできると思います。


序盤で $u$ は双曲線に近似するための値だというようなことを書きましたが、論文 (https://ieeexplore.ieee.org/document/5887299) の方には近似の形ではない式があります。

$$ P' = P + \frac{\alpha \lambda}{d + \lambda} M $$

ここで、$d$ は $|(P - B) \cdot N|$、$\alpha$ と $\lambda$ は $z$ や $c$ と同様、それぞれ歪みの大きさと鋭さを表すパラメータです。

$z$ ($\alpha$)、$c$ ($\lambda$) がともに 10 のとき、基準 $x_{\mathrm{L}}$ からの距離と $P' - P$ の大きさの関係を見てみると、近似あり (青線) が近似なし (橙線) に近い値となっていることがわかると思います。

というわけで、表現が違うだけでやりたいことは同じだと思いますので、どちらを使っても問題ないと思います。


渦の中心を $B$ として、渦から一定の距離にある $P$ は、

$$ P' = B + (P - B) \begin{pmatrix} \mathrm{cos} ~ a & -\mathrm{sin} ~ a \\ \mathrm{sin} ~ a & \mathrm{cos} ~ a \end{pmatrix} $$

で表される $P'$ に移動します。
ここで、$a$ は渦により $P$ が何度回転するかを表す値で、$\displaystyle a = \frac{zu^{||P - B| - r|}}{|P - B|}$ と表されます。

また $\begin{pmatrix} \mathrm{cos} ~ a & -\mathrm{sin} ~ a \\ \mathrm{sin} ~ a & \mathrm{cos} ~ a \end{pmatrix}$ は回転行列とよび、これが作用したベクトルは角度 $a$ だけ回転します。
(回転行列については、ベクトルの回転 ― Kinectで学ぶ数学 - Build Insider などに説明があります。)


$a$ について少し説明を加えます。

$P$ から $P'$ への移動量は、直線で歪めたときの関係式から、$zu^{d}$ と表すことができます。
$d$ はなぞった軌跡と $P$ の距離なので、渦の場合は $d = ||P - B| - r|$ となります。

また、軌跡を弧としたときの円の中心 $B$ と $P$ との距離は $|P - B|$ です。

弧度法から、$P$ が $B$ を中心とする円弧上を $zu^{d}$ だけ移動したときの中心角を $a$ とすると、$\displaystyle a = \frac{zu^{d}}{|P - B|}$ と表すことができるので、$\displaystyle a = \frac{zu^{||P - B| - r|}}{|P - B|}$ となります。


以上から、$B$ に $\overrightarrow{ BP }$ を角度 $a$ だけ回転させたベクトルを足すことにより、移動後の $P'$ を計算する、ということになります。


では、これをプログラムに組み込んでいきます。
Dropクラスに tineVortex() という関数を追加して処理を行うことにします。

回転行列はどう計算するのかというところですが、これは線形代数の基本である行列の積の計算を行えば良いです。
つまり、行列がかけられるベクトルを $V = (x, y)$ とすると、

$$ V \begin{pmatrix} \mathrm{cos} ~ a & -\mathrm{sin} ~ a \\ \mathrm{sin} ~ a & \mathrm{cos} ~ a \end{pmatrix} = (x \mathrm{cos} ~ a + y \mathrm{sin} ~ a, -x \mathrm{sin} ~ a + y \mathrm{cos} ~ a) $$

となります。

void tineVortex(PVector vcenter, float radius, float z, float c){
  for(PVector edge: edges){
    PVector p = edge;
    PVector b = vcenter;
    float r = radius;
      
    float u = 1 / pow(2, 1/c);
    PVector s = PVector.sub(p, b);
    float d = abs(s.copy().mag() - r);
    float a = z*pow(u, d) / s.copy().mag();
      
    edge.x = b.x + s.x*cos(a)+s.y*sin(a);
    edge.y = b.y - s.x*sin(a)+s.y*cos(a);
  }
}

軌跡の中心を vcenter、中心から軌跡までの距離を radius として受け取っています。
上で説明した計算をガリガリ行って、最後に edge の座標を更新して完了です。

これまでと同様、呼び出しは draw() 内で、キーを押したときに実行されるようにします。

if(keyPressed == true){
  for(Drop d: drops){
    PVector s = new PVector(mouseX, mouseY);
    d.tineVortex(s, 100, 1, 10);
  }
}


実行すると、

radius を 0 にすれば、The Mathematics of Marbling で Vortex として立項されているものになります。

混ぜすぎると描画がおかしくなってしまうので、描画方法は要検討かもしれません。


波型に乱すことによって、ある点 $P$ は、

$$ P' = P + A \mathrm{sin} ~ (\omega v + \phi) (\mathrm{cos} ~ t, \mathrm{sin} ~ t) $$

で表される $P'$ に移動します。
ここで、$A$ は振幅、$\omega$ は角振動数、$\phi$ は位相、$t$ は波の振動方向を表します。
また $v$ は $P \cdot (\mathrm{sin} ~ t, -\mathrm{cos} ~ t)$ という内積の計算結果が入ります。

$(\mathrm{sin} ~ t, -\mathrm{cos} ~ t)$ は $(\mathrm{cos} ~ t, \mathrm{sin} ~ t)$ を時計回りに90°回したベクトルに相当します。
このベクトルと $P$ との内積をとるので、波の振動方向と垂直な方向の成分を取り出すことができます。
これを変位としてサイン波を計算し、振動方向の単位ベクトルにかけることで波の形をつくりだします。
最後に、つくったベクトルを $P$ に足すことにより波型に変形させることができる、という感じだと思います。


では、これをプログラムに組み込んでいきます。
Dropクラスに tineWave() という関数を追加して処理を行うことにします。

void tineWave(float direction, float amplitude, float frequency, float phase){
  for(PVector edge: edges){
    PVector p = edge;
    float a = amplitude;
    float l = frequency;
    float o = phase;
    float t = direction;
      
    float dot = p.x*sin(t) - p.y*cos(t);
    float f = a*sin(l*dot+o);
      
    edge.x = p.x + f*cos(t);
    edge.y = p.y + f*sin(t);
  }
}

呼び出しは draw() 内で、キーを押したときに実行されるようにします。

if(keyPressed == true){
  for(Drop d: drops){
    d.tineWave(0, 1, 0.1, 0);
  }
}


実行すると、


まとめ

以上、processing でマーブリングを行えるようにしてみました。
このご時世、言語の壁は大きな障害にならないかもしれませんが、日本語版の解説記事的な感じで役立てていただければ幸いです。

乱し方として線分 (直線型をある程度の長さにとどめる) のものもあるようですが、ちらっと見た (https://people.csail.mit.edu/jaffer/Marbling/stroke.pdf) ところ難しそうだったので、追々ここに追記するか、新たに記事を書くかして触れてみようと思っています。
また、形が複雑になってくると描画に不満が出てきます。edges の点数を増やすと重くなってしまうので、別の方法で行う必要があるかなとも思っています。


新年度から割と忙しく、やっと時間が取れてリハビリがてらの更新のつもりがまあまあ長い記事になってしまいました。
最後まで読んでいただきありがとうございました。また次回。


参考