Xem dữ liệu thể tích 3D với Matplotlib

Phần lớn các bạn đã quen thuộc với dữ liệu ảnh chụp bằng máy ảnh thông thường (trong tài liệu khoa học thường gọi là “ảnh tự nhiên”), cũng như từ các thiết bị chuyên dụng như kính hiển vi hay kính thiên văn. Khi làm việc với ảnh trong Python, cách phổ biến nhất để hiển thị là dùng hàm imshow của Matplotlib, thư viện vẽ biểu đồ phổ biến nhất của Python.

Trong hướng dẫn này, chúng tôi sẽ chỉ cho bạn cách mở rộng hàm này để hiển thị dữ liệu thể tích 3D, bạn có thể hình dung như một chồng ảnh. Cùng nhau, chúng mô tả một cấu trúc 3 chiều. Ví dụ, chụp cộng hưởng từ (MRI) và chụp cắt lớp vi tính (CT) đo đạc cấu trúc 3D bên trong cơ thể người; X-ray microtomography đo cấu trúc 3D bên trong các vật liệu như thủy tinh hay hợp kim; và kính hiển vi light-sheet đo các hạt huỳnh quang trong mô sinh học.

Chúng ta sẽ minh họa cách tải một bộ dữ liệu MRI và hiển thị các lát cắt bằng matplotlib. Bạn sẽ học cách

• Lấy dữ liệu với sự trợ giúp của thư viện nibabel, vốn cung cấp bộ đọc cho các tệp định dạng NIfTI. Với những bạn muốn bỏ qua bước này và vẽ luôn, hướng dẫn cũng chỉ cách lấy dữ liệu bằng scikit-image.
• Tiếp theo, bạn sẽ thực hiện những bước đầu tiên để xây dựng một trình xem lát cắt hoàn chỉnh bằng API bộ xử lý sự kiện của matplotlib, và
• Bạn sẽ viết lại một phần mã để tận dụng công cụ mới.
• Cuối cùng, bạn sẽ xem cách dọn dẹp không gian làm việc sau khi hoàn tất.

Nhưng trước hết, hãy điểm qua một số kiến thức cơ bản: cách hiển thị ảnh bằng imshow của Matplotlib.

Bạn có thể bắt đầu bằng cách bật chế độ tương tác của matplotlib nếu đang làm việc với Jupyter Notebook:

%matplotlib notebookĐược hỗ trợ bởi

Bây giờ bạn có thể nhập matplotlib và hiển thị một số dữ liệu. Bạn sẽ tải một số dữ liệu ví dụ có sẵn trong mô-đun data của thư viện 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');Được hỗ trợ bởi

Lưu ý: Khi bạn chạy matplotlib ở chế độ notebook tương tác, hình đang mở sẽ là hình hoạt động duy nhất cho đến khi bạn tắt nó bằng biểu tượng nguồn ở góc trên bên phải của hình. Hãy chắc chắn làm điều đó trước khi chuyển sang mỗi biểu đồ.

Những ảnh này được gọi là ảnh 2 chiều hay 2D vì chúng sắp xếp theo 2 chiều: x và y, hoặc theo cách nói của NumPy là hàng và cột, hay r và c.

Một số ảnh là 3D, vì chúng có thêm chiều độ sâu (z, hay các mặt phẳng). Bao gồm MRI và hiển vi điện tử truyền qua theo lát nối tiếp (ssTEM), trong đó mẫu được cắt mỏng như xúc xích salami, và mỗi lát được chụp riêng.

Để xem các ảnh như vậy trong matplotlib, chúng ta phải chọn một lát và chỉ hiển thị lát đó. Hãy thử với một số dữ liệu MRI miễn phí trên mạng.

Gián đoạn: Lấy dữ liệu

Chúng ta sẽ tải một bộ dữ liệu được mô tả trong Buchel và Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Đầu tiên, chúng ta tạo một thư mục tạm để tải dữ liệu. Hãy nhớ xóa nó khi hoàn tất phân tích! Nếu bạn muốn giữ bộ dữ liệu này để dùng sau, hãy đổi d thành đường dẫn thư mục cố định theo ý bạn.

import tempfile

# Create a temporary directory
d = tempfile.mkdtemp()Được hỗ trợ bởi

Bây giờ, hãy tải dữ liệu:

import os

# Return the tail of the path
os.path.basename('http://google.com/attention.zip')Được hỗ trợ bởi

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'))Được hỗ trợ bởi

Và giải nén từ tệp zip vào thư mục tạm:

import zipfile

# Extract the contents into the temporary directory we created earlier
zipfile.ZipFile(fn).extractall(path=d)Được hỗ trợ bởi

Nếu bạn xem nội dung thực tế của tệp, bạn sẽ thấy một loạt tệp ‘.hdr’ và ‘.img’.

# List first 10 files
[f.filename for f in zipfile.ZipFile(fn).filelist[:10]]Được hỗ trợ bởi

Chúng ở định dạng tệp NIfTI, và chúng ta cần một bộ đọc cho định dạng này. May mắn là thư viện nibabel tuyệt vời cung cấp bộ đọc như vậy. Hãy cài nó bằng conda install -c conda-forge nibabel hoặc pip install nibabel, rồi:

import nibabelĐược hỗ trợ bởi

Giờ, chúng ta có thể đọc ảnh và dùng phương thức .get_data() để lấy một mảng NumPy để xem:

# 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()Được hỗ trợ bởi

Mẹo: nếu bạn muốn chuyển thẳng sang vẽ dữ liệu MRI, hãy chạy các dòng mã sau:

from skimage import io

struct_arr = io.imread("https://s3.amazonaws.com/assets.datacamp.com/blog_assets/attention-mri.tif")Được hỗ trợ bởi

Quay lại việc vẽ

Giờ hãy xem một lát trong mảng đó:

plt.imshow(struct_arr[75])Được hỗ trợ bởi

Ối! Trông khá bẹt! Đó là vì độ phân giải theo trục dọc trong nhiều ảnh MRI không giống với theo các trục ngang. Ta có thể khắc phục bằng cách truyền tham số aspect vào hàm imshow:

plt.imshow(struct_arr[75], aspect=0.5)Được hỗ trợ bởi

Nhưng để đơn giản, chúng ta sẽ chuyển vị dữ liệu và chỉ xem các lát ngang, vốn không cần chỉnh tỉ lệ như vậy.

struct_arr2 = struct_arr.T
plt.imshow(struct_arr2[34])Được hỗ trợ bởi

Đẹp! Dĩ nhiên, để xem lát khác, hoặc lát theo trục khác, ta cần gọi imshow lần nữa:

plt.imshow(struct_arr2[5])Được hỗ trợ bởi

Tất cả những lần gọi này khá nhanh trở nên nhàm chán. Lâu nay, tôi thường xem thể tích 3D bằng các công cụ ngoài Python, như ITK-SNAP. Nhưng hóa ra, rất dễ để thêm khả năng “cuộn” 3D vào trình xem matplotlib! Điều này cho phép chúng ta khám phá dữ liệu 3D ngay trong Python, giảm thiểu nhu cầu chuyển đổi ngữ cảnh giữa khám phá và phân tích dữ liệu.

Mấu chốt là dùng API bộ xử lý sự kiện của matplotlib, cho phép ta định nghĩa các hành động thực hiện trên biểu đồ — bao gồm thay đổi dữ liệu của biểu đồ! — để phản hồi các lần nhấn phím hoặc nhấp chuột nhất định.

Trong trường hợp của chúng ta, hãy gán phím J và K trên bàn phím cho “lát trước” và “lát sau”:

def previous_slice():
    pass

def next_slice():
    pass

def process_key(event):
    if event.key == 'j':
        previous_slice()
    elif event.key == 'k':
        next_slice()Được hỗ trợ bởi

Khá đơn giản! Tất nhiên, chúng ta cần tìm cách hiện thực các hành động này và cần cho hình biết rằng nó nên dùng hàm process_key để xử lý các lần nhấn phím! Việc sau thì đơn giản: chỉ cần dùng phương thức canvas của hình mpl_connect:

fig, ax = plt.subplots()
ax.imshow(struct_arr[..., 43])
fig.canvas.mpl_connect('key_press_event', process_key)Được hỗ trợ bởi

Bạn có thể tìm tài liệu đầy đủ cho mpl_connect tại đây, bao gồm các loại sự kiện khác có thể gán (như nhấp chuột).

Tôi chỉ cần khám phá một chút để phát hiện rằng imshow trả về một đối tượng AxesImage, đối tượng này “sống” bên trong đối tượng Axes của matplotlib nơi mọi việc vẽ diễn ra, trong thuộc tính .images. Và đối tượng này cung cấp phương thức tiện lợi set_array để thay thế dữ liệu ảnh đang hiển thị! Vậy, tất cả những gì ta cần làm là:

• vẽ một chỉ số bất kỳ và lưu chỉ số đó, có thể như một thuộc tính chạy thời gian bổ sung trên đối tượng Axes.
• cung cấp các hàm next_slice và previous_slice để thay đổi chỉ số và dùng set_array để đặt lát tương ứng của thể tích 3D.
• dùng phương thức draw của canvas hình để vẽ lại với dữ liệu mới.

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])Được hỗ trợ bởi

Hãy thử nhé!

multi_slice_viewer(struct_arr2)Được hỗ trợ bởi

Hoạt động rồi! Tuyệt! Nhưng nếu bạn thử tại nhà, bạn sẽ thấy cuộn lên bằng K cũng làm ép tỷ lệ ngang của biểu đồ. Hả? (Chỉ xảy ra khi chuột đang ở trên ảnh.)

Điều đang xảy ra là việc thêm các bộ xử lý sự kiện vào Matplotlib chỉ là chồng chất chúng lên nhau. Trong trường hợp này, K là phím tắt tích hợp để đổi trục x sang thang logarit. Nếu muốn dùng riêng K, ta phải loại nó khỏi sơ đồ phím mặc định của matplotlib. Chúng tồn tại dưới dạng danh sách trong từ điển plt.rcParams, vốn là kho thiết lập mặc định cấp hệ thống của matplotlib:

plt.rcParams['keymap.<command>'] = ['<key1>', '<key2>']Được hỗ trợ bởi

trong đó nhấn bất kỳ phím nào trong danh sách (tức <key1> hoặc <key2>) sẽ khiến <command> được thực thi.

Vì vậy, ta cần viết một hàm trợ giúp để xóa các phím ta muốn dùng ở bất cứ nơi nào chúng xuất hiện trong từ điển này. (Hàm này hiện chưa có trong matplotlib, nhưng có lẽ sẽ là một đóng góp đáng hoan nghênh!)

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)Được hỗ trợ bởi

Trình xem lát cắt hoàn chỉnh

Được rồi, hãy viết lại hàm để tận dụng công cụ mới này:

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])Được hỗ trợ bởi

Bây giờ, chúng ta có thể xem tất cả các lát trong thể tích MRI mà không bị làm phiền bởi sơ đồ phím mặc định nữa!

multi_slice_viewer(struct_arr2)Được hỗ trợ bởi

Một điểm hay của phương pháp này là nó hoạt động trên bất kỳ backend nào của matplotlib! Vì vậy, nếu bạn thử trong bảng điều khiển IPython terminal, bạn vẫn sẽ có tương tác giống như trên trình duyệt! Và điều tương tự đúng với ứng dụng Qt hay Tkinter nhúng biểu đồ matplotlib. Công cụ đơn giản này do đó cho phép bạn xây dựng các ứng dụng ngày càng phức tạp dựa trên khả năng trực quan hóa của matplotlib.

Trước khi rời đi

Đừng quên tự dọn dẹp và xóa thư mục tạm (nếu bạn đã tạo):

import shutil

# Remove the temporary directory
shutil.rmtree(d)Được hỗ trợ bởi