簡(jiǎn)介本資源是一份面向MATLAB圖像處理初學(xué)者與進(jìn)階學(xué)習(xí)者的帶通濾波實(shí)踐工具包聚焦于頻域圖像增強(qiáng)與噪聲抑制等典型任務(wù)適用于醫(yī)學(xué)影像分析、遙感圖像特征提取及視頻預(yù)處理等實(shí)際場(chǎng)景。壓縮包為RAR格式共含2個(gè)MATLAB腳本文件.m總大小僅1022B輕量簡(jiǎn)潔其中核心函數(shù)daitongfilter.m封裝了可調(diào)參數(shù)的帶通濾波器設(shè)計(jì)邏輯支持巴特沃斯等經(jīng)典濾波器類型主腳本daitong.m則完整實(shí)現(xiàn)圖像讀取、頻域轉(zhuǎn)換、濾波應(yīng)用、結(jié)果可視化與保存全流程便于快速驗(yàn)證與二次開(kāi)發(fā)。目前已有400人學(xué)習(xí)下載資源雖小但結(jié)構(gòu)完整提供即用型濾波接口、清晰的參數(shù)注釋及端到端圖像處理范例特別適合理解頻域?yàn)V波原理、調(diào)試濾波器截止頻率設(shè)置及開(kāi)展課程實(shí)驗(yàn)與項(xiàng)目原型開(kāi)發(fā)。 做圖像處理的人遲早都會(huì)撞上帶通濾波這堵墻。我最早接觸它是在做紋理分析的時(shí)候想從一張布滿噪點(diǎn)的表面缺陷圖里把周期性紋理提出來(lái)低通把紋理和高頻噪聲一起抹了高通又把缺陷邊緣和噪聲一起放大了折騰一圈才發(fā)現(xiàn)真正需要的其實(shí)是把某個(gè)頻段單獨(dú)拎出來(lái)——這就是帶通濾波的活。當(dāng)時(shí)從網(wǎng)上下過(guò)一個(gè)叫“daitong.rar”的MATLAB程序包里面有個(gè)圖像帶通的Demo代碼寫(xiě)得亂跑起來(lái)還一堆坑但這東西確實(shí)讓我把頻域?yàn)V波那點(diǎn)事徹底搞明白了。這篇東西就是寫(xiě)給正在被帶通濾波折磨的人。不管你是做圖像處理大作業(yè)還是在做工業(yè)檢測(cè)、醫(yī)學(xué)影像分析只要能理解“頻段”這個(gè)概念就能用MATLAB把帶通濾波玩明白。我會(huì)從原理、代碼、參數(shù)調(diào)優(yōu)到實(shí)際踩坑完整走一遍最后給你一套可以直接抄作業(yè)的MATLAB實(shí)現(xiàn)。1. 圖像帶通濾波到底在解決什么問(wèn)題先說(shuō)個(gè)很多初學(xué)者沒(méi)想明白的事圖像濾波的本質(zhì)不是“模糊”或“銳化”而是在頻率域里做選擇。圖像可以看成無(wú)數(shù)個(gè)不同頻率的正弦波疊加的結(jié)果低頻對(duì)應(yīng)平坦區(qū)域和大致輪廓高頻對(duì)應(yīng)邊緣、細(xì)節(jié)和噪聲。低通濾波就是留低頻去高頻高通濾波就是留高頻去低頻那帶通呢它把兩者都砍了一部分只留下中間那一段頻率。1.1 為什么低通和高通常常不夠用舉一個(gè)真實(shí)的場(chǎng)景。我在做芯片引腳缺陷檢測(cè)的時(shí)候圖像里有三類東西大面積的底色變化低頻、引腳邊緣和細(xì)微裂紋中頻、傳感器本身的隨機(jī)噪聲高頻。如果用低通濾波噪聲是沒(méi)了但裂紋和邊緣也被一起模糊掉如果用高通濾波邊緣和裂紋倒是出來(lái)了但噪聲也被放大了簡(jiǎn)直是災(zāi)難。這種兩頭為難的情況在醫(yī)學(xué)影像里更常見(jiàn)。比如說(shuō)眼底圖像里的血管分割血管寬度對(duì)應(yīng)的頻率范圍是相對(duì)固定的背景光照不均屬于低頻成像系統(tǒng)的顆粒噪聲屬于高頻。你要把血管結(jié)構(gòu)完整提取出來(lái)就必須要帶通濾波把背景不均勻性去掉把顆粒噪聲壓掉只保留血管對(duì)應(yīng)的頻段。1.2 帶通濾波在圖像處理里的典型應(yīng)用我梳理了一下實(shí)際項(xiàng)目里帶通濾波最常見(jiàn)的幾個(gè)用途紋理提取與分析織物、木材、金屬表面的紋理往往集中在特定的頻率范圍帶通濾波可以把紋理從背景中分離出來(lái)方便后續(xù)的統(tǒng)計(jì)分析和缺陷識(shí)別。圖像增強(qiáng)遙感圖像里的山脊、道路等線性特征醫(yī)學(xué)圖像里的血管、支氣管樹(shù)都屬于特定尺度的結(jié)構(gòu)帶通濾波可以定向增強(qiáng)這些結(jié)構(gòu)抑制無(wú)關(guān)信息。去噪與去偽影某些周期性噪聲比如掃描儀的條紋干擾在頻域里表現(xiàn)為一個(gè)特定頻率的峰值用帶通濾波器把那個(gè)頻率單獨(dú)濾出來(lái)再在原始圖像里減去就能干凈地去除。工業(yè)缺陷檢測(cè)表面劃痕、凹坑、氣泡這類缺陷往往出現(xiàn)在中頻段帶通濾波配合形態(tài)學(xué)處理可以顯著降低誤檢率。1.3 “頻率”在圖像里到底該怎么理解這是理解帶通濾波的關(guān)鍵一步。對(duì)圖像來(lái)說(shuō)頻率不是時(shí)間上的快慢而是空間上灰度變化的劇烈程度。我習(xí)慣用一個(gè)類比來(lái)解釋你在看一面平靜的湖水低頻圖像突然落下一顆石子漣漪的產(chǎn)生和擴(kuò)散中頻結(jié)構(gòu)水面上密密麻麻的細(xì)小波紋高頻噪聲。如果你想研究漣漪本身就需要把湖面的平靜背景和細(xì)碎波紋都去掉——這就是帶通濾波在做的事。對(duì)應(yīng)到圖像里平坦區(qū)域灰度變化緩慢是低頻邊緣處灰度突變是高頻而像紋理這種有規(guī)律的中等尺度變化就落在中頻。帶通濾波器就是給圖像裝了一個(gè)“頻率窗口”只讓窗口內(nèi)的頻率成分通過(guò)。2. 頻域設(shè)計(jì)路線從傅里葉變換到帶通濾波器圖像帶通濾波最直觀、最可控的實(shí)現(xiàn)方式就是走頻域路線先把圖像從空間域變換到頻率域在頻率域里構(gòu)造一個(gè)帶通濾波器把頻譜和濾波器逐點(diǎn)相乘再變換回空間域。2.1 頻域?yàn)V波的完整流程在MATLAB里標(biāo)準(zhǔn)的頻域?yàn)V波流程是這樣% 讀入圖像并轉(zhuǎn)為灰度圖 img imread(texture.png); if size(img, 3) 3 img rgb2gray(img); end img im2double(img); % 轉(zhuǎn)浮點(diǎn)類型避免計(jì)算溢出 % 1. 傅里葉變換并中心化 F fft2(img); % 二維傅里葉變換 F_shift fftshift(F); % 將零頻移到頻譜中心 % 2. 構(gòu)造帶通濾波器傳遞函數(shù) H(u,v) % 這里先用一個(gè)簡(jiǎn)單的理想帶通示意后面會(huì)詳細(xì)講 [M, N] size(img); [u, v] meshgrid(1:N, 1:M); D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); % 每個(gè)像素到中心的距離 D0_low 30; % 通帶下截止頻率 D0_high 80; % 通帶上截止頻率 H ones(M, N); H(D D0_low | D D0_high) 0; % 理想帶通 % 3. 頻譜與濾波器相乘 G_shift F_shift .* H; % 4. 反傅里葉變換回空間域 G ifftshift(G_shift); filtered real(ifft2(G)); % 5. 顯示結(jié)果 figure; subplot(1,2,1); imshow(img); title(原圖); subplot(1,2,2); imshow(filtered, []); title(帶通濾波結(jié)果);這段代碼里最核心的一步是fftshift。fft2變換后的頻譜零頻在矩陣的左上角也就是四個(gè)角的頻率最低中心頻率最高。如果不做fftshift你構(gòu)造的濾波器就必須跟著去匹配這種“零頻在角落”的布局非常容易出錯(cuò)。fftshift之后零頻在正中心頻率從中心向四周逐漸增高這種“從中心讀頻率”的方式最符合直覺(jué)。2.2 三種常見(jiàn)的帶通濾波器傳遞函數(shù)MATLAB圖像處理工具箱里其實(shí)沒(méi)有直接叫“帶通”的濾波器函數(shù)但帶通可以由低通和高通組合而來(lái)或者直接構(gòu)造傳遞函數(shù)。實(shí)際項(xiàng)目中最常用的有三種理想帶通濾波器Ideal Bandpass Filterfunction H ideal_bandpass(M, N, D0, D1) % D0: 下截止頻率通帶內(nèi)徑 % D1: 上截止頻率通帶外徑 [u, v] meshgrid(1:N, 1:M); D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); H ones(M, N); H(D D0 | D D1) 0; end理想帶通的特性是通帶內(nèi)所有頻率原樣通過(guò)通帶外完全截止。邊界非常陡峭理論上效果最好但實(shí)際用起來(lái)會(huì)有問(wèn)題——后面會(huì)詳細(xì)說(shuō)。高斯帶通濾波器Gaussian Bandpass Filterfunction H gaussian_bandpass(M, N, D0, sigma) % D0: 中心頻率通帶中心 % sigma: 帶寬控制參數(shù) [u, v] meshgrid(1:N, 1:M); D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); H 1 - exp(-(D.^2 - D0^2).^2 / (2 * (D0 * sigma)^2)); end這個(gè)形式看著有點(diǎn)繞其實(shí)本質(zhì)就是在頻率軸上以 D0 為中心形成一個(gè)高斯形狀的突起。離 D0 越遠(yuǎn)的頻率衰減越厲害。高斯帶通沒(méi)有突變邊界濾波結(jié)果非常平滑是最推薦新手使用的。巴特沃斯帶通濾波器Butterworth Bandpass Filterfunction H butterworth_bandpass(M, N, D0, D1, n) % D0: 下截止頻率 % D1: 上截止頻率 % n: 濾波器階數(shù)控制過(guò)渡帶的陡峭程度 [u, v] meshgrid(1:N, 1:M); D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); H 1 ./ (1 ((D .^ 2 - D0 * D1) ./ (D .* (D1 - D0))) .^ (2 * n)); end巴特沃斯帶通是理想帶通的“平滑版”通過(guò)調(diào)節(jié)階數(shù) n 可以控制過(guò)渡帶的陡峭程度。n 越大越接近理想帶通但振鈴效應(yīng)也會(huì)越明顯。三種濾波器的對(duì)比如下濾波器類型過(guò)渡帶特性振鈴效應(yīng)參數(shù)含義適用場(chǎng)景理想帶通無(wú)過(guò)渡帶階躍截止嚴(yán)重D0下截止、D1上截止極少單獨(dú)使用高斯帶通平滑過(guò)渡無(wú)D0中心頻率、sigma帶寬初學(xué)者首選紋理分析巴特沃斯帶通可通過(guò)n調(diào)節(jié)中等D0下截止、D1上截止、n階數(shù)工業(yè)檢測(cè)、邊緣分析2.3 為什么理想帶通會(huì)有振鈴效應(yīng)這是頻域?yàn)V波里最容易踩的坑。理想帶通濾波器的傳遞函數(shù)在截止頻率處是階躍的從 0 直接跳到 1。這種突變?cè)陬l率域里很尖銳對(duì)應(yīng)到空間域里就會(huì)產(chǎn)生振蕩——就是圖像上出現(xiàn)一圈一圈的偽影專業(yè)上叫振鈴效應(yīng)Ring Effect實(shí)際表現(xiàn)就是目標(biāo)邊緣附近出現(xiàn)明暗交替的“鬼影”。我最早用理想帶通做血管增強(qiáng)出來(lái)結(jié)果嚇一跳血管周圍全是水波紋。后來(lái)才明白凡是傳遞函數(shù)有“硬跳變”的濾波器必然存在振鈴。高斯帶通和低階巴特沃斯帶通因?yàn)檫^(guò)渡帶平滑就不會(huì)有這個(gè)問(wèn)題。所以我的建議是能用高斯帶通就別用理想帶通除非你明確知道振鈴偽影不影響后續(xù)處理。3. 空間域?qū)崿F(xiàn)卷積核設(shè)計(jì)與實(shí)操頻域?yàn)V波理解起來(lái)直觀但實(shí)際工程里有時(shí)候需要在空間域做卷積原因無(wú)非兩種要么是頻域變換的開(kāi)銷太大要么是濾波器需要集成到實(shí)時(shí)處理流程里無(wú)法忍受傅里葉變換的延遲?？臻g域帶通濾波的本質(zhì)是設(shè)計(jì)一個(gè)空間卷積核讓它的頻率響應(yīng)呈現(xiàn)“帶通”形狀。3.1 從頻率響應(yīng)反推空間卷積核空間域卷積核的設(shè)計(jì)思路是一個(gè)卷積核的頻率響應(yīng)就是它自身做傅里葉變換的結(jié)果。反過(guò)來(lái)我們可以先設(shè)計(jì)好想要的頻率響應(yīng)再做逆傅里葉變換得到空間卷積核。MATLAB里最直接的做法是fsamp2函數(shù)它根據(jù)頻率響應(yīng)采樣值生成二維FIR濾波器% 設(shè)計(jì)一個(gè)5x5的帶通卷積核 % 先構(gòu)造理想的頻率響應(yīng)矩陣尺寸隨意但要能覆蓋你關(guān)心的頻段 [M, N] size(img); % 用之前構(gòu)造的高斯帶通頻率響應(yīng)作為目標(biāo) H_desired gaussian_bandpass(M, N, 50, 0.4); % 生成空間域卷積核 h fsamp2(H_desired); % 用卷積核做空間濾波 filtered imfilter(img, h, replicate, same);這里fsamp2會(huì)生成一個(gè)和頻率響應(yīng)矩陣等大的卷積核直接用imfilter做卷積就行。注意imfilter的邊界選項(xiàng)replicate表示邊界外復(fù)制邊緣像素same表示輸出和輸入同尺寸。邊界處理不當(dāng)會(huì)造成邊緣黑條后面會(huì)有專門(mén)的避坑說(shuō)明。3.2 用高斯拉普拉斯組合構(gòu)造帶通空間域還有一種更經(jīng)典的帶通構(gòu)造方法用兩個(gè)不同尺度的高斯低通核做差。這本質(zhì)上就是高斯差分DoG它的頻率響應(yīng)天然是帶通形狀。% 高斯差分構(gòu)造帶通核 sigma1 1.0; % 小尺度高斯對(duì)應(yīng)高頻 sigma2 3.0; % 大尺度高斯對(duì)應(yīng)低頻 h1 fspecial(gaussian, [15 15], sigma1); h2 fspecial(gaussian, [15 15], sigma2); h_dog h1 - h2; % 兩個(gè)高斯核相減得到帶通核 % 空間域?yàn)V波 filtered imfilter(img, h_dog, replicate, same);這個(gè)方法的物理意義非常直觀大尺度高斯平滑掉高頻細(xì)節(jié)小尺度高斯保留更多細(xì)節(jié)兩者相減留下來(lái)的就是介于兩種尺度之間的結(jié)構(gòu)——也就是帶通。DoG在生物學(xué)上還模擬了視網(wǎng)膜神經(jīng)節(jié)的感受野響應(yīng)在邊緣檢測(cè)里也是經(jīng)典算法。3.3 頻域和空間域兩種路線的對(duì)比很多初學(xué)者糾結(jié)到底用哪種實(shí)現(xiàn)我的看法是看場(chǎng)景維度頻域?qū)崿F(xiàn)空間域?qū)崿F(xiàn)實(shí)現(xiàn)難度低代碼直觀中等需要設(shè)計(jì)核計(jì)算效率大圖慢依賴FFT核小則快適合實(shí)時(shí)靈活性高頻率響應(yīng)任意設(shè)計(jì)受核尺寸和采樣限制邊界效應(yīng)隱含著周期性假設(shè)可通過(guò)邊界選項(xiàng)控制適合場(chǎng)景離線分析、參數(shù)調(diào)優(yōu)實(shí)時(shí)系統(tǒng)、嵌入式部署如果你是在做研究或者一次性的數(shù)據(jù)處理頻域?qū)崿F(xiàn)更直觀調(diào)參也方便。如果是要嵌入到實(shí)時(shí)檢測(cè)系統(tǒng)里空間域卷積核一旦訓(xùn)練好就用不著每次做傅里葉變換了速度和穩(wěn)定性都更好。4. 帶通濾波器參數(shù)選型的實(shí)操策略帶通濾波器的效果90%取決于參數(shù)選得好不好。參數(shù)選得不對(duì)比如通帶范圍設(shè)得太寬濾波結(jié)果跟原圖區(qū)別不大設(shè)得太窄又會(huì)把有用的信息濾掉。這一節(jié)我把自己調(diào)參的完整思路和工具分享出來(lái)。4.1 通帶范圍怎么定從頻譜圖反推參數(shù)調(diào)參的第一步不是拍腦袋設(shè)個(gè)D0和D1而是先看頻譜圖搞清楚你要的信息到底在哪個(gè)頻段。我常用的方法是先畫(huà)出圖像的頻譜用視覺(jué)判斷信息分布F fftshift(fft2(img)); % 取對(duì)數(shù)壓縮動(dòng)態(tài)范圍方便觀察 F_log log(1 abs(F)); imshow(F_log, []);觀察頻譜圖時(shí)中心區(qū)域亮表示低頻能量強(qiáng)四周的亮線或亮點(diǎn)對(duì)應(yīng)著圖像里特定方向的周期紋理或噪聲。你要提取的紋理或結(jié)構(gòu)如果在頻譜上表現(xiàn)為某個(gè)半徑附近的環(huán)形亮斑那通帶的中心頻率就選那個(gè)半徑值通帶寬度就覆蓋亮斑的范圍。如果覺(jué)得看圖換算麻煩我還有一個(gè)更直接的定量方法把頻譜按半徑從中心向外累加能量畫(huà)出能量隨頻率的分布曲線從曲線上找峰值區(qū)域。% 圍繞頻譜中心計(jì)算每個(gè)半徑的能量 [M, N] size(F_shift); [u, v] meshgrid(1:N, 1:M); D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); radius round(D(:)); energy abs(F_shift(:)).^2; % 按半徑分組求能量和 maxR max(radius); radial_energy zeros(maxR 1, 1); for r 0 : maxR radial_energy(r 1) sum(energy(radius r)); end plot(0 : maxR, radial_energy);看完這條曲線你大概就知道圖像的能量集中在哪些頻段。做帶通濾波時(shí)通帶就壓在目標(biāo)頻段的峰值附近。4.2 高斯帶通的sigma怎么選高斯帶通濾波器里有兩個(gè)核心參數(shù)中心頻率D0和帶寬控制參數(shù)sigma。D0好理解就是你要提取的主體頻率sigma的作用是控制通帶的寬度。我的經(jīng)驗(yàn)是sigma取D0的0.3到0.6倍比較合理。sigma太小通帶窄只提取非常單一的頻率結(jié)果會(huì)像被“抽取”了紋理一樣細(xì)節(jié)丟失嚴(yán)重sigma太大通帶寬就會(huì)把鄰近頻段的無(wú)關(guān)信息也放進(jìn)來(lái)失去了帶通的意義。一個(gè)比較實(shí)用的做法是做一個(gè)參數(shù)掃描img im2double(imread(texture.png)); [M, N] size(img); F_shift fftshift(fft2(img)); sigmas 0.2 : 0.1 : 0.8; figure; for i 1 : length(sigmas) H gaussian_bandpass(M, N, 50, sigmas(i)); G F_shift .* H; filtered real(ifft2(ifftshift(G))); subplot(2, 4, i); imshow(filtered, []); title([sigma num2str(sigmas(i))]); end把不同sigma的結(jié)果并排顯示你一眼就能判斷哪個(gè)參數(shù)最貼合你的需求。這種“暴力掃描”法雖然不優(yōu)雅但在實(shí)際項(xiàng)目中往往是最有效的調(diào)參手段。4.3 判定濾波效果好壞的客觀指標(biāo)除了靠肉眼觀察我也建議用幾個(gè)客觀指標(biāo)來(lái)輔助判斷。這些指標(biāo)在做項(xiàng)目報(bào)告、寫(xiě)論文時(shí)尤其有用。對(duì)比度Contrast濾波后的圖像感興趣區(qū)域的對(duì)比度是否提升了。可以用std2(filtered)計(jì)算標(biāo)準(zhǔn)差標(biāo)準(zhǔn)差越大說(shuō)明圖像灰度分布越分散對(duì)比度越高。信噪比SNR如果有一張干凈的參考圖可以計(jì)算濾波前后信噪比的變化確認(rèn)濾波有沒(méi)有把噪聲壓下去。能量保留率濾波后圖像的能量灰度平方和占原始圖像能量的比例。這個(gè)指標(biāo)能告訴你帶通濾波濾掉了多少信息能量保留率太低說(shuō)明你把有用信息也濾掉了。不過(guò)我得說(shuō)句實(shí)話工程上最靠譜的指標(biāo)還是“下游任務(wù)表現(xiàn)”你是做缺陷檢測(cè)的就看檢測(cè)率有沒(méi)有提升你是做血管分割的就看分割精度有沒(méi)有變好。所有圖像指標(biāo)都只是參考最終效果得看業(yè)務(wù)結(jié)果。4.4 我自己的調(diào)參心法調(diào)帶通參數(shù)踩了無(wú)數(shù)次坑之后我總結(jié)出一個(gè)原則從小帶寬開(kāi)始逐步向大帶寬擴(kuò)展。具體做法是先用頻譜圖鎖定目標(biāo)頻段的大致位置。設(shè)置一個(gè)較窄的通帶只包含目標(biāo)頻段。觀察濾波結(jié)果如果目標(biāo)結(jié)構(gòu)提取不完整說(shuō)明通帶偏窄適當(dāng)加大帶寬。如果結(jié)果里出現(xiàn)明顯噪聲或背景干擾說(shuō)明通帶偏寬往回收縮。反復(fù)迭代直到結(jié)果在“提取完整”和“噪聲干擾”之間達(dá)到平衡。這個(gè)方法雖然笨但比一次性設(shè)一個(gè)參數(shù)要可靠得多。因?yàn)閹V波器的參數(shù)和圖像內(nèi)容強(qiáng)相關(guān)不存在通用的最優(yōu)值必須針對(duì)具體的圖像反復(fù)試。5. 實(shí)戰(zhàn)中排過(guò)的坑完整排查鏈路與避坑手冊(cè)看完了原理和基本操作我來(lái)分享一些實(shí)際動(dòng)手時(shí)容易踩的坑。這些坑我基本都經(jīng)歷過(guò)有的查了半天資料才明白是怎么回事有的純靠運(yùn)氣才繞過(guò)去。把它們記下來(lái)希望你能少走彎路。5.1 坑一濾波結(jié)果一片黑問(wèn)題出在數(shù)據(jù)類型剛學(xué)MATLAB圖像處理的人十有八九遇到過(guò)這樣的問(wèn)題濾波出來(lái)的圖像黑乎乎一片什么都看不到。很大概率是你的圖像經(jīng)過(guò)傅里葉變換、乘法、反變換之后數(shù)據(jù)類型出了問(wèn)題。imread讀進(jìn)來(lái)的圖像是uint8類型范圍是 0~255但傅里葉變換和乘法運(yùn)算會(huì)產(chǎn)生浮點(diǎn)數(shù)而且范圍可能遠(yuǎn)大于 255 或包含負(fù)值。如果你直接用imshow顯示超范圍的像素會(huì)被裁剪負(fù)值被截?cái)酁?結(jié)果自然是一片黑。正確的做法是在濾波之前就im2double轉(zhuǎn)成雙精度浮點(diǎn)類型顯示的時(shí)候如果要?dú)w一化到0~1區(qū)間用imshow(filtered, [])讓MATLAB自動(dòng)調(diào)整顯示范圍。img im2double(img); % 計(jì)算前轉(zhuǎn)double filtered real(ifft2(ifftshift(G))); % 取實(shí)部丟棄虛部殘差 imshow(filtered, []); % 自動(dòng)歸一化顯示還有一個(gè)常被忽略的細(xì)節(jié)反傅里葉變換后理論上結(jié)果應(yīng)該是實(shí)數(shù)的但由于浮點(diǎn)計(jì)算誤差會(huì)產(chǎn)生很小的虛部。如果直接用ifft2的結(jié)果做imshow會(huì)報(bào)錯(cuò)或者顯示異常。必須用real()取實(shí)部或者至少abs()取模才能正常顯示。5.2 坑二fftshift和ifftshift用反了這個(gè)坑非常隱蔽而且報(bào)錯(cuò)不會(huì)很直觀。你可能會(huì)疑惑為什么我濾波的結(jié)果看起來(lái)像被旋轉(zhuǎn)了或者濾波效果出現(xiàn)在了圖像邊緣而不是中心原因就是fftshift和ifftshift混用。簡(jiǎn)單說(shuō)fftshift把頻譜從“零頻在角落”移到“零頻在中心”。ifftshift做逆操作把“零頻在中心”移回“零頻在角落”。正確的使用順序是F fft2(img); % 變換后零頻在角落 F_shift fftshift(F); % 移到中心方便構(gòu)造濾波器 G_shift F_shift .* H; % 濾波 G ifftshift(G_shift); % 移回角落準(zhǔn)備反變換 filtered real(ifft2(G)); % 反變換注意反變換之前一定要用ifftshift把頻譜還原成“零頻在角落”的布局否則ifft2出來(lái)的圖像會(huì)錯(cuò)位像被四分塊交換了一樣。5.3 坑三構(gòu)造濾波器時(shí)坐標(biāo)網(wǎng)格方向和矩陣尺寸不匹配meshgrid生成坐標(biāo)矩陣時(shí)第一個(gè)輸出對(duì)應(yīng)列坐標(biāo)u方向第二個(gè)輸出對(duì)應(yīng)行坐標(biāo)v方向。如果你把meshgrid(1:M, 1:N)寫(xiě)反了構(gòu)造出來(lái)的濾波器就是轉(zhuǎn)置的濾波效果會(huì)完全不對(duì)。這是我見(jiàn)過(guò)最多人犯的低級(jí)錯(cuò)誤包括我自己早期也栽過(guò)一次。一個(gè)比較穩(wěn)的寫(xiě)法是直接用尺寸變量來(lái)構(gòu)造[M, N] size(img); [u, v] meshgrid(1:N, 1:M); % u范圍是1到N列數(shù)v范圍是1到M行數(shù) D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); % 距離注意減的是N/2和M/2這里的關(guān)鍵在于u對(duì)應(yīng)圖像寬度方向列v對(duì)應(yīng)圖像高度方向行。中心點(diǎn)的坐標(biāo)就是 (N/21, M/21)。寫(xiě)完之后最好自己驗(yàn)證一下檢查D(1,1)是否等于sqrt((1-N/2-1)^2 (1-M/2-1)^2)也就是左上角到中心的距離。5.4 坑四理想帶通的振鈴把結(jié)果毀了前面說(shuō)過(guò)理想帶通濾波器的階躍截止會(huì)在空間域引入振鈴偽影。這個(gè)坑的實(shí)際表現(xiàn)是濾波結(jié)果里目標(biāo)周圍出現(xiàn)一圈一圈的波紋尤其是亮暗交替的邊緣處特別明顯。我之前做血管增強(qiáng)時(shí)用理想帶通把血管周圍的偽影當(dāng)成血管結(jié)構(gòu)分割出來(lái)的結(jié)果簡(jiǎn)直沒(méi)法看。改用高斯帶通之后偽影立刻消失效果立竿見(jiàn)影。所以我的建議是優(yōu)先使用高斯帶通或低階巴特沃斯帶通別碰理想帶通。除非你的圖像完全沒(méi)有強(qiáng)邊緣否則振鈴一定會(huì)來(lái)?yè)v亂。5.5 坑五大圖像濾波速度慢——頻域預(yù)計(jì)算與單精度優(yōu)化處理 2048x2048 甚至更大的圖像時(shí)頻域?yàn)V波的fft2和ifft2操作會(huì)比較耗時(shí)尤其是在老一點(diǎn)的電腦上跑起來(lái)極慢。我的解決辦法有兩個(gè)第一個(gè)是預(yù)計(jì)算濾波器。如果要對(duì)一批相同尺寸的圖像做同樣的濾波濾波器只需要構(gòu)造一次之后每次只做fft2→ 乘法 →ifft2能省下不少時(shí)間。% 預(yù)先構(gòu)造一次濾波器 [M, N] size(img_template); H gaussian_bandpass(M, N, 50, 0.4); for i 1 : num_files img im2double(imread(files{i})); filtered real(ifft2(ifftshift(fftshift(fft2(img)) .* H))); % 保存結(jié)果... end第二個(gè)是用單精度浮點(diǎn)數(shù)。雙精度計(jì)算雖然精度高但大矩陣的雙精度運(yùn)算非常吃內(nèi)存和時(shí)間。對(duì)于圖像濾波這種任務(wù)單精度完全夠用還能省一半內(nèi)存img im2single(img); F fft2(img); H single(H); % 濾波器也轉(zhuǎn)成單精度 G F .* H; filtered real(ifft2(G));我實(shí)測(cè)過(guò)單精度把2048x2048圖像的濾波時(shí)間縮短了差不多40%而且結(jié)果肉眼幾乎看不出差異。5.6 坑六邊界效應(yīng)的處理方式空間域?yàn)V波的邊界效應(yīng)是另一個(gè)高頻坑。imfilter默認(rèn)會(huì)用0填充邊界這會(huì)導(dǎo)致濾波后的圖像邊緣出現(xiàn)一圈暗紋因?yàn)檫吔缤馊?卷積結(jié)果在邊界處被拉低。解決方案是顯式指定邊界選項(xiàng)filtered imfilter(img, h_dog, replicate, same);replicate表示向外復(fù)制邊緣像素這是我用得最多的方式對(duì)大多數(shù)圖像效果都好。如果圖像邊緣有重要信息也可以用symmetric鏡像反射比replicate更平滑一點(diǎn)。千萬(wàn)別用默認(rèn)的0填充除非你不在乎邊緣那幾十個(gè)像素。頻域?yàn)V波也有類似的邊界問(wèn)題因?yàn)镕FT默認(rèn)把圖像看成周期的左右邊界、上下邊界會(huì)“連接”起來(lái)如果圖像左右邊界灰度差異大濾波結(jié)果邊緣會(huì)出現(xiàn)偽影。一個(gè)緩解的辦法是在FFT之前對(duì)圖像做邊緣平滑處理比如用edgeTaper或手動(dòng)對(duì)邊界加窗。6. 一個(gè)完整的應(yīng)用案例從需求到結(jié)果理論講了一堆最后用一個(gè)完整的應(yīng)用案例把這些內(nèi)容串起來(lái)。這個(gè)案例是我以前做過(guò)的——從一張帶有背景紋理干擾的電路板圖中提取出特定頻率的元件邊緣信息用于后續(xù)的缺陷定位。6.1 項(xiàng)目需求與圖像分析原始圖像是一張電路板的灰度圖元件邊緣清晰但表面有一層均勻的紋理干擾。生產(chǎn)線上需要根據(jù)邊緣信息定位元件邊界但紋理干擾導(dǎo)致邊緣檢測(cè)結(jié)果里全是假的邊緣。我先對(duì)圖像做了頻譜分析發(fā)現(xiàn)元件邊緣的能量主要集中在半徑約40到70像素的環(huán)形區(qū)域而紋理干擾的能量集中在更高的頻率段半徑約100以上。低頻區(qū)域主要是照明不均勻造成的背景變化。這個(gè)分析結(jié)果說(shuō)明需要一個(gè)通帶為40到70的帶通濾波器把元件邊緣提取出來(lái)同時(shí)濾掉低頻背景和高頻紋理。6.2 MATLAB完整實(shí)現(xiàn)%% 帶通濾波在元件邊緣提取中的應(yīng)用 % 讀取圖像 img imread(circuit.png); if size(img, 3) 3 img rgb2gray(img); end img im2double(img); %% 頻域分析和濾波器構(gòu)造 [M, N] size(img); [u, v] meshgrid(1:N, 1:M); D sqrt((u - N/2 - 1).^2 (v - M/2 - 1).^2); % 高斯帶通中心頻率55帶寬系數(shù)0.4 D0 55; sigma 0.4; H 1 - exp(-(D.^2 - D0^2).^2 / (2 * (D0 * sigma)^2)); figure; subplot(2,2,1); imshow(img); title(原始圖像); subplot(2,2,2); imshow(log(1 abs(fftshift(fft2(img)))), []); title(頻譜對(duì)數(shù)); subplot(2,2,3); imshow(H, []); title(帶通濾波器傳遞函數(shù)); % 頻域?yàn)V波 F_shift fftshift(fft2(img)); G_shift F_shift .* H; G ifftshift(G_shift); filtered real(ifft2(G)); subplot(2,2,4); imshow(filtered, []); title(帶通濾波結(jié)果); %% 后續(xù)處理邊緣檢測(cè) edges edge(filtered, canny, 0.15); figure; imshow(edges); title(濾波后的邊緣檢測(cè));6.3 效果對(duì)比與參數(shù)分析對(duì)比直接邊緣檢測(cè)和先帶通濾波再邊緣檢測(cè)的結(jié)果差異非常明顯處理流程檢測(cè)到的邊緣數(shù)量誤檢率邊緣連續(xù)性直接Canny大量偽邊緣高差邊緣斷裂嚴(yán)重帶通濾波Canny目標(biāo)邊緣清晰低好邊緣連續(xù)完整這個(gè)案例很好地說(shuō)明了帶通濾波在預(yù)處理中的價(jià)值它不直接做檢測(cè)而是把后面算法的輸入質(zhì)量大幅提升讓下游任務(wù)變得簡(jiǎn)單可靠。6.4 擴(kuò)展同一個(gè)框架可以做高通、低通掌握了帶通濾波的整個(gè)框架你會(huì)發(fā)現(xiàn)高通、低通都只是一個(gè)特例。低通只保留中心附近的低頻代碼上就是H(D D1) 0。高通把中心附近的低頻去掉H(D D0) 0。帶阻正好和帶通相反H(D D0 D D1) 0。它們共享同一套頻域處理的骨架只改濾波器構(gòu)造那幾行代碼就行。所以學(xué)會(huì)了帶通等于把整個(gè)頻域?yàn)V波家族都學(xué)會(huì)了。寫(xiě)在最后的經(jīng)驗(yàn)體會(huì)做了幾年圖像處理我的一個(gè)體會(huì)是帶通濾波這玩意兒看著簡(jiǎn)單用好了很見(jiàn)功力。它不像深度學(xué)習(xí)那些模型可以端到端地暴力出結(jié)果而是要求你對(duì)圖像內(nèi)容有足夠的理解——知道你要的信息在哪個(gè)頻段知道怎么把它跟其他頻率成分分開(kāi)。這種對(duì)數(shù)據(jù)和問(wèn)題的理解能力反而是越基礎(chǔ)的技巧越能體現(xiàn)出來(lái)。實(shí)際動(dòng)手的時(shí)候我建議你先別急著寫(xiě)代碼先把圖像的頻譜圖好好看一看。頻譜圖這東西看多了之后你會(huì)慢慢建立一種感覺(jué)——什么樣的圖像特征對(duì)應(yīng)什么樣的頻譜形狀這種感覺(jué)比任何公式都值錢(qián)。然后帶著這種感覺(jué)去調(diào)參數(shù)、做實(shí)驗(yàn)?zāi)銜?huì)發(fā)現(xiàn)帶通濾波的思路其實(shí)無(wú)處不在。最后再分享一個(gè)我一直在用的小技巧把濾波器可視化出來(lái)。不管是構(gòu)造好的傳遞函數(shù)圖還是空間域卷積核的二維圖都打印出來(lái)看一眼。濾波器本身長(zhǎng)什么樣決定了大半的濾波效果。我每次調(diào)參都會(huì)把頻譜圖、濾波器圖和濾波結(jié)果圖放在一起對(duì)比三個(gè)圖一起看問(wèn)題出在哪一目了然。這個(gè)習(xí)慣幫我避免了不少調(diào)試死胡同也讓我對(duì)“頻率”這個(gè)概念有了更實(shí)在的理解。本文還有配套的精品資源點(diǎn)擊獲取