Chaotic Pearls

Lamunan dari seberang GKU lama …

Archive for the ‘Imaging’ Category

Pencitraan Sintesis (3)

Posted by suksmono on June 26, 2008

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 \omega  tertentu dari gelombang, maka fasor medan listrik yang bernilai kompleks ini adalah {E}_{\omega} (R)  . 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

E(r) = \int P(R,t) E(R) dS  (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

E(r) = \int \frac{E(R)exp(j2 \pi \omega |R-r|)/c)}{|R-r|} dS  (2)

dimana c adalah kecepatan rambat gelombang EM (3×108 m/s). Persamaan (2) menyatakan bentuk umum komponen kuasi-monokromatik dari medan listrik pada frekuensi \omega  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 r1, sedangkan yang kedua berada pada posisi r2. Korelasi antar keduanya dinyatakan sebagai visibility

V({r}_{1}, {r}_{2}) = \left< E({r}_{1}) {E}^{*} ({r}_{2}) \right>  (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 \left< E({R}_{1}) E({R}_{2}) \right>  akan berharga nol jika {R}_{1} \neq {R}_{2}  , sehingga persamaan (4) dapat disederhanakan menjadi

V({r}_{1},{r}_{2}) = \int \left<{|E(R)|}^{2} \right> |{R}^{2}| \frac{exp(j 2\pi \omega |R-{r}_{1}|)/c}{|R-{r}_{1}|}\frac{exp(-j2\pi \omega |R-{r}_{2}|)/c}{|R-{r}_{2}|} dS  (5)

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

|R-r| = \sqrt{{|R|}^{2} + {|r|}^{2} -2r.R} \approx |R|-\frac{r.R}{|R|} = |R|-r.s

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|>2 |R|2. Dengan menggantikan integrasi luas dari objek dS dengan sudut d \Omega  , maka persamaan (5) akan menjadi

V({r}_{1},{r}_{2}) = \int I(s) exp(-j2 \pi \omega s.({r}_{1}-{r}_{2})/c) d \Omega  (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 r1 maupun r2, 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= s0 + s, dimana s0 ortonormal dan tegaklurus terhadap s. Suatu sistem koordinat baru dimana s0 = (0,0,1) dapat dipakai, sehingga posisi relatif kedua stasiun penerima itu menjadi {r}_{1}-{r}_{2}=c(u,v,w)/ \omega  . Dengan demikian, vektor s akan menjadi s = (x,y,1) dan persamaan (6) dapat dinyatakan dalam koordinat citra berbentuk bidang datar (x,y)

V(u,v,w) = {e}^{-j2 \pi w} \int \int I(x,y) exp(-j2 \pi (ux+vy)) dx dy  (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

V(u,v) = \int \int I(x,y) exp(-j2 \pi (ux+vy)) dx dy  [8]

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

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

I(x,y) = \int \int V(u,v) exp(j2 \pi (ux+vy)) du dv  (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

S(u,v) = \sum_{k} {w}_{k} \delta (u-{u}_{k}) \delta (v-{v}_{k})  (10)

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

\hat{I} (x,y) = \sum_{k} {w}_{k} V({u}_{k},{v}_{k}) exp(j2 \pi (ux+vy))  (11)

Pembobot wk pada persamaan (11) dapat diatur agar terjadi perbaikan pada gambar kotor \hat{I} (x,y)  , 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.

Posted in Image Processing, Imaging, Nobel in Imaging, Uncategorized | 1 Comment »

Pencitraan Sintesis (2)

Posted by suksmono on June 23, 2008

Membuat teleskop radio dengan resolusi 100 kali dari teleskop radio Reber memerlukan aperture sebesar satu kilometer. Cara langsung dengan membuat antena sebesar ini bukanlah perkara yang mudah, tapi akhirnya dapat dicapai dengan kecerdikan. Alih-alih membuat satu buah antena raksasa, para Fisikawan eksperimental mempergunakan teknik interferometri.

1. Interferometer Optik
Secara singkat, interferensi adalah peristiwa penggabungan beberapa gelombang, sedangkan interferometri berarti pengukuran efek interferensi ini. Albert A. Michelson merupakan Fisikawan yang dianggap sebagai pelopor penggunaan efek interferensi didalam pengukuran.

Pada awal karirnya, riset yang paling sering dilakukan oleh Michelson adalah pengukuran kecepatan cahaya. Terobosan pertama dari hasil penelitian ini terjadi ketika pada tahun 1887, bersama dengan Morley, Michelson memperlihatkan sifat isotropik perambatan cahaya terhadap gerakan bumi. Ini berarti bahwa medium hipotetik bernama ether itu tidak ada dan fakta ini mendukung berlakunya Relativitas Khusus. Terobosan selanjunya terjadi ketika pada tahun 1920-1921 Michelson berhasil mengukur diameter bintang.

Bintang sebenarnya adalah matahari yang letaknya sangat jauh–sedemikian jauhnya hingga hanya nampak sebagai benda titik. Ini akibat mata manusia, atau teleskop optik saat itu, memiliki daya pisah yang rendah. Supaya diameter bintang dapat diukur, maka resolusi alat ukur harus ditingkatkan, misalnya dengan memperbesar ukuran lensa/cermin utama pada teleskop. Pembuatan cermin/lensa berukuran besar dengan permukaan halus bukanlah hal yang mudah.

Gambar 1. Teleskop dengan interferometer untuk memperbesar apertur

Michelson menemukan pemecahan masalah ini dengan cara memasang dua cermin berjauhan sebagai pemantul cahaya bintang sebelum masuk ke teleskop, seperti diperlihatkan pada Gb.1. Kedua berkas dari cermin ini selanjutnya dipadukan sehingga membentuk gambaran interferensi atau citra (kompleks) visibility. Dengan cara ini, apertur efektif dari teleskop dapat ditingkatkan dan Michelson-Pease dapat memperkirakan diamater sebuah bintang raksasa merah alpha-orionis (betelguese/ yad al-jawza) yang letaknya 427 tahun cahaya dari bumi, yaitu sebesar 0.047 arcseconds.

2. Interferometer Radio

Secara prinsipil, teknik interferometri radio mirip seperti pada interferometer optik Michelson. Fungsi dari cermin dan pemadu berkas digantikan oleh antena yang tersebar. Setiap antena seolah-olah merupakan cuplikan ruang (spasial) dari sebuah antena raksasa.

Gambar 2. Ekivalensi antenna (apertur) tunggal yang sangat besar dengan hasil sintesis apertur antena tersebar di cakupan yang sama (diadopsi dari [2]).

Pada dasarnya, setiap antena dalam teleskop interferometri mengukur koefisien transformasi Fourier dimensi dua pada kawasan frekuensi ruang, atau disebut juga kawasan k atau kawasan u-v, yang bernilai kompleks. Dengan demikian, setiap stasiun pengamat harus mampu mengukur baik magnitudo maupun fasa gelombang yang jatuh pada saat bersamaan di setiap titik. Gelombang yang datang di setiap antena selanjutnya disalurkan ke suatu tempat untuk dikorelasikan satu dengan yang lain. Untuk apertur yang sangat besar, cara ini tidak mudah dilakukan karena memerlukan saluran transmisi yang sangat panjang. Cara yang lain adalah mencatat cuplikan gelombang disetiap antena dan juga waktu kedatangannya. Cara seperti ini membutuhkan pencatat waktu yang sangat presisi sehingga digunakanlah jam atom.

Gambar 3. Rotasi bumi mengakibatkan sapuan berbentuk elips pada bidang u-v

Peningkatan apertur lebih lanjut dapat dilakukan, tidak dengan memperluas daerah pencuplikan, tetapi dengan mengikuti perputaran bumi. Dua buah antena yang diam dalam arah Utara-Selatan akan membentuk apertur sangat besar jika pengamatan dilakukan sehari penuh karena antena ini ”dibawa” bumi berputar pada porosnya, seperti diperlihatkan pada Gb.3. Dengan asumsi letak benda yang diamati sangat jauh, rotasi bumi akan memberikan sapuan berbentuk elips pada bidang u-v. Data visibility yang diperoleh selanjutnya dipakai untuk membentuk citra benda langit pemancar radio dengan inversi transform Fourier. Karena tidak seluruh bidang u-v terisi, harus ada pengolahan data dengan cara tertentu untuk melengkapi bidang ini, supaya hasil inversi memberikan citra yang berkualitas tinggi.

Daftar Pustaka

  1. A.A. Michelson and F.G. Pease, “Measurement of the diameter Alpha-Orionis by the interferometer,” Proc. Nat. Acad. of Sci., Vol.7, 1921, pp.143-146.
  2. K.S. Thorne, Black Holes and Time Warps: Einstein’s Outrageous Legacy, W.W. Norton & Company.

Posted in Imaging, Nobel in Imaging, Uncategorized | Leave a Comment »

Pencitraan Sintesis (1)

Posted by suksmono on June 22, 2008

Pada tahun 1930-an, hubungan telepon menggunakan radio telah mulai dipakai secara komersil. Bell Telephone adalah salah satu operator besar yang merintis bisnis dalam bidang ini, dengan menawarkan sambungan telepon radio antar benua.

Disamping kemampuannya menghubungkan dua orang yang ingin berkomunikasi yang letaknya sangat berjauhan, komunikasi radio menghadapi masalah teknis cukup penting akibat gangguan dari sumber yang tak diketahui. Bell Telephone Lab di Homdel (NJ) menugaskan Karl Jansky, salah seorang insinyur muda yang baru lulus dari perguruan tinggi untuk mencari sumber asal pengganggu ini.

Jansky terlebih dahulu ingin tahu persis dari mana arah gangguan ini berasal. Oleh karena itu, dibuatlah suatu susunan antena (antenna array) dari pipa-pipa logam. Pada tahun 1932, salah satu sumber utama gangguan ini akhirnya berhasil dikenali, yaitu badai petir. Namun demikian,bahkan ketika badai petir sudah mereda, gangguan ini tidak hilang sepenuhnya. Dengan mengarahkan susunan antenna ke berbagai penjuru, akhirnya Jansky menyimpulkan bahwa salah satu sumbernya ada di pusat galaksi kita, yaitu Bima Sakti.

Meskipun penemuan ini dipublikasikan secara luas, hanya sedikit peneliti yang menanggapi. Diantara yang sedikit ini adalah beberapa astronom dari Harvard. Mereka terheran-heran akan adanya sumber radiasi sekuat itu mengingat jaraknya yang begitu jauh. Mereka sama sekali tidak meragukan hasil pengamatan Jansky, masalahnya justru terletak pada teori Astrofisika saat itu.

Ditengah skeptisme kalangan ilmuwan; Grote Reber–seorang amatir radio dengan call sign W9GFZ tertarik pada bidang baru ini. Dengan menggunakan biaya pribadi tanpa bantuan teknis yang memadai, Reber berhasil membangun sebuah teleskop radio di belakang rumah ibunya. Teleskop radio pertama di dunia dengan diameter sekitar 9 meter ini berhasil digunakan untuk memetakan sumber radio ruang angkasa.

Peta Reber menunjukkan bahwa sumber radiasi bukan hanya dari pusat Bima Sakti, tetapi juga berasal dari arah konstelasi Cygnus dan Cassiopea. Hasil penelitiannya dilaporkannya kepada Chandrasekar, editor dari Astrophysical Journal. Meskipun semula banyak peneliti yang ragu-ragu, mengingat Reber bukanlah peneliti profesional, tulisan ini akhirnya diterbitkan setelah mereka melihat secara langsung teleskop radio buatan Reber.

Peta yang telah dibuat Reber masih kurang memadai untuk keperluan penentuan posisi benda langit pembangkit gelombang radio ini. Resolusi dari teleskop ini perlu ditingkatkan seratus kali lipat agar letak dan sumber radiasi dapat ditentukan dan dapat dipadukan dengan hasil pengamatan teleskop optik.

Penulis (jaket oranye) didepan teleskop radio 45 m di Nobeyama

Pada pencitraan, cara untuk  meningkatkan resolusi adalah dengan memperbesar ukuran apertur. Pada teleskop biasa (optik), hal ini dilakukan dengan memperbesar diameter lensa atau cermin-cekung. Panjang gelombang daerah optik sekitar sepersejuta meter, sedangkan panjang gelombang radio yang diamati berorde meter. Untuk dapat dipakai dalam penelitian astrofisika, teleskop radio memerlukan parabola dengan diameter sedikitnya satu kilometer. Ini tentu bukan hal yang mudah ditangani.

Daftar Pustaka

  1. K.S. Thorne, Black Holes and Time Warps-Einstein’s Outrageous Legacy, W.W. Norton & Company
  2. M. Ryle, “Radio telescopes of large resolving power,” Nobel Lecture, 1974.

Posted in Imaging, Nobel in Imaging, Uncategorized | 1 Comment »

Compressive SFCW-GPR/ Pencitraan Kompresif (3)

Posted by suksmono on June 13, 2008

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.

Posted in Applications, Compressive Imaging, Compressive Sensing, Image Processing, Imaging, Uncategorized | 3 Comments »

Radar Frekuensi Kontinyu Penembus Permukaan

Posted by suksmono on June 9, 2008

Pada dasarnya, radar bekerja dengan cara mengirimkan impuls gelombang elektromagnetik (EM) dan kemudian menangkap gema-nya. Pengiriman impuls dapat dilakukan langsung dalam kawasan waktu atau secara tak langsung dengan mensintesa gema radar pada kawasan frekuensi. Tulisan ini menjelaskan prinsip dasar dari radar penembus permukaan dengan menggunakan teknik sintesa frekuensi atau SFCW-GPR (Stepped-Frequency Continuous Wave-Ground Penetrating Radar). Hubungan impuls dalam kawasan waktu dan frekuensi melalui transformasi Fourier merupakan dasar dalam merancang radar jenis ini.

1. Pendahuluan

Radar penembus permukaan (GPR/Ground Penetrating Radar) adalah suatu alat pencitra gelombang elektromagnetik (EM) yang mampu melihat benda-benda di bawah permukaan tanah atau dibalik dinding. Desain GPR sangat bergantung pada tujuan aplikasi; GPR untuk mencitra benda dibawah permukaan dalam (mis. air) akan berbeda dengan yang akan dipakai untuk melihat benda di permukaan dangkal (mis. ranjau). Disamping itu, kebutuhan resolusi juga akan menentukan persyaratan desain dari radar.

Secara prinsipil ada dua macam teknologi pancaran radiasi EM yang bisa dipakai untuk membuat radar, yaitu pancaran impuls dan pancaran gelombang kontinyu. Tulisan ini akan menjelaskan GPR yang dibuat berdasarkan prinsip kedua, yakni radar dengan teknik SFCW (Stepped-Frequency Continuous Wave).

Untuk kepentingan praktis tertentu yang memerlukan resolusi tinggi, durasi impuls (efektif) yang diperlukan haruslah sangat singkat (sampai orde sub-nano detik). Sumber impuls yang demikian sulit untuk dibuat/ dicari, sangat mahal, memerlukan rangkaian RF yang tidak sederhana, dan perlu penguat daya yang tidak mudah dibeli. Masalah ini dapat diatasi dengan teknik SFCW, meskipun akan meningkatkan waktu akuisisi data dan pengolahan sinyal. Teknik multipleksing ruang-frekuensi, dimana beberapa frekuensi dapat dipancarkan secara serempak pada beberapa titik koordinat spasial sekaligus dapat dipakai sebagai pilihan dalam mengatasi masalah kecepatan pencitraan. Disamping itu, perkembangan mutakhir dari teknik pencitraan kompresif juga menjanjikan solusi bagi masalah yang sangat mendasar ini.

2. Prinsip Kerja dan Geometri Pencitraan GPR

Citra benda yang diperoleh sistem radar pada dasarnya merupakan sekumpulan pantulan gelombang EM sebagai fungsi dari posisi dan sifat benda pemantul. Untuk sebuah radar ideal, impuls radar dapat dianggap sebagai suatu fungsi delta Dirac \delta (t)  , seperti yang dilukiskan pada Gb.1: (a) skema pencitraan GPR dan (b) deretan gema impuls yang diterima sistem radar. Bagian pengirim (Tx) memancarkan impuls \delta (t)  , kemudian impuls akan mengenai objek dimana sebagian gelombang akan dipantulkan kembali dan akhirnya sampai ke sistem penerima (Rx). Sinyal yang diterima ini bisa dinyatakan sebagai \delta (t-\Delta t)  . Karena kecepatan gelombang EM dalam medium tertentu telah diketahui, maka jarak antara antenna ke benda dapat dihitung berdasarkan waktu tunda.

Gambar 1.a

(a)

Gambar 1.b

(b)

Gb.1 Prinsip radar: (a) geometri pencitraan GPR dan (b) A-scan ideal

Data pantulan untuk satu titik pencitraan akan berupa suatu fungsi waktu yang menyatakan letak dan kekuatan pemantul sepanjang perjalanan gelombang. Bentuk data radar yang paling mendasar ini disebut sebagai sapuan jenis-A atau A-scan. Untuk kasus impuls Dirac, pantulan ideal akan berupa impuls Dirac yang tertunda dan dilemahkan, seperti yang diperlihatkan pada Gb.1.b. Jika penyapuan dilakukan sepanjang suatu garis lurus, akan diperoleh sekumpulan A-scan yang menyatakan letak-letak pemantul pada kedalaman tertentu sepanjang garis. Untuk pemantuk berupa benda titik, profil pantulan akan berbentuk hiperbola tertelungkup, seperti diperlihatkan pada Gb.2. Hasil penyapuan yang demikian disebut sebagai B-scan. Jika B-Scan dilakukan berkali-kali sehingga meliputi suatu bidang datar, maka hasilnya adalah gambaran dimensi tiga yang disebut sebagai C-scan.

Gambar 2

Gb.2 Kurva hiperbola yang dihasilkan oleh B-scan

Pada kenyataanya, impuls delta Dirac tidak mungkin dapat diperoleh karena beberapa alasan. Yang pertama adalah keterbatasan rentang frekuensi kerja dari peralatan. Sebagai akibatnya akan terjadi pelebaran impuls sehingga membentuk suatu fungsi Gaussian. Disamping itu, antena juga berfungsi sebagai diferensiator. Ini bisa difahami karena komponen arus searah (DC) dari sinyal tidak mungkin disalurkan antenna tanpa kontak langsung dengan medium dan tanpa pembentukan rangkaian tertutup. Dengan demikian, keluaran dari antenna pemancar akan berupa turunan pertama dari fungsi Gaussian, yang biasa dikenal sebagai impuls monocycle. Sinyal monocycle dalam kawasan waktu dan frekuensi dilukiskan pada Gb.3.

Gambar 3

Gb.3 Hubungan waktu-frekuensi pada sinyal monocycle

Berdasarkan penjelasan sebelumnya dan mengacu pada Gb.3, terlihat bahwa sebuah impuls dapat dibangkitkan melalui dua cara, yaitu: (1) pembangkitan pada kawasan waktu, dan (2) pembangkitan pada kawasan frekuensi. Cara (1) disebut juga cara langsung, sedangkan cara (2) adalah cara yang tidak langsung dengan sintesa tanggapan frekuensi. Prinsip yang dipakai pada cara (2) ini adalah dualitas dari sinyal dalam kawasan waktu dan frekuensi melalui transformasi Fourier.

Tinjau suatu fungsi waktu atau sinyal kontinyu s(t). Penguraian sinyal ini kedalam komponen frekuensi dilakukan dengan transformasi Fourier sebagai berikut:

S(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} s(t) exp(-j \omega t) dt (1.a)

s(t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} S(\omega) exp(j \omega t) d\omega (1.b)

dimana \omega=2 \pi f  adalah frekuensi (angular), sedangkan j=\sqrt{-1}  adalah bilangan imajiner. Persamaan (1.a) disebut sebagai persamaan analisis, sedangkan persamaan (1.b) adalah persamaan sintesis.

Teknik SFCW berhubungan langsung dengan sintesis Fourier. Dalam hal ini, sinyal s diperoleh dengan terlebih dahulu mengukur nilai S sebagai fungsi frekuensi. Karena pada umumnya S bernilai kompleks, sintesis s memerlukan magnitudo maupun fasa dari S. Kebutuhan data yang demikian berakibat langsung pada sisi implementasi, yakni, sistem deteksi sinyal harus dibuat sedemikian hingga komponen magnitudo dan fasa, atau bagian riil dan imajiner bisa didapatkan. Perangkat yang mampu melakukan pengukuran sinyal kompleks ini tak lain adalah pen-demodulasi kuadratur (I/Q demodulator). Setelah S diperoleh, maka secara prinsipil sinyal s(t), yaitu sinyal A-scan, dapat ditentukan dari proses inversi Fourier dari persamaan (1.b).

3. Teknik Sintesa Frekuensi untuk Radar

Pada kenyataanya, nilai S tidak bisa diperoleh untuk seluruh frekuensi dalam kawasan waktu-kontinyu, melainkan hanya untuk sejumlah berhingga dari titik pengamatan. Ini berarti bahwa data S hanya dapat diperoleh untuk sekumpulan frekuensi yang berubah secara diskrit.

Dengan demikian, proses yang sebenarnya lebih sesuai untuk memodelkan sistem radar SFCW adalah analisis dan sintesis sinyal dengan transformasi Fourier diskrit (DFT) sbb:

{S}_{k}=\frac{1}{N}\sum_{n=0}^{N-1} {s}_{n} exp \left(- \frac{2\pi jn}{N}k \right)   (2.a)

{s}_{n}=\frac{1}{N}\sum_{k=0}^{N-1} {S}_{k} exp \left( \frac{2\pi jk}{N}n \right)   (2.b)

Sintesa untuk membentuk impuls dilakukan dengan persamaan (2.b), sedangkan koefisien Fourier Sk yang bernilai kompleks diperoleh dari pengukuran magnitudo dan fasa gelombang pantul. Indeks k menyatakan urutan frekuensi ke-k dari sinyal. Dengan demikian, pada saat akuisisi data dilakukan, sekumpulan sinyal dengan frekuensi tertentu dipancarkan kemudian tanggapannya diukur untuk mendapatkan estimasi dari koefisien Fourier .

Gambar 4

Gb.4 Konstruksi Dasar Sistem SFCW-GPR [1]

Diagram blok sederhana dari sistem SFCW-GPR diperlihatkan pada Gb.4. Pada gambar tersebut, sekumpulan gelombang dengan frekuensi tertentu {\omega}_{k}  dibangkitkan oleh frequency synthesizer dan dipancarkan secara berurutan melalui antena UWB (ultrawideband). Penerima akan menangkap pantulan gelombang melalui antena penerima untuk di-demodulasi dengan pen-demodulator kuadratur. Hasilnya adalah sinyal Ik (inphase) dan sinyal Qk (quadrature) yang secara bersama-sama membentuk koefisien Fourier kompleks

{S}_{k}=S({\omega}_{k}) = I({\omega}_{k}) + jQ({\omega}_{k})  (3)

Setelah dicuplik dengan ADC, pengolah sinyal akan merekonstruksi sinyal kawasan waktu sn dengan algoritma IDFT.

4. Pertimbangan Desain SFCW-GPR

Ada beberapa hal penting yang perlu dipertimbangkan dalam merancang sistem SFCW-GPR, diantaranya adalah: lama pancaran untuk satu frekuensi tertentu, frekuensi pusat, rentang frekuensi kerja, jumlah langkah-frekuensi yang digunakan (N), separasi antar frekuensi, tingkat kuantisasi ADC, dan ketersediaan komponen. Berikut ini penjelasan singkat dari masing-masing aspek:

  • Lama pancaran dan jumlah frekuensi N: kedua hal ini akan menentukan kecepatan pencitraan radar. Batas minimum-nya adalah settling-time dari pensintesis dan kecepatan ADC, sedangkan batas maksimumnya berhubungan dengan lama waktu pencitraan yang diinginkan. Jika Dt terlalu kecil, ada kemungkinan pembangkit frekuensi belum mencapai keadaan tunak. Disisi lain, jika nilainya terlalu besar, maka waktu untuk mendapatkan data akan terlalu lama. Hal terakhir harus diperhitungkan untuk menentukan kecepatan penyapuan radar.
  • Rentang frekuensi kerja dan pusat frekuensi: frekuensi kerja akan sangat mempengaruhi aplikasi dari GPR. Radar dengan frekuensi pusat yang rendah akan mampu menembus permukaan yang lebih dalam. Sebaliknya, radar dengan frekuensi pusat yang tinggi lebih cocok untuk pencitraan dangkal. Lebar pita (BW) berhubungan langsung dengan resolusi, semakin tinggi BW semakin besar pula resolusinya.
  • Tingkat kuantisasi ADC: hal ini akan menentukan lebar dinamika sinyal yang diperoleh. Sebagai gambaran, ADC 8 bit dapat dipakai untuk mengkuantisasi sinyal dengan lebar dinamika 41 dB, 12 bit bisa mencapai 65 dB, 14 bit mampu mengkuantisasi sinyal 77 dB, sedangkan untuk sinyal dengan lebar dinamik 89 dB memerlukan ADC 16 bit.
  • Ketersediaan komponen: meskipun dari perhitungan bisa diperoleh disain yang bagus, implementasi menjadi perangkat keras sangat ditentukan oleh ketersediaan komponen di pasaran.
  • Separasi antar frekuensi akan menentukan ambiguitas jangkauan (range ambiguity). Ini setara dengan PRF (pulse repetition frequency) pada radar impuls.

5. Eksperimen SFCW-GPR Dengan VNA (Vector Network Analyzer)

Pada dasarnya sebuah SFCW-GPR adalah pencitra kawasan Fourier. Sebuah citra GPR tersusun atas A-scan yang didapat dari inversi data kawasan frekuensi, yaitu koefisien Fourier kompleks untuk setiap frekuensi. Dengan demikian, perangkat SFCW-GPR bukan hanya dituntut untuk mampu mendeteksi magnitudo, tetapi juga fasa dari gelombang yang datang. Alat ukur yang bisa dipakai untuk keperluan ini adalah VNA (Vector Network Analyzer).

Konfigurasi bistatik seperti pada Gb.4 dapat dipakai untuk melakukan pencitraan, antena pertama berfungsi sebagai pemancar dan yang kedua sebagai penerima. Pada GPR-test range, antena dipasang pada penyapu elektro-mekanik tiga dimensi, tetapi penyapuan dapat juga dilakukan secara jika peralatan yang demikian tidak tersedia. VNA diatur untuk mengukur S21 dan mengeluarkan data kompleks. Sebuah komputer dapat dihubungkan ke VNA untuk mengambil dan merekonstruksi A-scan secara langsung. Jika antar-muka VNA dan komputer tidak memungkinkan pengambilan langsung, data bisa dipindahkan secara manual melalui floppy disk.

Sebelum scanning VNA dijalankan, rentang frekuensi kerja harus terlebih dahulu ditentukan. Hal ini bisa diperkirakan berdasarkan impuls atau resolusi yang diinginkan. Lebar pita yang dipakai harus berbanding terbalik dengan lebar impuls.

Setelah data S21 berhasil dikumpulkan, langkah selanjutnya adalah melakukan inversi Fourier dengan persamaan (2.b). Hal penting yang perlu dicatat, A-scan adalah sebuah sinyal riil. Dengan demikian hasil transformasi Fouriernya adalah sinyal kompleks yang bersifat konjugasi simetrik. Data dari S21 dari VNA hanya berisi setengah bagian dari koefisien Fourier, setengah bagian sisanya diisi dengan memanfaatkan sifat konjugasi simetrik ini. Contoh pengisian data bagian konjugasi simetrik diperlihatkan pada program Matlab diakhir tulisan.

Gambar 5

Gambar 5. A-scan hasil pencitraan SFCW-GPR dengan VNA

Contoh hasil rekonstruksi A-scan untuk pengujian di test-range diperlihatkan pada Gb.5. Disini terlihat adanya 2 buah pantulan dominan, yang pertama berasal dari pancaran langsung yang diterima antena loop yang terbenam dalam bak pasir, sedangkan yang kedua berasal dari pantulan benda. Konfigurasi pencitraan normal dengan kedua antena berada diatas permukaan akan menghasilkan A-scan yang bentuknya berbeda. Biasanya pantulan kuat terjadi di sekitar t=0 akibat pantulan langsung dari tanah atau kebocoran dari pemancar, sedangkan pantulan lain berasal dari benda didalam tanah. Karena jarak antena ke tanah maupun jarak pemancar ke penerima tetap, pantulan kuat ini dapat dihilangkan dengan penapisan sederhana.

Daftar Pustaka

  1. L.P. Ligthart, Course Handout: A 3-Days Short Course on GPR, IRCTR-Delft.
  2. A.B. Suksmono, A. Pramudita, E. Bharata, A.A. Lestari, and N. Rachmana, “A Novel Design of the Stepped Frequency Continuous Wave Radars Based on Non-Uniform Frequency Sampling Scheme”, Proc. of ICEEI 2007, Bandung, Indonesia, pp.270-273.

Kode MAtlab

% *****************************************************
% * display A-scan: menampilkan kurva A-scan dari data S21
% ****************************************************
path(path, ‘Data\’);
% baca data, format: FREQ. ## REAL-PART ## IMAGINARY-PART
data=load(‘VNA_TRACE.txt’);
N=length(data); NRows = 512; %extend seq. length
delta_t = 1/NRows; t=delta_t*(1:NRows);
dmy=zeros(NRows,1);
S21 = data(:,2) + data(:,3).*i; % Raw complex data
% filling the conjugate symmetric part of the data
trc_N = length(S21);
dmy(1) = 0; dmy(2:trc_N+1) = S21; dmy(trc_N+2:512) = 0;
dmy(NRows-trc_N+1:NRows) = conj(S21(trc_N:-1:1));
% inversi dan normalisasi
f=real(ifft(dmy)); f=f-mean(f); x=f/norm(f(:));
figure(1);plot(t,x); title(‘A-scan’); xlabel(‘time’);ylabel(‘amplitude’);

Posted in Imaging, Uncategorized | 2 Comments »

Pencitraan Kompresif (1.b)

Posted by suksmono on June 6, 2008

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.

haar rekonstruksi

(a)

Rekonstruksi Haar 50 Koefisien

(b)

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

Haar Coefficient H

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.

Original Checker Board

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.

Checker Coefficient by Haar

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 1.25 \sqrt{N}  . 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.

Observed Checker Board

(a)

Checker Reconstructed by Haar

(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;

Posted in Applications, Compressive Imaging, Compressive Sensing, Image Processing, Imaging, Uncategorized | Leave a Comment »

Pencitraan Kompresif (2)

Posted by suksmono on June 4, 2008

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

 \hat{s}= ( \Phi \Psi )  S  (1)

Setiap baris pada \Phi  hanya berisi satu buah angka bernilai satu, sedangkan sisanya nol. Perkalian dengan \Psi  mengakibatkan seleksi vektor baris pada indeks yang bersangkutan dengan letak kolom elemen aktif \Phi  tersebut. Sebagai akibatnya, vektor basis dalam \Psi  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.

diagram TV

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

Mengacu pada diagram ranah-jamak, sinyal s terlebih dahulu ditransformasikan ke kawasan frekuensi menjadi \zeta  . Pencuplikan kompresif kawasan ini dengan transformasi pengukur \Phi  akan menghasilkan cuplikan kompresif \hat{\zeta}  . 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 \hat{s}  . Selanjutnya yang akan diminimumkan adalah nilai variansi total:

|S|=|\nabla \hat{s}| = \sqrt{ \sum_{i} {({\nabla}_{i} s)}^{2}}  (2)

Basis pursuit untuk kasus ini dapat dirumuskan sebagai

(P1): min |S| s.t. \hat{s} = {F}^{-1} \hat{\zeta}  (3)

Optimisasi diatas mengatakan:

“pilih \hat{\zeta}  yang inverse transform Fouriernya, yaitu \hat{s}  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 l1-magic dari Caltech. Program utama sudah sedikit dimodifikasi supaya sesuai dengan keperluan, seperti diperlihatkan pada lampiran di akhir tulisan ini.

Original H-image

Gambar 2. Citra asli

H-observed-TV

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.

citra rekonstruks TV

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

Posted in Applications, Compressive Imaging, Compressive Sensing, Image Processing, Imaging, Uncategorized | Leave a Comment »

Pencitraan Kompresif (1)

Posted by suksmono on June 2, 2008

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.

one pixel camera

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 \Phi  . Jika N jumlah total dari cermin DMD, sedangkan M adalah cermin aktif dalam satu kali pengambilan gambar, maka \Phi  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:

\hat{s}= \Phi s = [ \Phi \Psi] S  (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 \psi  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 \Phi  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.

Original H-image

Gambar 2. Citra sederhana dalam pencitraan kompresif

measured image

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 \hat{s}  menyatakan hasil pengamatan, {\Psi}_{k}  vektor basis penjarang ke-k, maka pencarian basis dapat dirumuskan sebagai berikut:

P1: min ||S||1 s.t. \hat{s} = \sum_{k} {S}_{k} {\Psi}_{k}  (2)

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

DCT rekonstruksi

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 Sk yang terambil pada pencarian basis.

Ordered-magnitude of H

Gambar 5. Magnitudo koefisien terpilih yang diurutkan

Gambar 5 memperlihatkan kurva magnitudo Sk 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’);

Posted in Applications, Compressive Imaging, Compressive Sensing, Image Processing, Imaging, Uncategorized | Leave a Comment »