YFWu's Notes

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

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

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

  1. OpenCV 的 pointPolygonTest() 。StackOverflow 參考文章
  2. 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 也可以另外存起來當作機器學習使用。