Melihat Data Volumetrik 3D Dengan Matplotlib

Kebanyakan dari Anda sudah akrab dengan data citra yang diambil dengan kamera biasa (sering disebut “citra alami” dalam literatur ilmiah), tetapi juga dengan instrumen khusus seperti mikroskop atau teleskop. Saat bekerja dengan citra di Python, cara paling umum untuk menampilkannya adalah menggunakan fungsi imshow dari Matplotlib, pustaka plotting paling populer di Python.

Dalam tutorial ini, kami akan menunjukkan cara memperluas fungsi ini untuk menampilkan data volumetrik 3D, yang dapat Anda bayangkan sebagai tumpukan gambar. Bersama-sama, tumpukan tersebut menggambarkan sebuah struktur 3D. Misalnya, magnetic resonance imaging (MRI) dan computed tomography (CT) memetakan struktur 3D di dalam tubuh manusia; X-ray microtomography memetakan struktur 3D di dalam material seperti kaca atau paduan logam; dan light-sheet microscope memetakan partikel fluoresen di dalam jaringan biologis.

Kami akan mendemonstrasikan cara mengunduh dataset MRI dan menampilkan irisan menggunakan matplotlib. Anda akan belajar cara

• Mengambil data dengan bantuan pustaka nibabel, yang menyediakan pembaca untuk berkas berformat NIfTI. Bagi Anda yang ingin melewati langkah ini dan langsung mulai membuat plot, tutorial ini juga akan menunjukkan cara mendapatkan data dengan scikit-image.
• Selanjutnya, Anda akan mengambil langkah pertama menuju penampil irisan yang fungsional penuh dengan API event handler matplotlib, dan
• Anda akan menulis ulang sebagian kode untuk memanfaatkan alat baru Anda.
• Terakhir, Anda akan melihat cara merapikan workspace setelah selesai.

Tetapi pertama-tama, mari tinjau dasar-dasarnya: cara menampilkan gambar dengan imshow milik Matplotlib.

Anda bisa mulai dengan mengaktifkan mode interaktif matplotlib jika bekerja dengan Jupyter Notebook:

%matplotlib notebookDidukung oleh

Sekarang Anda dapat mengimpor matplotlib dan menampilkan sebagian data. Anda akan memuat beberapa data contoh yang disertakan dalam modul data pustaka scikit-image:

import matplotlib.pyplot as plt
from skimage import data

astronaut = data.astronaut()
ihc = data.immunohistochemistry()
hubble = data.hubble_deep_field()

# Initialize the subplot panels side by side
fig, ax = plt.subplots(nrows=1, ncols=3)

# Show an image in each subplot
ax[0].imshow(astronaut)
ax[0].set_title('Natural image')
ax[1].imshow(ihc)
ax[1].set_title('Microscopy image')
ax[2].imshow(hubble)
ax[2].set_title('Telescope image');Didukung oleh

Catatan: Saat Anda menjalankan matplotlib dalam mode notebook interaktif, figur yang terbuka akan tetap menjadi satu-satunya figur aktif hingga Anda menonaktifkannya menggunakan simbol daya di kanan atas figur. Pastikan Anda melakukannya sebelum beralih dari setiap plot.

Gambar-gambar ini disebut gambar 2 dimensi atau 2D karena disusun sepanjang 2 dimensi: x dan y, atau, dalam istilah NumPy, baris dan kolom atau r dan c.

Beberapa gambar bersifat 3D, karena memiliki dimensi kedalaman tambahan (z, atau bidang/plane). Ini termasuk magnetic resonance imaging (MRI) dan serial section transmission electron microscopy (ssTEM), di mana sampel diiris tipis seperti salami dan setiap irisan dipindai secara terpisah.

Untuk melihat gambar semacam itu di matplotlib, kita harus memilih satu irisan, dan hanya menampilkan irisan tersebut. Mari kita coba pada beberapa data MRI yang tersedia bebas secara online.

Selingan: Mengambil Data

Kita akan mengunduh dataset yang dijelaskan dalam Buchel dan Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Pertama, kita membuat direktori sementara untuk mengunduh data. Ingat untuk menghapusnya setelah selesai dengan analisis! Jika Anda ingin menyimpan dataset ini untuk digunakan nanti, ubah d ke lokasi direktori permanen pilihan Anda.

import tempfile

# Create a temporary directory
d = tempfile.mkdtemp()Didukung oleh

Sekarang, mari unduh datanya:

import os

# Return the tail of the path
os.path.basename('http://google.com/attention.zip')Didukung oleh

from urllib.request import urlretrieve

# Define URL
url = 'http://www.fil.ion.ucl.ac.uk/spm/download/data/attention/attention.zip'

# Retrieve the data
fn, info = urlretrieve(url, os.path.join(d, 'attention.zip'))Didukung oleh

Dan ekstrak dari berkas zip ke direktori sementara kita:

import zipfile

# Extract the contents into the temporary directory we created earlier
zipfile.ZipFile(fn).extractall(path=d)Didukung oleh

Jika Anda melihat isi sebenarnya dari berkas tersebut, Anda akan menemukan sejumlah berkas ‘.hdr’ dan ‘.img’.

# List first 10 files
[f.filename for f in zipfile.ZipFile(fn).filelist[:10]]Didukung oleh

Ini menggunakan format berkas NIfTI, dan kita memerlukan pembaca untuknya. Untungnya, pustaka nibabel yang sangat baik menyediakan pembaca tersebut. Pastikan Anda memasangnya dengan conda install -c conda-forge nibabel atau pip install nibabel, lalu:

import nibabelDidukung oleh

Sekarang, kita akhirnya bisa membaca gambar kita, dan menggunakan metode .get_data() untuk mendapatkan array NumPy yang akan ditampilkan:

# Read the image 
struct = nibabel.load(os.path.join(d, 'attention/structural/nsM00587_0002.hdr'))

# Get a plain NumPy array, without all the metadata
struct_arr = struct.get_data()Didukung oleh

Tip: jika Anda ingin langsung lanjut ke plotting data MRI, jalankan baris kode berikut:

from skimage import io

struct_arr = io.imread("https://s3.amazonaws.com/assets.datacamp.com/blog_assets/attention-mri.tif")Didukung oleh

Kembali ke Plotting

Sekarang mari lihat sebuah irisan dalam array tersebut:

plt.imshow(struct_arr[75])Didukung oleh

Wah! Terlihat cukup gepeng! Itu karena resolusi sepanjang sumbu vertikal pada banyak MRI tidak sama dengan sepanjang sumbu horizontal. Kita bisa memperbaikinya dengan meneruskan parameter aspect ke fungsi imshow:

plt.imshow(struct_arr[75], aspect=0.5)Didukung oleh

Tetapi, agar lebih mudah, kita akan cukup melakukan transpose pada data dan hanya melihat irisan horizontal, yang tidak memerlukan pengaturan seperti itu.

struct_arr2 = struct_arr.T
plt.imshow(struct_arr2[34])Didukung oleh

Bagus! Tentu saja, untuk melihat irisan lain, atau irisan sepanjang sumbu berbeda, kita memerlukan pemanggilan imshow lainnya:

plt.imshow(struct_arr2[5])Didukung oleh

Semua pemanggilan ini cepat terasa membosankan. Untuk waktu yang lama, saya melihat volume 3D menggunakan alat di luar Python, seperti ITK-SNAP. Tetapi, ternyata, cukup mudah menambahkan kemampuan “menggulir” 3D ke penampil matplotlib! Ini memungkinkan kita menjelajahi data 3D langsung di Python, meminimalkan kebutuhan berganti konteks antara eksplorasi data dan analisis data.

Kuncinya adalah menggunakan API event handler matplotlib, yang memungkinkan kita mendefinisikan aksi yang dilakukan pada plot — termasuk mengubah data plot! — sebagai respons terhadap penekanan tombol tertentu atau klik tombol mouse.

Dalam kasus kita, mari ikat tombol J dan K pada papan ketik ke “irisan sebelumnya” dan “irisan berikutnya”:

def previous_slice():
    pass

def next_slice():
    pass

def process_key(event):
    if event.key == 'j':
        previous_slice()
    elif event.key == 'k':
        next_slice()Didukung oleh

Cukup sederhana! Tentu, kita perlu mencari cara untuk benar-benar mengimplementasikan aksi ini dan kita perlu memberi tahu figur bahwa ia harus menggunakan fungsi process_key untuk memproses penekanan tombol! Yang terakhir ini sederhana: kita hanya perlu menggunakan metode kanvas figur mpl_connect:

fig, ax = plt.subplots()
ax.imshow(struct_arr[..., 43])
fig.canvas.mpl_connect('key_press_event', process_key)Didukung oleh

Anda dapat menemukan dokumentasi lengkap untuk mpl_connect di sini, termasuk jenis event lain yang bisa Anda ikat (seperti klik tombol mouse).

Perlu sedikit eksplorasi untuk mengetahui bahwa imshow mengembalikan objek AxesImage, yang hidup “di dalam” objek Axes matplotlib tempat semua penggambaran berlangsung, di atribut .images-nya. Dan objek ini menyediakan metode set_array yang praktis untuk menukar data gambar yang ditampilkan! Jadi, yang perlu kita lakukan hanyalah:

• mem-plot sebuah indeks sembarang, dan menyimpan indeks tersebut, mungkin sebagai atribut runtime tambahan pada objek Axes.
• menyediakan fungsi next_slice dan previous_slice yang mengubah indeks dan menggunakan set_array untuk menetapkan irisan yang sesuai dari volume 3D.
• menggunakan metode kanvas figur draw untuk menggambar ulang figur dengan data baru.

def multi_slice_viewer(volume):
    fig, ax = plt.subplots()
    ax.volume = volume
    ax.index = volume.shape[0] // 2
    ax.imshow(volume[ax.index])
    fig.canvas.mpl_connect('key_press_event', process_key)

def process_key(event):
    fig = event.canvas.figure
    ax = fig.axes[0]
    if event.key == 'j':
        previous_slice(ax)
    elif event.key == 'k':
        next_slice(ax)
    fig.canvas.draw()

def previous_slice(ax):
    """Go to the previous slice."""
    volume = ax.volume
    ax.index = (ax.index - 1) % volume.shape[0]  # wrap around using %
    ax.images[0].set_array(volume[ax.index])

def next_slice(ax):
    """Go to the next slice."""
    volume = ax.volume
    ax.index = (ax.index + 1) % volume.shape[0]
    ax.images[0].set_array(volume[ax.index])Didukung oleh

Mari kita coba!

multi_slice_viewer(struct_arr2)Didukung oleh

Berhasil! Bagus! Tetapi, jika Anda mencobanya sendiri, Anda akan melihat bahwa saat menggulir ke atas dengan K, skala horizontal plot juga jadi gepeng. Hah? (Ini hanya terjadi jika kursor mouse Anda berada di atas gambar.)

Yang terjadi adalah penambahan event handler ke Matplotlib hanya menumpuknya satu demi satu. Dalam kasus ini, K adalah pintasan bawaan papan ketik untuk mengubah sumbu x menjadi skala logaritmik. Jika kita ingin menggunakan K secara eksklusif, kita harus menghapusnya dari peta tombol bawaan matplotlib. Peta tersebut tersimpan sebagai daftar dalam kamus plt.rcParams, yang merupakan repositori matplotlib untuk pengaturan default sistem:

plt.rcParams['keymap.<command>'] = ['<key1>', '<key2>']Didukung oleh

di mana menekan salah satu tombol dalam daftar (yaitu <key1> atau <key2>) akan menyebabkan <command> dijalankan.

Jadi, kita perlu menulis fungsi pembantu untuk menghapus tombol yang ingin kita gunakan di mana pun tombol tersebut muncul dalam kamus ini. (Fungsi ini belum ada di matplotlib, tetapi mungkin akan menjadi kontribusi yang disambut baik!)

def remove_keymap_conflicts(new_keys_set):
    for prop in plt.rcParams:
        if prop.startswith('keymap.'):
            keys = plt.rcParams[prop]
            remove_list = set(keys) & new_keys_set
            for key in remove_list:
                keys.remove(key)Didukung oleh

Penampil Irisan yang Sepenuhnya Fungsional

Baik, mari kita tulis ulang fungsi kita untuk memanfaatkan alat baru ini:

def multi_slice_viewer(volume):
    remove_keymap_conflicts({'j', 'k'})
    fig, ax = plt.subplots()
    ax.volume = volume
    ax.index = volume.shape[0] // 2
    ax.imshow(volume[ax.index])
    fig.canvas.mpl_connect('key_press_event', process_key)

def process_key(event):
    fig = event.canvas.figure
    ax = fig.axes[0]
    if event.key == 'j':
        previous_slice(ax)
    elif event.key == 'k':
        next_slice(ax)
    fig.canvas.draw()

def previous_slice(ax):
    volume = ax.volume
    ax.index = (ax.index - 1) % volume.shape[0]  # wrap around using %
    ax.images[0].set_array(volume[ax.index])

def next_slice(ax):
    volume = ax.volume
    ax.index = (ax.index + 1) % volume.shape[0]
    ax.images[0].set_array(volume[ax.index])Didukung oleh

Sekarang, kita seharusnya bisa melihat semua irisan dalam volume MRI tanpa gangguan dari keymap bawaan!

multi_slice_viewer(struct_arr2)Didukung oleh

Satu fitur menarik dari metode ini adalah ia bekerja pada semua backend matplotlib! Jadi, jika Anda mencobanya di konsol terminal IPython, Anda tetap akan mendapatkan interaksi yang sama seperti di browser! Hal yang sama juga berlaku untuk aplikasi Qt atau Tkinter yang menyematkan plot matplotlib. Alat sederhana ini karenanya memungkinkan Anda membangun aplikasi yang semakin kompleks di sekitar kemampuan visualisasi matplotlib.

Sebelum Anda Pergi

Jangan lupa untuk merapikan setelahnya, dan hapus direktori sementara (jika Anda membuatnya):

import shutil

# Remove the temporary directory
shutil.rmtree(d)Didukung oleh