Poisson Image Editing (revised)

Tulisan ini melengkapi tulisan sebelumnya yang berjudul serupa. Waktu itu hasilnya masih kurang memuaskan dan sepertinya implementasinya masih ada yang salah. Akhirnya saya coba membuat implementasi ulang dengan Delphi, mulai dari pendefinisian persoalannya pada ruang berdimensi 1 secara diskret.

Sederhananya, misalkan Saya punya fungsi 1D yang terdefinisi oleh nilai di posisi awal f(a) dan akhir f(b) sementara nilai diantara selang (a,b) tidak diketahui. Kalau saya berasumsi bahwa turunan keduanya adalah konstan 0, maka interpolasi yang dibuat antara titik a dan b adalah sebuah garis lurus seperti gambar berikut.

interpolasi dengan turunan kedua = 0

interpolasi dengan turunan kedua = 0


Kalau fungsi turunan keduanya saya buat konstan positif atau negatif maka interpolasinya akan menjadi melengkung.

interpolasi dengan turunan kedua konstan positif ( = 2 )

interpolasi dengan turunan kedua konstan negatif ( = -2 )

Ide dari poisson image editing adalah dengan menggunakan informasi dari turunan kedua citra sumber untuk mengarahkan interpolasi ini pada citra tujuan.

interpolasi dengan turunan kedua pengarah berasal dari turunan kedua fungsi lain

Untuk menerapkan ide ini pada operasi pemindahan citra, maka perlu dibuat formulasi yang sesuai. Jika dilihat pada halaman wikipedia tersebut pada bagian example khususnya bagian formulasi untuk b pada tulisan saya sebelumnya agak berbeda. Pada tulisan wikipedia elemen dari b juga menambahkan kondisi batas (piksel di citra tujuan yang berbatasan/bertetangga dengan piksel yang ingin diinterpolasi) sedangkan pada tulisan saya sebelumnya tidak dimasukkan. Hal inilah yang menyebabkan pada gambar di tulisan saya sebelumnya seolah olah bagian batasnya masih terlihat gelap (karena intensitas di daerah batasnya dianggap 0 = hitam).

Ternyata setelah formulasinya saya sesuaikan, persamaan Ax = b bisa diselesaikan dengan cara langsung (spsolve, sebelumnya diselesaikan dengan cara relaksasi yaitu metode jacobi). Contohnya bisa dilihat dalam kode python berikut.

from scipy.ndimage.filters import laplace
from scipy.sparse import *
from scipy.sparse.linalg import *
from pyamg.relaxation import *

def poisson_blending(m, mask, src,iterations=5000,blendmode='paste'):
    h,w = mask.shape
    r,c = mask.nonzero()
    N = len(r)
    
    idx = np.zeros(mask.shape, dtype=np.uint32)
    for i in xrange(N):
        idx.itemset((r.item(i),c.item(i)), i+1)
    
    #laplace per channel
    slap = src.astype(np.float32)
    laplace(src[:h,:w,0], slap[:h,:w,0])
    laplace(src[:h,:w,1], slap[:h,:w,1])
    laplace(src[:h,:w,2], slap[:h,:w,2])
    
    mlap = m.astype(np.float32)
    laplace(m[:h,:w,0], mlap[:h,:w,0])
    laplace(m[:h,:w,1], mlap[:h,:w,1])
    laplace(m[:h,:w,2], mlap[:h,:w,2])
    
    #rhs per channel
    b_r = np.zeros((N, 1))
    b_g = np.zeros((N, 1))
    b_b = np.zeros((N, 1))
    
    #prepare for poisson matrix
    row = []
    col = []
    data = []
    for i in xrange(N):
        y,x = r.item(i),c.item(i)
        
        if blendmode=='blend':
            b_r.itemset(i, (-mlap.item((y,x,0))-slap.item((y,x,0)))*0.5)
            b_g.itemset(i, (-mlap.item((y,x,1))-slap.item((y,x,1)))*0.5)
            b_b.itemset(i, (-mlap.item((y,x,2))-slap.item((y,x,2)))*0.5)
        elif blendmode=='add':
            b_r.itemset(i, -mlap.item((y,x,0))-slap.item((y,x,0)))
            b_g.itemset(i, -mlap.item((y,x,1))-slap.item((y,x,1)))
            b_b.itemset(i, -mlap.item((y,x,2))-slap.item((y,x,2)))
        elif blendmode=='paste':
            b_r.itemset(i, -slap.item((y,x,0)))
            b_g.itemset(i, -slap.item((y,x,1)))
            b_b.itemset(i, -slap.item((y,x,2)))
        else: #smooth interpolation (laplace)
            b_r.itemset(i, 0)
            b_g.itemset(i, 0)
            b_b.itemset(i, 0)
            
        p = i
        if y>0: 
            if mask.item((y-1,x)):
                q = idx.item((y-1, x))-1
                row.append(p)
                col.append(q)
                data.append(-1.)
            else:
                b_r[i] += m.item((y-1,x,0))
                b_g[i] += m.item((y-1,x,1))
                b_b[i] += m.item((y-1,x,2))
        
        if x>0: 
            if mask.item((y,x-1)):
                q = idx.item((y, x-1))-1
                row.append(p)
                col.append(q)
                data.append(-1.)
            else:
                b_r[i] += m.item((y,x-1,0))
                b_g[i] += m.item((y,x-1,1))
                b_b[i] += m.item((y,x-1,2))
                
        if y<h-1: 
            if mask.item((y+1,x)):
                q = idx.item((y+1, x))-1
                row.append(p)
                col.append(q)
                data.append(-1.)
            else:
                b_r[i] += m.item((y+1,x,0))
                b_g[i] += m.item((y+1,x,1))
                b_b[i] += m.item((y+1,x,2))
                
        if x<w-1: 
            if mask.item((y,x+1)):
                q = idx.item((y, x+1))-1
                row.append(p)
                col.append(q)
                data.append(-1.)
            else:
                b_r[i] += m.item((y,x+1,0))
                b_g[i] += m.item((y,x+1,1))
                b_b[i] += m.item((y,x+1,2))
        row.append(p)
        col.append(p)
        data.append(4.)
        
    #reusable sparse matrix A
    A = csr_matrix((data, (row, col)), shape=(N,N))
    x = np.zeros((3,N))

    #rfn = sor
    #rfn = jacobi
    omega=1.0 #when omega=1.0 sor is equivalent to gauss-seidel
    n_iter = iterations
    #rfn(A, x[0,:], b_r, omega=omega, iterations=n_iter)
    #rfn(A, x[1,:], b_g, omega=omega, iterations=n_iter)
    #rfn(A, x[2,:], b_b, omega=omega, iterations=n_iter)
    
    #direct solver
    solver = spsolve
    x[0,:] = solver(A, b_r)
    x[1,:] = solver(A, b_g)
    x[2,:] = solver(A, b_b)
    
    #iterative solvers
    #solver = cg
    #solver = cgs
    #solver = gmres
    #solver = lgmres
    #solver = qmr
    #solver = minres #dark areas
    #solver = bicg
    #solver = bicgstab
    
    #least square
    #solver = lsqr
    #x[0,:] = solver(A, b_r)[0]
    #x[1,:] = solver(A, b_g)[0]
    #x[2,:] = solver(A, b_b)[0]
    
    
    x = x.transpose()
    for i in xrange(N):
        yy,xx = r.item(i), c.item(i)
        v = x[i]
        v[v<0] = 0
        v[v>255] = 255
        m[yy,xx] = v
    return m

Kode di atas melakukan penyelesaian persamaan poisson untuk citra berwarna pada masing-masing channel (merah, hijau, dan biru). Matriks A bersifat jarang (sparse) dan digunakan untuk tiap channel sehingga hanya perlu dibuat sekali. Kode di atas juga sengaja masih berisi bagian yang saya komentari untuk mencoba berbagai metode yang tersedia di numpy/scipy dan pyamg. Hasil akhir dari matriks kolom x perlu dibatasi pada rentang nilai [0,255] sesuai dengan representasi citra.

Selain itu kode di atas juga mengandung beberapa macam kombinasi yang bisa dilakukan dengan menggunakan turunan parsial kedua (laplacian) dari citra sumber.

        if blendmode=='blend':
            b_r.itemset(i, (-mlap.item((y,x,0))-slap.item((y,x,0)))*0.5)
            b_g.itemset(i, (-mlap.item((y,x,1))-slap.item((y,x,1)))*0.5)
            b_b.itemset(i, (-mlap.item((y,x,2))-slap.item((y,x,2)))*0.5)
        elif blendmode=='add':
            b_r.itemset(i, -mlap.item((y,x,0))-slap.item((y,x,0)))
            b_g.itemset(i, -mlap.item((y,x,1))-slap.item((y,x,1)))
            b_b.itemset(i, -mlap.item((y,x,2))-slap.item((y,x,2)))
        elif blendmode=='paste':
            b_r.itemset(i, -slap.item((y,x,0)))
            b_g.itemset(i, -slap.item((y,x,1)))
            b_b.itemset(i, -slap.item((y,x,2)))
        else: #smooth interpolation (laplace)
            b_r.itemset(i, 0)
            b_g.itemset(i, 0)
            b_b.itemset(i, 0)

Sebagai contoh, saya akan coba memindahkan wajah saya ke wajah Hawkeye (The Avengers)😛

citra sumber


citra mask


citra target


citra hasil

Di Tulisan berikutnya saya akan bahas bagaimana pemindahan wajah seperti pada contoh tersebut dilakukan secara otomatis, tunggu ya!😉

Tinggalkan Balasan

Isikan data di bawah atau klik salah satu ikon untuk log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout / Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout / Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout / Ubah )

Foto Google+

You are commenting using your Google+ account. Logout / Ubah )

Connecting to %s