MATLAB寻找三维曲面的有效波峰（滤除噪声产生的极大值）

最近在做课题研究的时候，在三维曲面寻找波峰方面遇到了些困难，虽然查找资料的时候有看到findpeaks和imregionalmax两个找极大值的函数，要么就是只适用于二位曲线，要么就是不能很好地滤除噪声引起的无效波峰，最后偶然看到一个大佬的思想，给了我很大的启发。

本文参考波峰识别算法介绍 | Zealseeker's Blog中的峰突prominence思想，对该识别算法进行说明改进，以利用MATLAB工具实现三维曲面中的有效波峰标定。

参考文章中对二维曲线进行举例说明，如下图所示，共有5个可见波峰，有效波峰用o表示，噪声引起的无效波峰用x表示，最大高度为5个单位长度。

假设对第3个x波峰进行检测，计算该点向左的波谷left_base，从波峰依次向左前进，若识别到该点比前一个点的高度更低，则继续向左前进，直至到达波谷，认为该点为left_base，同理向右前进可以得到right_base。

然后计算峰突强度（即峰突出的强度）prominence，假设该波峰x的高度为h，则

prominence=min(h-left_base,h-right_base)

可以得到，该波峰的left_base=3-2=1，right_base=3-1=2，则prominence=min(1，2)=1。

此时我们便可以根据实际情况设置阈值，滤除所有prominence小于阈值的波峰，例如上图中可以设阈值threshold=3，则三个x波峰都将被作为噪声滤除。

基于上述思想，同理可以寻找三维曲面中的有效波峰，首先用3*3窗口找出所有极大值（边界默认不是极值），然后计算每个极大值的prominence，由于是三维曲面，为了增加寻峰的有效性，除了左右两个波谷，还另加上前后两个波谷，最后取4个差值中的最小值作为prominence，再根据需求设定阈值便可以去除无效波峰啦~

 上图为直接使用imregionalmax函数所得到的波峰（红点标记）。

 上图为通过计算峰突强度prominence所得到的波峰（红点标记）。

比较两种方法，本文所述方法较好地滤除了噪声引起的无效波峰。

以下为本人手敲的函数代码，需要的自取，若有好的建议或是还有哪些更好的寻峰方法，还请各位大佬不吝赐教。

function peaks = findPeak_3D(Z,threshold)
% 基于峰突强度，找出二维数据中的波峰，Z为输入的二维数据（即曲面各点的高度），threshold为峰突强度阈值，
% 若对应极大值点的峰突强度大于所设阈值，则视为波峰。
% 返回大小与Z一致的二值矩阵SM，0代表不是波峰，1代表是波峰，边界默认为0。

m_size = size(Z);
width = m_size(1);
len = m_size(2);
peaks = zeros(width,len);

for i = 2:width-1
    for j = 2:len-1
       if Z(i,j)>=Z(i-1,j) && Z(i,j)>=Z(i+1,j) && Z(i,j)>=Z(i,j-1) &&...
               Z(i,j)>=Z(i,j+1) && Z(i,j)>=Z(i-1,j+1) && Z(i,j)>=Z(i+1,j+1) &&...
               Z(i,j)>=Z(i-1,j-1) && Z(i,j)>=Z(i+1,j-1) % 检测3*3窗口极大值
           peaks(i,j) = 1;
       end
    end
end

% 通过前后左右4个波谷计算峰突强度prominence
[x,y] = find(peaks == 1);  % 极大值坐标
for i = 1:length(x)
    
    count = 1;
    while Z(x(i),y(i)-count) < Z(x(i),y(i)-count+1)
        count = count + 1;
        if y(i)-count <= 0
            break;
        end
    end
    left_base = Z(x(i),y(i)-count+1);   % 左波谷
    count = 1;
    while Z(x(i),y(i)+count) < Z(x(i),y(i)+count-1)
        count = count + 1;
        if y(i)+count >= len-1
            break;
        end
    end
    right_base = Z(x(i),y(i)+count-1);   % 右波谷
    count = 1;
    while Z(x(i)-count,y(i)) < Z(x(i)-count+1,y(i))
        count = count + 1;
        if x(i)-count <= 0
            break;
        end
    end
    former_base = Z(x(i)-count+1,y(i));   % 前波谷
    count = 1;
    while Z(x(i)+count,y(i)) < Z(x(i)+count-1,y(i))
        count = count + 1;
        if x(i)+count >= width-1
            break;
        end
    end
    later_base =  Z(x(i)+count-1,y(i));   % 后波谷
    
    a = Z(x(i),y(i))-left_base;
    b = Z(x(i),y(i))-right_base;
    c = Z(x(i),y(i))-former_base;
    d = Z(x(i),y(i))-later_base;
    prominence = min([a b c d]);%峰突强度
    if prominence < threshold % 未达到阈值的极大值不视为波峰
        peaks(x(i),y(i)) = 0;
    end
end