本文參考博客：https://www.cnblogs.com/Imageshop/archive/2013/04/26/3045672.html
原生的中值濾波是基于排序演算法的，這樣的演算法復雜度基本在O(r2)左右，當濾波半徑較大時，排序演算法就顯得很慢，對此有多種改進演算法，這里介紹經典
的Huang演算法與O(1)演算法，兩者都是基于記憶性的演算法，只是后者記性更強，
排序演算法明顯的一個不足之處就是無記憶性，當核向右移動一列后，只是核的最左和最右列資料發生了變化，中間不變的資料應當被存盤起來，而排序演算法
并沒有做到這點，Huang演算法的思想是建立一個核直方圖，來統計核內的各灰度的像素數，當核向右移動時，就將新的一列所有資料加入到直方圖中，同時將最
左列的舊資料從直方圖中洗掉，如下圖所示，這樣做使得大部分資料能夠被記憶，減少重復操作，當直方圖更新完畢后，就可以通過從左到右累計像素來找到中值，
下面是演算法具體實作步驟與代碼：
1.每行開始都將直方圖、像素計數、中值變數清零，將核覆寫的所有像素加到直方圖中，
2.計算中值，sumcnt為小于中值灰度的像素數和，如果當前sumcnt大于等于閾值，則表明，實際中值比當前median小，則直方圖向左減去像素數，同時median
也減小，直到sumcnt小于閾值，則此時median為中值（sumcnt加上中值像素數將大于或等于閾值）；如果sumcnt加上median像素數仍小于閾值，表明，實
際中值比當前median大，則直方圖向右累加像素數，同時median增加，直到sumcnt加上median像素數大于等于閾值，則此時median為中值（sumcnt小于閾值），
3.核每向右移動一列，更新核直方圖，比較每一個新加入的灰度值與median，小于則median加1，比較每一個減去的灰度值與median，小于則median減1，執行
一次步驟2后繼續移動核，直到整幅圖都被濾波完畢，
%中值濾波黃法%輸入：8位深度影像、濾波半徑（3*3的半徑為1,中值百分比）function grayimg=medianfilterhuang(grayimg,radius,percentage) [sizex,sizey]=size(grayimg); thres=int16(((2*radius+1)^2*percentage)); %擴展邊界 extendimg=int16(horzcat(repmat(grayimg(:,1),1,radius),grayimg,repmat(grayimg(:,sizey),1,radius))); extendimg=int16(vertcat(repmat(extendimg(1,:),radius,1),extendimg,repmat(extendimg(sizex,:),radius,1))); for i=radius+1:sizex+radius %初始化直方圖，每行第一個將所有像素加入直方圖 histogram=int16(zeros(256,1)); sumcnt=int16(0); median=int16(0); for j=radius+1:sizey+radius %遍歷卷積覆寫行 for k=i-radius:i+radius if (j>radius+1) %刪去左列，加入新右列 histogram(extendimg(k,j-radius-1)+1)=histogram(extendimg(k,j-radius-1)+1)-1; histogram(extendimg(k,j+radius)+1)=histogram(extendimg(k,j+radius)+1)+1; %只關心中值灰度以左的像素數量 if (extendimg(k,j-radius-1)<median) sumcnt=sumcnt-1; end if (extendimg(k,j+radius)<median) sumcnt=sumcnt+1; end else %行首 for h=j-radius:j+radius histogram(extendimg(k,h)+1)=histogram(extendimg(k,h)+1)+1; end end end %sumcnt不將中值像素數累加進去 %舊中值偏大 while(sumcnt>=thres) median=median-1; sumcnt=sumcnt-histogram(median+1); end %舊中值偏小 while(sumcnt+histogram(median+1)<thres) sumcnt=sumcnt+histogram(median+1); median=median+1; end grayimg(i-radius,j-radius)=median; end endendHuang演算法
性能分析發現，Huang演算法的大部分時間花在了更新核直方圖上，并且由于它是O(r)的復雜度，因此當濾波半徑變大后，更新直方圖部分花費的時間也變多，而計算
中值部分所花時間一直都是固定的，因為它的復雜度是O(1)，與濾波半徑無關，那么為了使直方圖更新的復雜度能降到O(1)，就要增強演算法的記憶性，可以看到以相鄰行
相同列為中心的核之間也存在重合，這跟之前以相鄰列相同行為中心的核存在重合類似，于是，我們可以為每列都維護一個直方圖，而核直方圖可以由這些列直方圖相加
得到，更新核的時候，最右邊的一列直方圖向下移動一行，就更新了列直方圖，再減去最左邊一列直方圖，核就得到更新，這樣每次更新直方圖都只要移動一次列直方圖，
而其余的都是直方圖的加減操作，值得注意的是，這樣的實作很完美得利用了可并行的矩陣加減來代替順序執行的for迭代，因此，雖然更新程序看上去有許多加減操作，
實際上都是可快速并行執行的，
還有一點，由于并行化的矩陣加減，就不能用sumcnt來跟蹤閾值變化了，如果直接迭代累加直方圖，最壞的情況是要執行255加法，為了使中值計算更快，可以采用
兩層直方圖，高層用于保存高4位（直方圖大小為16*1），低層保存全位（直方圖大小為256*1），每次先累加高層進行范圍縮小，再從相應范圍累加低層找到中值，只是
這樣就犧牲了記憶體空間，
下面是演算法具體實作步驟與代碼：
1.擴展影像邊界，上邊界擴展r+1行，其他邊界擴展r行；為每列生成兩層直方圖，并為其初始化，添加各列的前2*r+1個像素；生成一個低層鏈接高層的索引矩陣，
2.每行開始清零核直方圖，更新核覆寫的各列直方圖（向下移動一行），將各列加到核直方圖中，
3.尋找中值，初始化sumcnt=0，這里它表示小于等于中值灰度的像素和，從左到右迭代高層直方圖，直到sumcnt大于等于閾值；再從高層確定的范圍內從右到左迭代
低層直方圖，直到sumcnt（減去了迭代灰度的像素數）小于閾值，則當前迭代值為中值，
4.核每向右移動，更新列直方圖與核直方圖，執行一次步驟3，然后繼續移動核，直到完成整個圖的濾波，
function grayimg=medianfilterdp(grayimg,radius,percentage) [sizex,sizey]=size(grayimg); thres=int16(((2*radius+1)^2*percentage)); %擴展邊界 extendimg=int16(horzcat(repmat(grayimg(:,1),1,radius),grayimg,repmat(grayimg(:,sizey),1,radius))); extendimg=int16(vertcat(repmat(extendimg(1,:),radius+1,1),extendimg,repmat(extendimg(sizex,:),radius,1))); %為每一列維護一個全位直方圖和一個高位直方圖 histogram_cols_low=int16(zeros(sizey+2*radius,256)); histogram_cols_high=int16(zeros(sizey+2*radius,16)); %中值計算所用的一個全位直方圖和一個高位直方圖 histogram_low=int16(zeros(1,256)); histogram_high=int16(zeros(1,16)); %全位轉高位的加速索引 grayrange=int16(zeros(1,256)); for i=1:16 grayrange(:,(i-1)*16+1:i*16)=i; end %初始化各列直方圖 for row=1:2*radius+1 for col=1:sizey+2*radius histogram_cols_low(col,extendimg(row,col)+1)=histogram_cols_low(col,extendimg(row,col)+1)+1; histogram_cols_high(col,grayrange(extendimg(row,col)+1))=histogram_cols_high(col,grayrange(extendimg(row,col)+1))+1; end end for row=int16(radius+2:sizex+radius+1) for col=int16(radius+1:sizey+radius) if(col>radius+1) %更新右邊新列，洗掉左邊舊列 histogram_cols_low(col+radius,extendimg(row+radius,col+radius)+1)=histogram_cols_low(col+radius,extendimg(row+radius,col+radius)+1)+1; histogram_cols_low(col+radius,extendimg(row-radius-1,col+radius)+1)=histogram_cols_low(col+radius,extendimg(row-radius-1,col+radius)+1)-1; histogram_cols_high(col+radius,grayrange(extendimg(row+radius,col+radius)+1))=histogram_cols_high(col+radius,grayrange(extendimg(row+radius,col+radius)+1))+1; histogram_cols_high(col+radius,grayrange(extendimg(row-radius-1,col+radius)+1))=histogram_cols_high(col+radius,grayrange(extendimg(row-radius-1,col+radius)+1))-1; histogram_low=histogram_low+histogram_cols_low(col+radius,:)-histogram_cols_low(col-radius-1,:); histogram_high=histogram_high+histogram_cols_high(col+radius,:)-histogram_cols_high(col-radius-1,:); else histogram_low(:,:)=0; histogram_high(:,:)=0; %行首手動初始化各列 for index=1:2*radius+1 histogram_cols_low(index,extendimg(row+radius,index)+1)=histogram_cols_low(index,extendimg(row+radius,index)+1)+1; histogram_cols_low(index,extendimg(row-radius-1,index)+1)=histogram_cols_low(index,extendimg(row-radius-1,index)+1)-1; histogram_cols_high(index,grayrange(extendimg(row+radius,index)+1))=histogram_cols_high(index,grayrange(extendimg(row+radius,index)+1))+1; histogram_cols_high(index,grayrange(extendimg(row-radius-1,index)+1))=histogram_cols_high(index,grayrange(extendimg(row-radius-1,index)+1))-1; histogram_low=histogram_low+histogram_cols_low(index,:); histogram_high=histogram_high+histogram_cols_high(index,:); end end %中值及之前的所有像素數 sumcnt=int16(0); for gray_high=1:16 sumcnt=sumcnt+histogram_high(gray_high); if (sumcnt>=thres) break; end end for gray_low=gray_high*16:-1:(gray_high-1)*16+1 sumcnt=sumcnt-histogram_low(gray_low); if (sumcnt<thres) break; end end grayimg(row-radius-1,col-radius)=gray_low-1; end end endO(1)中值濾波
下面是兩種演算法的性能對比，可以發現，在濾波半徑較小的時候，用Huang演算法更好；當濾波半徑>=8時，用后一種演算法更快，實際中我們一般用不到大半徑，但這兩種
演算法的思想很不錯，尤其是第二種的空間換時間，我們可能可以發現橫向存在重復性，但很容易忽略了縱向上的重復性，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/3805.html
標籤：其他
上一篇：CGAL的安裝與使用