一维曲线极值检测

clc;clear;close all;
BW1 = imread(‘BW1.BMP‘);
figure,imshow(BW1,[]);

% 提取图像中的曲线
[M,N] = size(BW1);
x=[];y=[];
j=1;
while j<=N,
    for i = 1 : M
     try
        if BW1(i, j) == 1
          x = [x, N - j]; %#ok<*AGROW>
          y = [y, M - i];
          j=j+1;
          break;
        end
     catch
        fprintf(1,‘%d, %d\n‘,i,j);
     end
    end
end


xx = max(x) : -0.1 : min(x);
yy = interp1(x, y, xx, ‘linear‘); % 一维数据插值cublic
x = xx;
y = yy;
N = size(y, 2);
y(N) = y(N - 1) + 1;
left = 1;
right = 0;
Sk = 0; % 0：斜率未知，1：上升，-1：下降
peak = [];

for j = 2 : N,
    if Sk == 0
        if (y(j) <= y(j - 1))
            Sk = 0;          % 下降及保持信号、初始信号
            continue;
        else
            Sk = 1;          % 上升信号
            left = j;
            continue;
        end;
    end;
    
    if (Sk == 1)
        if (y(j) < y(j - 1))
            Sk = -1;         % 下降信号
            continue;
        end       
    else if (Sk == -1)
            if (y(j) < y(j - 1))
                right = j;
            end;
            if (y(j) > y(j - 1))
                Sk = 0;
                disp([left, right]);
                peak = [peak; left, right];
                continue;
            end           
        end     
    end
    
    if Sk == 0
        if (y(j) <= y(j - 1))
            Sk = 0;
            continue;
        end;
        if (y(j) > y(j - 1))
            Sk = 1;
            left = j;
            continue;
        end;
    end;    
end;

% 画出峰值
figure,
plot(y);
hold on;
[M,N] = size(peak);
for j = 1 : M
    minP = min(peak(j, 1), peak(j, 2));
    maxP = max(peak(j, 1), peak(j, 2));
    [C, I] = max(y(minP : maxP));
    plot([minP + I, minP + I], [0, C], ‘g‘);
end;