以前の記事で利用したPoisson Surface Reconstructionのアルゴリズムを調べる前に2次元でポアソン方程式を解く手法について調べたところ、Poisson Image Editingを見つけた。この手法を3次元に拡張すると3次元メッシュも変形できる。
概要
2つの画像を境界値を指定したポアソン方程式を解くことで合成する方法。画像を人が見るときに注目される特徴として勾配があり、これを保存することが重要。f*をターゲットとなる画像、fを埋め込まれたあとの画像として、境界で$f|_{\partial\Omega} = f^{\ast}|_{\partial\Omega}$を満たしながら$\int\!\!\int_{\Omega} |\nabla f - v|^{2}$を最小とするfを求めればいい。これは$\triangle f = \mathrm{div}v$に等しい。証明はオイラー・ラグランジュ方程式をつかって行われる(参考pdf)。vを変化させることで様々な画像の合成に対応できる。
離散化
$\frac{\partial^{2} f}{\partial x^{2}} \simeq (f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1}-4f_{i,j})/{\Delta x}^{2}$を利用して、$|N_{p}|f_{p}-\sum_{q\in N_{p}\cap\Omega}f_{q} = \sum_{q \in N_{p}\cap\partial\Omega}f^{\ast}_{q} + \sum_{q\in N_{p}}v_{pq}$と書ける。Ax = bの形の行列なのでそのまま解くこともできるが、$f_{p}^{(k+1)} = (\sum_{q\in N_{p}\cap\Omega}f_{q}^{(k)} + \sum_{q \in N_{p}\cap\partial\Omega}f^{\ast}_{q} + \sum_{q\in N_{p}}v_{pq}) / |N_{p}|$と変形してガウス＝ザイデル法やヤコビ法で高速に解くことができる。
- 参考
実装(バグ取れてないしまともに動かない)
あと少しいじったら動くはず。判断ミスだったのはヤコビ法をpythonで実装してしまったこと。pythonで書くなら実行速度も考えてpoisson solverライブラリを拾って来るべき。
# coding: utf-8 import numpy as np from PIL import Image import sys from time import time import matplotlib.pyplot as plt # ヤコビ法の収束許容範囲 EPS = 0.0001 def __seamless_cloning(target, mask, source, v): #im = Image.fromarray(mask.astype('uint8')*255) #Image.merge('RGB', (im, im, im)).save("mask_clip.png") # target、mask、sourceの形は揃えてから渡す w, h = mask.shape # 初期値 #f = source[mask == 1] f = np.empty(np.sum(mask)) + 10 f = f.flatten() # ベクトル化できるためヤコビ法を利用する # |Np|==4と仮定して、4f(n+1) = Af(n) + bとなるようにA、bを作る # TODO: Aを疎行列に置き換える A = np.zeros((f.size, f.size)) pos = (np.cumsum(mask) - 1).reshape(mask.shape) * mask print(w,h) for i in range(1, w - 1): for j in range(1, h - 1): if mask[i,j] == 1: A[pos[i,j], pos[i + 1, j]] = mask[i + 1, j] A[pos[i,j], pos[i - 1, j]] = mask[i - 1, j] A[pos[i,j], pos[i, j + 1]] = mask[i, j + 1] A[pos[i,j], pos[i, j - 1]] = mask[i, j - 1] b = np.zeros(f.size) for i in range(1, w - 1): for j in range(1, h - 1): if mask[i,j] == 1: # mask == 0のとき近傍が境界に含まれるのでfq*は0でない fq_ast = np.sum([target[m, n] * (1 - mask[m, n]) for m, n in ((i + 1, j), (i - 1,j), (i,j + 1), (i,j - 1))]) b[pos[i, j]] == fq_ast - v[i,j] b_out = np.zeros(mask.shape) for i in range(w): for j in range(h): if mask[i,j] == 1: b_out[i,j] = f[pos[i,j]] #plt.matshow(b_out) #plt.savefig(str(time())+".png") # ヤコビ法 while (True): f, f_old = (np.dot(A,f) + b) / 4, f print(f) error = np.sum((f - f_old)**2) #/ np.sum(f**2) # out = np.zeros(mask.shape) for i in range(w): for j in range(h): if mask[i,j] == 1: out[i,j] = f[pos[i,j]] #plt.matshow(out) #plt.savefig(str(time())+".png") #Image.fromarray(out).convert('RGB').show() # print(error) if error < EPS: break out = np.array(target) for i in range(w): for j in range(h): if mask[i,j] == 1: out[i,j] = f[pos[i,j]] return out def clip_where(target, mask, source, offset): # mask.shape == source.shapeを前提とする t = (max(0, offset[0]-1), min(target.shape[0], mask.shape[0] + offset[0]+1), max(0, offset[1]-1), min(target.shape[1], mask.shape[1] + offset[1]+1)) s = (max(0, -offset[0]-1), min(mask.shape[0], target.shape[0] - offset[0]+1), max(0, -offset[1]-1), min(mask.shape[1], target.shape[1] - offset[1]+1)) return t, s def seamless_cloning(target, mask, source, offset = (0, 0)): t, s = clip_where(target, mask, source, offset) ds_x = source[1:,1:] - source[:-1,1:] ds_y = source[1:,1:] - source[1:,:-1] div_s = (ds_x + ds_y)[s[0]:s[1], s[2]:s[3]] #Image.fromarray(div_s).show() ct, cm, cs = (target[t[0]:t[1], t[2]:t[3]], mask[s[0]:s[1], s[2]:s[3]], source[s[0]:s[1], s[2]:s[3]]) out = np.array(target) for i in range(3): blend = __seamless_cloning(ct[:,:,i], cm, cs[:,:,i], v = div_s) out[t[0]:t[1], t[2]:t[3], i] = blend plt.matshow(blend) plt.show() Image.fromarray(out).show() return out def mixed_seamless_cloning(target, mask, source, offset = (0, 0)): ds_x = source[1:,1:] - source[:-1,1:] ds_y = source[1:,1:] - source[1:,:-1] div_s = (ds_x + ds_y)[s[0]:s[1], s[2]:s[3]] dt_x = target[1:,1:] - target[:-1,1:] dt_y = target[1:,1:] - target[1:,:-1] div_t = (dt_x + dt_y)[t[0]:t[1], t[2]:t[3]] v = np.array([[max(div_s[i,j], div_t[i,j]) for j in range(div_s.shape[0])] for i in range(div_s.shape[1])]) ct, cm, cs = (target[t[0]:t[1], t[2]:t[3]], mask[s[0]:s[1], s[2]:s[3]], source[s[0]:s[1], s[2]:s[3]]) out = np.array(target) for i in range(3): blend = __seamless_cloning(target[t[0]:t[1], t[2]:t[3]], mask[s[0]:s[1], s[2]:s[3]], source[s[0]:s[1], s[2]:s[3]], v) out[t[0]:t[1], t[2]:t[3]] = blend return __seamless_cloning(target, mask, source, v) def local_illumination_change(): raise NotImplementedError def local_color_change(): raise NotImplementedError def seamless_tiling(): raise NotImplementedError def main(): target = np.asarray(Image.open("test/target.png", 'r')) mask = np.asarray(Image.open("test/mask.png", 'r')) mask = (mask != 0) * 1 source = np.asarray(Image.open("test/source.png", 'r')) seamless = seamless_cloning(target, mask, source, offset = (40, -30)) Image.fromarray(seamless).save("test/seamless.png") #mixed = mixed_seamless_cloning(target, mask, source, offset = (40, -30)) #Image.fromarray(mixed).save("test/mixed.png") if __name__ == "__main__": main()
