光条中心提取代码

        本篇文章主要给出极大值法，灰度重心法及steger方法光条中心提取的代码，代码及方法未做任何优化处理，鉴于自身水平所限，如有错误，欢迎批评指正。（欢迎进Q群交流：874653199）

clc;
clear;

srcImg=imread('laser.bmp');

if(numel(size(srcImg))>2)
    srcImg=rgb2gray(srcImg);
end

thresh= graythresh(srcImg)*255-50;

coor1=grayCentroidFunction(srcImg,thresh,0) ;
coor2=maxValueFunction(srcImg,thresh,0) ;
coor3=stegerFunction(srcImg);

figure(1);
imshow(uint8(srcImg));
hold on
plot(coor1(:,1),coor1(:,2),'r.','MarkerSize',1)
title('极大值法','color','r')
hold off

figure(2);
imshow(uint8(srcImg));
hold on
plot(coor2(:,1),coor2(:,2),'r.','MarkerSize',1)
title('灰度重心法','color','r')
hold off

figure(3);
imshow(uint8(srcImg));
hold on
plot(coor3(:,1),coor3(:,2),'r.','MarkerSize',1)
title('steger方法','color','r')
hold off


function centerCoor=maxValueFunction(srcImg,thresh,flag) 
  [row,col]=size(srcImg);
  centerCoor=zeros(row*col,2);
  sigma=3;
  kernel_size=5;
  kernel=fspecial('gaussian',kernel_size,sigma);
  srcImg=imfilter(srcImg,kernel);
  srcImg=double(srcImg); 
  if(flag==0)%水平
      max=0;
      y=0;
      count=1;
      for i=1:row
          for j=1:col
              current=srcImg(i,j);
              if(current>max)
                  max=current;
                  y=j;
              end
          end
          if(max>=thresh)  
             centerCoor(count,:)=[y,i] ;
             count=count+1;
          end
          max=0;
          y=0;
      end
  else%垂直
      max=0;
      y=0;
      count=1;
      for i=1:col
          for j=1:row
              current=srcImg(j,i);
              if(current>max)
                  max=current;
                  y=j;
              end
          end
          if(max>=thresh)  
             centerCoor(count,:)=[i,y] ;
             count=count+1;
          end
          max=0;
          y=0;
      end
      
  end    
centerCoor(centerCoor(:,1)==0,:)=[];
end

function centerCoor=grayCentroidFunction(srcImg,thresh,flag)
[row,col]=size(srcImg);
  centerCoor=zeros(row*col,2);
  sigma=3;
  kernel_size=5;
  kernel=fspecial('gaussian',kernel_size,sigma);
  srcImg=imfilter(srcImg,kernel);
  srcImg=double(srcImg); 
  if(flag==0)%水平
      count=1;
      for i=1:col
        value=zeros(row,1);
        coor=zeros(row,1);
        for j=1:row
            current=srcImg(j,i);
            if(current>thresh)
               value(j)=current;
               coor(j)=j;
            end
        end
        coor(coor(:,1)==0,:)=[];
        value(value(:,1)==0,:)=[];
        num=size(value,1);
        if(num>=3)
            sum_valuecoor=sum(value.*coor);
            sum_value=sum(value);
            x=sum_valuecoor/sum_value;
            centerCoor(count,:)=[i,x];
            count=count+1;
        end        
      end      
      centerCoor(centerCoor(:,1)==0,:)=[];      
  else  %垂直
      count=1;
      for i=1:row
        value=zeros(col,1);
        coor=zeros(col,1);
        for j=1:col
            current=srcImg(i,j);
            if(current>thresh)
               value(j)=current;
               coor(j)=j;
            end
        end
        coor(coor(:,1)==0,:)=[];
        value(value(:,1)==0,:)=[];
        num=size(value,1);
        if(num>=3)
            sum_valuecoor=sum(value.*coor);
            sum_value=sum(value);
            x=sum_valuecoor/sum_value;
            centerCoor(count,:)=[x,i];
            count=count+1;
        end
        
      end    
      centerCoor(centerCoor(:,1)==0,:)=[]; 
  end
end

function centerCoor=stegerFunction(srcImg)
sigma=3;
thresh= graythresh(srcImg);%otsu  
srcImg=double(srcImg);
[m,n]=size(srcImg);
ky=[-1,1];
kx=[-1;1];
kyy=[1,-2,1];
kxx=[1;-2;1];
kxy=[1,-1;-1,1];
gausFilter = fspecial('gaussian',17,sigma);   %高斯滤波
dstImg=imfilter(srcImg,gausFilter,'replicate');  
dx=imfilter(dstImg,kx);
dy=imfilter(dstImg,ky);
dxx=imfilter(dstImg,kxx);
dyy=imfilter(dstImg,kyy);
dxy=imfilter(dstImg,kxy);
hessian=zeros(2,2);
points=zeros(m*n,2);
for i=1:m
    for j=1:n
        if(srcImg(i,j)/255>thresh)
            hessian(1,1)=dxx(i,j);
            hessian(1,2)=dxy(i,j);
            hessian(2,1)=dxy(i,j);
            hessian(2,2)=dyy(i,j);
            [eigenvectors,eigenvalues]=eig(hessian);           
            if(abs(eigenvalues(1,1))>= abs(eigenvalues(2,2)))
            nx=eigenvectors(1,1);
            ny=eigenvectors(2,1);
            fmax_dist=eigenvalues(1,1);
            else
            nx=eigenvectors(1,2);
            ny=eigenvectors(2,2);
            fmax_dist=eigenvalues(2,2);                
            end            
            t=-(nx*dx(i,j)+ny*dy(i,j))/(nx*nx*dxx(i,j)+2 * nx*ny*dxy(i,j)+ny*ny*dyy(i,j));          
            if(abs(t*nx) <= 0.5 && abs(t*ny) <= 0.5)
                points((i-1)*m+j,:)=[ j+t*ny,i+t*nx];
            end                 
        end      
    end
end

index = find(points(:,1)==0);
points(index,:) = [];
index = find(points(:,2)==0);
points(index,:) = [];

centerCoor=points;

end

光条中心提取结果：