“Nerusin Yogi”: Menghitung Porositas pada Sayatan Tipis Karbonat menggunakan OpenCV di Python
(PS: bisa digunakan untuk menghitung komponen lain selain porositas!)
Pesan sponsor: seluruh konten dan kode yang saya tulis akan selalu gratis dan dapat digunakan oleh siapa saja. Apabila Anda merasa terbantu, dukung saya melalui laman Saweria berikut. 10% dari dukungan Anda akan saya sumbangkan melalui Kitabisa. Terima kasih :)
Selamat datang di tulisan yang mulai telat terbit ini! Ingin rasanya menyalahkan wabah COVID-19, tapi ‘kan wabahnya malah memberikan banyak waktu luang ya? Ya, pada akhirnya memang hanya bisa menyalahkan diri sendiri…
Kali ini, saya membuka kembali foto-foto sayatan tipis petrografi dari penelitian Tugas Akhir (TA) saya di tahun 2013 (tuaaaa!), dan berusaha menghitung porositas secara kuantitatif menggunakan sebuah library untuk computer vision yang ngetop yaitu OpenCV.
Untuk pembaca yang memiliki latar belakang non-geosains, “porositas” adalah rasio antara volume ruang antar butir dalam batuan terhadap keseluruhan volume batuan, umumnya dinyatakan dalam %. Ruang ini menyimpan fluida, termasuk minyak atau gas bumi. Sehingga, porositas batuan sangat penting dalam eksplorasi migas.
OpenCV ditulis dalam C++, tetapi API nya tersedia untuk Python. Kebetulan sekali, saya sedang memulai belajar Python, dan membaca buku Think Python oleh Allen Downey (terima kasih Aji Samudra atas pinjaman bukunya!). Di tulisan ini saya mencoba mempraktekkan penggunaan struktur data dasar Python yaitu list, tuple, dan dictionary.
Tulisan ini sangat terinspirasi oleh tulisan suhu Yogi Pamadya tentang hal yang sama. Saya berusaha membuat fungsi-fungsi yang lebih generalized dan dapat dipakai untuk mencari fraksi komponen-komponen lain dalam sayatan tipis berdasarkan warnanya.
Saya juga melakukan tes terhadap foto sayatan tipis dari laman online learning milik Pak Jerry Lucia (BEG UT Austin) yang terkalibrasi oleh core plug porosity, dan hasilnya tidak jauh berbeda. Sepertinya, hasil akhirnya bisa digunakan untuk membantu kawan-kawan mahasiswa dalam mengerjakan TA-nya. Berikut library yang saya gunakan dalam tulisan ini.
Tugas Akhir dan Sayatan Tipis
Tahun 2013. Harga minyak dunia mencapai lebih dari $100 per barel. “We Are Getting Back Together”-nya Taylor Swift merajai tangga lagu. Saya masih kuliah di tingkat akhir. Pastinya lebih baik daripada hari ini, dimana COVID-19 merajalela!
Di bawah bimbingan Prof. Benyamin Sapiie, Ph.D, saat itu saya melakukan pemetaan tugas akhir ke lapangan. Pemetaan tersebut adalah sebuah proyek dari salah satu perusahan migas multinasional, sebagai studi geologi awal untuk eksplorasi. Bersama dengan rekan saya Azmi dan Fajar, serta tim Geodynamics Research Group ITB, saya diterjunkan ke Pulau Biak-Supiori, Papua, untuk “berburu” komponen-komponen petroleum system di pulau tersebut.
Saya kebagian melakukan pemetaan di Pulau Supiori, tepatnya di Desa Sabarmiokre dan sekitarnya. Berdasarkan studi literatur, daerah pemetaan saya terdiri dari formasi-formasi karbonat. Karena itu, saya memfokuskan studi khusus diagenesis karbonat untuk meneliti potensi salah satu formasi karbonat yang ada di daerah pemetaan, yaitu Formasi Wainukendi, sebagai reservoir.
Singkat cerita, saya pulang dengan selamat, membawa banyak cerita, bab studi khusus diagenesis karbonat bisa dibaca di sini, serta skripsi beserta peta bisa ditemukan di sini.
Gambar di bawah menunjukkan foto sayatan tipis dari TA saya, mewakili seluruh fasies yang saya temukan di lapangan. Secara keseluruhan, terdapat 14 foto sayatan tipis yang saya ambil dari sampel batuan yang saya bawa pulang. Warna biru adalah porositas.
Sejak 2013, saya sudah berkali-kali berganti laptop, dan sulit untuk menemukan sumber gambar asli dari semua foto sayatan tipis saya. Foto yang saya gunakan dalam tulisan ini adalah hasil screenshoot dari tulisan di skripsi saya. Untuk keperluan tulisan ini, rupanya kualitas screenshoot sudah cukup baik.
HSV Filtering
Sayatan tipis yang saya miliki diberikan blue dye untuk menandai porositas, serta alizarin red untuk membedakan mineral dolomit dan kalsit. Sehingga, untuk melakukan penghitungan porositas, saya harus mem-filter warna biru blue dye.
OpenCV memiliki cv2.inRange() untuk melakukan hal tersebut. Dalam colorspace HSV merepresentasikan warna untuk keperluan filtering akan menjadi lebih mudah dibanding dalam BGR, karena warna hanya direpresentasikan oleh satu komponen yaitu Hue (H). Sedangkan parameter Saturation (S) mengatur saturasi warna, dan Value (V) mengatur kecerahan tiap warna.
Untuk memfilter warna, kita hanya perlu mengatur satu parameter (H). Sedangkan dalam colorspace BGR (Blue-Green-Red, default colorspace dari OpenCV) kita perlu mengatur semua parameter untuk mendapatkan suatu warna.
Untuk melakukan hal tersebut, saya mem-wrap cv2.inrange() ke dalam fungsi hsv_filter(). Fungsi ini me-return tuple berisi objek gambar asli dan mask (img dan mask) hasil dari filtering.
hsv_filter() juga memerlukan array berisi batas bawah dan batas atas dari tiap parameter HSV (lower dan upper). Apabila pixel dalam gambar asli masuk ke dalam ambang batas yang kita tentukan, maka akan di-return sebagai pixel putih (True) dalam objek mask.
Berarti, kita harus memainkan batas bawah dan atas ini hingga mendapatkan mask yang menseleksi warna yang kita inginkan. Dengan cara ini, kita bisa melakukan masking terhadap komponen apapun, asalkan dapat dibedakan berdasarkan warna. Misalkan warna biru blue dye untuk porositas, atau warna kemerahan untuk butiran kalsit pada sayatan yang telah di-staining dengan alizarin red.
Lalu bagaimana cara kita mengecek apakah masking yang kita lakukan sudah representatif? Tentunya harus kita amati secara kualitatif dengan menggunakan plot.
Membuat Plot Tes
Untuk menilai apakah masking telah menyeleksi warna yang diinginkan dengan benar, perlu dilakukan plot. Saya menulis fungsi threshold_test_plot() untuk melakukan hal tersebut.
threshold_test_plot() membutuhkan objek gambar asli dan mask sebagai input, kemudian melakukan operasi logika and antara kedua gambar tersebut dan mengeplotnya dengan matplotlib. Untuk menghasilkan True, operator and memerlukan kondisi True di kedua ruas. Dalam OpenCV, operasi ini diimplementasikan dalam operator cv2.bitwise_and().
Objek mask adalah “binary image”: sekumpulan nilai True yang direpresentasikan oleh warna putih, dan False yang direpresentasikan oleh warna hitam. Apabila cv2.bitwise_and() me-return True, maka gambar asli akan terlihat, begitu pula sebaliknya. Dengan kata lain, operasi ini akan menghasilkan gambar asli yang telah di-masking.
Pengeplotan dilakukan menggunakan matplotlib yang hanya memahami colorspace RGB, sehingga saya menuliskan fungsi untuk mengkonversi gambar menjadi RGB (para pembaca, tolong koreksi apabila saya salah dan ada cara tanpa konversi ya!).
Tiba saatnya untuk menggunakan kedua fungsi yang telah saya tulis pada sayatan tipis TA saya. Saya juga berlatih menggunakan dictionary dan tuple dalam proses ini. Hasil proses saya simpan dalam suatu dictionary (img_dict), sebuah struktur data berupa key-value pair yang menggunakan:
- string nama file gambar asli sebagai key
- tuple berisi objek gambar asli dan mask sebagai value
Kemudian, saya bisa melakukan looping terhadap img_dict menggunakan fungsi hsv_filter() dan threshold_test_plot().
Saya menggunakan nilai HSV [70, 0, 80] sebagai batas bawah dan [179, 255, 255] sebagai batas atas filter untuk memfilter warna blue dye yang mengisi porositas. Menurut saya, hasilnya cukup memuaskan di sebagian besar gambar. Berikut dua contoh plot dari 14 sayatan tipis yang saya punya.
Tetapi ada satu sayatan yang hasilnya kurang memuaskan untuk saya. Tentunya hal ini wajar: tidak mungkin ada satu nilai threshold yang bisa menyelesaikan seluruh permasalahan anda. Sehingga, untuk satu sayatan di bawah ini, saya menggunakan threshold khusus.
Dalam Python, dictionary adalah objek yang mutable. Kita bisa mengubah isinya tanpa perlu mengubah identitasnya. Sehingga, untuk mengubah salah satu elemen img_dict, kita dapat menggunakan cara berikut:
There you go. Terlihat lebih baik!
Setelah yakin dengan semua masking yang saya lakukan, saatnya melakukan penghitungan porositas…
Menghitung Porositas
Cara penghitungan porositas batuan cukup mudah. Mengacu kepada definisi porositas batuan, kita tinggal menghitung % pixel berwarna putih dalam objek mask. Saya menulis compute_fraction() untuk melakukan hal tersebut.
Pro tip: operasi mean terhadap sekumpulan nilai boolean menghasilkan %
Truedari keseluruhan kumpulan tersebut. Berlaku di mana saja!
Lalu bagaimana caranya yakin bahwa tahapan yang kita lakukan sudah benar? Saya mencari sayatan tipis yang juga memiliki porositas dari sampel core batuan, dan menemukan materi online dari Pak Jerry Lucia yang membahas mengenai kaitan antara rock fabric batuan karbonat dengan hubungan porositas-permeabilitas batuan tersebut.
Saya memilih salah satu foto batugamping grainstone untuk mengetes fungsi yang saya buat. Dalam keterangannya, sampel core plug dari grainstone ini memiliki porositas sebesar 20,1%. Sedangkan, metode yang saya gunakan menghasilkan porositas sebesar 20,94%. Pretty good! Dengan pengaturan threshold lebih lanjut, mungkin kita bisa mendapatkan akurasi yang lebih baik.
Setelah yakin dengan fungsi tersebut, saya menggunakan compute_fraction() ke seluruh objek mask dari sayatan tipis TA saya. Sekali lagi, saya menggunakan dictionary (yang dinamai por_dict) untuk menyimpan string nama file dan hasil hitungan porositas sebagai key-value pair.
Hasil
Tabel di bawah menunjukkan perbedaan antara penghitungan menggunakan kode di atas, dengan porositas hasil pengamatan kualitatif yang saya tulis dalam skripsi. And, ha… what a joke! Tidak ada konsistensi dalam pengamatan kualitatif yang saya lakukan. Saya juga cenderung melakukan over-estimate, serta — yang paling fatal dan baru saya sadari sekarang — mengira mineral berwarna kebiruan di XPL sebagai porositas!
//==============[]============[]================================\\
|| Filename || OpenCV (%) || Qualitative (%) by yours truly ||
|]==============[]============[]================================[|
|| Wmd1.03.png || 0 || 10 ||
|| Wmd1.05.png || 0 || 5 ||
|| Rsd1.02.png || 4.97 || 10 ||
|| Bnr1.07.png || 0.98 || 15 ||
|| Bnr1.04.png || 0.02 || 5 ||
|| Wpr1.11.png || 0 || 5 ||
|| Rsd1.05.png || 0 || 5 ||
|| Bnr1.01.png || 0.01 || 2 ||
|| Rwrt1.04.png || 0.36 || 5 ||
|| Imk1.01.png || 0.06 || 5 ||
|| Wbr1.04.png || 1.35 || 7.5 ||
|| Rsd1.05A.png || 0.5 || 5 ||
|| Wbr1.07.png || 2.94 || 9 ||
|| Rsd1.05B.png || 0.44 || 3 ||
\\==============[]============[]================================//Mengingat cara saya dalam mengerjakan tugas akhir, ditambah hasil kalibrasi dengan foto yang ada core plug porosity-nya, tentunya saya lebih yakin dengan hasil penghitungan hari ini sih…
Kesimpulan dan Saran
Tulisan ini adalah salah satu hal yang paling seru yang pernah saya lakukan dengan Python (karena masih jarang juga sih). Saya belajar penggunaan list, tuple, dan dictionary, serta mencoba ide-ide dasar computer vision. Untuk itu saya berterima kasih kepada Yogi Pamadya yang mengenalkan ide ini!
Selanjutnya, berikut kesimpulan dan saran yang dapat saya sampaikan:
- Metode penghitungan porositas menggunakan HSV filtering dari foto sayatan tipis dapat menghasilkan angka porositas yang cukup akurat, namun pengetesan dengan foto yang telah dikalibrasi dengan core plug porosity tetap diperlukan
- Pengamatan kualitatif pada sayatan tipis cenderung menghasilkan over-estimate porositas
- Untuk pengembangan selanjutnya, dapat dilakukan identifikasi tipe porositas menggunakan OpenCV, yang kemudian dapat digunakan dalam estimasi permeabilitas seperti yang dicontohkan Pak Jerry Lucia
- Pengembangan lebih jauh: penggunaan OpenCV dan convolutional neural network untuk mengenali spesies foraminifera kecil pada sayatan tipis
- Metode yang mirip dapat digunakan untuk menyelesaikan permasalahan lain dalam geologi, seperti picking kelurusan dari citra satelit/drone, digitasi peta geologi lama, dan lain-lain. At scale!
Terima kasih telah membaca tulisan kali ini, semoga bermanfaat. Seperti biasa, notebook Python yang berisi kode yang saya gunakan dapat ditemukan di repo GitHub berikut. Diskusi, kritik, dan saran saya tunggu di kolom komentar. Apabila anda mencoba proses ini pada dataset yang anda miliki, kabari saya soal hasilnya, ya!
