function [peakindex2D, peakval2D] = findpeak2D(data2D, th)

[Nrow, Ncol] = size(data2D);

peakindex2D = zeros(Nrow, Ncol);
peakval2D = zeros(Nrow, Ncol);

if isempty(th)
    for rowidx = 1 : Nrow
        for colidx = 1 : Ncol
            indexset = [rowidx-1 colidx; rowidx+1 colidx; rowidx colidx-1; rowidx colidx+1];
            det_indexset = (indexset(:,1) > 0) .* (indexset(:,1) < Nrow+1) .* (indexset(:,2) > 0) .* (indexset(:,2) < Ncol+1);
            det_indexset = find(det_indexset > 0);
            cond_peak = 0;
            index_iter_len = length(det_indexset);
            for indexset_idx = 1 : length(det_indexset)
                if data2D(rowidx, colidx) >= data2D(indexset(det_indexset(indexset_idx),1), indexset(det_indexset(indexset_idx),2))
                    cond_peak = cond_peak + 1;
                end            
            end

            if cond_peak == index_iter_len
                peakindex2D(rowidx, colidx) = 1;
                peakval2D(rowidx, colidx) = data2D(rowidx, colidx);
            end
        end
    end
else
    for rowidx = 1 : Nrow
        for colidx = 1 : Ncol
            indexset = [rowidx-1 colidx; rowidx+1 colidx; rowidx colidx-1; rowidx colidx+1];
            det_indexset = (indexset(:,1) > 0) .* (indexset(:,1) < Nrow+1) .* (indexset(:,2) > 0) .* (indexset(:,2) < Ncol+1);
            det_indexset = find(det_indexset > 0);
            cond_peak = 0;
            index_iter_len = length(det_indexset);
            for indexset_idx = 1 : length(det_indexset)
                if data2D(rowidx, colidx) >= data2D(indexset(det_indexset(indexset_idx),1), indexset(det_indexset(indexset_idx),2))
                    cond_peak = cond_peak + 1;
                end            
            end

            if cond_peak == index_iter_len && data2D(rowidx, colidx) >= th
                peakindex2D(rowidx, colidx) = 1;
                peakval2D(rowidx, colidx) = data2D(rowidx, colidx);
            end
        end
    end
end

end