Seamless Cloning using Poisson Image Blending in Python

04Apr12

Beberapa waktu lalu ketika saya sedang mempersiapkan tulisan sebelumnya tentang face tracking, Saya terdampar di sebuah situs berbahasa Jepang yang cukup membuat penasaran (sampai sekarang saya belum ketemu lagi dengan situs tsb). Judul yang saya ingat waktu itu adalah “Poisson Image Editing”. Kalau dari judulnya saja saya tidak terlalu tertarik, tapi saya justru tertarik karena melihat gambar yang ditampilkan.

image cloning

pemindahan citra kepala beruang ke citra bukit dengan panduan mask. (gambar bukan milik saya)

Singkat kata, “Poisson Image Editing” ini merupakan judul makalah yang ditulis oleh Perez et al. Isinya kurang lebih sebagai memaparkan cara menggabungkan dua buah citra (sumber + mask + target = hasil). mask merupakan citra biner yang menandakan bagian citra sumber yang boleh dipindahkan ke citra hasil. Contohnya terlihat pada gambar di atas. Pada gambar tersebut terdapat 3 citra yaitu citra sumber, mask, dan citra hasil. citra target adalah citra hasil yang sebelumnya tidak ada gambar kepala beruang yang merupakan bagian dari citra sumber yang dipindahkan menggunakan panduan dari citra mask. citra mask dapat dianggap sebagai cetakan atau penyaring bagian yang dipindahkan dan yang tidak.

Proses pemindahan pada umumnya menggunakan mask sebagai pemandu dengan cara menganggap citra berwarna putih bernilai True (piksel yang bersesuaian di citra sumber akan dipindahkan) atau False. Pada proses matting citra mask dapat berperan sebagai alpha mask yaitu citra mask dapat berisi piksel bernilai [0..1] dan piksel pada citra mask digunakan sebagai faktor kombinasi linier untuk menginterpolasi nilai piksel pada citra sumber dengan nilai piksel pada citra target. Biasanya citra mask ini merupakan hasil segmentasi (biner) dan kemudian diperhalus baik secara manual maupun otomatis (misal. menggunakan fasilitas Feather pada program pengolah citra). Untuk mendapatkan hasil penggabungan yang baik maka citra mask diasumsikan sudah sesuai dengan objek yang mau dipindahkan sehingga citra mask tidak lagi mencakup piksel yang merupakan objek yang tidak diinginkan untuk dipindahkan. Oleh sebab itu peran inspeksi manual tidak jarang menjadi diperlukan. Kelebihan dari penggabungan dengan menggunakan teknik Poisson Blending ini adalah area mask tidak harus melingkupi objek yang ingin dipindahkan saja tetapi bisa juga sedikit berlebih akibat penyeleksian yang tidak presisi. Setidaknya itulah yang diklaim di makalah yang saya baca walaupun setelah saya coba ternyata hasilnya tidak demikian (bisa jadi kode yang saya implementasi pun masih belum benar).

Dasar Teori Poisson Image Editing

Nama Poisson pada teknik ini berasal dari persamaan diferensial yang disebut persamaan Poisson (ini berbeda dengan fungsi distribusi peluang Poisson). Persamaan Poisson adalah persamaan diferensial parsial orde dua. Persamaan Poisson pada kasus pemindahan citra didefinisikan sebagai persoalan nilai batas (boundary value problem) dengan kondisi pembatas Dirichlet (nilai fungsi di daerah perbatasan diketahui).

Sederhananya, teknik pemindahan Poisson ini mencari nilai piksel-piksel yang dilingkupi oleh citra mask pada citra target/hasil dengan mengandalkan informasi berupa nilai piksel di perbatasan dan citra laplacian dari citra sumber.

Nilai piksel pada citra hasil/tujuan yang dilingkupi mask dianggap tidak diketahui. Informasi nilai piksel di perbatasan merupakan piksel-piksel pada citra tujuan yang tidak dilingkupi oleh citra mask tetapi bertetangga dengan piksel yang dilingkupi mask. Citra laplacian merupakan citra hasil penerapan filter laplacian atau turunan parsial kedua dari citra sumber.

Persamaan diferensial ini diselesaikan dengan metode numerik dengan cara mentransformasi persamaan Poisson menjadi diskret dan dianggap sebagai sistem persamaan linier Ax=b. Matriks A adalah operator laplacian. Vektor x merupakan nilai-nilai piksel yang dilingkupi oleh mask. Vektor b adalah nilai laplacian pada citra sumber (b = Ay; y adalah piksel piksel di citra sumber yang juga dilingkupi oleh mask). Matriks A merupakan matriks jarang (sparse) yang seluruh elemen diagonalnya bernilai positif dan sisanya bernilai -1 atau 0 (kebanyakan bernilai 0). Perez menganjurkan penyelesaian sistem persamaan linier ini dengan menggunakan metode Gauss-seidel atau SOR (Successive over relaxation). Namun karena saya sendiri tidak sempat mempelajari ini sewaktu kuliah metode numerik dulu, akhirnya saya menggunakan metode Jacobi yang menurut saya paling sederhana untuk diimplementasi.

Implementasi

Saya membuat tiga macam implementasi dari teknik ini. Implementasi pertama hanya dilakukan menggunakan python murni yang berjalan sangat lambat karena ada banyak pengulangan. Implementasi kedua dipercepat dengan menggunakan inline C++ melalui scipy.weave yang berjalan paling cepat (diantara dua lainnya). Implementasi terakhir saya buat karena penasaran dengan hasil yang tidak sesuai harapan. Saya mencari library lain yang mengimplementasi metode SOR/Gauss-Seidel seperti yang disarankan oleh Perez. Akhirnya saya menemukan libaray pyamg (Python Algebraic Multi Grid) yang di dalamnya mengandung berbagai metode relaksasi sistem persamaan linier Ax=b termasuk jacobi, sor, dan gauss_seidel.

Implementasi dengan Python I : murni python

def poisson_blending(tgt, mask, src, iterations=5):
    h,w = mask.shape
    dg = src.astype(float32)
    laplace(src[:h,:w,0], dg[:h,:w,0])
    laplace(src[:h,:w,1], dg[:h,:w,1])
    laplace(src[:h,:w,2], dg[:h,:w,2])
    r,c = np.nonzero(mask)
    cnt = len(r)

    g1 = tgt.astype(float)
    for iter in xrange(iterations):
        for i in xrange(cnt):
            y,x = r.item(i), c.item(i)
            Np = 0
            val = np.zeros(3)
            if y>0:# and mask.item((y-1,x)):
                val += g1[y-1, x]
                Np += 1
            if x>0:# and mask.item((y,x-1)):
                val += g1[y, x-1]
                Np += 1
            if y<h:# and mask.item((y+1,x)):
                val += g1[y+1, x]
                Np += 1
            if x<w:# and mask.item((y,x+1)):
                val += g1[y, x+1]
                Np += 1
            val -= dg[y,x]
            val /= Np
            val[val<0] = 0
            val[val>255] = 255
            g1[y,x] = val
    return g1.astype(uint8)

Implementasi pada Python II : Inline C++ dengan scipy.weave

def poisson_blending_weave(tgt, mask, src):
    h, w = mask.shape
    r,c = np.nonzero(mask)
    N = len(r)
    dg = src.astype(int32)
    g = tgt.astype(float32)
    val = np.array([0,0,0], dtype=int32)
    vd = np.array([0,0,0], dtype=int32)
    g[mask] = np.array([255.,255.,255.])
    
    MAX_ITER = 1000
    code = r"""
        for (int y=0; y<h-1; ++y){
            for (int x=0; x<w-1; ++x){
                if (MASK2(y,x)){
                    int count = 0;
                    float f_val = 0.0;
                    float f_vd = 0.0;
                    
                    for(int c=0; c<3; ++c){
                        VAL1(c) = 0.0;
                        VD1(c) = 0.0;
                    }
                    
                    if(y>0 && MASK2(y-1,x)){//top
                        count++;
                        for(int c=0; c<3; ++c){
                            VAL1(c) += SRC3(y-1,x,c);
                            VD1(c) += TGT3(y-1,x,c);
                        }
                    }
                    if(x>0 && MASK2(y,x-1)){//left
                        count++;
                        for(int c=0; c<3; ++c){
                            VAL1(c) += SRC3(y,x-1,c);
                            VD1(c) += TGT3(y,x-1,c);
                        }
                    }
                    if(y<h && MASK2(y+1,x)){//bottom
                        count++;
                        for(int c=0; c<3; ++c){
                            VAL1(c) += SRC3(y+1,x,c);
                            VD1(c) += TGT3(y+1,x,c);
                        }
                    }
                    if(x<w && MASK2(y,x+1)){//right
                        count++;
                        for(int c=0; c<3; ++c){
                            VAL1(c) += SRC3(y,x+1,c);
                            VD1(c) += TGT3(y,x+1,c);
                        }
                    }
                    for(int c=0; c<3; ++c){//center
                        VAL1(c) -= count*SRC3(y,x,c);
                        VD1(c) -= count*TGT3(y,x,c);
                        
                        f_val += abs(VAL1(c));
                        f_vd += abs(VD1(c));
                    }
                    
                    if (f_vd > f_val){
                        for(int c=0; c<3; ++c){
                            DG3(y,x,c) = VD1(c);
                        }
                    }
                    else{
                        for(int c=0; c<3; ++c){
                            DG3(y,x,c) = VAL1(c);
                        }
                    }
                    
                }
            }
        }
        
        float delta;
        //Jacobi method
        int y,x,count;
        for (int iter=0; iter<MAX_ITER; iter++){
            delta = 0.0;
            for (int i=0; i<N-1; ++i){
                y = R1(i);
                x = C1(i);
                count = 0;
                for(int c=0; c<3; ++c){
                    VAL1(c) = 0.0;
                }
                
                if(y>0){//top
                    count++;
                    if (!MASK2(y-1,x)){
                        for(int c=0; c<3; ++c){
                            VAL1(c) += TGT3(y-1,x,c);
                        }
                    }else
                    for(int c=0; c<3; ++c){
                        VAL1(c) += G3(y-1,x,c);
                    }
                }
                if(x>0){//left
                    count++;
                    if (!MASK2(y,x-1)){
                        for(int c=0; c<3; ++c){
                            VAL1(c) += TGT3(y,x-1,c);
                        }
                    }else
                    for(int c=0; c<3; ++c){
                        VAL1(c) += G3(y,x-1,c);
                    }
                }
                if(y<h){//bottom
                    count++;
                    if (!MASK2(y+1,x)){
                        for(int c=0; c<3; ++c){
                            VAL1(c) += TGT3(y+1,x,c);
                        }
                    }else
                    for(int c=0; c<3; ++c){
                        VAL1(c) += G3(y+1,x,c);
                    }
                }
                if(x<w){//right
                    count++;
                    if (!MASK2(y,x+1)){
                        for(int c=0; c<3; ++c){
                            VAL1(c) += TGT3(y,x+1,c);
                        }
                    }else
                    for(int c=0; c<3; ++c){
                        VAL1(c) += G3(y,x+1,c);
                    }
                }
                
                //updating
                for(int c=0; c<3; ++c){
                    VAL1(c) -= DG3(y,x,c);
                    float newval = VAL1(c)/(float)count;
                    if (newval < 0.) newval = 0.;
                    if (newval > 255.) newval = 255.;
                    delta += abs(G3(y,x,c)-newval);
                    G3(y,x,c) = newval;
                }
            }
            if (delta<1.0) break;
        }
    """
    inline(code, ['src','mask','tgt','dg','g','val', 'vd','r','c','N','w','h','MAX_ITER'], compiler='gcc')
    return g

Implementasi pada Python III : menggunakan library pyamg

Berbeda dengan implementasi sebelumnya yang dilakukan di tempat, implementasi ini tiap piksel diubah dulu menjadi vektor dan kemudian operator laplacian dibuat sebagai matriks jarang. Setelah sistem persamaan liniernya diselesaikan, nilai vektor x dikembalikan ke piksel di citra tujuan.

def poisson_blending_pyamg(m, mask, src,iterations=50):
    h,w = mask.shape
    r,c = mask.nonzero()
    N = len(r)
    idx = np.zeros(mask.shape, dtype=uint32)
    for i in xrange(N):
        idx.itemset((r.item(i),c.item(i)), i+1)
    
    #laplace per channel
    slap = src.astype(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])
    
    #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)
        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)))
        p = i
        Np = 0
        if y>0 and mask.item((y-1,x)):
            q = idx.item((y-1, x))-1
            row.append(p)
            col.append(q)
            data.append(-1.)
            Np += 1
        if x>0 and mask.item((y,x-1)):
            q = idx.item((y, x-1))-1
            row.append(p)
            col.append(q)
            data.append(-1.)
            Np += 1
        if y<h and mask.item((y+1,x)):
            q = idx.item((y+1, x))-1
            row.append(p)
            col.append(q)
            data.append(-1.)
            Np += 1
        if x<w and mask.item((y,x+1)):
            q = idx.item((y, x+1))-1
            row.append(p)
            col.append(q)
            data.append(-1.)
            Np += 1
        row.append(p)
        col.append(p)
        data.append(Np*1.)
    
    #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
    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)
    
    x = x.transpose()
    for i in xrange(N):
        yy,xx = r.item(i), c.item(i)
        v = m[yy,xx] - x[i]
        v[v<0] = 0
        v[v>255] = 255
        m[yy,xx] = v
    return m

Kode pemanggil

import numpy as np
from scipy.sparse import *
from scipy.ndimage.filters import laplace
from scipy import *
from scipy.weave import inline

from pyamg import solve, smoothed_aggregation_solver
from pyamg.gallery import poisson
from pyamg.relaxation import *
from pyamg.util.linalg import norm

import cv2
import cv2.cv as cv

def combine(a, amask, b, x, y):
    c = b.copy()
    h,w = a.shape[:2]
    aa = c[y:y+h,x:x+w]
    
    #pure python solving (jacobi)
    #aa = poisson_blending(aa, amask, a)
    
    #uses relaxation module in pyamg solver
    #poisson_blending_pyamg(aa, amask, a)
    
    #jacobi solver using weave
    aa = poisson_blending_weave(aa, amask, a)
    c[y:y+h,x:x+w] = aa
    
    cv2.imshow("blended", c)
    return c
    
    
files = ["bear.jpg", "bear-mask.jpg", "mountain.jpg", 350,180, 0]

im1 = cv2.imread(files[0])
im1mask = cv2.imread(files[1])
im2 = cv2.imread(files[2])

graymask = cv2.cvtColor(im1mask, cv.CV_BGR2GRAY)
mask = graymask>250

cv2.imshow("im1", im1)
cv2.imshow("im1mask", im1mask)
combine(im1[files[5]:,:,:], mask[files[5]:,:], im2, files[3],files[4])

cv2.waitKey(0)

Penutup

Beberapa hal masih mengganjal pada implementasi yang sudah saya buat. Isu pertama adalah penentuan banyaknya iterasi pada implementasi python murni dan pyamg. Hal ini disebabkan implementasi python murni terlalu lambat dan implementasi menggunakan pyamg sulit diprediksi secara heuristik kondisi yang cukup untuk berhenti. Pada implementasi menggunakan inline c++ sudah ada mekanisme berhenti otomatis. Isu berikutnya adalah pemetaan kembali nilai hasil vektor X pada penyelesaian yang menggunakan metode relaksasi. Hal ini disebabkan nilai vektor x memiliki rentang negatif ke positif dan rentang ini akan semakin besar jika nilai parameter banyak iterasi juga diperbesar. Meskipun demikian tetapi gambar yang dihasilkan tidak semakin baik. Kualitas hasil penggabungan citra biasanya bergantung pada kondisi piksel di citra sumber yang tercakup area mask. Keberhasilan penggabungan ini biasanya tidak akan tampak aneh jika memang nilai piksel-piksel ini tidak jauh berbeda dengan piksel perbatasan di citra tujuan khususnya warna dan gradiennya.

About these ads


One Response to “Seamless Cloning using Poisson Image Blending in Python”


  1. 1 Poisson Image Editing (revised) « GAIBlog

Tinggalkan Balasan

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

WordPress.com Logo

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

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s


Ikuti

Get every new post delivered to your Inbox.

Bergabunglah dengan 2.049 pengikut lainnya.

%d bloggers like this: