Julia Advent Calendar 2020/12/01

Julia好きの読者の皆様、おはようございます。堀川です。

突然ですが、私の名刺をご覧ください。

実はこの名刺、中に遊星歯車機構が入っていて、回るようになっています。(すごいでしょう！)

こちらの動画は、ストップモーションと呼ばれる手法を使って2年前に製作したものです。ストップモーションというのは、被写体やカメラ等の位置を変えて撮影し、また被写体やカメラ等の位置を変えて撮影し、またまた被写体やカメラ等の位置を変えて撮影し、、という作業を繰り返して動画を作る手法のことです。

さきほどの動画では、歯車の角度を少しずつ変えて、合計108枚の画像を撮影していました。

動画の問題点

2年前のストップモーションの動画ではちゃんとカメラや名刺の位置を調整せずに撮影していました。
そのため、このようなカクついた動画になっていたのでした。
「あとで手動でトリミングすればええやろ！」と思ってたら面倒になったケース…。

ふたたび画像を108枚撮り直すの面倒だし、画像処理でいい感じにスムーズに作り直せないかなあ…。

今回の成果物

そこで、Juliaでいい感じに画像処理してカクつきを減らしました！
それがこちらです↓

以降では今回の画像処理について解説していきます。

画像処理やっていき

方針

以下の方針でやっていけそうですね。

画像読み込み
エッジ・コーナー検出
コーナー情報をもとに画像を変形
色の調整
連番で画像を保存
連番画像から動画生成

① 画像読み込み

とりあえず画像を読み込みましょう。話はそれからです。

ここで必要なパッケージはImages.jlです。あらかじめusing Imagesで読み込んでおいてください。

# Load image
img = load("source_images/005.jpg")
h,w = size(img)

実行すれば分かると思いますが、こちらのjpg画像はArray{RGB{N0f8},2}として読み込まれましたね。つまり、画像とは2次元配列のことで、その要素はRGB{N0f8}です。N0f8は8bit固定小数点数ですね。

後の処理のためにグレースケールに変換しておきましょう。Gray.(img)はimgの各要素に関数Grayを適用するという意味です。GrayはColorのsubtypeであり、ここでは各要素のRGB(<:Color)をGrayへ変換しているという訳です。

# Grayscale
img_gray = Gray.(img)

② エッジ・コーナー検出

エッジ検出のためのフィルターは何種類かありますが、ここではLoGフィルターを使いました。
ImageFiltering.jlパッケージが必要です。

"""
Return a value between Gray(0.0) and Gray(1.0).
"""
function clp(c)
    return clamp(c,Gray(0.0),Gray(1.0))
end

# LoG filter
img_LoG = clp.(imfilter(img_gray, Kernel.LoG(3))*50)

ここで関数clpは戻り値をGray(0.0)からGray(1.0)に制限するための関数です。

上辺と左辺はいい感じにエッジが検出できそうですね。ただ、右辺に関してはあまり鮮明には得られませんでした。
(これは撮影時の影の影響ですが、正しくは照明を工夫すべきでしたね…。)

ちょっと分かりにくいですが、「影」で発生する白色より、「名刺右辺」で発生する白色の方がパキッとしてる感じに見えますね。この性質を使えば、別のフィルタリングで「名刺右辺」のエッジを強調できそうな気がします。

上下方向のモーションフィルタを適用してみましょうか。

# Vertical filter
n = 200
kernel_vrt = zeros(2n+1,2n+1)
kernel_vrt[:,n-1:n+1] .= 1/n
img_vrt = clp.(imfilter(img_LoG, kernel_vrt))

名刺の右辺がこれで強調できましたね！

同様に左右方向のモーションフィルタを適用しましょう。

# Horizontal filter
n = 200
kernel_hrz = zeros(2n+1,2n+1)
kernel_hrz[n-1:n+1,:] .= 1/n
img_hrz = clp.(imfilter(img_LoG, kernel_hrz))

じゃあimg_vrtとimg_hrzを掛けるとどうでしょうか..？

# Multiply
img_mul = img_hrz.*img_vrt

そうです！
上下の白線と左右の白線の交点が強調されるようになります。ピクセルの黒色は0で白色は1なので、乗算すると交点のみが1として現れて強調される訳ですね！

四隅の200x200を切り出して確認してみましょう。

# Corner image
corner = 200
img_corner00 = img_mul[end-corner+1:end,1:corner]
img_corner10 = img_mul[end-corner+1:end,end-corner+1:end]
img_corner01 = img_mul[1:corner,1:corner]
img_corner11 = img_mul[1:corner,end-corner+1:end]
mosaicview([img_corner01,img_corner00,img_corner11,img_corner11], npad=20, nrow=2, fillvalue=colorant"white")

やはりコーナーの頂点部分が一番明るいですね。

これらの200x200領域内で最も輝度の高い場所をコーナーとして推定して良いでしょう。配列の最大値とそのインデックスを取得するにはfindmaxを使います。

# Get coordinates of corners
corner_index00 = findmax(img_corner00)[2]+CartesianIndex(h-corner,0)
corner_index10 = findmax(img_corner10)[2]+CartesianIndex(h-corner,w-corner)
corner_index01 = findmax(img_corner01)[2]+CartesianIndex(0,0)
corner_index11 = findmax(img_corner11)[2]+CartesianIndex(0,w-corner)
corner_coordinates00 = [corner_index00[1],corner_index00[2]]
corner_coordinates10 = [corner_index10[1],corner_index10[2]]
corner_coordinates01 = [corner_index01[1],corner_index01[2]]
corner_coordinates11 = [corner_index11[1],corner_index11[2]]

…実は上記の「最も輝度の高い場所をコーナーとして推定して良い」というのは少し間違った主張です。詳しくは後述します。

③ コーナー情報をもとに変形

やりたいことは、こういうことです。
コーナーとして得られた4点をいい感じに動かして… $\textcolor{blue}{青四角形}$ を $\textcolor{red}{赤長方形}$ に重ねる。

では、どのような変換が使えるでしょうか？

平行性が崩れる変形が必要なので、Affine変換は使えない。
しかし、変換の前後で直線は保たれて欲しい。

ここで射影変換ですよ！
射影変換とは同次座標(定数倍を無視する同値関係を入れたもの)を使って以下のように定義される変換のことです。

\begin{aligned} \begin{pmatrix} x' \\ y' \\ 1 \end{pmatrix} \sim \begin{pmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{pmatrix} \begin{pmatrix} x \\ y \\ 1 \end{pmatrix} \end{aligned}

あるいは以下でも同じ定義です。

\begin{aligned} \begin{pmatrix} x' \\ y' \end{pmatrix} =\frac{\begin{pmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \end{pmatrix} \begin{pmatrix} x \\ y \\ 1 \end{pmatrix}}{\begin{pmatrix} h_{31} & h_{32} & h_{33} \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}} \end{aligned}

しかしながら、JuliaのパッケージImageTransformations.jlの周辺には射影変換は見つかりませんでした。(探した日付: 2020/11/22)

なので以下のように実装しました。

# Definition of projective transformation
"""
Projective transformation.
"""
struct Projective
    H::SArray{Tuple{3,3},Float64,2,9}
end

"""
Definition of projective transformation.
"""
function (p::Projective)(x)
    x_ = SA[x[1],x[2],1.0]
    x′_1, x′_2, x′_3 = p.H*x_
    return SA[x′_1/x′_3, x′_2/x′_3]
end

"""
Inverse of projective transformation.
"""
function Base.inv(p::Projective)
    Projective(inv(p.H))
end

"""
Generate matrix for projective transformation,
which satisfy (0,0)↦p00, (1,0)↦p10, (0,1)↦p01, (1,1)↦p11.
"""
function _projectivematrix(p00,p01,p10,p11)
    k1,k2 = hcat(p11-p10,p11-p01)\(p10+p01-p00-p11)
    v1 = (k1+1)*p10-p00
    v2 = (k2+1)*p01-p00
    h = hcat(v1,v2,p00)
    H = vcat(h,[k1,k2,1]')
    return H
end

"""
Generate projective transformation,
which satisfy (0,0)↦p00, (1,0)↦p10, (0,1)↦p01, (1,1)↦p11.
"""
function Projective(p00,p01,p10,p11)
    H = _projectivematrix(p00,p01,p10,p11)
    p = Projective(H)
    return p
end

"""
Generate projective transformation,
which satisfy p00↦q00, p10↦q10, p01↦q01, p11↦q11.
"""
function Projective(p00,p01,p10,p11,q00,q01,q10,q11)
    Hp = _projectivematrix(p00,p01,p10,p11)
    Hq = _projectivematrix(q00,q01,q10,q11)
    H = Hq*inv(Hp)
    p = Projective(H)
    return p
end

詳細な数学的解説はここでは避けようと思います。重要なことは、最後の関数Projective(p00,p01,p10,p11,q00,q01,q10,q11)は任意の4点 $(\bm{p}_{00},\bm{p}_{01},\bm{p}_{10},\bm{p}_{11})$ を別の任意の4点 $(\bm{q}_{00},\bm{q}_{01},\bm{q}_{10},\bm{q}_{11})$ に移す射影変換(のインスタンス)を生成する関数だということです。

これを使って名刺の4隅を長方形の頂点に移動させましょう。このような画像の変換にはwarp関数を使います。この関数を使うにはImageTransformations.jlパッケージが必要になります。ピクセル単位の補間とかもいい感じにやってくれます。

# Transformation
xmin = 92
xmax = 2950
ymin = 90
ymax = 1810

mark00 = [ymax,xmin]
mark10 = [ymax,xmax]
mark01 = [ymin,xmin]
mark11 = [ymin,xmax]

p = Projective(mark00,mark10,mark01,mark11,corner_coordinates00,corner_coordinates10,corner_coordinates01,corner_coordinates11)
img_positioned = warp(img, p)

出てきた画像はOffsetArrayになってますね。原点の場所を覚えるよう配慮してくれてるのでOffsetArrayが使われているという訳です。変形後に画像の当たらなかった場所(つまり背景)は黒(RGB{N0f8}(0.0,0.0,0.0))になってます。

トリミングしましょう。

# Crop
img_cropped = img_positioned[ymin-40:ymax+40,xmin-40:xmax+40]

これで位置合わせが完了しました！

④ サイズ・色の調整

さて、トリミング後の画像は少し大きすぎ(2939x1801)ますね。適当に1/5くらいに小さくしましょう。

# Resize
img_resized = imresize(img_cropped, ratio=1/5)

実は、照明やカメラ撮影時の自動露光調整の影響などで、画像ごとに明るさが違う問題もありました。

ここでは001番の画像の明るさに合わせることにしましょう。

色味の調節については深く考えるのはやめて、平均値を合わせるように適当に加減することにしました。meanを使うにはStatisticsパッケージが必要です。

img_target = load("source_images/001.jpg")
meancolor = mean(img_target)
img_coloraligned = img_resized.-mean(img_resized).+meancolor

ちょっと明るくなりましたね。

⑤ 連番で画像を保存

これまでやってきた画像処理をすべての写真(001.jpg,....,108.jpg)に対して適用しましょう。

for i in 1:108
    img = load("source_images/$(lpad(i,3,'0')).jpg")
    # (画像処理いろいろ)
    save("output_images/$(lpad(i,3,'0')).png",img_coloraligned)
end

ここでlpadは左側に文字を埋めてくれる関数です。つまりlpad(i,3,'0')だと左側に0を埋めるようにしてiを3桁の数字にします。例えばi==42なら"042"の文字列を返してくれます。

⑥ 連番画像から動画生成

ffmpeg氏に頑張ってもらいましょう。

ffmpeg -r 15 -i output_images/%03d.png -vcodec libx264 -r 30 namecard.mp4

あるいはgif動画を作るためにimagemagick氏にお願いしましょう。

convert -delay 10 output_images/*.png namecard.gif

おわりに

今回のコードのリポジトリこちらです↓

感想

画像処理を真面目(?)にやったのは初めてでしたが、ImageFiltering.jlやImageTransformations.jlが優秀だったので思ったより苦労せずにできました。他の言語、例えばpythonのcv2だと画像は3次元配列として読み込まれ、さらにBGRの順番なのでとても使いにくかった記憶です。Images.jl (が内部で使っているColors.jl) パッケージではRGBもBGRも用意されていますが、インターフェースが共通化されいるので使いやすいです。さらに「ラスタ画像とは色(Colorant)を2次元的に並べたもの(Array{<:Colorant,2})である」というのはかなり合理的だと思います。

エッジ検出の他のアルゴリズムなど

また、今回はLoGフイルタ+モーションフィルタで頂点を検出しましたが、これには「名刺の回転方向がある程度揃っていること」が必要でした。もっと回転(>10°くらい?)している図形に対しても処理したい場合は、Hough変換などを使う必要があります。

LoGフィルタについて

今回の実装ではLoG(Laplacian of Gaussian)フィルタで最も輝度の高い場所をもとに、コーナー位置を推定しました。しかし、Laplacianフィルタでの最大値は正確にはエッジ位置に一致していません。これはLaplacianが2階微分を参照しているためです。なので、より正確にエッジ位置を検出するためにはLaplacianの結果がゼロになる場所(ゼロ交差)を探す必要があります。今回は「ストップモーションの動画を滑らかに動かすこと」だけが目的だったので、エッジ位置がずれていたとしても、すべてのフレームで同程度ずれるなら問題ないと判断してゼロ交差を求めませんでした。

参考文献など

画像処理100本ノック
- 画像処理はこれで入門しています。
- いまJuliaで解いている途中です。
- https://yoyoyo-yo.github.io/Gasyori100knock/
ディジタル画像処理
- 網羅的に書かれているので辞書的に見ています。
- 射影変換・フィルタ処理の辺りを参考にしました。
- https://www.cgarts.or.jp/book/img_engineer/
Fablab北加賀屋
- 名刺はここのレーザー加工機で作りました。
- 会員になれば色々な機械が使えるのでとても有り難いです。
- https://fablabkitakagaya.org/

明日(2020/12/02)のJuliaアドベントカレンダーは、、、空いてました！記事誰か書いてくれると嬉しいです！！

Discussion