読者です 読者をやめる 読者になる 読者になる

tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

Poisson Image Editing

以前の記事で利用した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()

参考