Rでインボリュート曲線を波打たせて人の顔っぽいのを描く

タイトルどう書けばいいのかわからなかったんですが、要はこれをRでやりたいという話です。途中までやったのでメモ。

Mathematicaを使って，シュレーディンガーの顔をこのようなアニメーションにされたユーザの方がいらっしゃいます。コードも掲載されています。https://t.co/IDIzM8Xfy2 pic.twitter.com/ZTIbjtmXBm

— Wolfram Japan (@WolframJapan) 2017年6月15日

Mathematicaのコード読めないので適当に想像でやってみます。

インボリュート曲線を描く

こういう渦巻ってなんて言うんだったか思い出せずググったら、インボリュート曲線というらしいです。

インボリュート曲線 - Wikipedia

これをRで描くとこんな感じです。↑に載っている数式を、Rのコードに置き換えるだけです。

library(tidyverse)

steps <- 10000L
coil_turns <- 3L

d <- tibble(
  phi = 2 * pi * seq(0, steps) / steps * coil_turns,
  x = cos(phi) + phi * sin(phi),
  y = sin(phi) - phi * cos(phi)
)

ggplot(d) +
  geom_path(aes(x, y)) +
  theme_minimal() +
  coord_equal()

インボリュート曲線を回転させる

次に、これを回転させてみます。rotationという変数を追加して、ここに指定した値だけ時計回りに回転するようにします。

rotation <- pi/4

d2 <- tibble(
  radius = 2 * pi * seq(0, steps) / steps * coil_turns,
  phi = radius - rotation,
  x = cos(phi) + radius * sin(phi),
  y = sin(phi) - radius * cos(phi)
)

ggplot(d2) +
  geom_path(aes(x, y)) +
  theme_minimal() +
  coord_equal()

インボリュート曲線を波立たせる

次に、インボリュート曲線を波立たせます。より周期の短い波をつくって、それにsin()、cos()をそれぞれかけることで回転させます。インボリュート曲線の何分の一の周期にするかというcoil_wavinessというパラメータと、波立つ大きさwave_heightというパラメータを導入します。

coil_waviness <- 100
wave_height <- 1.5

d3 <- tibble(
  radius = 2 * pi * seq(0, steps) / steps * coil_turns,
  phi = radius - rotation,
  x = cos(phi) + (radius + wave_height * sin(phi * coil_waviness)) * sin(phi),
  y = sin(phi) - (radius + wave_height * sin(phi * coil_waviness)) * cos(phi)
)

ggplot(d3) +
  geom_path(aes(x, y)) +
  theme_minimal() +
  coord_equal()

インボリュート曲線に強弱をつける

これは、さっきのwave_heightに当たる部分を変化させればいいので、適当にsin()でつくったまた別の波をかけて強弱をつけてみます。

wave_height_base <- 1.5

d4 <- tibble(
  radius = 2 * pi * seq(0, steps) / steps * coil_turns,
  phi = radius - rotation,
  wave_height = wave_height_base * sin(phi * 1.76),
  x = cos(phi) + (radius + wave_height * sin(phi * coil_waviness)) * sin(phi),
  y = sin(phi) - (radius + wave_height * sin(phi * coil_waviness)) * cos(phi)
)

ggplot(d4) +
  geom_path(aes(x, y)) +
  theme_minimal() +
  coord_equal()

そう、これを画像の濃淡にあわせてやればいいわけです。

画像の濃淡に合わせてインボリュート曲線に強弱をつける

画像を取得

まず、httrパッケージを使って適当な画像を取ってきて、content()で読み込みましょう。

library(httr)
res <- GET("http://hoxo-m.com/img/team/makiyama.png")
img <- content(res)
img
#> , , 1
#> 
#>              [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
#>   [1,] 0.81960784 0.82352941 0.81960784 0.82352941 0.81568627 0.82745098
#>   [2,] 0.81960784 0.81960784 0.81960784 0.81960784 0.82352941 0.81960784
#>   [3,] 0.81960784 0.81960784 0.81960784 0.81960784 0.81960784 0.81960784
#>   [4,] 0.81960784 0.82352941 0.82352941 0.82352941 0.82352941 0.82352941
#> ...

dim(img)
#> [1] 138 138   3

これは3次元の配列になっていて、3つ目のインデックスがRGBそれぞれに対応しているらしいです。（content()がPNG画像をパースするときに使うのはpngパッケージのreadPNG()ですが、これは色情報によって結果の形は異なります。今回はめんどくさいので他の形式は想定しないことにします）

グレースケールに変換

ここを参考にしました。

色々あるみたいですが、NTSC方式というやつがRGBそれぞれに係数をかけて足して割るだけで、計算しやすそうだったのでこの簡易版にします。

img_grey <- (2 * img[,,1] + 4 * img[,,2] + img[,,3]) / 7

ただし、これは1に近いほど白いことを示します。知りたいのは黒さなので、1 - img_grayを計算します。さらに、細かい波が出ないように、小さい値は足切りするなど適当に加工します。今回は2乗して、ゼロに近い値がより小さくなるようにしてみます。

img_grey <- (1 - img_grey) ^ 2

インボリュート曲線の各点に対応する画像のピクセルを見つける

このあと、この画像とインボリュート曲線を重ねて、各点が重なっているピクセルの黒さに応じて波立たせ方を強くします。このために、まずはxとyの範囲が画像のピクセル数より少し小さい範囲に収まるように調整します。これで、対応する画像中のピクセルのインデックスを求めたことになります。（ここでxとyが逆になるあたりとかいまいち理解できてません。。）

dim(img_grey)
#> [1] 138 138

max_radius <- max(d2$radius)
pixels_x <- dim(img_grey)[1] - 1
pixels_y <- dim(img_grey)[2] - 1

d5_ <- d2 %>%
  mutate(
    x_index = ceiling(pixels_x - (1 + x / max_radius) / 2 * pixels_x),
    y_index = ceiling(pixels_y - (1 + y / max_radius) / 2 * pixels_y)
  )

インボリュート曲線を画像の濃さに合わせて波立たせる

インデックスが分かったので、対応する画像のピクセルを抜き出し、その値を波にかけ合わせます。

注意すべき点は、行列に複数のインデックスを指定してベクトルとして要素を抜き出したいときは、インデックスも行列にしなければいけないという点です。インデックスにベクトルを別々に指定すると、行列として抜き出されてしまいます。

m <- matrix(1:9, ncol = 3)

m[c(1, 1, 2), c(1, 2, 1)]
#>      [,1] [,2] [,3]
#> [1,]    1    4    1
#> [2,]    1    4    1
#> [3,]    2    5    2

こうならないようにcbind()を使います。

m[cbind(c(1, 1, 2), c(1, 2, 1))]
#> [1] 1 4 2

つまり、こうです。

d5 <- d5_ %>%
  mutate(
    blackness = img_grey[cbind(y_index, x_index)],
    x = x + blackness * wave_height * sin(phi * coil_waviness) * sin(phi),
    y = y - blackness * wave_height * sin(phi * coil_waviness) * cos(phi)
  )

ggplot(d5) +
  geom_path(aes(x, y)) +
  theme_minimal() +
  coord_equal()

ちょっとわかりづらいですね。もうちょっと巻き数を増やしてみるとこんな感じです。

# せっかくなので関数にする
draw_wavy_involute <- function(steps = 10000L,
                               coil_turns = 3L,
                               rotation = pi/2,
                               coil_waviness = 100,
                               wave_height = 1.5,
                               img_grey = NULL) {
  pixels_x <- dim(img_grey)[1] - 1
  pixels_y <- dim(img_grey)[2] - 1
  
  d <- tibble(
    radius = 2 * pi * seq(0, steps) / steps * coil_turns,
    phi = radius - rotation,
    x = cos(phi) + radius * sin(phi),
    y = sin(phi) - radius * cos(phi)
  )
  
  max_radius <- max(d$radius)
  
  d <- d %>%
    mutate(
      x_index = ceiling(pixels_x - (1 + x / max_radius) / 2 * pixels_x),
      y_index = ceiling(pixels_y - (1 + y / max_radius) / 2 * pixels_y)
    )
  
  d <- d %>%
    mutate(
      blackness = img_grey[cbind(y_index, x_index)],
      x = x + blackness * wave_height * sin(phi * coil_waviness) * cos(phi),
      y = y - blackness * wave_height * sin(phi * coil_waviness) * sin(phi)
    )
  
  ggplot(d) +
    geom_path(aes(x, y)) +
    theme_minimal() +
    coord_equal()
}

draw_wavy_involute(coil_turns = 30, wave_height = 6, img_grey = img_grey)

ちょっとそれっぽくなりました。

TODO

• アニメーション
• 画像の拡大縮小
• ノイズの除去