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’);