昨天的攝影留下的待辦是檢查 Condat 與 Hirakawa 的排名在雜訊約 3.7 的地方會不會翻轉。作者的 2×3 六色程式找不到,我改成自己重寫:Bayer、Hirakawa 2×4、Condat 2×3 三種圖樣,每個相位各訓練一個線性最小平方重建器,在 Kodak 24 張上掃雜訊標準差 σ(0 到 255 的尺度)從 0 到 20。結果是翻轉真的出現,而且視窗取 13×13 時,差距的三個點和論文對得上:σ=0 時 Hirakawa 領先 0.24 dB(論文 Table I 是 0.25),σ=4 時兩者打平到 0.03 dB 以內(論文寫 37.60 對 37.58),σ=10 時 Condat 領先 0.13 dB(論文 0.16)。不過視窗縮到 5×5,翻轉就消失了,Condat 從 σ=0 就領先。昨天另一個待辦,取得 RR-08-04,也完成了:它是頻率結構的表示法報告,沒有雜訊模型,也沒有最佳化目標,所以無法回答昨天的疑問。
📖 學
一、兩個待辦的處理結果
RR-08-04 的標題是〈Color Filter Arrays: Representation and Analysis〉,作者 Yan Li、Pengwei Hao、Zhouchen Lin,2008 年 5 月,Queen Mary 電腦系技術報告,共 20 頁,我用 Wayback 的存檔下載。文字層是亂碼(字母被系統性替換),但摘要與表格的符號還讀得出來。摘要的內容是:提出用「頻率結構」表示任意矩形週期的 CFA 圖樣,頻率結構記錄亮度與色度成分被 CFA 濾過後在各頻率點的多工關係,可由圖樣的符號離散傅立葉變換求得,並用幾個例子說明 CFA 的特性看頻率結構就能懂。我看到的一張表,列出 Bawolek(1999)、Hamilton(2001)、Gindele(2002)、Bean(2003)、Roddy(2006)、Hirakawa(2007)、Kodak(2007)、FillFactory 等圖樣和各自的頻率結構。我用關鍵字數過:「noise」只出現在一條參考文獻標題(藍雜訊圖樣)裡,沒有雜訊模型、沒有實驗、沒有 的最佳化。所以昨天問的「TIP 2011 的目標是不是問題 (15)」,RR-08-04 不是答案的所在,我讀到的只有摘要、表格與關鍵字,不是逐頁讀完。
那張表有一個我之前沒注意的地方:Hirakawa(2007)被寫成 2×4 的 M Y C G 接 G C Y M(青、洋紅、黃、綠),每格總透光量為 1,和我 10-03 從 Condat 論文的增益反推出的那一組四色(每格總和 1.5,顏色如 (½, 1, 0))不是同一個圖樣。後者是 Condat 拿來當「Hirakawa 圖樣」比較的版本。今天一併把表中的版本也跑了模擬,見第四節。
Condat 的程式碼:他的軟體頁 lcondat.github.io/software.html 今天打得開(昨天寫作者頁面 404,我不確定當時試的是哪個網址),頁面列的是全變分與 Tikhonov 正則化的去馬賽克、去雜訊程式,以及隨機 RGB 圖樣的產生程式,沒有 2×3 六色圖樣的程式。所以重跑只能靠自己寫。另外要更正一件事:10-05 寫的「翻轉點約 σ=3.7」,現在回頭核對,它就是用論文 σ=0 的差距(Hirakawa 領先 0.25)與 σ=4 的差距(Condat 領先 0.02)做線性內插,;論文自己只寫了「優勢在 σ 低到 4 仍成立」,沒有給翻轉點,3.7 是內插,不是論文的數字。
二、模擬怎麼做
三種圖樣的每格透光率(R, G, B)如下:Bayer 是純色;Hirakawa 2×4 是 (½,1,0)、(0,1,½)、(1,0,½)、(½,0,1) 一列,下一列為 (½,0,1)、(1,0,½)、(0,1,½)、(½,1,0);Condat 2×3 由論文附錄的式 (26) 在 時產生,每格為 ,,實算每格 R+G+B 都是 1.5、範圍在 0 到 1 之內。Hirakawa 與 Condat 的總透光量是 Bayer 的 1.5 倍,這在後面的分解會用到。
資料與估計器:Kodak 24 張,把每張當真值,逐格乘透光率相加成單通道感測值,再加零均值高斯雜訊(σ 以 0 到 255 為尺度,三種圖樣用同一組雜訊,所以是配對比較)。重建器對每個相位各解一組嶺迴歸:輸入是以該像素為中心的視窗內所有感測值加常數項,輸出是該像素的 R、G、B。訓練和測試用兩折交叉(前 12 張訓練、後 12 張測試,反過來再一次),每個相位每張訓練 3000 個樣本;訓練時的 σ 與測試相同,所以每個 σ 都是各自最佳的線性重建器。這個設定與昨天 2×2 實驗相同,驗證:σ=0 時 Bayer 是 37.38 dB,昨天是 37.36;σ=10 時 30.06,昨天 30.05。
下面是精簡的單檔版本(貼上就能跑,影像先從 r0k.us 的 Kodak 頁抓到 kodak 資料夾;兩折都跑,結果與我的完整版差距在 0.07 dB 以內,應是抽樣隨機性造成):
import numpy as np
from PIL import Image
S3 = 3 ** .5
BAYER = [[(0,1,0),(1,0,0)], [(0,0,1),(0,1,0)]]
HIRA = [[(.5,1,0),(0,1,.5),(1,0,.5),(.5,0,1)], [(.5,0,1),(1,0,.5),(0,1,.5),(.5,1,0)]]
ph = 2*np.pi/3/4 - np.pi/2
CONDAT = [[tuple(.5 + (-1)**k1*np.sin(2*np.pi/3*(1-i) + 2*np.pi/3*k2 - ph)/S3
for i in range(3)) for k2 in range(3)] for k1 in range(2)]
imgs = [np.asarray(Image.open(f'kodak/kodim{i:02d}.png').convert('RGB'), float) for i in range(1, 25)]
def cpsnr(cfa, sigma, half=6, seed=0): # half=6 -> 13x13 視窗
h = np.array(cfa, float); ny, nx, _ = h.shape; r = half
dy, dx = np.mgrid[-r:r+1, -r:r+1]
rng = np.random.default_rng(seed)
def mosaic(im):
t = np.tile(h, (im.shape[0]//ny+1, im.shape[1]//nx+1, 1))[:im.shape[0], :im.shape[1]]
return (im*t).sum(2) + rng.normal(0, sigma, im.shape[:2])
def pick(M, im, py, px, n):
ys = np.arange(r+(py-r) % ny, M.shape[0]-r, ny); xs = np.arange(r+(px-r) % nx, M.shape[1]-r, nx)
yy, xx = [a.ravel() for a in np.meshgrid(ys, xs, indexing='ij')]
s = rng.choice(len(yy), n, replace=False); yy, xx = yy[s], xx[s]
P = M[yy[:,None,None]+dy, xx[:,None,None]+dx].reshape(n, -1)/255
return np.hstack([P, np.ones((n,1))]), im[yy, xx]/255
Ms = [mosaic(im) for im in imgs]; mse = []
for tr, te in [(range(0,12), range(12,24)), (range(12,24), range(0,12))]:
W = {}
for py in range(ny):
for px in range(nx):
XY = [pick(Ms[i], imgs[i], py, px, 3000) for i in tr]
X = np.vstack([a for a, b in XY]); Y = np.vstack([b for a, b in XY])
W[py, px] = np.linalg.solve(X.T@X + 1e-3*len(X)/1000*np.eye(X.shape[1]), X.T@Y)
for i in te:
e = 0
for py in range(ny):
for px in range(nx):
X, Y = pick(Ms[i], imgs[i], py, px, 20000//(ny*nx))
e += (((X@W[py, px]-Y)*255)**2).mean()/(ny*nx)
mse.append(e)
return 10*np.log10(255**2/np.mean(mse)) # 10*log10(平均 MSE)我跑這個版本的結果:σ=0 時 Bayer 38.03、Hirakawa 39.42、Condat 39.19;σ=4 時 34.59、36.50、36.52。
三、7×7 視窗的完整掃描
主表用 7×7 視窗(49 個感測值加常數項),σ 取 0 到 20。「差距」是 Condat 減 Hirakawa,以 10·log10(平均 MSE) 計;「逐張」是 24 張各自 dB 差的平均(也就是論文 Table 平均列的算法)與標準誤,最後是 Condat 逐張勝出的張數。
| σ | Bayer | Hirakawa | Condat | 差距(平均 MSE) | 逐張差(標準誤) | Condat 勝 |
|---|---|---|---|---|---|---|
| 0 | 37.38 | 38.96 | 38.85 | −0.115 | −0.086(0.045) | 6/24 |
| 2 | 36.30 | 38.05 | 38.01 | −0.045 | −0.012(0.037) | 10/24 |
| 3 | 35.32 | 37.18 | 37.18 | +0.005 | +0.039(0.034) | 15/24 |
| 4 | 34.33 | 36.23 | 36.27 | +0.048 | +0.075(0.030) | 18/24 |
| 6 | 32.56 | 34.43 | 34.56 | +0.134 | +0.154(0.022) | 23/24 |
| 10 | 30.06 | 31.68 | 31.90 | +0.211 | +0.223(0.018) | 24/24 |
| 20 | 26.67 | 27.77 | 28.07 | +0.301 | +0.319(0.019) | 24/24 |
兩個圖樣的差距隨 σ 單調地從負變正。翻轉點的位置取決於怎麼平均:用平均 MSE 的 dB,在 σ≈2.92 穿過零;用逐張 dB 的平均,在 σ≈2.27。逐張勝負張數(12 張是對半)穿過 12 的位置在 σ≈2.3。三個數都落在 2 到 3 之間,但互相差 0.6,所以「翻轉點」這個詞要連同平均方式一起報。低雜訊那一端的「Hirakawa 領先」在 7×7 下只有 0.09 dB(標準誤 0.045),σ=2 與 3 的逐張差距都在一個標準誤上下,證據不算強;σ 到 6 以上才有 23 到 24 張同向,那一段是穩的。
四、視窗大小會改變結論
同一組實驗改視窗,只看 Condat 對 Hirakawa 的逐張 dB 差:
| 視窗 | σ=0 | σ=3 | σ=4 | σ=10 | 翻轉 |
|---|---|---|---|---|---|
| 5×5 | +0.182 | +0.224 | +0.233 | +0.333 | 沒有,全程 Condat 領先 |
| 7×7 | −0.086 | +0.039 | +0.075 | +0.223 | σ≈2.3(逐張)/2.9(平均 MSE) |
| 9×9 | −0.132 | +0.029 | +0.074 | +0.184 | σ≈2.9(平均 MSE) |
| 13×13 | −0.240 | −0.028 | +0.030 | +0.130 | σ≈3.5(逐張)/4.0(平均 MSE) |
| 論文(13×13 LS 濾波器,Table I/II) | −0.25 | +0.02 | +0.16 | 內插 3.7 |
視窗越大,Hirakawa 在無雜訊端的領先越大,翻轉點也越往高 σ 移。13×13 的三個點(−0.24、+0.03、+0.13)和論文的三個點(−0.25、+0.02、+0.16)只差 0.03 dB 以內,看起來吻合得很好,但我不想把它當成「重現」:論文的 σ=4 與 σ=10 是用 BM3D 去雜訊的聯合流程,σ=0 是頻率選擇的 13×13 LS 濾波器;我的是純線性重建器,沒有任何去雜訊步驟,絕對水準也低(例如 σ=0 時 Hirakawa 是 39.4,論文是 40.36)。差距吻合,可以說兩個不同的重建管線給出同一個方向與相近量級,不能說我的管線等於論文的。
5×5 視窗的反轉有一個我不確定的解釋:Hirakawa 圖樣水平週期是 4,色度載波在 ,視窗只有 5 格寬,幾乎只裝得下一個週期,載波分離得不乾淨;Condat 的週期是 3,同樣的視窗寬度相對較夠用。這是我的推測,我沒有做能區分它與其他原因的實驗。確定的是:「翻轉」不是圖樣的內在性質,而是「圖樣加重建器」這一對的性質,論文自己也提醒過比較結果可能依重建流程而定。
另外把 RR-08-04 表中的 Hirakawa(2007)也跑了(M Y C G 接 G C Y M,透光量總和 1,7×7 視窗):σ=0 時 38.44 dB(比 Bayer 高 1.07),σ=3 時 35.55(高 0.23),σ=10 時 29.43(比 Bayer 低 0.63)。同樣叫 Hirakawa 的圖樣,總透光量 1.5 的版本在 σ=10 時比 Bayer 高 1.63 dB,總透光量 1 的版本卻低 0.63 dB。這個圖樣我是從亂碼文字層的表格抄出來的,抄錯顏色的風險存在,所以只當一個提醒:文獻裡的「Hirakawa 圖樣」不止一個,比較前要確認是哪一個。
五、光量與結構的分解
Hirakawa 與 Condat 的總透光量是 1.5,Bayer 是 1。在加性雜訊模型裡(雜訊大小與訊號無關),把訊號乘 1.5 等於把雜訊除以 1.5,所以我另外跑了一個對照組「等光量 Bayer」:Bayer 圖樣,雜訊改成 。它不是真實存在的濾片,而是「Bayer 多曝光 級」的加性模型等價物,忽略飽和與散粒雜訊。用它可以把 Condat 對 Bayer 的差距拆成兩項:
| σ | Condat − Bayer | 等光量 Bayer − Bayer(光量的貢獻) | Condat − 等光量 Bayer(結構的淨貢獻) |
|---|---|---|---|
| 0 | +1.47 | 0.00 | +1.47 |
| 3 | +1.86 | +0.98 | +0.88 |
| 6 | +2.00 | +1.76 | +0.24 |
| 8 | +1.91 | +1.94 | −0.03 |
| 10 | +1.84 | +2.01 | −0.17 |
| 20 | +1.40 | +1.94 | −0.54 |
結構的淨貢獻在 σ≈7.8 穿過零。σ 小的時候,Condat 的好處幾乎全來自取樣結構(σ=0 時光量沒有用處,我推測是混疊較少,但沒有另外驗證);σ 超過 7.8,把同樣多的光量「直接給 Bayer」比換成六色圖樣更划算。完全的光量優勢是 dB,但等光量 Bayer 最多只多了 2 dB 左右,因為 MSE 裡還含有不受雜訊影響的重建誤差。這張表告訴我昨天說的「範數越小畫質越好」在這裡需要加一個前提:範數只刻畫雜訊被放大多少,沒有算「這個圖樣多收了多少光」,而後者在高雜訊端是主要的好處。這個分解的限制要講清楚:對照組是我自己設計的假想物,雜訊加在 gamma 編碼後的 0 到 255 值上(真實感光元件的雜訊在線性域),所以 σ=2.9 或 7.8 不能直接換算成某個 ISO。
六、仍未解決
- 論文 Table II 的聯合去雜訊流程我沒有重現(BM3D 沒有跑),σ=4 與 σ=10 的吻合有一部分可能是巧合。
- 2×3 圖樣的作者程式找不到,我用的是自己從論文式 (26) 寫的圖樣與自己的重建器。
- RR-08-04 的全文我只讀了摘要、表格與關鍵字;文字層亂碼,表中各圖樣的頻率結構矩陣沒有逐一核對。
- 昨天的 Table II「Sony」欄與 CFA23 的係數和兩個疑問,今天都沒有進展。
🧠 記
翻轉點要連同平均方式與重建器一起報。 同一組模擬,用平均 MSE 的 dB 是 σ≈2.9,用逐張 dB 的平均是 σ≈2.3(7×7 視窗);視窗換成 13×13,這兩個數字變成 4.0 與 3.5。論文的 σ=4 是「優勢仍成立」的最低測試點,不是翻轉點。
線性重建器下,視窗越寬 Hirakawa 越強。 5×5 時 Condat 全程領先,13×13 時 σ=0 領先差距是 −0.24 dB,與論文的 −0.25 對得上。圖樣的優劣要在足夠寬的重建視窗下比較。
「加性雜訊下的好處」要把光量拆出來看。 Condat 比 Bayer 多的 1.5 倍光量,在 σ>7.8 時就是它全部的好處,結構本身反而淨虧;σ 小的時候是結構在賺。
同名圖樣可能是不同圖樣。 「Hirakawa」在 RR-08-04 是總透光量 1 的 M Y C G,在 Condat 的比較裡是總透光量 1.5 的四色,σ=10 時對 Bayer 的差距是 −0.63 與 +1.63 dB。
邊界:這些結論只在我的設定內成立:加性高斯雜訊(不含散粒雜訊)、加在 gamma 編碼值上、每相位一個線性最小平方重建器(沒有邊緣自適應、沒有去雜訊)、Kodak 24 張(768×512,影像內容偏向自然與建築)、兩折交叉。換成非線性解馬賽克、真實的 RAW 雜訊(散粒與讀出並存)或不同影像集,翻轉點與排序都可能移動。等光量 Bayer 是假想對照,不是任何市售濾片。
✍️ 實踐
一、用自己的相機量「曝光多 0.585 級」對訊噪比的影響(約 20 分鐘)。 這是把第五節的光量貢獻搬到實拍上。
- 器材:三腳架、18% 灰卡或一面均勻的牆、可輸出 RAW 的相機。光源用窗光或白熾燈(台灣市電 60 Hz,燈光亮度每秒閃 120 次,1/60 秒是 2 個週期、1/40 秒是 3 個週期,兩個快門都是整數個週期,不會拍到閃爍條紋;若用 LED 燈,先連拍兩張看平均亮度是否一致,不一致就換成窗光)。
- 設定:手動曝光,光圈 f/4,ISO 3200,關閉高 ISO 與長曝光降噪,鏡頭略失焦讓灰卡均勻。做兩組:暗組先調燈光或距離,使 1/60 秒時灰卡的 RAW 平均約為滿載值的 5%;亮組使 1/60 秒時約為滿載值的 40%。每組都是先用 1/60 秒拍兩張,再用 1/40 秒拍兩張(相機的 1/40 實際比 1/60 長多少,以量到的平均值比為準)。
- 計算:讀 RAW(例如 Python 的 rawpy),每個色彩通道各取畫面中央 200×200 的平面,先扣掉黑電平。訊號 取兩張的平均,雜訊取 (兩張相減可抵消固定圖樣與暈影),訊噪比 。
- 預測與判準:兩組的 比應該落在 1.45 到 1.55 之間;訊噪比增加量落在 1.76 dB(散粒雜訊主導,理論值 )與 3.52 dB(讀出雜訊主導,理論值 )之間。把變異數對訊號做兩點直線 :,若 就是讀出雜訊主導,增加量應靠近 3.52。我預期暗組靠近 3.52、亮組靠近 1.76,如果兩組的增加量差不到 0.5 dB,表示 ISO 3200 的設定下兩組落在同一個雜訊區間,把暗組再調暗或亮組再調亮重拍。證偽條件:訊噪比增加量小於 1.5 dB 或大於 3.7 dB,先查有沒有忘了扣黑電平、有沒有達到飽和、相機是不是在 RAW 裡偷做了降噪; 比不在 1.45 到 1.55 之內,就是光源在閃或相機自動調了東西,這組數據作廢重拍。
二、重現翻轉(約 10 分鐘)。 貼上第二節的程式,把視窗半徑 half 設成 3(7×7)與 6(13×13),各算 σ=0、3、4 三個點。預測:7×7 時 σ=3 的 Condat 減 Hirakawa 是正的(約 +0.005 dB,逐張約 +0.04),σ=0 是負的(約 −0.11 dB);13×13 時 σ=3 仍是負的(約 −0.06 dB),σ=4 約為 0。證偽條件:若你的 7×7 在 σ=0 就是正的,先檢查圖樣的列順序有沒有照第二節輸入。差距小於 0.05 dB 的格子,換一個隨機種子(seed 參數)可能就變號,請各跑三個種子再下判斷。
🔗 延伸學習
- Condat,〈A new color filter array with optimal properties for noiseless and noisy color image acquisition〉(作者網站預印本 PDF) — 今天對照的 Table I、Table II 與「優勢在 σ 低到 4 仍成立,37.60 對 37.58」那句話的出處。
- Condat 的軟體頁 — 列有全變分與 Tikhonov 去馬賽克/去雜訊的 Matlab 程式與隨機 CFA 程式;沒有 2×3 六色圖樣的程式,今天確認過。
- RR-08-04〈Color Filter Arrays: Representation and Analysis〉(Wayback 存檔) — 頻率結構的表示法報告,含各家 CFA 的頻率結構對照表;文字層亂碼,表格可讀。
- Kodak 無損影像集(24 張) — 模擬用的資料,要加
-L追蹤轉址才能下載。
💬 問 AI
我用線性重建器比較三種 CFA 在加性高斯雜訊(σ 以 0 到 255 為尺度)下的 CPSNR:Bayer、Hirakawa 2×4(四色,每格總透光量 1.5)、Condat 2×3(六色,每格總透光量 1.5)。重建器是每個相位各解一組嶺迴歸(輸入為視窗內所有感測值),Kodak 24 張,兩折交叉。我的結果:視窗 13×13 時 Condat 減 Hirakawa 在 σ=0 是 −0.24 dB、σ=4 約 0、σ=10 是 +0.13 dB;視窗 5×5 時全程為正;另外若把 Bayer 的雜訊改成 σ/1.5(等光量對照),它在 σ≈7.8 以上勝過 Condat。
(1) 請檢查「5×5 視窗讓 Hirakawa 變差」是否可以用其色度載波在 (π, ±π/2)、水平週期 4 解釋,並說明如何設計實驗區分這個解釋與其他可能(例如嶺迴歸係數、訓練樣本數)。
(2) 等光量 Bayer 的對照把「訊號乘 1.5」等同於「雜訊除以 1.5」,在散粒雜訊主導、或有飽和時這個等價怎麼失效?請給出會讓結論反轉的條件。
(3) 論文 Table II 的聯合去雜訊(BM3D)與我的純線性重建,各自的翻轉點為什麼可能不同?
(4) 請把回答分成三類逐條標示:「文獻明確寫的」「你查證到的(附出處)」「你自己的推論」;查不到就說查不到,不要補數字;不要給任何買賣器材的建議。