Algoritma Papoulis-Gerchberg

Algoritma Papoulis-Gerchberg adalah salah satu algoritma yang dapat digunakan untuk melakukan rekonstruksi gambar. Contoh yang dibahas kali ini adalah mengenai rekonstruksi gambar dengan resolusi yang lebih tinggi.
Dalam kasus ini, algoritma Papoulis-Gerchberg akan dibandingkan dengan algoritma Bicubic Interpolation sebagai salah satu teknik dasar dalam melakukan rekonstruksi gambar.

* Tentukan gambar utama yang akan digunakan dalam perhitungan
Diasumsikan gambar utama yang digunakan adalah

Sebelum masuk kedalam langkah-langkah pembahasan algoritma, ada beberapa konstanta atau parameter yang harus diketahui, yaitu:
* Tentukan estimasi rotasi, estimasi pergeseran, dan faktor rekonstruksi
Nilai parameter ini bisa diganti sesuai kebutuhan
Diasumsikan dalam kasus ini semua nilai yang digunakan adalah 2

estimasiRotasi = 2;
estimasiPergeseran = [2 2];
faktorRekonstruksi = 2;

Langkah-langkah penggunaan algoritma ini adalah

1. dapatkan nilai pixel dari masing-masing komponen warna merah, hijau, dan biru dari gambar tersebut

imgKomponenMerah = img(:,:,1); 
imgKomponenMerah = mat2cell(double(imgKomponenMerah), size(imgKomponenMerah,1));
imgKomponenHijau = img(:,:,2); 
imgKomponenHijau = mat2cell(double(imgKomponenHijau), size(imgKomponenHijau,1));
imgKomponenBiru = img(:,:,3); 
imgKomponenBiru = mat2cell(double(imgKomponenBiru), size(imgKomponenBiru,1));

2. Lakukan perhitungan rekonstruksi gambar menggunakan Bicubic Interpolation untuk masing masing komponen warna

interpolationKomponenMerah = interpolation(imgKomponenMerah,estimasiPergeseran,estimasiRotasi,faktorRekonstruksi); 
interpolationKomponenHijau = interpolation(imgKomponenHijau,estimasiPergeseran,estimasiRotasi,faktorRekonstruksi); 
interpolationKomponenBiru = interpolation(imgKomponenBiru,estimasiPergeseran,estimasiRotasi,faktorRekonstruksi); 

Memasuki perhitungan pada fungsi interpolation

2a. Dapatkan jumlah gambar dan ukuran gambar
Kemudian dapatkan titik pusat dan hitung sudut rotasi dalam radian

jumlahGambar=length(s);
ukuranGambar = size(s{1});
if (length(ukuranGambar)==2) ukuranGambar=[ukuranGambar 1]; end
pusat = (ukuranGambar+1)/2;
sudutRotasiRadian = estimasiRotasi*pi/180;

2b. Lakukan perhitungan pada masing-masing gambar (poin 2b1 – 2b8)

for k=1:ukuranGambar(3)
    for i=1:jumlahGambar
	. . .

2b1. Dapatkan vektor dari input data

vektorGambar{i}=s{i}(:,:,k);
vektorGambar{i} = vektorGambar{i}(:);

2b2. Buat matriks dengan isian indeks baris dan indeks kolom dari matriks input

idxBarisMatriks{i} = [1:faktorRekonstruksi:faktorRekonstruksi*ukuranGambar(1)]'*ones(1,ukuranGambar(2));
idxKolomMatriks{i} = ones(ukuranGambar(1),1)*[1:faktorRekonstruksi:faktorRekonstruksi*ukuranGambar(2)];

2b3. Lakukan pergeseran indeks baris dan indeks kolom menuju titik pusat

idxBarisMatriks{i} = idxBarisMatriks{i}-faktorRekonstruksi*pusat(1);
idxKolomMatriks{i} = idxKolomMatriks{i}-faktorRekonstruksi*pusat(2);

2b4. Lakukan rotasi menggunakan sudut rotasi yang sudah ditentukan sebelumnya

koordinatVektor{i} = [idxKolomMatriks{i}(:) idxBarisMatriks{i}(:)]*[cos(sudutRotasiRadian(i)) sin(sudutRotasiRadian(i)); -sin(sudutRotasiRadian(i)) cos(sudutRotasiRadian(i))];

2b5. Kembalikan pergeseran data yang sudah dilakukan sebelumnya sebesar parameter estimasi pergeseran

idxBarisMatriks{i} = koordinatVektor{i}(:,2)+faktorRekonstruksi*pusat(1)+faktorRekonstruksi*estimasiPergeseran(i,1);
idxKolomMatriks{i} = koordinatVektor{i}(:,1)+faktorRekonstruksi*pusat(2)+faktorRekonstruksi*estimasiPergeseran(i,2);

2b6. Dapatkan indeks baris, indeks kolom, beserta nilai data menggunakan indeks yang baru

indeksBaris{i} = idxBarisMatriks{i}((idxBarisMatriks{i}>0)&(idxBarisMatriks{i}<=faktorRekonstruksi*ukuranGambar(1))&(idxKolomMatriks{i}>0)&(idxKolomMatriks{i}<=faktorRekonstruksi*ukuranGambar(2)));
indeksKolom{i} = idxKolomMatriks{i}((idxBarisMatriks{i}>0)&(idxBarisMatriks{i}<=faktorRekonstruksi*ukuranGambar(1))&(idxKolomMatriks{i}>0)&(idxKolomMatriks{i}<=faktorRekonstruksi*ukuranGambar(2)));
nilai{i} = vektorGambar{i}((idxBarisMatriks{i}>0)&(idxBarisMatriks{i}<=faktorRekonstruksi*ukuranGambar(1))&(idxKolomMatriks{i}>0)&(idxKolomMatriks{i}<=faktorRekonstruksi*ukuranGambar(2)));

2b7. Masukkan semua data yang ditemukan kedalam variabel

dataNilai = []; dataIndeksBaris = []; dataindeksKolom = [];
for i=1:jumlahGambar
	dataIndeksBaris = [dataIndeksBaris; indeksBaris{i}];
	dataindeksKolom = [dataindeksKolom; indeksKolom{i}];
	dataNilai = [dataNilai; nilai{i}];
end

2b8. Lakukan perhitungan interpolation menggunakan vektor indeks baris, indeks kolom dan nilai data yang baru
Kemudian lakukan reshape untuk mengembalikan gambar menjadi ukuran semula

outputInterpolation = griddata(dataindeksKolom,dataIndeksBaris,dataNilai,[1:ukuranGambar(2)*faktorRekonstruksi],[1:ukuranGambar(1)*faktorRekonstruksi]','cubic');
rekonstruksi(:,:,k) = reshape(outputInterpolation,ukuranGambar(1)*faktorRekonstruksi,ukuranGambar(2)*faktorRekonstruksi);

3. Lakukan perhitungan rekonstruksi gambar menggunakan Papoulis Gerchberg untuk masing masing komponen warna

pgKomponenMerah = papoulisgerchberg(imgKomponenMerah,estimasiPergeseran,faktorRekonstruksi); 
pgKomponenHijau = papoulisgerchberg(imgKomponenHijau,estimasiPergeseran,faktorRekonstruksi); 
pgKomponenBiru = papoulisgerchberg(imgKomponenBiru,estimasiPergeseran,faktorRekonstruksi); 

3a. Dapatkan vektor x dan y sebagai indeks baris dan indeks kolom
x: vektor kolom yang berisi sampel data
sebagai contoh: jika ukuranY = 2, ukuranX = 3, maka x' = 1 1 2 2 3 3; dan y' = 1 2 1 2 1 2;

[ukuranY, ukuranX] = size(s{1});
[X,Y] = meshgrid(1:ukuranX,1:ukuranY);
x = X(:);
y = Y(:);

3b. Hitung estimasi pergeseran setelah memperhitungan faktor rekonstruksi

estimasiRotasi = round(estimasiRotasi*faktorRekonstruksi);

3c. Lakukan inisialisasi variable HRimg sebagai hasil gambar setelah dilakukan perkiraan
sebagai contoh:
jika s{1} = [1 1 1]
faktorRekonstruksi = 4
maka N0 = 1+faktorRekonstruksi*(size(s{1})-1);
= 3*(3-1)+1=9;
HRimg = N0+(maksEstimasiRotasi-minEstimasiRotasi);
= [1 0 0 0 1 0 0 0 1];

minEstimasiRotasi = min(estimasiRotasi,[],1);
maksEstimasiRotasi = max(estimasiRotasi,[],1);
Delta = maksEstimasiRotasi-minEstimasiRotasi;
HRimg = zeros(1+faktorRekonstruksi*(size(s{1})-1)+Delta);
ukuranHR = size(HRimg);

3d. Masukkan data input dengan ukuran yang baru sesuai dengan indeks baris dan kolom yang sudah ditentukan sebelumnya
indeksBaris: koordinat y untuk menyusun gambar menggunakan vektor x
indeksKolom: koordinat x untuk menyusun gambar menggunakan vektor x

indeksBaris = [];indeksKolom = [];
for i = 1:length(s),
    indeksBaris_temp=1-minEstimasiRotasi(1)+faktorRekonstruksi*(y-1)+estimasiRotasi(i,1);
    indeksKolom_temp=1-minEstimasiRotasi(2)+faktorRekonstruksi*(x-1)+estimasiRotasi(i,2);
    indeksBaris = [indeksBaris; indeksBaris_temp];
    indeksKolom = [indeksKolom; indeksKolom_temp];
    HRimg(indeksBaris_temp,indeksKolom_temp) = s{i}(y,x);
end

3e. Lakukan inisialisasi variable idxKolomX, NLR, barisFFT0, kolomFFT0
barisFFT0: data indeksBaris saat nilai spektrum FFT = 0
kolomFFT0: data indeksKolom saat nilai spektrum FFT = 0

idxKolomX = HRimg((indeksKolom-1)*ukuranHR(1)+indeksBaris);
NLR = floor(sqrt(length(s))*size(s{1})/2)*2;
barisFFT0 = [(1+NLR(1)/2+1):(ukuranHR(1)-NLR(1)/2-1)];
kolomFFT0 = [(1+NLR(2)/2+1):(ukuranHR(2)-NLR(2)/2-1)];

3f. Lakukan inisialisasi variabel y dan salin ke dalam variabel lain
Kemudian beri nilai awal dengan nilai indeks kolom x

y=zeros(ukuranHR);
ySebelumnya=y;
y((indeksKolom-1)*ukuranHR(1)+indeksBaris)=idxKolomX;

%3g. Lakukan perhitungan sebanyak jumlah iterasi (poin 3g1 -3g4)

3g1. Lakukan transformasi data menggunakan teknik FFT (Fast Fourier Transform)
Kemudian kosongkan semua data pada posisi dengan nilai spektrum = 0

Y=fft2(y);
Y(barisFFT0,:)=0;
Y(:,kolomFFT0)=0;

3g2. Lakukan invers dari transformasi yang sudah dilakukan sebelumnya
Kemudian tambahkan nilai awal pada posisi dengan nilai spektrum = 0

y=ifft2(Y);
y((indeksKolom-1)*ukuranHR(1)+indeksBaris)=idxKolomX;

3g3. Hitung delta sebagai selisih antara hasil perhitungan pada iterasi ini dengan iterasi sebelumnya
Jika nilai delta sudah mendekati 0 maka hentikan perhitungan

delta=norm(y-ySebelumnya)/norm(y);
Epoch=[Epoch; iterasi delta];
iterasi= iterasi+1;
if iterasi>3 
	if abs(Epoch(iterasi-3,2)-delta) < 1e-8
		break  
	end
end

3g4. Simpan nilai perhitungan untuk digunakan pada perulangan selanjutnya

ySebelumnya=y;

Hasil akhir adalah: (klik untuk perbesar gambar)

Contoh modul / source code menggunakan Matlab dapat didownload disini:

Jika membutuhkan jasa kami dalam pembuatan program, keterangan selanjutnya dapat dilihat di Fasilitas dan Harga
Jika ada yang kurang paham dengan langkah-langkah algoritma diatas, silahkan berikan komentar Anda.
Selamat mencoba.