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