2021-12-29 雙月會筆記

放射醫學會每兩個月會有個雙月會，由輪值醫院的放射科報告，本質上是一個分散型的年會。這次輪到台大醫院，台大影像科決定報告他們科關於 AI 方面的進展。底下是筆記和感想。

筆記部分

創業（TFDA 的 SaMD 審核介紹）
- 講師非來自台大，而是任職於雲象科技的林佩蓁醫師
- 基於機器學習的醫療器材術語上稱為 SaMD，software as a medical device，歸於 TFDA 管轄。官方文件是《醫療器材分類分級管理辦法》。講者據此介紹一二級軟體與三四級軟體的差別及申請 TFDA 審核所需要文件。
- 所謂四級 SaMD，即提供危急情況的治療或診斷且直接 bypass 醫師；如果不 bypass 醫師，只是單純輔助就是一二級
確效機制
- 討論醫院對於橋接商業軟體到院內系統的審核機制。其中，醫院若要採用未經 FDA / TFDA 許可的三四軟體，需要申請臨床試驗。
- 確效機制是確保軟體在本地環境下（病人族群、軟體系統、放射機器）是否能按照廠商所宣稱的運作。另外也有資安及病安的議題。
- 確效機制包含哪些項目？例如：資料傳輸過程是否保留、資料 ground truth 最終由誰判定、系統資訊安全維護。
案例分享
- COVID-19 detection based on chest radiography
- CV calcium score auto-calculation based on HeaortaNet
- Integration of LunitAI（chest radiography、mammography）
- ICH detection

COVID-19 detection based on chest radiography

台大的 COVID-19 偵測模型使用了 3‑step hierarchical CNN。

先取出 lung field
判斷是否異常
再判斷異常是否為 COVID-19

Model

不過其實模型讀取的是壓縮過的 512 x 512 影像。究竟機器對這樣的影像其判別能力為何？如果用其他 atypical pneumonia 的影像做測試，表現能力如何？似乎沒有特別提到。報告者用幾個範例介紹過去這個項目。

感想

隔天我們科的 AI 負責人阿吉主任也有來跟我聊感想。

整合與確效：這部分我們科也有出台自己的一套方案，跟台大的差不多，包括資料管理以及旁路安全性測試等。在整合部分，目前確定接入的是 ClearRead 這套去除 CT 血管的演算法，可以一鍵讀取，很方便。相較之下，台大的 Lunit 接口似乎是下載再上傳（？），沒有那麼方便。
硬體資源：我覺得台大額外的強項是來自大學方面的軟體技術支援，比我們強很多。
住院醫師參與：台大這次雙月會出動了多名住院醫師報告。我很好奇他們的 R 對 AI 的感興趣程度。基本上，不做介入就一定得做 AI，這是放射科未來的宿命。不過，這些題目似乎都是是由其他科大佬所推動（例如以 HeaortaNet 為例，它是 TWCVAI 計畫所衍生的產物）。當然，我們小團隊無法跟醫學中心聯盟其幾百萬張的標註影像和來自大學資工科系的博士生們抗衡，所以必須尋找較冷門的題目發展一整套環環相扣的研究題目，才是生存之道（例如我選擇主動脈疾病、支架放置到血流動力學模擬等）。

在 Windows 上安裝多個 Python

今天幫學長在 Windows 上裝一套叫 ASReview 的軟體，而這個軟體不曉得為什麼特別指定，需要 Python 3.8.1 版本。學長弄了半天還是沒搞定，所以來著找我幫忙看看。我自己 Windows 基本只用於打電玩；主機另外有一套 Arch Linux，而同為 Unix 衍生的 macOS 配置多版本 Python 也不是什麼困難事情，花了一點時間才搞懂如何在 Windows 上面設定。本來想說建議學長乾脆直接用 Conda，不過看學長已經很忙了，還是用比較簡單的方案。整個步驟如下：

電腦裡面原本裝了 Python 3.10.0 在家目錄下，有設定 PATH
可以在 CMD 下面呼叫到 pip，於是安裝 virtualenvwrapper-win。需要注意的是，virtualenvwrapper 其實是一套 shell script，不適用於 Windows 環境，需要裝 -win 版本。
接著安裝所依賴的 Python 版本，同樣是裝在家目錄，不過要取消所有預設項目（例如：關聯副檔名為 .py、創建環境變數）
最後，使用 virtualenvwrapper-win 的指令（如同 virtualenvwrapper），新建環境，並通過 p 參數指定 Python 解釋器的位置。接著就可以切換進去並通過 pip 安裝需要的軟體了。

Marp 使用細節

用 Marp 做投影片有些額外的需求，上網讀了官方文件有了如下筆記。

圖片大小及全圖背景

Marp 利用了 markdown 圖片語法中原本用來放 alternative text 的地方來放圖片描述詞。可以用的選項包括：

寬度 width:100px、高度 height:100px
模糊 blur、灰階 grayscale、亮度 brightness

例如：![width:100px](image.png)。另外，在替代文字的開頭設置 bg 則可以讓圖片變成背景圖，例如：![bg](image.png)。繼續添加 right 或 left 則可以變成半幅背景。完整設定參考官方文件。

HTML 元素細節調整

可以利用 HTML 註解來影響單獨的頁面，例如：

￰YFWU00000000￱

或者在 YAML header 使用 style 來描述，例如：

---
style: |
  table {
    font-size: 25px;
  }
  section {
    font-family: "Roboto Slab";
  }
  code {
    font-family: 'Roboto Mono';
  }
---

額外的 YAML header 設定

投影片 footer / header：在開頭的 YAML 添加 footer 項，例如：footer: Author information。Header 部分也是比照使用。
不再用 --- 作為分頁標誌，而是使用 header level 作為表示，需在 YAML 設定 headingDivider: 2 就會讓 Marp 自動按標題等級切割。
可以設定 math 來改變 TeX 渲染引擎，選項可以是 MathJax 或 KaTeX。MathJax 支援的符號類型較完整，但是 KaTeX 速度較快。

模糊字串比對

小學妹的研究中有個要從病歷中獲取關鍵字的需求。例如：pneumonia，不過考慮到我們的病歷都是疲累的住院醫師打出來的，所以有不少的錯字（例如：多了一個字 pneumoniaa / 字序有誤 peunomonia），需要模糊比對。稍微研究了一下，了解到這個區塊的正式知識叫做 approximate string matching

字符串匹配是很著名的經典計算機科學主題。例如 Knuth-Morris-Pratt 算法、Trie 字典樹、編輯距離後綴自動機（Suffix Levenshtein automata）、Soundex。在 bioinformatics 也有很多特化演算法。我這次簡單使用了 Classical Levenshtein distance 和 Damerau–Levenshtein distance 來解決這個題目。

預計步驟：

資料清理：去除空格、小寫化
- 使用正則 \[A-Za-z]+ 撈出單字
使用 Levenshtein distance 的算法來篩選

import pandas as pd
import re

notes = pd.read_excel("notes.xlsx",
                      index_col=0,
                      engine='openpyxl')

pattern = r'[A-Za-z]'
texts = [[t.lower() for t in re.findall(pattern, p)]
         for p in list(notes['diagnosis'])]

上面程式碼將每個病人的診斷從 Excel 中讀出，通過巢狀的 list comprehension 轉換成單字列表。

from Levenshtein import distance
Text = list[str]

def compare(word: str, texts: Text, n: int):
    for text in texts:
        result = distance(text, word)
        if result <= n:
            return(result, text)
    return None

上面程式碼通過 Levenshtein.distance() 這個函數來比對每個單字與目標字的距離，小於閾值 n 的單字則返回。
由於目標是「偵測」，所以用的是非貪婪搜尋模式（也就是找到第一個即返回，剩下的字忽略）。
另用 type hint 標明型別。
如果改用 Damerau-Levenshtein distance 的話，則可以如下改寫 compare（把 import 部分替換掉即可）。使用的是叫 fastDamerauLevenshtein 的模組 link

from fastDamerauLevenshtein import damerauLevenshtein as distance

最後就是比較每個病人的診斷資料所構築的字串列表與關鍵字群。當然，這操作也可以簡化成 list comprehension，方便將資料寫回 Pandas dataframe。

keywords = {"pneumonia": 3,
            "dyspnea": 3}

for index, text in enumerate(texts):
    for keyword in keywords.keys():
        result = compare(text, keyword, keywords[keyword])
        if result is not None:
            print(index, result)

Enthought 生態系

最近回老家，發現了一本 2013 年買的《用 Python 做科學計算》（現在有 Gitbook 的版本。這本書其實是一個大自然的搬運工，因為看起來是翻譯或改寫工具的介紹文件）。前半段介紹的是 Python 科學計算的底層 - NumPy 和 SciPy，以及重要的兩個核心工具 - SymPy（符號運算）和 Matplotlib（繪圖）。

書中後半部分則介紹了很多庫（lib） … 仔細一看，它們其實都來自於 Enthought 這間公司。他們著名的產品是 Python(X,Y)（現在改名 Enthought Tool Suite），這是一套類似 Conda + Spyder 的環境，給他們的客戶一個架構好的資料處理環境而不用去煩惱相依性或是功能的問題，本質上是 Matlab 和 Rstudio 的競爭對手。

作者介紹了如下的幾個由 Enthought 開發的庫：

Python 類型工具 - Traits
使用者介面 - TraitsUI
互動式資料呈現 - Chaco
3D 資料呈現 - Mayavi
VTK 包裝 - Traited VTK

不過，其實我自己現在做的也不是科學計算：以後若開始接觸流體相關的才是。總之，這篇文章就權充記錄。Enthought 的這些工具看起來都是開源，商業模式應該是仿照紅帽，即：開源以擴大使用群眾、然而優先滿足付費使用者的環境建構及開發需求。

開源統計工具

前幾天因緣際會稍微詢問了推友統計工作站。商業的軟體除了老牌的 SPSS 之外，還有最近熱門的 Graphpad Prism 以及新思維蔡校長獨愛的 MedCalc。開源部分，搜集了幾個建議如下：

不過我覺得軟體「好不好操作」是一回事，因為重點還是理解各種分佈的適用情境。用 R 來跑統計其實沒有比拖拉這些元件慢多少，而且可維護性更好。

ROI 備用匯出方案：mask 製備及 ROI 獲取

得到 ROI 座標後，下一步就是取得封閉多邊形內的每個座標點其 CT 密度（HU value）；在我原始的 pyOsiriX 腳本內，這件事情直接呼叫 OsiriX 給的。現在當然就得自幹了。（不知道偉大的 OpenCV 有沒有實作好了的！）。如果簡化成算法，就是要先回答「如何快速決定一個點在封閉多邊形內部或外部」，然後把所有在多邊形內的點的座標拿去存取原始圖片，就能得到 ROI 區域影像了。這邊搜尋了一些資料後，找到幾種經典教科書解法，例如 ray-tracing，而 Randolph Franklin 的 PNPOLY 算法應該是效率最好的。不過，這些都是求解數學問題的，實際應用上，對於像素圖型，會有各種共線問題。

繼續搜尋，發現了其他的圖學解決方案：

OpenCV 的 pointPolygonTest() 。StackOverflow 參考文章
Shapely 的 Polygon.contains() 。StackOverflow 參考文章

我其實想多用 OpenCV，熟悉 OpenCV 的 C++ 風格處理方式。然而，因為 pointPolygonTest 需要先把座標轉換成它某個資料格式 Contour，比較麻煩。Shapely 的解法直截了當。

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

point = Point(256, 256)
polygon = Polygon(ROIs[0][21])
polygon.contains(point)
# True

有了這個方法，我就可以構建一個函數，把 512 x 512 的矩陣傳入，就能拿到該 ROI 對應的 mask matrix。實際操作上是建構一個 512 x 512 個點類別（shapely.geometry.Point）的序列，用 polygon.contains 遍歷，再把處理後的序列 reshape 來達成 vectorization 的目的。接著是套用到原始影像取得 ROI。步驟異常簡單，因為前面搞定的 mask 是一個 512 x 512 的矩陣，跟原始影像一樣大小，所以直接用基本的矩陣乘法就可以了。

值得注意，這裡還要考慮到 DICOM 檔案存檔時候的平移（為了降低檔案大小做的像素數值調整），需要通過 DCM tag 改回來，否則 HU 值會有很大的誤差。

from pydicom import dcmread
ds = dcmread(dicomfile)
real_HU = ds.pixel_array * ds.RescaleSlope + ds.RescaleIntercept
roi = real_HU * mask

當然，因為 ROI 有時候會有兩個（我們的研究題材是大腿的表淺跟深層組織在某些疾病下的差異），所以需要大圈剪掉小圈。判斷哪個 ROI 序列代表外圈的步驟只需要一次。之後就用 subtraction 來取得內圈。

if sum(ROIs_a[0].flatten()) > sum(ROIs_b[0].flatten()):
    ROI_out = ROIs_a
    ROI_int = ROIs_a - ROIs_b

後續就可以接上之前寫好的程式進行 composition analysis 相關研究了。所產生的 ROIs 也可以另外存起來當作機器學習使用。

ROI 備用匯出方案：存取輸出檔

首先通過 Horos 公開的程式碼，了解 OrisiX MD 內的 ROI 曲線是通過 Spline interpolation 生成；然後可從官方匯出功能得到 ROI 標注點的座標。所以，我們用 Numpy 照樣刻一套。StackOverflow 參考文章

import numpy as np
from scipy.interpolate import splprep, splev

nodes = np.array(point_pairs)
x = nodes[:,0]
y = nodes[:,1]

tck, u = splprep([x, y], s=0, per=True)
xi, yi = splev(np.linspace(0, 1, 1200), tck)

接著用 Shoelace formula 計算面積。StackOverflow 參考文章

S1 = np.sum(xi * np.roll(yi, -1))
S2 = np.sum(yi * np.roll(xi, -1))
area = .5 * np.absolute(S1 - S2)

面積部分是計算來說服老闆，home-made software 的計算效果跟 OsiriX MD 針對 ROI 區域的面積一樣。不過，實際上誤差大約 0.1% 到 0.05% 之間，原因我還沒找到。效果大概如圖所示：

ROI

藍色是我生成的。紅色的點是研究助理圈選 ROI 的時候標注的。綠色的線是 OrisirX MD 自己生成的。

Reeder 5 and Goodlinks

升級了 Reeder 到 5。作爲 RSS reader 愛用者，Reeder 這種擠牙膏式的升級模式還真是令人討厭卻又不得不買（因爲它是目前 macOS / iOS 上能找到最好的閱讀器）。這次升級的亮點是終於能在 App 之間互通訂閱清單，總算可以離開 Inoreader 以及其嵌入式的網頁廣告了。

Goodlinks 則是一個稍後閱讀的服務。兩者的共同點在於：所有的資訊都是通過 iCloud 同步，也就是說，除了 Apple 的核心雲建設，基本上我不需要其他工具。

這兩個工具目前構成了我網路閱讀的主力。通常零散的網頁都是存到 Goodlinks，而部落格則是通過 Reeder 5 來追蹤。

Rust 寫的 Leetcode-cli 工具

其實比較熱門的是 NodeJS 寫的同名工具，也有豐富的插件，但是不知爲何我總是無法登入，所以改用 Rust 寫成的 Leetcode-cli。注意：具體如何存取 Cookie 請參考官方頁面的介紹。

search 搭配參數（例如 easy）搜尋問題列表。通過 pick 觀看問題敘述。
edit 來編輯問題。可以設定要用什麼編輯軟體。我是會開啟一個 Emacs 作爲編輯。（本來想用 Emacsclient 但是不知道爲何一直失敗）
test 以官方案例測試（不計入記錄）
submit 提交（會列入提交記錄）。最後會返回成績；如果失敗的話也會返回失敗的其中一個結果。

整體來說，我覺得這個工具很棒！分幾個面向：

程式碼統一保存，可以在其他筆記中引入這些程式碼。我是把 leetcode-cli 的 code 移到 Dropbox 然後通過 soft-link 放回去原本的 .leetcode。
可通過其他測試工具添加其他測試條件（例如：第 7 題需要檢測 int32 轉換的 overflow，但是官方的測試案例並沒有包含這一塊）。
版本管理。

使用 `unittest`

補充上述第二點。將 .leetcode/leetcode.toml 中檔案名字改成單純的 ${slug}，除去檔案開頭的 ID 以及分格檔名用的 .。之後使用下列的格式匯入（有點醜陋，但是這是檔名中有 dash 的解法）（以第 7 題 reverse-integer 舉例）：

import unittest
Solution = __import__("reverse-integer").Solution()

class SolutionTestCase(unittest.TestCase):
    def test_official(self): # 用來測試官方結果
        self.assertEqual(21, Solution.reverse(120))

    def test_personal(self): # 測試自己添加的案例
        self.assertEqual(0, Solution.reverse(231))

然後可以通過執行器 TextTestRunner 執行測試，例如：

suite = (unittest.TestLoader()
                 .loadTestsFromTestCase(SolutionTestCase))
unittest.TextTestRunner(verbosity=2).run(suite)