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.
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.
GAIBlog 
terima kasih, bagus peneranganya