MATLAB實戰(zhàn)：OMP、CoSaMP與IRLS算法實現(xiàn)與調(diào)優(yōu))
簡介這份資源面向數(shù)字圖像處理、計算機視覺與醫(yī)學(xué)成像方向的學(xué)習(xí)者和研究者提供二維圖像壓縮感知重構(gòu)算法的完整MATLAB實現(xiàn)幫助讀者從代碼層面理解CS理論如何用遠少于像素數(shù)的測量值重建圖像。壓縮包共26個文件以18個m腳本為核心覆蓋測量矩陣生成、圖像采樣、ISTA/FISTA等重構(gòu)算法與性能評估流程另含7個fig結(jié)果圖和1個bmp測試圖像整體約1.59MB結(jié)構(gòu)緊湊便于逐文件調(diào)試。已有4295人學(xué)習(xí)下載說明其在教學(xué)與實驗場景中具備一定參考價值。讀者可據(jù)此對比高斯、伯努利等不同測量矩陣的重構(gòu)質(zhì)量調(diào)整松弛因子與閾值參數(shù)并借助MSE、PSNR指標和可視化結(jié)果驗證算法效果也可作為向視頻壓縮、醫(yī)學(xué)影像、遙感分析及深度學(xué)習(xí)結(jié)合方向拓展的實驗起點。1. 二維圖像壓縮感知重構(gòu)一份能直接跑通的 MATLAB 代碼包做圖像處理的人遲早會撞上壓縮感知這個坎。采樣率明明壓到了 30% 甚至更低重構(gòu)出來的圖卻要跟原圖肉眼難辨這中間的落差就是重構(gòu)算法要填的坑。這份資源給的不是論文里的偽代碼而是一套完整的 MATLAB 程序代碼專門針對二維圖像做壓縮感知重構(gòu)。它把測量矩陣生成、稀疏基構(gòu)造、重構(gòu)求解、PSNR 對比這幾步全串起來了拿到手改個圖片路徑就能出結(jié)果。適合正在做圖像處理大作業(yè)、想驗證 OMP 或 CoSaMP 這類算法實際效果的人也適合需要快速搭一個重構(gòu) baseline 再往上改的研究場景。MATLAB 版本不用追新R2016b 之后的都能跑核心函數(shù)沒用到太偏的工具箱。2. 壓縮感知重構(gòu)的數(shù)學(xué)骨架從測量矩陣到稀疏基的選型邏輯2.1 為什么二維圖像不能直接套一維壓縮感知一維壓縮感知的模型是 y Φxx 是長度 N 的列向量Φ 是 M×N 的測量矩陣M 遠小于 N。但二維圖像是 N×N 的矩陣如果直接按列拉伸成 N2×1 的向量測量矩陣的尺寸會膨脹到 M×N2內(nèi)存和計算量都吃不消。常見做法是分塊處理把圖像切成 32×32 或 64×64 的小塊每塊單獨做壓縮感知最后拼回去。另一種做法是用可分離的測量矩陣行方向和列方向各乘一個矩陣把二維運算拆成兩次一維運算。這份代碼走的是分塊路線塊大小在參數(shù)區(qū)可以改默認給的是 32×32兼顧了重構(gòu)質(zhì)量和速度。分塊之后有個細節(jié)容易被忽略塊與塊之間的邊界效應(yīng)。如果每塊獨立重構(gòu)再硬拼塊邊緣會出現(xiàn)明顯的方格感。代碼里用了一個簡單的重疊分塊策略塊間重疊 4 個像素重構(gòu)時對重疊區(qū)域做加權(quán)平均。這個處理不復(fù)雜但有沒有它視覺效果差很多。2.2 測量矩陣的選取高斯矩陣還是伯努利矩陣測量矩陣要滿足 RIP 條件實際中沒法直接驗證所以選的都是大概率滿足 RIP 的隨機矩陣。高斯矩陣是最常用的每個元素服從 N(0, 1/M)跟絕大多數(shù)稀疏基都不相關(guān)。伯努利矩陣元素取 ±1/√M硬件實現(xiàn)時更友好但重構(gòu)質(zhì)量在高壓縮比下會略遜于高斯矩陣。代碼里默認用的是高斯矩陣生成方式如下% 生成 M x N 的高斯測量矩陣M 為測量數(shù)N 為信號長度 function Phi generate_gaussian_matrix(M, N) % 每個元素服從 N(0, 1/M)保證列向量近似正交 Phi randn(M, N) / sqrt(M); % 對每一列做歸一化避免數(shù)值量級差異過大 Phi Phi ./ repmat(sqrt(sum(Phi.^2, 1)), M, 1); end這里做了列歸一化原因是后續(xù)重構(gòu)算法對矩陣列向量的能量敏感。如果不歸一化某些列的范數(shù)偏大會導(dǎo)致求解時權(quán)重失衡重構(gòu)出來的圖像局部過亮或過暗。歸一化之后每列能量為 1求解穩(wěn)定性明顯提升。參數(shù) M 就是測量數(shù)M/N 就是采樣率代碼里通過采樣率反算 M比如采樣率 0.3、塊大小 32×32 時N1024M307。2.3 稀疏基的選擇DCT 還是小波自然圖像在 DCT 基下是近似稀疏的在小波基下更稀疏。DCT 基構(gòu)造簡單就是余弦函數(shù)采樣計算快小波基需要濾波器組多一層分解但壓縮感知重構(gòu)質(zhì)量通常更好。代碼里兩種基都給了通過一個 flag 切換。DCT 基的構(gòu)造代碼% 構(gòu)造 N x N 的 DCT 正交基 function Psi generate_dct_basis(N) Psi zeros(N, N); for k 0:N-1 for n 0:N-1 if k 0 Psi(k1, n1) 1 / sqrt(N); else Psi(k1, n1) sqrt(2/N) * cos(pi * k * (2*n1) / (2*N)); end end end end這個雙重循環(huán)在 N1024 時會比較慢實際跑的時候可以用 dctmtx 函數(shù)替代MATLAB 自帶一行就能生成。代碼里保留手寫版本是為了讓讀者看清 DCT 基的每個元素怎么來的實際使用時建議換成 dctmtx。小波基用的是 sym4 小波通過 wavedec2 做兩層分解再把系數(shù)拉成列向量。注意小波分解后的系數(shù)長度跟原信號長度一致但排列順序會影響稀疏性代碼里做了按能量排序的處理把大系數(shù)排在前面這樣重構(gòu)時閾值收縮更有效。3. 重構(gòu)算法實戰(zhàn)OMP、CoSaMP 與 IRLS 的代碼實現(xiàn)與參數(shù)調(diào)優(yōu)3.1 OMP 算法的 MATLAB 實現(xiàn)與停止條件設(shè)置OMP 是正交匹配追蹤思路是每次迭代選一個跟殘差最相關(guān)的原子然后對已選原子做最小二乘更新。它的優(yōu)點是實現(xiàn)簡單、速度快缺點是需要知道稀疏度 KK 估不準就重構(gòu)不好。代碼里的 OMP 核心循環(huán)function x_rec omp_reconstruct(y, Phi, K) % y: 測量向量 M x 1 % Phi: 測量矩陣 M x N % K: 稀疏度 [M, N] size(Phi); residual y; support_set []; x_rec zeros(N, 1); for iter 1:K % 計算殘差與各原子的內(nèi)積找最大相關(guān)原子 correlations abs(Phi * residual); [~, idx] max(correlations); support_set [support_set, idx]; % 在已選原子集上做最小二乘 Phi_sub Phi(:, support_set); x_sub pinv(Phi_sub) * y; % 更新殘差 residual y - Phi_sub * x_sub; % 殘差能量低于閾值時提前退出 if norm(residual) 1e-6 break; end end x_rec(support_set) x_sub; end停止條件有兩個迭代次數(shù)達到 K或者殘差能量低于 1e-6。實際跑的時候如果 K 設(shè)得偏大OMP 會把噪聲也選進來重構(gòu)圖像出現(xiàn)顆粒感。K 設(shè)得偏小細節(jié)丟失。經(jīng)驗值是 K 取 N 的 10% 到 20%具體看圖像紋理復(fù)雜度。代碼里 K 是通過稀疏度估計函數(shù)自動算的用 DCT 系數(shù)排序后取累計能量 95% 的位置作為 K比固定值靠譜。3.2 CoSaMP 算法的回溯機制與參數(shù)影響CoSaMP 跟 OMP 的區(qū)別在于每次迭代選 2K 個原子然后回溯剪枝只保留 K 個。這個回溯機制讓它在稀疏度估計不準時更魯棒代價是每次迭代計算量翻倍。代碼里 CoSaMP 的關(guān)鍵步驟function x_rec cosamp_reconstruct(y, Phi, K) [M, N] size(Phi); residual y; support_set []; x_rec zeros(N, 1); max_iter 50; for iter 1:max_iter % 步驟1形成代理信號找 2K 個最大相關(guān)原子 proxy Phi * residual; [~, sort_idx] sort(abs(proxy), descend); new_support sort_idx(1:min(2*K, N)); % 步驟2合并支撐集 support_set union(support_set, new_support); % 步驟3在合并支撐集上做最小二乘 Phi_sub Phi(:, support_set); x_sub pinv(Phi_sub) * y; % 步驟4剪枝只保留 K 個最大系數(shù) [~, keep_idx] sort(abs(x_sub), descend); keep_idx keep_idx(1:min(K, length(keep_idx))); support_set support_set(keep_idx); x_sub x_sub(keep_idx); % 步驟5更新殘差 residual y - Phi(:, support_set) * x_sub; if norm(residual) 1e-6 break; end end x_rec(support_set) x_sub; endmax_iter 默認 50實際跑的時候大多數(shù)圖像在 20 次以內(nèi)就收斂了。如果跑到 50 次還沒收斂通常是測量矩陣條件數(shù)太差或者稀疏基選得不合適。回溯剪枝那一步是 CoSaMP 的核心它保證了支撐集大小始終不超過 K避免 OMP 那種選錯原子就一路錯到底的問題。代價是每次迭代要做兩次排序N 大時排序開銷不可忽略。3.3 IRLS 算法用迭代重加權(quán)最小二乘逼近 L0 范數(shù)IRLS 不走貪婪路線它用一系列加權(quán) L2 問題去逼近 L0 范數(shù)。每次迭代根據(jù)當前解的幅度調(diào)整權(quán)重幅度小的系數(shù)權(quán)重加大逼著它們往零走。function x_rec irls_reconstruct(y, Phi, max_iter, epsilon) [M, N] size(Phi); x_rec Phi * y; % 初始解用匹配濾波 weights ones(N, 1); for iter 1:max_iter % 構(gòu)造加權(quán)矩陣 W diag(weights); % 求解加權(quán)最小二乘min ||Wx||_2 s.t. Phi*x y % 用拉格朗日乘子法轉(zhuǎn)為無約束問題 x_rec W \ (Phi * ((Phi / W) * Phi \ y)); % 更新權(quán)重epsilon 防止除零 weights 1 ./ (abs(x_rec) epsilon); % 權(quán)重變化很小時退出 if norm(weights - 1./(abs(x_rec)epsilon)) 1e-4 break; end end endepsilon 是個關(guān)鍵參數(shù)默認給 1e-3。太大則稀疏性不夠太小則數(shù)值不穩(wěn)定。IRLS 的優(yōu)點是重構(gòu)精度通常比 OMP 高缺點是要解線性方程組每次迭代計算量比 OMP 大一個量級。圖像塊小的時候32×32還能接受塊大到 128×128 就明顯慢了。代碼里給 IRLS 設(shè)了最大迭代 30 次實際跑下來 15 次左右權(quán)重就基本穩(wěn)定了。4. 避坑與排查重構(gòu)圖像質(zhì)量不達標時先查這五處4.1 重構(gòu)圖像全黑或全白現(xiàn)象跑完代碼imshow 出來的圖要么全黑要么全白PSNR 低到個位數(shù)。原因測量矩陣沒有歸一化或者稀疏基的逆變換用錯了。全黑通常是重構(gòu)系數(shù)全為零全白是系數(shù)爆炸。檢查 generate_gaussian_matrix 里有沒有做列歸一化以及 DCT 基的轉(zhuǎn)置有沒有用對。解決在重構(gòu)完成后加一句 x_rec x_rec / max(abs(x_rec(:)))先歸一化再顯示。如果還是不對把測量矩陣換成 dctmtx 的前 M 行這個是確定性的能排除隨機矩陣的問題。4.2 PSNR 卡在 20dB 上不去現(xiàn)象采樣率提到 0.5 了PSNR 還是 20 出頭跟論文里的 35dB 差很遠。原因稀疏度 K 設(shè)得不對或者分塊大小太小。32×32 的塊在采樣率 0.5 時M512對紋理復(fù)雜的塊來說測量數(shù)還是不夠。另外 OMP 的 K 如果按固定值給遇到平滑塊會過擬合遇到紋理塊會欠擬合。解決把塊大小改成 64×64采樣率不變的情況下每塊測量數(shù)翻倍。K 改成自適應(yīng)用 DCT 系數(shù)累計能量 95% 來定。如果還不行換 CoSaMP它對 K 的估計誤差容忍度更高。4.3 重構(gòu)時間隨圖像尺寸指數(shù)增長現(xiàn)象512×512 的圖跑了十分鐘還沒出結(jié)果內(nèi)存占用飆升。原因沒有分塊直接把整圖拉成 262144×1 的向量測量矩陣尺寸 131072×262144內(nèi)存直接爆?；蛘叻謮K了但每塊都重新生成測量矩陣重復(fù)計算。解決確認代碼走的是分塊路徑塊大小不超過 64×64。測量矩陣生成一次就存下來所有塊共用。如果內(nèi)存還是緊張把測量矩陣改成稀疏存儲MATLAB 的 sparse 函數(shù)能省不少內(nèi)存。4.4 塊邊界出現(xiàn)方格狀偽影現(xiàn)象重構(gòu)圖像放大看塊與塊之間有明顯的亮度跳變像馬賽克。原因分塊時沒有重疊每塊獨立重構(gòu)后硬拼塊邊緣的系數(shù)不連續(xù)。解決啟用重疊分塊重疊像素數(shù)設(shè)為塊大小的 1/8 到 1/4。重構(gòu)后對重疊區(qū)域做加權(quán)平均權(quán)重用余弦窗中間高兩邊低。代碼里 overlap 參數(shù)默認是 4塊大小 32 時相當于 12.5% 重疊夠用了。4.5 換圖片后重構(gòu)完全失效現(xiàn)象用代碼自帶的測試圖跑得好好的換成自己的圖就重構(gòu)不出來。原因圖片尺寸不是塊大小的整數(shù)倍分塊時邊界處理出錯?；蛘邎D片是彩色三通道代碼只處理了灰度。解決讀圖后先做尺寸裁剪把長寬都裁成塊大小的整數(shù)倍。彩色圖先用 rgb2gray 轉(zhuǎn)灰度或者對三個通道分別做壓縮感知再合并。代碼里加一句 size_check 函數(shù)自動裁剪并給出提示。5. 從重構(gòu)到評估PSNR、SSIM 與采樣率曲線的自動化腳本跑通單張圖的重構(gòu)只是第一步真正做實驗的時候需要批量跑不同采樣率、不同算法然后畫對比曲線。這份代碼里帶了一個評估腳本把 PSNR 和 SSIM 的計算封裝好了直接調(diào)用就行。PSNR 的計算function psnr_val calc_psnr(original, reconstructed) % 確保數(shù)據(jù)類型一致 original double(original); reconstructed double(reconstructed); mse mean((original(:) - reconstructed(:)).^2); if mse 0 psnr_val Inf; return; end max_val max(original(:)); psnr_val 10 * log10(max_val^2 / mse); end注意 max_val 取的是原圖最大值不是固定的 255。如果圖像歸一化到 [0,1] 了max_val 就是 1PSNR 數(shù)值會跟 [0,255] 范圍差 20dB 左右對比論文數(shù)據(jù)時要統(tǒng)一量綱。SSIM 的計算用 MATLAB 自帶的 ssim 函數(shù)但要注意它默認是彩色圖模式灰度圖要加參數(shù)ssim_val ssim(reconstructed, original, DynamicRange, max(original(:)));批量實驗的腳本結(jié)構(gòu)是三層循環(huán)外層采樣率從 0.1 到 0.5 步進 0.1中層算法遍歷 OMP、CoSaMP、IRLS內(nèi)層遍歷測試圖片。每次跑完把 PSNR 和 SSIM 存到矩陣里最后用 plot 畫曲線。我一般會加一個計時器記錄每種算法的平均重構(gòu)時間。因為實際選算法的時候不只看質(zhì)量還要看能不能在可接受的時間內(nèi)跑完。OMP 在采樣率 0.3 時單塊重構(gòu)大概 0.05 秒CoSaMP 是 0.12 秒IRLS 要 0.8 秒。整圖 512×512 分 256 個塊OMP 總共 13 秒左右IRLS 要三分鐘以上。如果做實時性要求高的場景IRLS 基本不用考慮。還有一個技巧是固定隨機種子。測量矩陣是隨機生成的不同種子跑出來的 PSNR 會有 0.5dB 左右的波動。做對比實驗時把 rng(42) 寫在循環(huán)外面保證所有算法用的是同一個測量矩陣這樣對比才公平。這個細節(jié)論文里經(jīng)常不寫但實際復(fù)現(xiàn)時不知道就會覺得結(jié)果對不上。從那以后我每次跑壓縮感知實驗都先把隨機種子固定、量綱統(tǒng)一、塊大小和重疊參數(shù)確認一遍這三步走完再開始跑數(shù)據(jù)。希望幫到你。本文還有配套的精品資源點擊獲取 二維圖像壓縮感知重構(gòu)MATLAB實戰(zhàn)：OMP、CoSaMP與IRLS算法實現(xiàn)與調(diào)優(yōu)
簡介這份資源面向數(shù)字圖像處理、計算機視覺與醫(yī)學(xué)成像方向的學(xué)習(xí)者和研究者提供二維圖像壓縮感知重構(gòu)算法的完整MATLAB實現(xiàn)幫助讀者從代碼層面理解CS理論如何用遠少于像素數(shù)的測量值重建圖像。壓縮包共26個文件以18個m腳本為核心覆蓋測量矩陣生成、圖像采樣、ISTA/FISTA等重構(gòu)算法與性能評估流程另含7個fig結(jié)果圖和1個bmp測試圖像整體約1.59MB結(jié)構(gòu)緊湊便于逐文件調(diào)試。已有4295人學(xué)習(xí)下載說明其在教學(xué)與實驗場景中具備一定參考價值。讀者可據(jù)此對比高斯、伯努利等不同測量矩陣的重構(gòu)質(zhì)量調(diào)整松弛因子與閾值參數(shù)并借助MSE、PSNR指標和可視化結(jié)果驗證算法效果也可作為向視頻壓縮、醫(yī)學(xué)影像、遙感分析及深度學(xué)習(xí)結(jié)合方向拓展的實驗起點。1. 二維圖像壓縮感知重構(gòu)一份能直接跑通的 MATLAB 代碼包做圖像處理的人遲早會撞上壓縮感知這個坎。采樣率明明壓到了 30% 甚至更低重構(gòu)出來的圖卻要跟原圖肉眼難辨這中間的落差就是重構(gòu)算法要填的坑。這份資源給的不是論文里的偽代碼而是一套完整的 MATLAB 程序代碼專門針對二維圖像做壓縮感知重構(gòu)。它把測量矩陣生成、稀疏基構(gòu)造、重構(gòu)求解、PSNR 對比這幾步全串起來了拿到手改個圖片路徑就能出結(jié)果。適合正在做圖像處理大作業(yè)、想驗證 OMP 或 CoSaMP 這類算法實際效果的人也適合需要快速搭一個重構(gòu) baseline 再往上改的研究場景。MATLAB 版本不用追新R2016b 之后的都能跑核心函數(shù)沒用到太偏的工具箱。2. 壓縮感知重構(gòu)的數(shù)學(xué)骨架從測量矩陣到稀疏基的選型邏輯2.1 為什么二維圖像不能直接套一維壓縮感知一維壓縮感知的模型是 y Φxx 是長度 N 的列向量Φ 是 M×N 的測量矩陣M 遠小于 N。但二維圖像是 N×N 的矩陣如果直接按列拉伸成 N2×1 的向量測量矩陣的尺寸會膨脹到 M×N2內(nèi)存和計算量都吃不消。常見做法是分塊處理把圖像切成 32×32 或 64×64 的小塊每塊單獨做壓縮感知最后拼回去。另一種做法是用可分離的測量矩陣行方向和列方向各乘一個矩陣把二維運算拆成兩次一維運算。這份代碼走的是分塊路線塊大小在參數(shù)區(qū)可以改默認給的是 32×32兼顧了重構(gòu)質(zhì)量和速度。分塊之后有個細節(jié)容易被忽略塊與塊之間的邊界效應(yīng)。如果每塊獨立重構(gòu)再硬拼塊邊緣會出現(xiàn)明顯的方格感。代碼里用了一個簡單的重疊分塊策略塊間重疊 4 個像素重構(gòu)時對重疊區(qū)域做加權(quán)平均。這個處理不復(fù)雜但有沒有它視覺效果差很多。2.2 測量矩陣的選取高斯矩陣還是伯努利矩陣測量矩陣要滿足 RIP 條件實際中沒法直接驗證所以選的都是大概率滿足 RIP 的隨機矩陣。高斯矩陣是最常用的每個元素服從 N(0, 1/M)跟絕大多數(shù)稀疏基都不相關(guān)。伯努利矩陣元素取 ±1/√M硬件實現(xiàn)時更友好但重構(gòu)質(zhì)量在高壓縮比下會略遜于高斯矩陣。代碼里默認用的是高斯矩陣生成方式如下% 生成 M x N 的高斯測量矩陣M 為測量數(shù)N 為信號長度 function Phi generate_gaussian_matrix(M, N) % 每個元素服從 N(0, 1/M)保證列向量近似正交 Phi randn(M, N) / sqrt(M); % 對每一列做歸一化避免數(shù)值量級差異過大 Phi Phi ./ repmat(sqrt(sum(Phi.^2, 1)), M, 1); end這里做了列歸一化原因是后續(xù)重構(gòu)算法對矩陣列向量的能量敏感。如果不歸一化某些列的范數(shù)偏大會導(dǎo)致求解時權(quán)重失衡重構(gòu)出來的圖像局部過亮或過暗。歸一化之后每列能量為 1求解穩(wěn)定性明顯提升。參數(shù) M 就是測量數(shù)M/N 就是采樣率代碼里通過采樣率反算 M比如采樣率 0.3、塊大小 32×32 時N1024M307。2.3 稀疏基的選擇DCT 還是小波自然圖像在 DCT 基下是近似稀疏的在小波基下更稀疏。DCT 基構(gòu)造簡單就是余弦函數(shù)采樣計算快小波基需要濾波器組多一層分解但壓縮感知重構(gòu)質(zhì)量通常更好。代碼里兩種基都給了通過一個 flag 切換。DCT 基的構(gòu)造代碼% 構(gòu)造 N x N 的 DCT 正交基 function Psi generate_dct_basis(N) Psi zeros(N, N); for k 0:N-1 for n 0:N-1 if k 0 Psi(k1, n1) 1 / sqrt(N); else Psi(k1, n1) sqrt(2/N) * cos(pi * k * (2*n1) / (2*N)); end end end end這個雙重循環(huán)在 N1024 時會比較慢實際跑的時候可以用 dctmtx 函數(shù)替代MATLAB 自帶一行就能生成。代碼里保留手寫版本是為了讓讀者看清 DCT 基的每個元素怎么來的實際使用時建議換成 dctmtx。小波基用的是 sym4 小波通過 wavedec2 做兩層分解再把系數(shù)拉成列向量。注意小波分解后的系數(shù)長度跟原信號長度一致但排列順序會影響稀疏性代碼里做了按能量排序的處理把大系數(shù)排在前面這樣重構(gòu)時閾值收縮更有效。3. 重構(gòu)算法實戰(zhàn)OMP、CoSaMP 與 IRLS 的代碼實現(xiàn)與參數(shù)調(diào)優(yōu)3.1 OMP 算法的 MATLAB 實現(xiàn)與停止條件設(shè)置OMP 是正交匹配追蹤思路是每次迭代選一個跟殘差最相關(guān)的原子然后對已選原子做最小二乘更新。它的優(yōu)點是實現(xiàn)簡單、速度快缺點是需要知道稀疏度 KK 估不準就重構(gòu)不好。代碼里的 OMP 核心循環(huán)function x_rec omp_reconstruct(y, Phi, K) % y: 測量向量 M x 1 % Phi: 測量矩陣 M x N % K: 稀疏度 [M, N] size(Phi); residual y; support_set []; x_rec zeros(N, 1); for iter 1:K % 計算殘差與各原子的內(nèi)積找最大相關(guān)原子 correlations abs(Phi * residual); [~, idx] max(correlations); support_set [support_set, idx]; % 在已選原子集上做最小二乘 Phi_sub Phi(:, support_set); x_sub pinv(Phi_sub) * y; % 更新殘差 residual y - Phi_sub * x_sub; % 殘差能量低于閾值時提前退出 if norm(residual) 1e-6 break; end end x_rec(support_set) x_sub; end停止條件有兩個迭代次數(shù)達到 K或者殘差能量低于 1e-6。實際跑的時候如果 K 設(shè)得偏大OMP 會把噪聲也選進來重構(gòu)圖像出現(xiàn)顆粒感。K 設(shè)得偏小細節(jié)丟失。經(jīng)驗值是 K 取 N 的 10% 到 20%具體看圖像紋理復(fù)雜度。代碼里 K 是通過稀疏度估計函數(shù)自動算的用 DCT 系數(shù)排序后取累計能量 95% 的位置作為 K比固定值靠譜。3.2 CoSaMP 算法的回溯機制與參數(shù)影響CoSaMP 跟 OMP 的區(qū)別在于每次迭代選 2K 個原子然后回溯剪枝只保留 K 個。這個回溯機制讓它在稀疏度估計不準時更魯棒代價是每次迭代計算量翻倍。代碼里 CoSaMP 的關(guān)鍵步驟function x_rec cosamp_reconstruct(y, Phi, K) [M, N] size(Phi); residual y; support_set []; x_rec zeros(N, 1); max_iter 50; for iter 1:max_iter % 步驟1形成代理信號找 2K 個最大相關(guān)原子 proxy Phi * residual; [~, sort_idx] sort(abs(proxy), descend); new_support sort_idx(1:min(2*K, N)); % 步驟2合并支撐集 support_set union(support_set, new_support); % 步驟3在合并支撐集上做最小二乘 Phi_sub Phi(:, support_set); x_sub pinv(Phi_sub) * y; % 步驟4剪枝只保留 K 個最大系數(shù) [~, keep_idx] sort(abs(x_sub), descend); keep_idx keep_idx(1:min(K, length(keep_idx))); support_set support_set(keep_idx); x_sub x_sub(keep_idx); % 步驟5更新殘差 residual y - Phi(:, support_set) * x_sub; if norm(residual) 1e-6 break; end end x_rec(support_set) x_sub; endmax_iter 默認 50實際跑的時候大多數(shù)圖像在 20 次以內(nèi)就收斂了。如果跑到 50 次還沒收斂通常是測量矩陣條件數(shù)太差或者稀疏基選得不合適。回溯剪枝那一步是 CoSaMP 的核心它保證了支撐集大小始終不超過 K避免 OMP 那種選錯原子就一路錯到底的問題。代價是每次迭代要做兩次排序N 大時排序開銷不可忽略。3.3 IRLS 算法用迭代重加權(quán)最小二乘逼近 L0 范數(shù)IRLS 不走貪婪路線它用一系列加權(quán) L2 問題去逼近 L0 范數(shù)。每次迭代根據(jù)當前解的幅度調(diào)整權(quán)重幅度小的系數(shù)權(quán)重加大逼著它們往零走。function x_rec irls_reconstruct(y, Phi, max_iter, epsilon) [M, N] size(Phi); x_rec Phi * y; % 初始解用匹配濾波 weights ones(N, 1); for iter 1:max_iter % 構(gòu)造加權(quán)矩陣 W diag(weights); % 求解加權(quán)最小二乘min ||Wx||_2 s.t. Phi*x y % 用拉格朗日乘子法轉(zhuǎn)為無約束問題 x_rec W \ (Phi * ((Phi / W) * Phi \ y)); % 更新權(quán)重epsilon 防止除零 weights 1 ./ (abs(x_rec) epsilon); % 權(quán)重變化很小時退出 if norm(weights - 1./(abs(x_rec)epsilon)) 1e-4 break; end end endepsilon 是個關(guān)鍵參數(shù)默認給 1e-3。太大則稀疏性不夠太小則數(shù)值不穩(wěn)定。IRLS 的優(yōu)點是重構(gòu)精度通常比 OMP 高缺點是要解線性方程組每次迭代計算量比 OMP 大一個量級。圖像塊小的時候32×32還能接受塊大到 128×128 就明顯慢了。代碼里給 IRLS 設(shè)了最大迭代 30 次實際跑下來 15 次左右權(quán)重就基本穩(wěn)定了。4. 避坑與排查重構(gòu)圖像質(zhì)量不達標時先查這五處4.1 重構(gòu)圖像全黑或全白現(xiàn)象跑完代碼imshow 出來的圖要么全黑要么全白PSNR 低到個位數(shù)。原因測量矩陣沒有歸一化或者稀疏基的逆變換用錯了。全黑通常是重構(gòu)系數(shù)全為零全白是系數(shù)爆炸。檢查 generate_gaussian_matrix 里有沒有做列歸一化以及 DCT 基的轉(zhuǎn)置有沒有用對。解決在重構(gòu)完成后加一句 x_rec x_rec / max(abs(x_rec(:)))先歸一化再顯示。如果還是不對把測量矩陣換成 dctmtx 的前 M 行這個是確定性的能排除隨機矩陣的問題。4.2 PSNR 卡在 20dB 上不去現(xiàn)象采樣率提到 0.5 了PSNR 還是 20 出頭跟論文里的 35dB 差很遠。原因稀疏度 K 設(shè)得不對或者分塊大小太小。32×32 的塊在采樣率 0.5 時M512對紋理復(fù)雜的塊來說測量數(shù)還是不夠。另外 OMP 的 K 如果按固定值給遇到平滑塊會過擬合遇到紋理塊會欠擬合。解決把塊大小改成 64×64采樣率不變的情況下每塊測量數(shù)翻倍。K 改成自適應(yīng)用 DCT 系數(shù)累計能量 95% 來定。如果還不行換 CoSaMP它對 K 的估計誤差容忍度更高。4.3 重構(gòu)時間隨圖像尺寸指數(shù)增長現(xiàn)象512×512 的圖跑了十分鐘還沒出結(jié)果內(nèi)存占用飆升。原因沒有分塊直接把整圖拉成 262144×1 的向量測量矩陣尺寸 131072×262144內(nèi)存直接爆?；蛘叻謮K了但每塊都重新生成測量矩陣重復(fù)計算。解決確認代碼走的是分塊路徑塊大小不超過 64×64。測量矩陣生成一次就存下來所有塊共用。如果內(nèi)存還是緊張把測量矩陣改成稀疏存儲MATLAB 的 sparse 函數(shù)能省不少內(nèi)存。4.4 塊邊界出現(xiàn)方格狀偽影現(xiàn)象重構(gòu)圖像放大看塊與塊之間有明顯的亮度跳變像馬賽克。原因分塊時沒有重疊每塊獨立重構(gòu)后硬拼塊邊緣的系數(shù)不連續(xù)。解決啟用重疊分塊重疊像素數(shù)設(shè)為塊大小的 1/8 到 1/4。重構(gòu)后對重疊區(qū)域做加權(quán)平均權(quán)重用余弦窗中間高兩邊低。代碼里 overlap 參數(shù)默認是 4塊大小 32 時相當于 12.5% 重疊夠用了。4.5 換圖片后重構(gòu)完全失效現(xiàn)象用代碼自帶的測試圖跑得好好的換成自己的圖就重構(gòu)不出來。原因圖片尺寸不是塊大小的整數(shù)倍分塊時邊界處理出錯?；蛘邎D片是彩色三通道代碼只處理了灰度。解決讀圖后先做尺寸裁剪把長寬都裁成塊大小的整數(shù)倍。彩色圖先用 rgb2gray 轉(zhuǎn)灰度或者對三個通道分別做壓縮感知再合并。代碼里加一句 size_check 函數(shù)自動裁剪并給出提示。5. 從重構(gòu)到評估PSNR、SSIM 與采樣率曲線的自動化腳本跑通單張圖的重構(gòu)只是第一步真正做實驗的時候需要批量跑不同采樣率、不同算法然后畫對比曲線。這份代碼里帶了一個評估腳本把 PSNR 和 SSIM 的計算封裝好了直接調(diào)用就行。PSNR 的計算function psnr_val calc_psnr(original, reconstructed) % 確保數(shù)據(jù)類型一致 original double(original); reconstructed double(reconstructed); mse mean((original(:) - reconstructed(:)).^2); if mse 0 psnr_val Inf; return; end max_val max(original(:)); psnr_val 10 * log10(max_val^2 / mse); end注意 max_val 取的是原圖最大值不是固定的 255。如果圖像歸一化到 [0,1] 了max_val 就是 1PSNR 數(shù)值會跟 [0,255] 范圍差 20dB 左右對比論文數(shù)據(jù)時要統(tǒng)一量綱。SSIM 的計算用 MATLAB 自帶的 ssim 函數(shù)但要注意它默認是彩色圖模式灰度圖要加參數(shù)ssim_val ssim(reconstructed, original, DynamicRange, max(original(:)));批量實驗的腳本結(jié)構(gòu)是三層循環(huán)外層采樣率從 0.1 到 0.5 步進 0.1中層算法遍歷 OMP、CoSaMP、IRLS內(nèi)層遍歷測試圖片。每次跑完把 PSNR 和 SSIM 存到矩陣里最后用 plot 畫曲線。我一般會加一個計時器記錄每種算法的平均重構(gòu)時間。因為實際選算法的時候不只看質(zhì)量還要看能不能在可接受的時間內(nèi)跑完。OMP 在采樣率 0.3 時單塊重構(gòu)大概 0.05 秒CoSaMP 是 0.12 秒IRLS 要 0.8 秒。整圖 512×512 分 256 個塊OMP 總共 13 秒左右IRLS 要三分鐘以上。如果做實時性要求高的場景IRLS 基本不用考慮。還有一個技巧是固定隨機種子。測量矩陣是隨機生成的不同種子跑出來的 PSNR 會有 0.5dB 左右的波動。做對比實驗時把 rng(42) 寫在循環(huán)外面保證所有算法用的是同一個測量矩陣這樣對比才公平。這個細節(jié)論文里經(jīng)常不寫但實際復(fù)現(xiàn)時不知道就會覺得結(jié)果對不上。從那以后我每次跑壓縮感知實驗都先把隨機種子固定、量綱統(tǒng)一、塊大小和重疊參數(shù)確認一遍這三步走完再開始跑數(shù)據(jù)。希望幫到你。本文還有配套的精品資源點擊獲取