簡介Curvelet MATLAB工具箱是一套基于Curvelet變換的MATLAB實現(xiàn)庫面向圖像處理與信號分析領(lǐng)域的科研人員和工程師用于圖像去噪、壓縮、增強及邊緣特征提取等任務(wù)相比傳統(tǒng)小波變換更能捕捉圖像中的曲線與邊緣結(jié)構(gòu)適用于科研教學(xué)及工程實踐。該壓縮包共151個文件包含96個.m函數(shù)腳本、16個.hpp頭文件與16個.cpp源碼、9張示例圖像以及readme、makefile等配套文檔整體大小僅737KB便于快速部署和二次開發(fā)。目前已有578人學(xué)習(xí)下載是理解和應(yīng)用Curvelet變換的實用參考。工具箱基于CurveLab 1.0實現(xiàn)提供USFFT和Wrapping兩種算法的完整MATLAB接口用戶可直接調(diào)用函數(shù)完成正反變換也可結(jié)合C源碼理解底層原理并可通過options文件調(diào)整參數(shù)。包內(nèi)還附帶說明文檔和示例圖片適合開展閾值去噪、圖像壓縮重建等實驗源碼結(jié)構(gòu)清晰便于按需修改與擴展。 上個月在跑圖像去噪實驗時被手頭的 curvelet matlab 工具箱折騰得夠嗆——問題不是找不到代碼而是裝上之后第一條fdct_usfft就報錯查遍論壇才發(fā)現(xiàn)是 mex 編譯器和版本不對。后來把整個流程理順又花了兩個晚上。今天把整個使用體驗、踩坑記錄和核心細節(jié)都整理出來給準(zhǔn)備用 curvelet 做多尺度幾何分析的同行們一個可以照著走的參考。如果你正在用 Matlab 做圖像處理、地震信號分析、偏微分方程數(shù)值解或者想找一個比小波更擅長“捕捉曲線狀結(jié)構(gòu)”的工具這篇文章應(yīng)該能幫到你。curvelet 這個工具箱并不復(fù)雜但它和普通 matlab 函數(shù)包不一樣很多文件需要編譯而且參數(shù)設(shè)置直接決定結(jié)果好壞。我盡量用實際實驗的角度來講不堆理論。1. 為什么還要用 curvelet小波之后的幾何多尺度分析1.1 小波、脊波與 curvelet 的關(guān)系先說清楚 curvelet 是什么。小波變換大家熟它能做多尺度分析但在二維圖像里小波基是“點狀”的也就是說它在表示邊緣、曲線這類奇異特征時系數(shù)衰減較慢需要很多系數(shù)才能近似一條曲線。脊波Ridgelet通過沿直線方向積分能較好處理直線狀特征但真實圖像里直線太少曲線太多。curvelet 的思路是在脊波基礎(chǔ)上加了一個“尺度-方向”的參數(shù)化讓基函數(shù)像一個個“小條帶”能夠沿著曲線方向自適應(yīng)變形。通俗地講curvelet 可以把圖像分解為不同尺度、不同方向上的“線元素”就像用小波處理點用 curvelet 處理線。這個工具箱正是基于這種理論實現(xiàn)的。它不是簡單封裝一個函數(shù)而是包含多個 mex 編譯的底層 C/Fortran 程序以及大量 Matlab 腳本。核心功能就是把二維數(shù)組做 curvelet 正變換、逆變換并提取系數(shù)。1.2 適用場景與不適合的場景curvelet 在以下場景里表現(xiàn)非常突出圖像去噪、增強尤其是醫(yī)學(xué)圖像、遙感圖像中細膩的邊緣保護地震數(shù)據(jù)中反射波同相軸的分離和重建圖像壓縮、融合以及反卷積問題中的正則化項設(shè)計。但如果你的數(shù)據(jù)本身就是各向同性的紋理或者你只需要全局統(tǒng)計特征curvelet 的優(yōu)勢就不明顯了。更關(guān)鍵的一點是它的變換速度比小波慢一個數(shù)量級內(nèi)存開銷也大。如果只是做一個小小的實驗可能感覺不到但處理大量高分辨率遙感圖時一定要先測試數(shù)據(jù)量是否可行。2. 工具箱選型CurveLab 還是官方 Curvelet Toolbox2.1 兩個版本的區(qū)別Matlab 相關(guān)討論中常出現(xiàn)的 curvelet 工具箱主要有兩個來源。首先是 CurveLab這是 curvelet 變換發(fā)明人 Emmanuel Candès 等人發(fā)布的官方庫包含 USFFT 和 Wrapping 兩種實現(xiàn)也包含在 Matlab 、C 和 Fortran 中的調(diào)用接口。大多數(shù)論文、開源代碼里的fdct_usfft、fdct_wrapping就是來自這個庫。它的穩(wěn)定性好功能完整但需要自己編譯。另一個是 MathWorks File Exchange 上一些第三方封裝版本比如名為 Curvelet Toolbox 的工具箱。這類版本往往只包含 Wrapping 實現(xiàn)接口更簡單有的甚至不需要編譯直接調(diào)用。但問題在于這些版本多年未更新對 R2020 之后的 Matlab 兼容性很差且缺少很多底層參數(shù)的控制能力。我自己下載過兩個版本對比發(fā)現(xiàn) CurveLab 更適合做研究。雖然編譯麻煩點但函數(shù)命名統(tǒng)一、文檔完整、測試用例多出了問題還能查源碼。第三方版本雖然“開箱即用”但在處理任意尺寸圖像時經(jīng)常報維度錯誤而且因為代碼結(jié)構(gòu)被簡化很難定位問題。2.2 我最終選型的原因最終我選擇 CurveLab版本 2.1.1。三個理由它是原作者的實現(xiàn)數(shù)學(xué)定義和論文一致做實驗時與理論對照不怕有偏差同時提供fdct_usfft和fdct_wrapping可以在同一套環(huán)境下對比兩種快速算法的差異它內(nèi)置了去噪、插值、壓縮的示例腳本對快速上手非常有幫助。所以下面的所有操作都以 CurveLab 2.1.1 為基礎(chǔ)。如果你手頭是其他版本函數(shù)名可能有所不同但核心流程類似。3. 安裝與配置Windows/Linux 下的實操記錄3.1 下載與目錄結(jié)構(gòu)打開 CurveLab 的官網(wǎng)下載對應(yīng)的壓縮包解壓后目錄結(jié)構(gòu)大致是fdct_usfft_matlab/ fdct_wrapping_matlab/ fdct3d_matlab/ mecv/ ...其中fdct_usfft_matlab和fdct_wrapping_matlab是我們最常用的兩個文件夾。每個文件夾內(nèi)都有MakefileUnix/Linux 下用mex后綴的 C 源文件一些.m文件主函數(shù)和輔助函數(shù)test或demo腳本建議不要直接解壓到 Matlab 的安裝目錄下而是放到自己的工作區(qū)比如D:\tools\curvelet\或者/home/user/tools/curvelet方便后續(xù)維護。3.2 編譯命令解析這個工具箱的 C 源碼需要用 mex 編譯成.mexw64Windows或.mexa64Linux文件。編譯過程是安裝中坑最多的一步。在 Windows 下直接進入fdct_wrapping_matlab目錄在 Matlab 命令行運行cd(D:\tools\curvelet\fdct_wrapping_matlab) make但請確保你已經(jīng)通過mex -setup配置了 C 編譯器。在 R2020a 之后Matlab 默認推薦 MinGW-w64。如果電腦上裝的是 Visual Studio也可以選它但注意版本要匹配。我實測下來R2021b 用 MinGW 10.3 編譯 CurveLab 2.1.1 是沒有問題的R2022b 也可以。如果運行make時報找不到mex命令說明當(dāng)前目錄不對或者沒進入 Matlab 的 Command Window。建議直接在 Matlab 中cd到目標(biāo)目錄再運行。在 Linux 下進入目錄后直接運行make即可前提是你已經(jīng)安裝了 gcc 和 Matlab 的 mex 支持。如果報gcc: error: unrecognized command-line option -mno-cygwin這是老版本源代碼里針對 Cygwin 的選項在 Linux 下需要手動刪掉 Makefile 里的這一項。這個坑在較新版 gcc 上很常見。編譯成功后目錄下會出現(xiàn)一系列.mexw64或.mexa64文件同時提示不需要額外配置。為了驗證我建議立刻運行fdct_wrapping(rand(64,64))如果能返回一個系數(shù)結(jié)構(gòu)體說明安裝成功。如果只是輸出錯誤提示那多半是矩陣尺寸不匹配或編譯未完成。3.3 環(huán)境變量設(shè)置雖然編譯成功后可以直接使用但為了其他腳本也能調(diào)用需要在 Matlab 中把工具箱路徑加入搜索路徑。在 Matlab 命令行運行addpath(genpath(D:\tools\curvelet\fdct_usfft_matlab)) addpath(genpath(D:\tools\curvelet\fdct_wrapping_matlab)) savepath注意savepath會保存到當(dāng)前用戶的 pathdef.m如果因為權(quán)限問題保存失敗可以在MATLABPATH環(huán)境變量里加入這兩個目錄。不過我更推薦直接建一個startup.m把 addpath 寫進去這樣每次啟動 Matlab 時都會自動加載不用每天重新 set path。4. 核心函數(shù)使用與參數(shù)選擇fdct_usfft 和 fdct_wrapping 的差異4.1 兩個核心變換函數(shù)CurveLab 提供了兩組接口fdct_usfft基于非均勻快速傅里葉變換USFFT實現(xiàn)比較接近原始論文里的頻域采樣方式fdct_wrapping基于“wrapping”技巧將頻域中的楔形區(qū)域像素“包裹”回矩形用普通 FFT 實現(xiàn)速度更快。兩者的數(shù)學(xué)結(jié)果幾乎一致但在邊界效應(yīng)、參數(shù)選擇和計算復(fù)雜性上略有差別。fdct_usfft對尺寸要求更嚴(yán)通常要求圖像尺寸是某些數(shù)的倍數(shù)fdct_wrapping更靈活幾乎適用于任意大小。4.2 參數(shù)詳解基本調(diào)用方式是C fdct_wrapping(X, is_real, finest, nbscales, nbangles_coarse);其中X是輸入二維矩陣is_real為1表示處理實數(shù)圖像通常選 10表示復(fù)數(shù)域finest控制是否在最細尺度上也進行分解取1或0nbscales是分解尺度數(shù)一般取ceil(log2(min(M,N)))附近nbangles_coarse是最粗尺度上的方向數(shù)通常設(shè)定為16或32后續(xù)尺度方向數(shù)會按倍數(shù)增加。這些參數(shù)的物理意義可以這樣理解nbscales控制把圖像分解成幾個層次類似于小波分解的層數(shù)。nbangles_coarse控制最粗尺度上劃分方向的數(shù)量數(shù)量越多對曲線方向的捕捉越細致但計算量也越大。默認nbangles_coarse 16對大多數(shù)圖像已經(jīng)足夠。如果你處理的是紋理極其豐富的圖像比如遙感中的城市區(qū)域建議增加到 32但同時要注意內(nèi)存消耗會顯著上升。4.3 一個完整的重構(gòu)演示安裝完成后可以先做一個最簡單的重構(gòu)測試驗證工具鏈?zhǔn)欠裢暾? 生成 256x256 的測試圖像 X zeros(256,256); X(64:192,64:192) 1; X(80:120,80:180) 0; % 正向 curvelet 變換 C fdct_wrapping(X, 1, 0, 5, 16); % 逆變換 Y ifdct_wrapping(C, 1, 0, 5, 16); % 檢查誤差 disp(max(abs(X(:) - Y(:))));如果一切正常誤差應(yīng)該是0或者是1e-10級別的浮點誤差。如果誤差大大概率是參數(shù)沒對齊后面會詳細講。注意這里的逆變換參數(shù)必須和正變換完全一致is_real、finest、nbscales、nbangles_coarse缺一不可。我在最初測試時只改了尺度數(shù)忘了改方向數(shù)導(dǎo)致重建出來的圖像被嚴(yán)重扭曲那種情況還很難察覺是參數(shù)不匹配導(dǎo)致的。5. 常見問題與排查技巧實錄5.1 編譯報錯與 mex 版本這是我遇到最多的一類問題?？偨Y(jié)起來主要有幾種mex 未找到支持的編譯器這是 R2019b 之后最常見的。解決方式是安裝 MinGW-w64然后在 Matlab 中執(zhí)行mex -setup C選中 MinGW。錯誤使用 mex未找到匹配的函數(shù)大概率是.m文件里有同名函數(shù)干擾。檢查當(dāng)前文件夾是否存在其他fdct_wrapping.m如果有把路徑放到最前面。error LNK2019: unresolved external symbol這是 Windows 下老版本 CurveLab 的鏈接問題通常是因為 Matlab 版本太新而 C 源碼里用的舊接口已經(jīng)不兼容?？梢試L試在源碼中修改 C 源文件將mxCreateDoubleMatrix等函數(shù)調(diào)用全部加上mx前綴其實本來就是但更保險的是安裝一個 2017 之后的版本并且將編譯器換成較新的 MinGW。我個人遇到最多的是“unrecognized command-line option”主要是舊 Makefile 與新版編譯器不兼容。Linux 下刪掉-mno-cygwin和-fopenmp如果你不需要并行計算就能解決。5.2 坐標(biāo)系與索引的坑curvelet 變換后的系數(shù)按尺度、方向、位置儲存在一個大 cell 數(shù)組里。很多人第一次使用時不知道方向數(shù)的編號方式正不正確。我用一個經(jīng)驗規(guī)則fdct_wrapping的輸出C{n}{m}中n是尺度編號n1是哪個尺度取決于finest參數(shù)。如果finest為 0第一層是最粗尺度最后一層是最細細節(jié)如果finest為 1最后一層既包含細節(jié)也包含粗尺度信息。方向編號m的范圍取決于該尺度上的方向數(shù)不同尺度方向數(shù)不同所以不能用一個固定的m值。如果操作時越界會報Index exceeds matrix dimensions。我的建議是在動手處理系數(shù)前先用size(C{n}{m})打印所有尺度的方向數(shù)量并畫一張系數(shù)分布圖確認你的索引習(xí)慣和輸出一致。這個工作只要做一次以后就不會再被繞暈。5.3 計算速度與內(nèi)存問題一個 512x512 的圖像nbscales6、nbangles_coarse32正變換大約耗時 0.5~1 秒逆變換類似。如果換成 4096x4096耗時就會飆到幾十秒而且內(nèi)存很容易爆掉尤其是同時保存正變換和逆變換結(jié)果時。如果是大批量處理有幾個實測有效的辦法把nbangles_coarse降到 16通常能減少 30% 左右的計算量圖像尺寸盡量是 2 的冪次既方便尺度劃分也能提升 FFT 效率如果只做去噪不用保存所有尺度系數(shù)可以直接在變換域處理然后重建避免多次復(fù)制 cell 結(jié)構(gòu)。另外在 Matlab R2021b 之后fdct_wrapping底層還是默認單線程多核優(yōu)化不明顯。如果數(shù)據(jù)量特別大建議考慮用 C 版本或 Python 的包裝器但那就超出了 Matlab 工具箱的范疇了。6. 實際應(yīng)用案例圖像去噪與增強6.1 基于 curvelet 的去噪流程curvelet 最經(jīng)典的應(yīng)用是去噪。原理很簡單圖像中噪聲是高頻、各向同性的而邊緣和紋理通常是曲線狀、方向性強的特征。curvelet 變換后邊緣對應(yīng)的系數(shù)幅值大而集中噪聲系數(shù)幅值小而分散通過軟閾值或硬閾值處理就能在保護邊緣的同時去掉噪聲。具體腳本可以參考% 讀入圖像并添加高斯噪聲 X im2double(imread(cameraman.tif)); Xnoisy X 0.05*randn(size(X)); % 參數(shù)設(shè)置 nbscales 5; nbangles_coarse 16; % 正變換 C fdct_wrapping(Xnoisy, 1, 0, nbscales, nbangles_coarse); % 對每一尺度、每一方向做閾值收縮 Cthresh C; sigma 0.05; % 噪聲標(biāo)準(zhǔn)差 for s 1:nbscales for w 1:length(C{s}) % 每個系數(shù)的閾值尺度越細閾值越大 thr 3*sigma*sqrt(2*(s1)); Cthresh{s}{w} wthresh(C{s}{w}, s, thr); end end % 逆變換重建 Y ifdct_wrapping(Cthresh, 1, 0, nbscales, nbangles_coarse);這里的thr公式是一個經(jīng)驗值并不是數(shù)學(xué)上的最優(yōu)值。實際使用中如果去噪太強會導(dǎo)致圖像過于平滑可以降低倍數(shù)系數(shù)。我通常會在2.5到3.5之間調(diào)整。6.2 參數(shù)調(diào)優(yōu)經(jīng)驗在去噪實驗中nbscales決定了分解層次。如果噪聲特別嚴(yán)重比如噪聲標(biāo)準(zhǔn)差超過 0.1建議用 6 個尺度讓高頻細節(jié)與噪聲更好地分離如果噪聲較輕5 個尺度就夠太多反而容易把微小細節(jié)也當(dāng)作噪聲濾掉。nbangles_coarse則直接影響對曲線邊緣的保護效果。我用實驗做過對比從 8 增加到 16PSNR峰值信噪比大約能提升 0.3~0.5 dB從 16 增加到 32PSNR 提升很小但邊緣保留的主觀視覺質(zhì)量稍微好一些。從計算效率考慮16 是最平衡的選擇。還有一個很關(guān)鍵的細節(jié)去噪前最好對圖像做邊緣延拓避免傅里葉變換的周期邊界假象。CurveLab 本身不自動處理邊界所以如果你處理的圖像尺寸不是偶數(shù)或者有強烈邊緣建議先對圖像做鏡像延拓變換后再裁掉。我都是用padarray(X, [32 32], symmetric)做延拓去噪完成后用crop裁回原尺寸效果會好很多。最后再分享一個小技巧。在mex編譯時如果你用了較新的 Matlab建議在編譯命令中加上-largeArrayDims標(biāo)志這樣生成的 mex 文件能處理大于 2GB 的內(nèi)存訪問在處理高分辨率圖像時更安全。具體操作是在 Makefile 的mex命令里加上這個參數(shù)或者在 Matlab 中手動調(diào)用mex(-largeArrayDims, fdct_wrapping_mex.c);這個方法是從一個老外的工程博客里看到的實測確實能在處理 4096×4096 的遙感圖像時避免部分“內(nèi)存不足”的報錯。如果實驗結(jié)果經(jīng)常不理想優(yōu)先檢查參數(shù)匹配和邊界處理不要一上來就懷疑算法本身。curvelet 工具箱雖然老了點但在邊緣保護這個方向上現(xiàn)在的很多深度學(xué)習(xí)模型依然會拿它做對比基準(zhǔn)它依然是值得一用的經(jīng)典工具。本文還有配套的精品資源點擊獲取