% imregionalmax のプログラム
%
%   使い方: IMREGION_MAX_TEST を実行
%
%   使用するファイル: imregin_max_test.m
%
% See also IMREGIONALMAX, MAX, IMREGIONALMIN, MIN

clear;
close all;

%
% 入力画像の作成
input_image = [ 0 0 0 1 2 3;
                0 1 0 2 3 3;
                0 0 0 2 3 4;
                1 3 1 5 4 5;
                0 0 0 3 2 4;
                2 1 0 1 2 3];

% input_image = [ 0 2 3 3 1 0;
%                 0 0 0 0 0 0;
%                 0 0 5 0 0 0;
%                 0 0 0 0 0 0;
%                 4 0 0 0 0 0;
%                 4 5 0 0 0 0];

size_input_x = 6; size_input_y = 6; 

output_image = ones(size_input_x, size_input_y);

% ここから開始
input = input_image;
output = output_image;
output_cp = zeros(size_input_x, size_input_y);
nloop = 0;
%
% 最大値の可能性の行列が変化しなくなったら終了とする．
while sum(sum(output_cp ~= output)) ~= 0
    output_cp = output;
    nloop = nloop + 1;

    %% Middle Section
    for m=1+1:size_input_y-1
        for n=1+1:size_input_x-1
            %
            % 周りの8近傍について，最大値を求める．
            max_val = max([...
                input(m-1,n-1), input(m-1,n), input(m-1,n+1),...
                input(m  ,n-1),               input(m  ,n+1),...
                input(m+1,n-1), input(m+1,n), input(m+1,n+1)]);
            %
            % 周りの最大値より注目画素の値が小さい場合，
            % 注目画素は最大値にならないため，
            % 最大値の可能性を保持する行列にFalse（0）を代入
            if max_val > input(m,n)
                output(m,n) = 0;
            end
            %
            % 周りの最大値と注目画素の値が同値の場合，
            % 注目画素は最大値となり得る可能性があるため，
            % 最大値の可能性を保存してる配列を確認する。
            if max_val == input(m,n)
                %
                % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
                % 注目画素も最大値になり得ないため、False（0）を代入　
                % 以下，8近傍について，最大値の可能性があるかどうかを見ている．
                if max_val == input(m-1,n-1)
                    if output(m-1,n-1) == 0
                        output(m,n) = 0;
                    end
                end
                if max_val == input(m-1,n)
                    if output(m-1,n) == 0
                        output(m,n) = 0;
                    end
                end
                if max_val == input(m-1,n+1)
                    if output(m-1,n+1) == 0
                        output(m,n) = 0;
                    end
                end

                if max_val == input(m,n-1)
                    if output(m,n-1) == 0
                        output(m,n) = 0;
                    end
                end
                if max_val == input(m,n+1)
                    if output(m,n+1) == 0
                        output(m,n) = 0;
                    end
                end

                if max_val == input(m+1,n-1)
                    if output(m+1,n-1) == 0
                        output(m,n) = 0;
                    end
                end
                if max_val == input(m+1,n)
                    if output(m+1,n) == 0
                        output(m,n) = 0;
                    end
                end
                if max_val == input(m+1,n+1)
                    if output(m+1,n+1) == 0
                        output(m,n) = 0;
                    end
                end
            end
        end 
    end

    %% Upper Edge
    m = 1;
    %%% Upper Left
    for n=1:1
    % 周りの3近傍について，最大値を求める．
        max_val = max([...
                                              input(m  ,n+1),...
                                input(m+1,n), input(m+1,n+1)]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．    
            if max_val == input(m,n+1)
                if output(m,n+1) == 0
                    output(m,n) = 0;
                end
            end
    
            if max_val == input(m+1,n)
                if output(m+1,n) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m+1,n+1)
                if output(m+1,n+1) == 0
                    output(m,n) = 0;
                end
            end
        end
    end
    %%% Upper Edge
    for n=1+1:size_input_x-1
        %
        % 周りの5近傍について，最大値を求める．
        max_val = max([...
                input(m  ,n-1),               input(m  ,n+1),...
                input(m+1,n-1), input(m+1,n), input(m+1,n+1)]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．    
            if max_val == input(m,n-1)
                if output(m,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m,n+1)
                if output(m,n+1) == 0
                    output(m,n) = 0;
                end
            end
    
            if max_val == input(m+1,n-1)
                if output(m+1,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m+1,n)
                if output(m+1,n) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m+1,n+1)
                if output(m+1,n+1) == 0
                    output(m,n) = 0;
                end
            end
        end
    end
    %%% Upper Right
    for n=size_input_x:size_input_x
        % 周りの3近傍について，最大値を求める．
        n = size_input_x;
        max_val = max([...
                input(m  ,n-1),               ...
                input(m+1,n-1), input(m+1,n)]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．    
            if max_val == input(m,n-1)
                if output(m,n-1) == 0
                    output(m,n) = 0;
                end
            end
    
            if max_val == input(m+1,n-1)
                if output(m+1,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m+1,n)
                if output(m+1,n) == 0
                    output(m,n) = 0;
                end
            end
        end
    end

    %% Lower Edge
    m = size_input_y;
    %%% Lower Left
    for n=1:1
    % 周りの3近傍について，最大値を求める．
        max_val = max([...
            input(m-1,n), input(m-1,n+1),...
                           input(m  ,n+1),...
            ]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．    
            if max_val == input(m-1,n)
                if output(m-1,n) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m-1,n+1)
                if output(m-1,n+1) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m,n+1)
                if output(m,n+1) == 0
                    output(m,n) = 0;
                end
            end
        end
    end
    %%% Lower Edge
    for n=1+1:size_input_x-1
        %
        % 周りの5近傍について，最大値を求める．
        max_val = max([...
                input(m-1,n-1), input(m-1,n), input(m-1,n+1),...
                input(m  ,n-1),               input(m  ,n+1),...
                ]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．    
            if max_val == input(m-1,n-1)
                if output(m-1,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m-1,n)
                if output(m-1,n) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m-1,n+1)
                if output(m-1,n+1) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m,n-1)
                if output(m,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m,n+1)
                if output(m,n+1) == 0
                    output(m,n) = 0;
                end
            end
        end
    end
    %%% Lower Right
    for n=size_input_x:size_input_x
        % 周りの3近傍について，最大値を求める．
        n = size_input_x;
        max_val = max([...
            input(m-1,n-1), input(m-1,n), ...
            input(m  ,n-1),               ...
            ]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．
            if max_val == input(m-1,n-1)
                if output(m-1,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m-1,n)
                if output(m-1,n) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m,n-1)
                if output(m,n-1) == 0
                    output(m,n) = 0;
                end
            end
        end
    end

    %% Left Edge
    n = 1;
    for m=1+1:size_input_y-1
        %
        % 周りの5近傍について，最大値を求める．
        max_val = max([...
             input(m-1,n), input(m-1,n+1),...
                           input(m  ,n+1),...
             input(m+1,n), input(m+1,n+1)]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．
            if max_val == input(m-1,n)
                if output(m-1,n) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m-1,n+1)
                if output(m-1,n+1) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m,n+1)
                if output(m,n+1) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m+1,n)
                if output(m+1,n) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m+1,n+1)
                if output(m+1,n+1) == 0
                    output(m,n) = 0;
                end
            end
        end    
    end

    %% Right Edge
    n = size_input_x;
    for m=1+1:size_input_y-1
        % 周りの5近傍について，最大値を求める．
        max_val = max([...
            input(m-1,n-1), input(m-1,n), ...
            input(m  ,n-1),               ...
            input(m+1,n-1), input(m+1,n)]);
        %
        % 周りの最大値より注目画素の値が小さい場合，
        % 注目画素は最大値にならないため，
        % 最大値の可能性を保持する行列にFalse（0）を代入
        if max_val > input(m,n)
            output(m,n) = 0;
        end
        %
        % 周りの最大値と注目画素の値が同値の場合，
        % 注目画素は最大値となり得る可能性があるため，
        % 最大値の可能性を保存してる配列を確認する。
        if max_val == input(m,n)
            %
            % 周りの画素の最大値について，最大値の可能性がない場合（False 0），
            % 注目画素も最大値になり得ないため、False（0）を代入　
            % 以下，8近傍について，最大値の可能性があるかどうかを見ている．
            if max_val == input(m-1,n-1)
                if output(m-1,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m-1,n)
                if output(m-1,n) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m,n-1)
                if output(m,n-1) == 0
                    output(m,n) = 0;
                end
            end

            if max_val == input(m+1,n-1)
                if output(m+1,n-1) == 0
                    output(m,n) = 0;
                end
            end
            if max_val == input(m+1,n)
                if output(m+1,n) == 0
                    output(m,n) = 0;
                end
            end
        end
    end
    
    fprintf('Loop %d\n',nloop);
    disp(output);
end

% matlab_func_output_image = imregionalmax(input);

fprintf('Iutput\n');
disp(input);
fprintf('Output\n');
disp(output);

%% 
% MATLABの関数と比較したい場合は，以下をコメントアウト
% (但し，Signal Processing Toolboxが必要)
%
% fprintf('Matlab function : imregionalmax\n');
% disp(matlab_func_output_image);