killer.m

关于偏微分方程进行图像处理的电子书,对于学习动态轮廓实现图像的分割非常有用

function killer( image_name,threshold,out_file_name )
% 最大最小算子
if nargin<3
    disp('请按参数说明输入参数!');
    return;
end

image_I=imread(image_name);
subplot(1,2,1);
imshow(image_I);
[size_r,size_c]=size(image_I);   % 图像大小
flag_image=zeros(size_r+2,size_c+2); % 标志矩阵

image_ex=matrix_expand(image_I);  % 边界处理
% －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
for row=2:size_r+1
    for col=2:size_c+1        
        if flag_image(row,col)==0 & image_ex(row,col)>=image_ex(row-1,col) & image_ex(row,col)>=image_ex(row+1,col) & image_ex(row,col)>=image_ex(row,col-1) & image_ex(row,col)>=image_ex(row,col+1)
               stack_r=[ row ];  % 行列堆栈初始化
               stack_c=[ col ];
               flag_image(row,col)=1;
               alpha=image_ex(row,col);
            
               while 1
                     [max_row,max_col]=seek_max_elem(stack_r,stack_c,image_ex );  % 取堆栈内元素的所有邻域的最大值元素 
                   
                     stack_r=[ stack_r max_row ];
                     stack_c=[ stack_c max_col ];                
                     flag_image(max_row,max_col )=1;
                   
                     if image_ex(max_row,max_col)<=alpha
                          alpha=image_ex(max_row,max_col);
                     end
                     if image_ex(max_row,max_col)>alpha | length(stack_r)>threshold 
                          for index=1:length(stack_r)
                              index_r=stack_r(index);
                              index_c=stack_c(index);
                              image_ex(index_r,index_c)=alpha;
                          end
                          break;
                     end   
               end
        elseif flag_image(row,col)==0 & image_ex(row,col)<=image_ex(row-1,col) & image_ex(row,col)<=image_ex(row+1,col) & image_ex(row,col)<=image_ex(row,col-1) & image_ex(row,col)<=image_ex(row,col+1)
               stack_r=[ row ];  % 行列堆栈初始化
               stack_c=[ col ];
               alpha=image_ex(row,col);
               flag_image(row,col )=1;
            
               while 1
                     [min_row,min_col]=seek_min_elem( stack_r,stack_c,image_ex );  % 取堆栈内元素的所有邻域的最大值元素 
                   
                     stack_r=[ stack_r min_row ];
                     stack_c=[ stack_c min_col ];                
                     flag_image(min_row,min_col )=1;
                   
                     if image_ex(min_row,min_col)>=alpha  % 取最大
                         alpha=image_ex(min_row,min_col);
                     end                   
                     if image_ex(min_row,min_col)<alpha | length(stack_r)>threshold 
                          for index=1:length(stack_r)
                              index_r=stack_r(index);
                              index_c=stack_c(index);
                              image_ex(index_r,index_c)=alpha;
                          end
                          break;
                     end                
               end
        end
    end
end

for row=2:size_r+1
    for col=2:size_c+1 
        image_I(row-1,col-1)=image_ex(row,col);
    end
end
subplot(1,2,2);
imshow(image_I);
imwrite(image_I,out_file_name);
% －－－－－－－－－－－－－－－－－－－

function [max_r,max_c]=seek_max_elem( array_r,array_c,matrix ) 
% 求堆栈内元素的邻域内的最大值及其坐标
    max=-1;
   [size_rr,size_cc]=size(matrix); 
   flag_matrix=zeros(size_rr,size_cc);   % 标志栈内元素
    
    for index=1:length(array_r)
        index_r=array_r(index);
        index_c=array_c(index);   
        flag_matrix(index_r,index_c)=1;       
    end
    
    for index=1:length(array_r)
        index_r=array_r(index);
        index_c=array_c(index);         
          
        if index_r-1>1 & matrix((index_r-1),index_c)>max & flag_matrix(index_r-1,index_c)==0 % 上
            max=matrix(index_r-1,index_c);
            max_r=index_r-1;
            max_c=index_c;          
        end    
        
        if index_r+1<size_rr & matrix(index_r+1,index_c)>max & flag_matrix(index_r+1,index_c)==0 % 下
            max=matrix(index_r+1,index_c);
            max_r=index_r+1;
            max_c=index_c;     
        end     
        
        if index_c+1<size_cc & matrix(index_r,index_c+1)>max & flag_matrix(index_r,index_c+1)==0
            max=matrix(index_r,index_c+1);
            max_r=index_r;
            max_c=index_c+1;
        end   
        
        if index_c-1>1 & matrix(index_r,index_c-1)>max & flag_matrix(index_r,index_c-1)==0
            max=matrix(index_r,index_c-1);
            max_r=index_r;
            max_c=index_c-1;
        end
    end
% －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
    
function [min_r,min_c]=seek_min_elem( array_r,array_c,matrix ) 
% 求堆栈内元素的领域内的最小值及其坐标
   min=256;
   [size_rr,size_cc]=size(matrix); 
   flag_matrix=zeros(size_rr,size_cc);   % 标志栈内元素
    
    for index=1:length(array_r)
        index_r=array_r(index);
        index_c=array_c(index);   
        flag_matrix(index_r,index_c)=1;       
    end
    
    for index=1:length(array_r)
        index_r=array_r(index);
        index_c=array_c(index);         
          
        if index_r-1>1 & matrix((index_r-1),index_c)<min & flag_matrix(index_r-1,index_c)==0 % 上
            min=matrix(index_r-1,index_c);
            min_r=index_r-1;
            min_c=index_c;          
        end    
        
        if index_r+1<size_rr & matrix(index_r+1,index_c)<min & flag_matrix(index_r+1,index_c)==0 % 下
            min=matrix(index_r+1,index_c);
            min_r=index_r+1;
            min_c=index_c;     
        end     
        
        if index_c+1<size_cc & matrix(index_r,index_c+1)<min & flag_matrix(index_r,index_c+1)==0
            min=matrix(index_r,index_c+1);
            min_r=index_r;
            min_c=index_c+1;
        end   
        
        if index_c-1>1 & matrix(index_r,index_c-1)<min & flag_matrix(index_r,index_c-1)==0
            min=matrix(index_r,index_c-1);
            min_r=index_r;
            min_c=index_c-1;
        end
    end
% －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function matrix_ex=matrix_expand(matrix)
% 对矩阵进行延拓   
   [size_r,size_c]=size(matrix);
   matrix_ex=zeros(size_r+2,size_c+2);
   
   % 中
   for row=1:size_r
       for col=1:size_c
           matrix_ex(row+1,col+1)=matrix(row,col);
       end
  end
  
  % 左，右
  for row=1:size_r
      matrix_ex(row+1,1)=matrix(row,1);
      matrix_ex(row+1,size_c+2)=matrix(row,size_c);
  end
  
  % 上，下
  for col=1:size_c
      matrix_ex(1,col+1)=matrix(1,col);
      matrix_ex(size_r+2,col+1)=matrix(size_r,col);
  end

  matrix_ex(1,1)=matrix(1,1);
  matrix_ex(1,size_c+2)=matrix(1,size_c);

  matrix_ex(size_r+2,1)=matrix(size_r,1);
  matrix_ex(size_r+2,size_c+2)=matrix(size_r,size_c);
  % －－－－－－－－－－－－－－－－－－－－－－－－－－