NIR 高光譜影像的化學計量學實作

1. 認識資料集：Art Paint Pigment Concentrations

這套資料由巴塞隆納大學 José F. García 教授與 Clarimma Sessa 提供（IASIM-10 工作坊資料集）： 24 個塗在畫布上的油畫顏料樣品，各含已知比例的三種藍色色素—— 普魯士藍（Prussian blue）、酞菁藍（Heliogen blue）、群青（Ultramarine blue）——再加上油性展色劑（oil binder）。

2. 互動探索：先「看見」資料

以下圖表全部由真實資料計算後內嵌於本頁（Cube B 的 24 個樣品平均光譜與影像），可直接操作。

3. 化學計量學原理：從光譜到濃度

不管你用 Python、R 還是 Orange，底層的化學計量學流程都一樣。四個工具只是同一條路的四種走法。

3.1 前處理：為什麼要 log(1/R)？

Beer–Lambert 定律說吸收度 A 與濃度 c 成正比（A = εlc）。反射式量測拿到的是反射率 R， 常用 A = log(1/R) 當「擬吸收度（apparent absorbance）」，讓訊號與濃度的關係更接近線性，利於之後的線性模型（PCA/PLS）。 其他常見前處理還有 SNV（標準常態變數，校正散射造成的乘法效應）、Savitzky–Golay 微分（去基線飄移、突顯峰形）等。

3.2 非監督探索：PCA

主成分分析（PCA）把 207 維光譜投影到少數幾個「最大變異方向」。分數圖（score plot）能一眼看出樣品分群、離群值與趨勢—— 但 PC 軸是「變異最大的方向」，不等於某個色素的濃度軸。

3.3 監督式建模：PLS 迴歸

偏最小平方法（PLS regression）在壓縮 X（光譜）時同時考慮 y（濃度）， 以少數幾個潛在變數（latent variables, LV）建立線性預測模型，是 NIR 定量的黃金標準。 LV 太少會欠擬合、太多會把雜訊學進去（過擬合），要用交叉驗證（cross-validation）選。

3.4 驗證設計：為什麼 A 建模、B/C 驗證？

同一樣品的相鄰像素高度相關。如果把所有像素混在一起「隨機切」訓練/測試集， 測試集裡的像素跟訓練集裡的像素根本來自同一塊顏料——這是資料洩漏（data leakage），會嚴重高估模型表現。 正確做法是照資料集設計：用一次拍攝（Cube A 或 B）建模，用另一次獨立拍攝（Cube C）驗證，才能反映真實的預測能力。 本頁示範：Cube B 訓練 → Cube C 逐像素預測，24 樣品平均濃度的 RMSEP ≈ 3.6%。

4. Python / Jupyter Notebook

最完整彈性的路線。需要 numpy、scipy、scikit-learn、matplotlib。 完整可執行的 notebook 在 examples/artpaint_pls.ipynb。

4.1 讀取 DSO .mat 與還原影像

# pip install numpy scipy scikit-learn matplotlib
import numpy as np, scipy.io as sio
import matplotlib.pyplot as plt

m    = sio.loadmat("PaintCubeB.mat")
dso  = m["PaintCubeB"][0, 0]          # DSO struct（structured array）
X    = dso["data"].astype(float) / 65536.0   # V = R*65536 → 反射率 R，(57600, 207)
wl   = np.ravel(dso["axisscale"][1, 0])      # 207 個波長 988.9–1674.7 nm
mask = m["PaintMask"].astype(int)            # (240, 240)，值 1–24

# MATLAB 是 column-major：像素列 ↔ 影像要用 order='F' 互轉
mean_img = X.mean(axis=1).reshape(240, 240, order="F")
plt.imshow(mean_img, cmap="viridis"); plt.title("mean reflectance"); plt.colorbar()

4.2 前處理與各樣品平均光譜

A = np.log10(1.0 / np.clip(X, 1e-6, None))     # log(1/R) 擬吸收度
mask_f = mask.reshape(-1, order="F")            # 與 X 的列對齊

S = np.array([A[mask_f == k].mean(axis=0) for k in range(1, 25)])  # (24, 207)
for s in S: plt.plot(wl, s)
plt.xlabel("wavelength (nm)"); plt.ylabel("log(1/R)")

4.3 y 值、PCA 與 PLS（留一交叉驗證）

from sklearn.decomposition import PCA
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import LeaveOneOut

ym   = sio.loadmat("ArtImageY.mat")
yTot = ym["yTotalPercent"][0, 0]["data"].astype(float)   # (57600, 4) 逐像素
y24  = np.array([yTot[mask_f == k].mean(axis=0) for k in range(1, 25)])
yP   = y24[:, 0]                                 # 第 0 欄 = Prussian %

scores = PCA(n_components=2).fit_transform(S - S.mean(axis=0))
plt.scatter(scores[:, 0], scores[:, 1], c=yP); plt.colorbar(label="Prussian %")

pred = np.zeros(24)
for tr, te in LeaveOneOut().split(S):
    pls = PLSRegression(n_components=5).fit(S[tr], yP[tr])
    pred[te] = pls.predict(S[te]).ravel()
rmsecv = np.sqrt(np.mean((pred - yP) ** 2))
print(f"RMSECV = {rmsecv:.2f} %")               # ≈ 4.9 %

4.4 逐像素預測 → 濃度分布圖（chemical image）

# 用 Cube B 的像素訓練，預測獨立拍攝的 Cube C
mC   = sio.loadmat("PaintCubeC.mat")
XC   = np.log10(65536.0 / np.clip(mC["PaintCubeC"][0, 0]["data"].astype(float), 1, None))

rng  = np.random.default_rng(0)
idx  = rng.choice(A.shape[0], 6000, replace=False)   # 抽 6000 個像素當訓練集就夠
pls  = PLSRegression(n_components=5).fit(A[idx], yTot[idx, 0])

pred_map = pls.predict(XC).ravel().reshape(240, 240, order="F")
plt.imshow(pred_map, cmap="inferno", vmin=0, vmax=35)
plt.colorbar(label="predicted Prussian %")

5. Google Colab：免安裝、雲端跑

程式碼與第 4 節完全相同（Colab 已內建 numpy/scipy/sklearn/matplotlib），差別只在「檔案怎麼進來」。.mat 檔共約 65 MB，建議放 Google Drive。

# 方法一（推薦）：掛載 Google Drive——檔案大、session 重啟不必重傳
from google.colab import drive
drive.mount("/content/drive")
DATA = "/content/drive/MyDrive/ArtImageDataA/"
m = sio.loadmat(DATA + "PaintCubeB.mat")

# 方法二：直接上傳（每次 session 重啟都要重傳，大檔不建議）
from google.colab import files
uploaded = files.upload()

6. R：R.matlab + pls 套件

R 生態的化學計量學主力是 pls（PLS/PCR）與 prospectr（光譜前處理）。讀 DSO .mat 用 R.matlab。

# install.packages(c("R.matlab", "pls", "prospectr"))
library(R.matlab); library(pls)

m    <- readMat("PaintCubeB.mat")
dso  <- m$PaintCubeB[, , 1]              # DSO struct → list
X    <- dso$data / 65536                 # (57600, 207) 反射率
wl   <- as.numeric(dso$axisscale[[3]])   # 波長（cell {2,1}；用 str(dso$axisscale) 確認位置）
mask <- m$PaintMask                      # 240×240
A    <- log10(1 / pmax(X, 1e-6))
maskf <- as.vector(mask)                 # R 本身就是 column-major，天生和 MATLAB 對齊！

# 各樣品平均光譜與 y
S   <- t(sapply(1:24, function(k) colMeans(A[maskf == k, ])))
ym  <- readMat("ArtImageY.mat")
yT  <- ym$yTotalPercent[, , 1]$data      # (57600, 4)
y24 <- t(sapply(1:24, function(k) colMeans(yT[maskf == k, ])))
yP  <- y24[, 1]                          # Prussian %

matplot(wl, t(S), type = "l", xlab = "wavelength (nm)", ylab = "log(1/R)")

# PCA 與 PLS（留一交叉驗證內建在 plsr）
pc  <- prcomp(S, center = TRUE)
plot(pc$x[, 1:2], col = "steelblue", pch = 19)

fit <- plsr(yP ~ S, ncomp = 8, validation = "LOO")
plot(RMSEP(fit))                          # 看 RMSECV vs LV 數，選 ~5 個 LV
plot(fit, ncomp = 5, line = TRUE)         # 預測 vs 實際

# 濃度分布圖：逐像素預測 Cube C 再排回 240×240
mC  <- readMat("PaintCubeC.mat")
AC  <- log10(65536 / pmax(mC$PaintCubeC[, , 1]$data, 1))
fitPix <- plsr(yT[idx, 1] ~ A[idx, ], ncomp = 5)   # idx = sample(57600, 6000)
predC  <- predict(fitPix, AC, ncomp = 5)
image(matrix(predC, 240, 240), col = hcl.colors(64, "Inferno"))

7. Orange Data Mining：能用 HyperSpectra widget 嗎？

7.1 安裝

1. Orange 選單 Options → Add-ons → 勾選 Spectroscopy → Install → 重啟 Orange；
2. 或 pip 環境：pip install orange-spectroscopy；
3. 或直接安裝 Quasar（已把 Orange + Spectroscopy 打包好，光譜課首選）。

7.2 資料進 Orange：先轉檔

Orange 讀不了 PLS_Toolbox 的 DSO .mat（它的 .mat 支援僅限單純矩陣），所以先用我們提供的轉檔腳本 examples/artpaint_to_orange.py 把 hypercube 轉成 Orange 原生 .tab 表格——每列一個像素、207 個波長欄，外加 map_x/map_y/sample/Prussian/Heliogen/Ultramarine/Oil meta 欄：

# 在 ArtImageDataA 資料夾內執行；--stride 2 = 每 2 像素取 1（120×120，檔案小、載入快）
python artpaint_to_orange.py PaintCubeB.mat --stride 2 -o artpaint_B.tab
python artpaint_to_orange.py PaintCubeC.mat --stride 2 -o artpaint_C.tab

7.3 建議 workflow

1. File 載入 artpaint_B.tab（轉檔腳本已把 Prussian 設為 target、map_x/map_y 設為 meta）；
2. HyperSpectra（Spectroscopy 分類下）：自動用 map_x/map_y 排出 120×120 影像，左邊光譜、右邊影像，滑鼠框選區域即可比較不同樣品的光譜；
3. Preprocess Spectra：試 Savitzky–Golay、Vector Normalization——記得用 Test and Score 比較前處理前後的 RMSE，呼應第 3.1 節「前處理用驗證選」；
4. PLS（Model 分類下，Orange 3.32+ 內建）：設定成分數（LV）；Test and Score 看 cross-validation 的 RMSE/R²；
5. 把 artpaint_C.tab 接到 Predictions，再把輸出接回 HyperSpectra，選擇以預測欄位上色——濃度分布圖完成，跟 Python/R 的結果對照。

8. 互動測驗：你掌握了嗎？

10 題單選，涵蓋化學計量學原理與 工具操作兩類。 作答完按「交卷評分」，每題都會給解說；可重複作答。