快速径向对称 只检测暗对称 代码（2）

<span style="font-size:18px;">%RST  径向对称变换 
% function [FilterImage,Image_out] = RST(I);


I = imread('lena.jpg');
% I=zeros(80,80) ;
% for x=1:1:80
%     for y=1:1:80
%         I(y,x) =255;
%     end;
% end;
% 
% r2 =8*8 ;
% for x=1:1:80
%     for y=1:1:80
%         tmp1 =(x-40)*(x-40)+(y-40)*(y-40) ;
%         if(tmp1 <= r2)
%             I(y,x) = 0 ;
%         end;
%     end;
% end;

%%
% 计算图像梯度
[h w dim]  = size(I);
if(dim == 3)
    I = rgb2gray(I);
end;


%如果均值小于100就开始做log变换
mean_v = mean(I(:));
if(mean_v <60)
I_log =log(double(1+I));
Min =min(I_log(:));
Max =max(I_log(:));
I_log = (I_log-Min)*255/(Max-Min);
%I =    I_log ; 
end;


 [FX,FY] = gradient(double(I)) ;
%  FX = textread('D:\RST\XG.txt');
%  [t1 t2]= size(FX);
%  FX =FX(:,1:t2-1);
%  FY = textread('D:\RST\YG.txt');
%   [t1 t2]= size(FY);
%  FY =FY(:,1:t2-1);
%计算梯度方向图和幅度图
% theta = atan2( FY,FX);
Mag = sqrt(FY.^2 +FX.^2);%图像的幅度图



O_map_pos=zeros(h,w) ;%方向投影图   正反向 就是梯度方向
O_map_neg=zeros(h,w) ;
M_map_pos=zeros(h,w) ;   %图像的幅度图
M_map_neg=zeros(h,w) ;

M_map_Total =zeros(h,w) ;

S_map=zeros(h,w) ;




%输出高斯模板

Moban = fspecial('gaussian',[9 9],2); %高斯模板对图像进行处理

for r =3:1:6
    %计算RST   正反方向的  方向投影和幅度图
  O_map_pos=zeros(h,w) ;
  M_map_pos=zeros(h,w) ;
  O_map_neg=zeros(h,w) ;
  M_map_neg=zeros(h,w) ;
  a =2 ;
  %对图像进行处理
  for y=(r+1):1:(h-r-1)
      for x=(r+1):1:(w-r-1)
          
          if(Mag(y,x) == 0) continue; end;
          %就是原文中的正向影响像素点坐标偏差值
          tmp_y = FY(y,x)*r/Mag(y,x) ;%在x和y方向上的r*单位增量  (g(p)/||g（p）||)*n
          tmp_x = FX(y,x)*r/Mag(y,x) ;

     %     if( ((tmp_y>=0) &&(tmp_x>=0) )||((tmp_y<0) &&(tmp_x<0) )) end;
          y_to =floor(y - tmp_y  +0.5) ;%在当前像素点的x和y方向上的增量
          x_to =floor(x - tmp_x  +0.5) ;
%          if(I(y,x) >200) continue; end;
%            if( ((tmp_y>=0) &&(tmp_x>=0) )||((tmp_y<0) &&(tmp_x<0) )) 
%               O_map_pos(y_to,x_to) = O_map_pos(y_to,x_to)+1 ;
%                M_map_pos(y_to,x_to) = M_map_pos(y_to,x_to) +Mag(y,x);
%            else
               O_map_neg(y_to,x_to) = O_map_neg(y_to,x_to)+1 ;
               M_map_neg(y_to,x_to) = M_map_neg(y_to,x_to) +Mag(y,x);
%         end;
      end;  %x
  end; %y
end

 
 K_1 =max(abs(O_map_neg(:))) ;
 K_2 =max(abs(M_map_neg(:))) ;  
 K_1=1 ;
 K_2=1 ;
 S_map = (((O_map_neg/K_1).^2 ).* M_map_neg/K_2) ;
 M_map_Total = M_map_Total + S_map ;
  
% end; %r

S_map_Filter = conv2( M_map_Total,Moban ,'same');
Image_out =M_map_Total ;
FilterImage =S_map_Filter ;
imshow(FilterImage ,[]); 
imshow(Image_out,[]) 





</span>