Pencitraan Sintesis (3)

Pencitraan sintesis memanfaatkan prinsip interferometri untuk membentuk sebuah apertur “maya” raksasa dari beberapa aperture kecil yang tersebar dalam ruang. Teknik ini disebut sebagai sintesis apertur. Tulisan ini menjelaskan geometri pencitraan sintesis, formulasi matematika, dan proses sintesis apertur.

Formulasi Matematika Pencitraan Sintesis

Tinjau benda langit yang memancarkan gelombang radio yang terletak pada posisi R dari titik acuan di bumi. Medan listrik berubah waktu sebesar E(R,t) dari sumber akan dipancarkan ke segala penjuru, hingga tiba pada pengamat di Bumi yang terletak pada posisi r. Analisis dapat dilakukan dengan lebih sederhana dalam notasi fasor, dimana perubahan sinusoidal terhadap waktu diserap dan tidak berpengaruh langsung pada manipulasi matematis. Jika ditinjau komponen frekuensi tertentu dari gelombang, maka fasor medan listrik yang bernilai kompleks ini adalah . Subscript penanda frekuensi akan dihilangkan, namun tetap diingat bahwa gelombang hanya ditinjau pada suatu nilai frekuensi. Kondisi ini disebut kuasi monokromatik dan selanjutnya fasor medan listrik hanya akan dituliskan sebagai E(R).

Gambar 1. Geometri pencitraan interferometri radio

Komponen medan listrik dari gelombang radio/elektromagnetik (EM) yang tiba pada titik pengamat adalah hasil superposisi linier dari seluruh gelombang yang berasal dari sumber radio

(1)

Fungsi P(R, r) disebut juga propagator; suatu faktor yang menyatakan bagaimana perubahan medan listrik pada sumber di R mempengaruhi perubahan medan di sisi penerima pada posisi r. Jika ruang antara pengamat dan sumber dapat dianggap kosong, propagator akan memiliki bentuk sederhana dan persamaan (1) akan menjadi

(2)

dimana c adalah kecepatan rambat gelombang EM (3×10⁸ m/s). Persamaan (2) menyatakan bentuk umum komponen kuasi-monokromatik dari medan listrik pada frekuensi dari sumber radio kosmik. Teleskop radio dapat mengamati kuat medan listrik ini. Namun demikian, besaran yang lebih penting lagi adalah korelasi antar medan listrik yang teramati di dua stasiun pengamat atau teleskop radio yang berbeda; stasiun pertama terletak pada posisi r₁, sedangkan yang kedua berada pada posisi r₂. Korelasi antar keduanya dinyatakan sebagai visibility

(3)

dimana tanda ”*” menyatakan konjugasi kompleks. Substitusi persamaan (2) ke persamaan (3) akan menghasilkan

(4)

Dengan asumsi yang cukup realistis bahwa radiasi benda langit secara spasial bersifat tak-koheren, maka akan berharga nol jika , sehingga persamaan (4) dapat disederhanakan menjadi

(5)

Penyederhanaan lebih lanjut dapat diperoleh mengingat R>>r₁ dan R>>r₂, karena jarak antara pengamat (di Bumi) ke sumber ke radio berada pada orde tahun cahaya, sedangkan posisi pengamat atau baseline r₁ maupun r₂ dalam orde kilometer. Akibatnya pendekatan berikut ini dapat dipakai

dengan s menyatakan vektor satuan pada arah R, atau s = R/|R|.

Nilai magnitudo kuadrat dari medan listrik menyatakan intensitas atau tingkat kecerahan sumber radio. Dapat didefinisikan intensitas radiasi sebagai I= <|E|>² |R|². Dengan menggantikan integrasi luas dari objek dS dengan sudut , maka persamaan (5) akan menjadi

(6)

Inilah rumus dasar (governing equation) dalam pencitraan sintesis radio-interferometri, suatu formula yang menghubungkan visibility V(r) yang diukur oleh (jaringan) radio-teleskop dengan (sebaran) intensitas pemancar I(s) yang ingin ditentukan. Persamaan diatas menyatakan bahwa nilai visibility yang terukur tidak bergantung pada letak absolut r₁ maupun r₂, melainkan letak relatifnya, sehingga pengukuran dapat dilakukan dengan cara membuat satu stasiun pengamat tetap dan stasiun yang lain bergerak melakukan penyapuan ruang. Visibility juga disebut sebagai fungsi koherensi spasial atau fungsi autokorelasi spasial. Gambar 2 memperlihatkan diagram blok interferometer teleskop radio.

Gambar 2. Diagram blok interferometer teleskop radio (diadopsi dari [1])

Proses Inversi Citra Interferometri Radio

Jika sumber radio cukup kecil, dengan mengabaikan orde tinggi maka vektor s dapat dijabarkan sebagai s= s₀ + s, dimana s₀ ortonormal dan tegaklurus terhadap s. Suatu sistem koordinat baru dimana s₀ = (0,0,1) dapat dipakai, sehingga posisi relatif kedua stasiun penerima itu menjadi . Dengan demikian, vektor s akan menjadi s = (x,y,1) dan persamaan (6) dapat dinyatakan dalam koordinat citra berbentuk bidang datar (x,y)

(7)

Faktor eksponensial kompleks didepan integral pada ruas kanan biasanya diserap ke ruas kiri, lalu lambang yang sama untuk visibility dipakai lagi. Dengan demikian, persamaan (7) dapat disederhanakan lagi menjadi

[8]

dimana V (u,v) adalah fungsi koherensi relatif terhadap pusat fasa s₀.

Persamaan [8] adalah ekspresi transformasi Fourier dari citra I(x,y), dengan demikian citra yang dicari dapat diperoleh dari visibility melalui inversi Fourier

(9)

Pada prakteknya, tidak semua nilai fungsi visibility V di ruang u-v dapat diketahui, kecuali di sebagian kecil titik cuplikan di bidang Fourier dengan kualitas yang kurang baik. Dengan demikian, nilai visibilitas pada ruas kanan digantikan dengan nilai cuplikannya S(u,v)V(u,v) dimana

(10)

dan persamaan inversi pada persamaan (9) menghasilkan citra kotor (dirty image)

(11)

Pembobot w_k pada persamaan (11) dapat diatur agar terjadi perbaikan pada gambar kotor , misalnya dengan memasukkan waktu integrasi atau faktor kerapatan lokal dari cuplikan.

Daftar Pustaka

1. B.G. Clark, “Chapter 1: Interferometers and Coherence Theory,” in J.A. Zensus, P.J. Diamond, and, P.J. Napier (eds.), Very Long Baseline Interferometry and VLBA, ASP Conf. Series, Vol.82, 1995.

2. T. Cornwell, “Chapter 3: Imaging Concepts,” in J.A. Zensus, P.J. Diamond, and, P.J. Napier (eds.), Very Long Baseline Interferometry and VLBA, ASP Conf. Series, Vol.82, 1995.

Compressive SFCW-GPR/ Pencitraan Kompresif (3)

Radar pensintesa frekuensi SFCW-GPR memiliki kelemahan dalam hal kecepatan pengumpulan data. Inilah salah satu alasan mengapa sistem radar ini kurang berkembang. Berbagai upaya telah dilakukan untuk mengatasi masalah ini, misalnya dengan mengirimkan beberapa frekuensi sekaligus secara serempak. Munculnya teknik pencitraan kompresif memungkinkan masalah ini dapat diatasi dengan pendekatan lain. Tulisan ini menjelaskan prinsip kerja SFCW-GPR kompresif dan hasil-hasil penelitian yang telah dieproleh hingga saat ini.

Cerita lengkap dapat dibaca dari presentasi di sini atau di situ.

In English

Data acquisition speed is an inherent problem of the stepped-frequency continuous wave (SFCW) radars, which discouraging further usage and development of this technology. We propose an emerging paradigm called the compressed sensing (CS), to manage this problem. In the CS, a signal can be reconstructed exactly based on only a few samples below the Nyquist rate. Accordingly, the data acquisition speed can be increased significantly. A novel design of SFCW ground penetrating radar (GPR) with a high acquisition speed capability is proposed and evaluated. Simulation by a monocycle waveform and actual measurement by a Vector Network Analyzer in a GPR test-range confirm the implementability of the proposed system.

You can read the full story here or there.

Pencitraan Kompresif (1.b)

Pencitraan kompresif pada kamera piksel-tunggal dapat dilakukan dalam kawasan spasio-temporal maupun kawasan frekuensi. Pencitraan kawasan spasio-temporal merupakan tafsiran langsung dari pencuplikan kompresif; lebih sederhana, dan intuitif. Namun demikian, penggunaan DCT sebagai transformasi penjarang tidak memberikan hasil yang memuaskan. Tulisan ini memberikan alternatif transformasi penjarang berbentuk wavelet paling sederhana, yaitu transformasi Haar.

Transformasi Haar untuk Pencuplikan Kompresif

Ketidaksesuaian DCT sebagai transformasi penjarang terlihat dari peluruhan koefisien dari citra berbentuk ”H” yang relatif kurang cepat. Hal ini telah dapat diperbaiki dengan meminimumkan nilai TV (total variance). Cara ini cukup ampuh, namun kurang intuitif karena tidak secara langsung memberikan transformasi pada kawasan basis penjarang-nya. Tulisan ini menjelaskan penggunaan transformasi Haar, salah satu bentuk transformasi wavelet paling tua dan paling sederhana, sebagai transformasi penjarang-nya.

Konstruksi matriks penjarang PSI pada kode sebelumnya diubah dari transformasi DCT menjadi transformasi Haar dengan mengubah baris terkait sbb :

PSI = haar_T(N)’;%(eye(N,N)); % ** buat matriks transformasi Haar **

dimana ”haar_T” diatas adalah fungsi untuk membentuk matriks transformasi Haar yang kode Matlab-nya diperlihatkan di akhir tulisan.

Hasil dari transformasi Haar dengan citra dan parameter kompresi sebesar 1.4 (seperti pada tulisan pertama) diberikan pada Gambar 1 berikut ini.

(a)

(b)

Gambar 1. Hasil rekonstruksi pencitraan kompresif dengan transformasi Haar: (a) K-buah koefisien pertama, (b) 50 buah koefisien pertama.

Gambar 2. Plot koefisien terurut dari citra ”H” terhadap basis Haar

Hasil rekonstruksi terlihat masih kurang baik, meskipun bagian utama dari citra sudah terambil. Ini bisa dijelaskan dengan melihat kembali plot koefisien penyusun citra hasil optimasi yang ditampilkan pada Gambar 2. Terlihat bahwa meskipun koefisien meluruh jauh lebih cepat dari DCT, nilai koefisien nol tercapai pada indeks urutan ke sekitar 50. Dengan demikian, asumsi derajat sparsity sebesar 32 masih kurang mencukupi. Bahkan, melibatkan 50 buah koefisien hasil optimisasi-pun masih belum menghasilkan rekonstruksi yang sempurna, seperti diperlihatkan pada Gambar 1.b.

Gambar 3. Citra asli papan catur

Sebagai perbandingan, berikut ini disajikan hasil pencitraan kompresif dengan basis Haar untuk citra papan catur (checker board) yang diperlihatkan pada Gambar 3. Relatif terhadap basis Haar, citra ini sangat sparse, seperti diperlihatkan pada Gambar 4.

Gambar 4. Koefisien terurut citra papan-catur hasil pemilihan basis

Kali ini tingkat kompresi diperbesar menjadi 2.3 kalinya, dengan cara memperkecil K menjadi sekitar . Gambar 5 menunjukkan citra teramati, dimana piksel abu-abu menyatakan piksel yang tak terlihat oleh pengamatan DMD. Hasil rekonstruksi diperlihatkan pada Gambar 5. Terlihat bahwa rekonstruksi eksak sudah tercapai.

(a)

(b)

Gambar 5. (a) Citra teramati dan (b) citra rekonstruksi

Kode Matlab

% ——————————————————————

% ** fungsi membentuk matriks transformasi Haar **

% ——————————————————————

function [H]=haar_T(h_size);

n=ceil(log2(h_size));N=2^n;

H=zeros(N,N);

norm_factor=1/sqrt(N);

%difine z-axis

z=(0:1:N-1)/N;

H(1,:)=1;

for p=0:n-1;

for q=1:2^p;

k=2^p+q-1;

%disp(sprintf(‘k=%d , p=%d, q= %d \n’,k,p,q)); % debug p and q

hk=zeros(N,1);

hk_1=2^(p/2); hk_2=-(2^(p/2));

lb_1=ceil(N*(q-1)/2^p);ub_1=ceil(N*(q-0.5)/2^p);

lb_2=ceil(N*(q-0.5)/2^p);ub_2=ceil(N*q/2^p);

% *debug z-bound

% disp(sprintf(‘lb_1=%d , ub_1=%d, lb_2= %d, ub_2= %d\n’,lb_1,ub_1,lb_2,ub_2));

H(k+1,lb_1+1:ub_1)=hk_1;

H(k+1,lb_2+1:ub_2)=hk_2;

end

end;

H=H*norm_factor;

Pencitraan Kompresif (2)

Pemilihan transformasi penjarang yang tepat sangat berpengaruh pada hasil rekonstruksi pencitraan kompresif. Tulisan ini menjelaskan pemakaian teknik minimisasi variansi total (TV-total variance), salah satu estimasi kokoh (robust) yang sebenarnya sudah sering dipakai namun baru disadari hubungannya dengan pencuplikan kompresif. Karena parameter dari TV minimum saja tidak bisa dipakai untuk melakukan inversi, maka harus ada basis lain yang terhubung dengan seleksi TV minimum. Diagram ranah-jamak dapat dipakai untuk menjelaskan keterkaitan antar transformasi dan parameter yang digunakan dalam pencitraan kompresif kawasan frekuensi.

Pencuplikan Kompresif Kawasan Frekuensi

Operasi pengamatan pada kawasan spasial yang diberikan pada persamaan (1) didalam tulisan sebelumnya dapat ditafsirkan sebagai pencuplikan terhadap vektor basis penjarang jika dituliskan sebagai

(1)

Setiap baris pada hanya berisi satu buah angka bernilai satu, sedangkan sisanya nol. Perkalian dengan mengakibatkan seleksi vektor baris pada indeks yang bersangkutan dengan letak kolom elemen aktif tersebut. Sebagai akibatnya, vektor basis dalam akan berkurang dan letaknya diacak.

Kinerja rendah yang dihasilkan akibat penggunaan DCT sebagai transformasi penjarang dapat ditingkatkan dengan transformasi bentuk lain. Salah satu yang telah sering dipakai tetapi baru disadari kaitannya dengan pencuplikan kompresif adalah variansi total (TV-total variance). Pendekatan dengan cara ini cukup berbeda jika dibandingkan dengan metoda sebelumnya yang merupakan tafsiran langsung pencuplikan spasial. Kaitan berbagai transformasi dan proses recovery sinyal dilukiskan dalam diagram ranah jamak berikut ini.

Gambar 1. Diagram ranah-jamak [1] pencuplikan kompresif kawasan frekuensi

Mengacu pada diagram ranah-jamak, sinyal s terlebih dahulu ditransformasikan ke kawasan frekuensi menjadi . Pencuplikan kompresif kawasan ini dengan transformasi pengukur akan menghasilkan cuplikan kompresif . Karena basis penjarang ada dalam kawasan lain yang perhitungannya didasarkan pada sinyal kawasan spasio-temporal, maka hasil cuplikan kompresif terlebih dahulu dikembalikan ke kawasan semula menjadi sinyal . Selanjutnya yang akan diminimumkan adalah nilai variansi total:

(2)

Basis pursuit untuk kasus ini dapat dirumuskan sebagai

(P1): min |S| s.t. (3)

Optimisasi diatas mengatakan:

“pilih yang inverse transform Fouriernya, yaitu memiliki variansi total |S| terkecil dan bisa menjelaskan hasil pengamatan .

Pada diagram, pengaruh hasil pemilihan |S| terkecil pada koefisien Fourier digambarkan sebagai panah bergaris ganda pada bagian kiri diagram. Operasi pengambilan nilai variansi total digambarkan sebagai transformasi merugi karena pengetahuan nilai variansi total saja tidak dapat dipakai untuk merekonstruksi sinyal semula.

Eksperimen dan Analisis

Pada eksperimen dipakai citra yang sama seperti sebelumnya, tetapi mekanisme rekonstruksi-nya sangat berbeda; yang dapat disimpulkan dengan membandingkan diagram ranah jamak pada Gambar 1 dengan yang ada pada tulisan sebelumnya. Fungsi optimisasi yang digunakan ada dalam paket program l₁-magic dari Caltech. Program utama sudah sedikit dimodifikasi supaya sesuai dengan keperluan, seperti diperlihatkan pada lampiran di akhir tulisan ini.

Gambar 2. Citra asli

Gambar 3. Citra teramati

Gambar 2 menunjukkan citra asli, sedangkan Gambar 3 adalah citra teramati dilihat dari kawasan spasial. Piksel berharna abu-abu menyatakan ketiadaan data pada posisi bersangkutan, sedangkan piksel hitam dan putih menyatakan hasil pengamatan terkompresi sebanyak 4 kali. Tidak seperti pada tulisan sebelumnya, sangat sulit menebak citra sebenarnya dari citra teramati ini karena tingkat kompresi yang lebih tinggi.

Gambar 4. Citra hasil rekonstruksi berdasarkan TV

Hasil rekonstruksi yang diperlihatkan pada Gambar 4 jauh lebih bagus dari yang sebelumnya—bahkan tepat sama dengan citra asalnya. Disini terlihat betapa pentingnya pemilihan basis penjarang yang tepat didalam pencuplikan kompresif. Nilai TV menggambarkan smoothness karena berhubungan dengan variansi atau keragaman citra. Meminimumkan nilai ini berarti memilih citra dengan variansi rendah, atau smoothness yang tinggi. Ini sangat cocok untuk citra-citra tomografi dalam biomedika, karena organ tubuh manusia sudah semestinya smooth.

Daftar Pustaka

1. AB. Suksmono, ”A graphical representation of multi-domain signal processing,” Proc. of SICE-ICASE 2006, Busan, Korea.

2. EJ. Candes, J. Romberg, and T. Tao, ”Robust uncertainty principles: Exact signal recovery from highly incomplete frequency information,” IEEE Trans. Information Theory, Vol.52, no.2, Feb.2006, pp.489-509.

3. MF. Duarte, MA. Davenport, D. Takhar, JN. Laska, T. Sun, KF. Kelly, and RG. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Proc. Magazine (83), March 2008, pp.83-90.

4. AB. Suksmono, “A brief review on compressive imaging,” Doc. RGQ14-2/2/ -E/Jul 2008, ITU-D, Rapp. Group Meeting on Q.14-2/2, Tokyo, 3-4 July 2008.

Kode Matlab

% *******************************************

% kamera piksel tunggal: cs_camera.m

% ditulis kembali berdasarkan kode L1-magic

% fungsi-fungsi lain dapat di-download dari Caltech

% *******************************************

% use implicit, matrix-free algorithms ?

largescale = 1;

path(path, ‘./Optimization’); path(path, ‘./Measurements’); path(path, ‘./Data’);

% read original image

I=double(imread(‘IMG\H.bmp’,'bmp’));

%tampilkan citra asal;

figure(1); imagesc(I); colormap(gray); title(‘citra asli’);

[n_row,n_col] = size(I); n = n_row; % panjang baris/kolom

N = n_row*n_col; % total cuplikan asal/ jumlah piksel

I = I/norm(I(:)); %normalisasi energi citra

I = I – mean(I(:)); %mean dibuat nol

x = reshape(I,N,1); %ubah menjadi larik 1-D /vektor

% jumlah observasi yang diinginkan

f_comp=4;%faktor kompresi

K = N/f_comp;%jumlah pengamatan

% Simulasi DMD kawasan frekuensi

P = randperm(N)’;

q = randperm(N/2-1)+1;

OMEGA = q(1:K/2)’;

% measurement matrix

if (largescale)

A = @(z) A_f(z, OMEGA, P);

At = @(z) At_f(z, N, OMEGA, P);

% obsevations

b = A(x);

% initial point

x0 = At(b);

else

FT = 1/sqrt(N)*fft(eye(N));

A = sqrt(2)*[real(FT(OMEGA,:)); imag(FT(OMEGA,:))];

A = [1/sqrt(N)*ones(1,N); A];

At = [];

% observations

b = A*x;

% initial point

x0 = A’*b;

end

I0=reshape(x0, n, n);

figure(2);imagesc(I0);colormap(gray);title(‘citra teramati’);

tvI = sum(sum(sqrt([diff(I,1,2) zeros(n,1)].^2 + [diff(I,1,1); zeros(1,n)].^2 )));

disp(sprintf(‘Original TV = %.3f’, tvI)); %TV asal

time0 = clock;

xp = tveq_logbarrier(x0, A, At, b, 1e-3, 5, 1e-8, 200);

Ip = reshape(xp, n, n);

disp(sprintf(‘Total elapsed time = %f secs\n’, etime(clock,time0)));

figure(3);imagesc(Ip);colormap(gray);title(‘citra rekonstruksi’);

Pencitraan Kompresif (1)

Sebuah kamera dijital terdiri dari suatu susunan optik dan sebuah CCD (Charged-Coupled Device) yang berfungsi untuk menangkap bayangan benda, lalu mengubahnya menjadi citra dijital. Jika tidak ada CCD, satu buah fotodioda-pun dapat dipakai untuk membuat kamera dijital. Bagaimana caranya?

Kamera Piksel-Tunggal

Technology Review adalah sebuah majalah bulanan yang sering dijadikan acuan untuk melihat perkembangan teknologi terkini. Seperti pada tahun-tahun sebelumnya, sekitar bulan Maret majalah ini menampilkan sepuluh buah teknologi terkini yang akan menjadi trend di masa datang. Pada tahun 2006, satu dari 10 teknologi ini menarik perhatian khalayak luas karena dianggap sangat tidak biasa, yaitu Kamera Piksel-Tunggal. Kamera ini tidak memakai CCD sebagai pengubah bayangan benda menjadi gambar dijital, melainkan hanya satu buah fotodioda saja.

Sebuah kamera piksel-tunggal, seperti yang diperlihatkan pada Gambar 1, memanfaatkan susunan cermin renik (DMD-Digital Micromirror Device) yang dapat dikendalikan komputer untuk melakukan penyapuan bayangan benda secara acak. Bagian bayangan yang jatuh pada elemen cermin yang sedang ”aktif” diteruskan kesebuah fotodioda untuk diubah menjadi sinyal listrik. Perangkat DMD biasanya dipakai pada televisi proyeksi.

Gambar 1. Kamera piksel-tunggal

Informasi elemen DMD aktif dan intensitas cahaya yang sudah didijitasi selanjutnya disampaikan ke sebuah pengolah data. Karena tidak seluruh cermin diaktifkan untuk menangkap bayangan benda, jumlah piksel yang direkam lebih kecil jika dibandingkan dengan kamera dijital biasa. Selanjutnya rekonstruksi dilakukan menurut prinsip pencuplikan kompresif.

Dari sudut pandang pencuplikan kompresif, operasi pengambilan elemen bayangan secara acak oleh DMD adalah proses pengamatan yang dapat dinyatakan sebagai proyeksi s dengan matriks pengukur . Jika N jumlah total dari cermin DMD, sedangkan M adalah cermin aktif dalam satu kali pengambilan gambar, maka akan berupa sebuah matriks berukuran MxN. Setiap baris dari matriks ini akan berisi satu buah angka satu yang letaknya acak, sedangkan elemen sisanya pada baris ini bernilai nol. Operasi pengukuran ini dapat dituliskan sebagai:

(1)

Berdasarkan persamaan diatas, proses pengamatan dapat dipandang sebagai pencuplikan kompresif dengan dua cara yang berbeda, yaitu:

1. Pencuplikan kawasan spasial terhadap sinyal s secara acak
2. Pencuplikan kawasan transform sebagian dari basis penjarang secara acak

Dua cara pandang ini dapat memberikan dua pendekatan berbeda dalam mencari solusi atau merekonstruksi citra hasil pengamatan. Tulisan ini akan menjelaskan terlebih dahulu cara pertama, sedangkan cara kedua akan dijelaskan pada tulisan berikutnya.

Pencuplikan kompresif kawasan spasial

Sebelumnya telah diperkenalkan formulasi teknik pencuplikan kompresif dalam kawasan ruang-waktu. Secara prinsipil, teknik pencitraan yang pertama ini tidak berbeda dari yang telah dijelaskan dalam pencuplikan kompresif. Dengan demikian, diagram ranah-jamaknya masih tetap sama seperti sebelumnya.

Operasi proyeksi dilakukan dengan mengalikan sinyal asli s dengan operator proyeksi yang dibuat dengan ”mengaktifkan” satu elemen saja dari setiap vektor baris matriks ini. Setelah mengubah ”bayangan” benda dari larik dua dimensi menjadi satu dimensi, instruksi untuk mensimulasikan proses pengukuran dengan Matlab adalah:

% *** Simulasikan Digital Micromirror Device ***

rp=randperm(N); act_cell(1:M_sub)=rp(1:M_sub); %randomize mirror in DMD

PHI=zeros(M_sub,N);

for k=1:M_sub; PHI(k,act_cell(k))=1;end;

PHI=orth(randn(M_sub,N)’)'>0;

%*** Lakukan pengukuran x_sub=PHI*x0 ***

x_sub=PHI*x0;

Pada contoh ini akan diamati sebuah objek sederhana, yaitu citra yang menyerupai huruf “H” seperti yang diperlihatkan pada Gambar 2. Hasil pengukuran akan merupakan citra tak-sempurna karena piksel yang terambil hanya sedikit dan letaknya acak. Hasil ini diperlihatkan pada Gambar 3.

Gambar 2. Citra sederhana dalam pencitraan kompresif

Gambar 3. Citra hasil pengukuran

Berdasarkan hasil pengukuran, gambar dapat direkonstruksi dengan melalui algoritma basis pursuit. Sebagai contoh bisa diambil DCT sebagai transformasi penjarang. Jika menyatakan hasil pengamatan, vektor basis penjarang ke-k, maka pencarian basis dapat dirumuskan sebagai berikut:

P₁: min ||S||₁ s.t. (2)

Setelah seluruh vektor basis ditemukan, citra bisa direkonstruksi dengan menjumlahkan vektor-vektor basis ini yang diboboti dengan S_k yang sesuai. Hasil rekonstruksi ditunjukkan pada Gambar 4.

Gambar 4. Citra rekonstruksi

Meskipun citra huruf ”H” masih terlihat pada hasil rekonstruksi, Gambar 4 memperlihatkan bahwa secara visual hasilnya masih kurang baik. Hal ini bisa dijelaskan dengan melihat seluruh komponen S_k yang terambil pada pencarian basis.

Gambar 5. Magnitudo koefisien terpilih yang diurutkan

Gambar 5 memperlihatkan kurva magnitudo S_k yang telah diurutkan. Meskipun koefisien ini meluruh, terlihat bahwa peluruhan ini kurang cepat. Pengambilan koefisien sebanyak derajat sparsity yang diperkirakan masih akan menyisakan informasi yang cukup penting. Dengan demikian, masalah utamanya adalah pemilihan transformasi penjarang yang kurang tepat. Bagaimana cara mengatasinya?

Daftar Pustaka

1. Michael Wakin, Jason Laska, Marco Duarte, Dror Baron, Shriram Sarvotham, Dharmpal Takhar, Kevin Kelly, and Richard Baraniuk, “An Architecture for Compressive Imaging,” Proc. International Conference on Image Processing — ICIP 2006, Atlanta, GA, Oct. 2006.
2. Dharmpal Takhar, Jason Laska, Michael Wakin, Marco Duarte, Dror Baron, Shriram Sarvotham, Kevin Kelly and Richard Baraniuk, “A New Compressive Imaging Camera Architecture using Optical-Domain Compression,” Proc. of Computational Imaging IV at SPIE Electronic Imaging, San Jose, CA, Jan. 2006.
3. Marco F. Duarte, Mark A. Davenport, Dharmpal Takhar, Jason N. Laska, Ting Sun, Kevin F. Kelly and Richard G. Baraniuk, “Single Pixel Imaging via Compressive Sampling”, IEEE Signal Processing Magazine, March 2008.

Kode Matlab

%———————————————————-

% Pencitraan Kompresif: kamera dijital dng 1 foto-dioda

%———————————————————-

clear;%close all;clc;

% read original image

%cam01=double(imread(‘cameraman.tif’,'tif’));

%I=cam01(140:155,35:50); %load cam01; I=cam01;

%baca citra

I=double(imread(‘IMG\H.bmp’,'bmp’));

%tampilkan citra asal;

figure(1); imagesc(I); colormap(gray); title(‘citra asli’);

[N_col,N_row]=size(I);

N=N_col*N_row; %panjang sinyal

%Pembentukan basis sparsity /DCT

PSI=dct(eye(N,N))’; %buat matriks DCT

K_sparse=2*round(sqrt(N));%tingkat sparsity

x = reshape(I,N,1); %ubah menjadi vektor

x0=(x-mean(x))/(max(x)-min(x));%normalisasi

%proses pengamatan/pengukuran

%estimasikan minimal jumlah cuplikan M>CKlog(N)

mu=1*(1/sqrt(N)); %tentukan koherensi->perkiraan

C=sqrt(N)*mu; %koherensi tak ternormalisasi=sqrt(N).koh

M_sub=ceil(C*K_sparse*log(N));

% *** Simulasikan Digital Micromirror Device ***

rp=randperm(N); act_cell(1:M_sub)=rp(1:M_sub); %randomize mirror in DMD

PHI=zeros(M_sub,N);

for k=1:M_sub; PHI(k,act_cell(k))=1;end;

PHI=orth(randn(M_sub,N)’)'>0;

% *** Lakukan pengukuran x_sub=PHI*x0 ***

x_sub=PHI*x0;

x_obs=zeros(N,1);x_obs(act_cell)=x0(act_cell);

%Lakukan recovery

%Bentuk kamus basis D=PHI*PSI

D=PHI*PSI;

%x_sub=D*x0;

figure(4);imagesc(reshape(x_obs,N_row,N_col));colormap(gray);title(‘DMD aktif’);

[D_row,D_col]=size(D);

D_pos=-D;

%ubah L1 menjadi LP

A=zeros(D_row,2*D_col);

A(:,1:D_col)=D;

A(:,D_col+1:2*D_col)=D_pos;

f=ones(2*D_col,1);

lb=zeros(2*D_col,1);

f1=ones(2*D_col,1);

lb1=zeros(2*D_col,1);

A1=-A;

x_sub1=x_sub;

%*******************************

%* Jalankan algoritma simpleks *

%*******************************

alpha0=linprog(f,[],[],A,x_sub,lb); %get positive part

alpha1=linprog(f1,[],[],A1,x_sub1,lb1);%get negative part

alpha=alpha0-alpha1;

% Urutkan

[z,imax]=sort(abs(alpha(1:D_col)),’descend’);

x_abs=1:D_col;

figure(2); plot(x_abs, z,’b-’);

%Reconstrusi Citra

%===

x_hat=zeros(N,1);

for k=1:K_sparse;

x_hat=x_hat+alpha(imax(k))*PSI(:,imax(k));

end;

Ip = reshape(x_hat, N_row, N_col);

figure(3);imagesc(Ip);colormap(gray);title(‘citra rekonstruksi’);

Prinsip Ketidakpastian dan Kompresi Data (1)

Selama ini teknik kompresi data (sinyal) telah banyak dipakai untuk mengefisienkan penyimpanan dan pengiriman suara atau gambar pada peralatan komputer maupun ponsel. Biasanya orang mengenal kompresi sebagai metoda yang berkembang dalam bidang pengolahan sinyal dijital. Tanpa disadari, sebenarnya ada kaitan erat antara paradigma kompresi sinyal ini dengan prinsip fundamental yang dikenal dalam bidang Fisika Kuantum yaitu prinsip ketidakpastian.

Siapapun yang pernah belajar Fisika Modern tentu mengenal prinsip ketidakpastian Hiesenberg atau HUP (Heisenberg Uncertainty Principle). Prinsip ini menyatakan bahwa pengukuran posisi dan momentum suatu partikel tidak mungkin dibuat teliti kedua-duanya sekaligus, yakni

HUP adalah konsekuensi langsung dari sifat dualitas gelombang-partikel. Sebelum Teori Kuantum muncul, orang mengira bahwa ketelitian pengukuran posisi suatu partikel dan momentum (pada saat bersamaan) hanya akan dibatasi oleh ketelitian alat ukur. Prinsip diatas menyatakan larangan alam untuk mengetahui keduanya sekaligus secara teliti, sepresisi apapun alat ukur yang digunakan.

Selain relasi ketidakpastian dari momentum dengan posisi, HUP juga dapat dinyatakan sebagai hubungan ketidakpastian antara energi dengan waktu, yaitu

Menurut Max Planck, suatu foton yang memiliki frekuensi , akan memiliki energi sebesar dimana adalah konstanta Planck. Dengan demikian HUP untuk energi-waktu akan setara dengan

Pertidaksamaan diatas dinamakan juga prinsip ketidakpastian Weyl-Heisenberg (WHUP) yang berbunyi:

“a continuous-time signal cannot be simultaneously well-localized in both time and frequency“

Apakah bedanya dengan prinsip ketidakpastian yang sebelumnya? Jika HUP berbicara mengenai partikel dan gelombang sebagai objek fisik, WHUP menyatakan sifat umum dari prinsip ketidakpastian akan berlaku untuk sebarang sinyal atau fungsi kontinyu. Dalam bahasa pengolahan sinyal bisa ditafsirkan bahwa jika suatu sinyal s(t) terlokalisir dalam kawasan waktu maka transform Fourier dari sinyal ini,

akan tersebar dikawasan frekuensi dan demikian pula sebaliknya, suatu sinyal yang terlokalisir dikawasan frekuensi akan memiliki transform Fourier yang tersebar dikawasan waktu.

Gambar 1. Pasangan transform Fourier dari sinyal delta Dirac

Contoh dari sinyal kontinyu terlokalisir waktu adalah sinyal delta Dirac. Transform Fourier dari sinyal ini akan tersebar keseluruh rentang frekuensi, seperti yang diperlihatkan oleh Gambar 1. Karena delta Dirac adalah fungsi riil, maka haruslah transform Fouriernya bernilai kompleks. Gambar 1 sebelah kanan hanya menunjukkan nilai magnitude dari koefisien Fourier-nya untuk memperlihatkan sebaran komponen dominan dan terlihat disini bahwa semua komponen sama dominan-nya.

Fungsi delta dapat dilihat sebagai fungsi (distribusi) Gaussian dengan limit variansi mendekati nol. Kita akan memakai fungsi Gaussian ini untuk memperlihatkan (bukan membuktikan) bahwa implikasi tidak langsung dari WHUP adalah,

“sinyal kontinyu yang tersebar dalam kawasan waktu, akan terlokalisir pada kawasan frekuensi dan begitu pula sebaliknya.”

Pasangan Fourier dari fungsi Gaussian dengan mean nol dapat dituliskan sebagai berikut

Dari relasi diatas terlihat bahwa peningkatan variansi di kawasan waktu atau pelebaran sinyal Gaussian akan mengakibatkan penyempitan spektrum dikawasan frekuensi dan demikian pula sebaliknya, penyempitan sinyal di kawasan waktu akan melebarkan spektrum-nya. Dengan demikian, implikasi tak langsung dari WHUP tersebut menjadi jelas, khususnya untuk fungsi Gaussian. Fungsi-fungsi atau sinyal-sinyal waktu-kontinyu yang lain juga akan memiliki sifat yang demikian. Matlab Script berikut ini dapat dipakai untuk memperlihatkan fenomena ini.

%—————————————————–———-
% Script Matlab untuk memperlihatkan WHUP
%—————————————————————–
clear;clc;
t=-5:0.1:5; w=t; %koordinat waktu dan frekuensi.
sigma1=1; sigma2=0.5; sigma3=0.25;

%Hitung Fungsi Gaussian
g1=normal_baku(sigma1,t); g2=normal_baku(sigma2,t); g3=normal_baku(sigma3,t);

%Hitung transform Fourier-nya
G1=F_normal_baku(sigma1,w);
G2=F_normal_baku(sigma2,w);
G3=F_normal_baku(sigma3,w);

figure(1);plot(t,g1,t,g2,t,g3);
legend(’sigma’,’0.5sigma’,’0.25sigma’)
figure(2);plot(w,G1,w,G2,w,G3);
legend(’sigma’,’0.5sigma’,’0.25sigma’)

%—————————————————–
% Fungsi Gaussian
%—————————————————–
function [g]= normal_baku(sigma,t)
g=(1/sqrt(2*pi*sigma*sigma))*exp(-t.*t/(sigma*sigma));

%————————————————————————–
% Transform Fourier dari fungsi Gaussian
%————————————————————————–
function [G]= F_normal_baku(sigma,w)

G=exp(-w.*w*(sigma*sigma));

Keluaran dari program Matlab diatas diperlihatkan pada Gambar 2 dan Gambar 3. Disini jelas terlihat penyempitan fungsi Gauss pada kawasan waktu akan mengakibatkan pelebaran pada kawasan frekuensi.

Gambar 2. Fungsi Gaussian dengan berbagai nilai variansi

Gambar 3. Transform Fourier dari fungsi Gaussian pada Gambar 2

Sebelum diskusi dilanjutkan, akan terlebih dahulu disinggung sekilas mengenai suatu teorema yang sangat penting dalam transformasi sinyal, yaitu Teorema Parseval. Teorema ini pada prinsipnya adalah Hukum Kekekalan Energi dari sinyal, yang menyatakan bahwa pengukuran energi dalam kawasan waktu (misalnya dengan bantuan Osciloscope) dan pengukuran energi pada kawasan frekuensi (misalnya dengan Spectrum Analyzer) akan memberikan hasil yang sama. Teorema ini juga berimplikasi bahwa penghilangan komponen frekuensi tertentu dari sinyal akan mendistorsi sinyal kawasan waktunya, sepadan dengan magnitudo koefisien Fourier tersebut. Dengan demikian, hilangnya komponen dengan magnitudo rendah tidak akan mengubah bentuk sinyal semula terlalu banyak.

Untuk aplikasi pengolahan sinyal atau citra dijital, perhitungan dilakukan terhadap data dijital yang tak lain adalah fungsi-fungsi diskrit. Fungsi ini diperoleh dari fungsi kontinyu dengan melalui proses pencuplikan dan kuantisasi yang dapat dilakukan dengan sebuah ADC (Analog to Digital Converter).

• Apakah WHUP juga berlaku untuk fungsi diskrit?
• Apakah hubungan antara prinsip ketidakpastian, Teorema Parseval, dan kompresi data?
• Selain transformasi Fourier, adakah transformasi lain yang punya sifat mirip?
• Mengapa justru DCT dan bukannya DFT yang dipilih untuk kompresi JPEG?

Jawaban dan penjelasan lebih lanjut dari pertanyaan-pertanyaan diatas akan disajikan pada tulisan mendatang di blog ini.

Pipis Kecoa versus LCT (Limit Central Theorem)

Saya sering menyebut pengolahan citra dijital (Digital Image Processing) sebagai ilmu sihir. Betapa tidak, sebuah gambar yang tadinya penuh dengan kotoran—akibat dimakan waktu atau di-pipis-in kecoa, dapat kembalikan bersih (hampir) seperti sediakala melalui proses yang sederhana.

Cerita tadi menggambarkan salah satu teknik perbaikan citra (image restoration). Cara yang paling sederhana adalah dengan teknik perata-rataan dengan jendela yang berjalan (sliding window). Berdasarkan formulasi sistem linier, setiap proses pengolahan sinyal dapat dilihat dari dua kawasan, yaitu kawasan ruang-waktu dan kawasan frekuensi.

Dalam kawasan frekuensi, perata-rataan — yang sebenarnya adalah konvolusi citra dengan sebuah kernel, merupakan perkalian titik-demi-titik antara citra dengan (tanggapan impuls dari) tapis lolos rendah. Dengan demikian, proses ini akan menghilangkan komponen frekuensi menengah dan tinggi. Karena derau tersebar ke seluruh rentang frekuensi, penghilangan komponen frekuensi ini sekaligus menghilangkan derau. Penalaran dari kawasan frekuensi dapat secara mudah dipakai untuk menjelaskan hilangnya sebagian besar derau dari citra yang terkotori.

Bagi saya, yang lebih sulit adalah menjelaskan proses ini dari kawasan ruang-waktu. Ini sangat penting untuk menjelaskan proses hilangnya derau kepada mahasiswa. Setelah cukup lama bersarang di alam bawah sadar, teka-teki ini muncul kembali sewaktu memberikan kuliah dasar-dasar statistik dalam pengolahan citra. Sub topik kuliah hari itu adalah sebaran normal dan LCT (Limit Central Theorem).

LCT adalah teorema fundamental kedua dalam teori peluang setelah LLN (the Law of Large Numbers). Dalam kalimat yang agak bebas, LCT mengatakan: jika X₁, X₂, …, X_n adalah n-buah peubah acak saling bebas dengan sebaran identik (i.i.d.-independent and identically distributed) dengan mean m dan variansi s² yang berhingga, ketika jumlah titik cuplikan meningkat, maka nilai rata-rata cuplikan akan semakin tersebar secara Gaussian dengan variansi s²/N.

Pada penapisan dengan jendela berjalan, kita memilih satu piksel untuk diubah nilainya. Kemudian, piksel-piksel tetangganya yang masuk dalam jendela dihitung nilai rata-ratanya. Nilai rata-rata inilah yang selanjutnya dipakai untuk memperbarui nilai piksel yang ditinjau.

Hasil terakhir dalam LCT, yakni mengecilnya variansi, menjelaskan mengapa nilai derajat keabuan dari satu piksel ke piksel tetangganya hampir seragam atau tidak berubah banyak. Ini akibat turunnya variansi dari s² menjadi s²/N. Artinya gambar menjadi halus, dan derau berkurang. Dengan demikian, teka-teki hilangnya derau terjawab sudah.

Kekurangan dari teknik perata-rataan ini adalah citra akan menjadi kabur, detail dari citra bisa lenyap. Ada berbagai cara lain yang dapat dipakai untuk perbaikan citra, misalnya menggantikan mean dengan median. Cara yang kedua ini malahan lebih kokoh terhadap derau impuls. Namun demikian, penapisan median tidak dapat dilakukan dalam kawasan frekuensi. Jika dibutuhkan proses yang cepat, penapisan kawasan fekuensi menjadi pilihan karena tersedianya algoritma cepat FFT (Fast Fourier Transform). Lahirnya beberapa teknik alternatif yang tidak bisa dijelaskan dengan teori sistim linier, seperti penapisan median, dapat mengubah pendekatan pengolahan sinyal menjadi proses komputasi. Setelah masa DSP (Digital Signal Processing) berlalu, barangkali akan lahir bidang ilmu baru– Algorithmic Signal Processing.