多分辨率峰谷探测(阈值选取)

利用图像的直方图帮助选择阈值是常用的方法，其中的关键是确定峰点和谷点。由于场景的复杂性，图像成像过程中各种干扰因素的存在等原因，峰点和谷点的有无检测和位置确定常常比较困难。峰点和谷点的检测和直方图的尺度有密切的联系。一般在较大尺度下常能较可靠地检测到真正的峰点和谷点，但在大尺度下对峰点和谷点的定位不易准确。相反，在较小尺度下对真正峰点和谷点的定位常常比较准确，但在小尺度下误检或漏检的比例会增加。因此可以考虑进行多分辨率峰谷探测，在较大尺度下检测峰谷的有无，再在精细尺度下确定峰谷的位置。

具体来说，多尺度是通过调节高斯滤波器的sigma实现的。滤波平滑后信号（直方图）的二阶导的过零点对应着信号的拐点，即波峰或波谷的起始或终止点（波峰的起始点可视为波谷的终止点）。其中，从负值到正值的零交叉点对应峰的终点，从正值到负值的零交叉点对应峰的起点。在波峰的起点与终点之间的最大值点对应该尺度下峰的位置，在波峰的终点与下一个波峰的起点之间的最小值点对应该分辨率下波谷的位置。原理参见下面示意图（其中黄线代表放大后的二阶导）：

 

在确定了一个尺度的峰谷位置后，通过向上一尺度（更低尺度）搜索与当前尺度下波峰/波谷距离最近的波峰/波谷位置来替代当前尺度下的波峰/波谷位置来进行更精确的定位。该过程从最大尺度不断向上迭代，直到到达设定的起始尺度（可以是原始分辨率）为止。

下面给出本人编写的matlab参考代码：

clear all;close all;clc;
x=(pi/50)*(1:256);
hist=50*cos(x)+30*cos(2*x)+randi([1,10],[1,256]);
figure,plot(1:256,hist);
%%
[top,bottom]=multiScaleThSelect(hist,3,1);

function [top,bottom]=multiScaleThSelect(hist,S,s0)
if ~exist('s0','var')
    s0=0;
end
% 多尺度阈值选取
% S为总的分解层数
% s0为起始尺度，0则对应原始分辨率
% 各层分辨率分别为2^(s-1), s=1,2,...S
hist=reshape(hist,[1,length(hist)]);
multiScaleTop=cell(S+1-s0,1);
multiScaleBottom=cell(S+1-s0,1);
for s=S:-1:s0
    if s>0
        sigma=2^(s-1);
        f=fspecial('gaussian',[1,6*sigma],sigma);
        fHist=imfilter(hist,f,'symmetric');
    else
        fHist=hist;
    end
    
    d2f_fHist=diff(diff(fHist));
    d2f_fHist=[d2f_fHist(1) d2f_fHist(1)  d2f_fHist];
    
    zeroPt=[]; %二阶导过零点，为nx2矩阵，第一列标记位置，第二列标记过零点类型 
    for i=2:length(hist)-1
        if d2f_fHist(i-1)<0 && d2f_fHist(i+1)>0
            zeroPt=[zeroPt;i,-1]; %峰的终点
        else
            if d2f_fHist(i-1)>0 && d2f_fHist(i+1)<0
                zeroPt=[zeroPt;i,1]; %峰的起点
            end
        end
    end
    
    top=[];
    bottom=[];
    for i=2:length(zeroPt)
        if zeroPt(i-1,2)==-1 && zeroPt(i,2)==1
            interval=fHist(zeroPt(i-1,1):zeroPt(i,1));
            index=find(interval==min(interval));
            bottom=[bottom,zeroPt(i-1,1)-1+index(1)];
        elseif zeroPt(i-1,2)==1 && zeroPt(i,2)==-1
            interval=fHist(zeroPt(i-1,1):zeroPt(i,1));
            index=find(interval==max(interval));
            top=[top,zeroPt(i-1,1)-1+index(1)];
        end
    end
    
    if s0>0
        multiScaleTop(s)={top};
        multiScaleBottom(s)={bottom};
    else
        multiScaleTop(s+1)={top};
        multiScaleBottom(s+1)={bottom};
    end
end

top=multiScaleTop{S+1-s0};
bottom=multiScaleBottom{S+1-s0};
for i=S-s0:-1:1
    for k=1:length(top)
        dis=abs(multiScaleTop{i}-top(k));
        index=find(dis==min(dis(:)));
        top(k)=multiScaleTop{i}(index(1));
    end
    for k=1:length(bottom)
        dis=abs(multiScaleBottom{i}-bottom(k));
        index=find(dis==min(dis(:)));
        bottom(k)=multiScaleBottom{i}(index(1));
    end
end
top=unique(top);
bottom=unique(bottom);
end