%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     水果大小检测    2013年10月10日  zhengrui  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%读取水果图像
IClr = imread('大小样本\sideviewclr.png', 'png');

%转换为灰度图
IGry = rgb2gray(IClr);

%用中值滤波去除椒盐噪声
IFltM = medfilt2(IGry,[3 3]);

%创建高斯滤波模板
H = fspecial('gaussian',[3,3],0.5);

%用高斯滤波器对原图像进行平滑
IFltG = filter2(H,IFltM);

%获得Sobel边缘图像
ISbl = edge(IFltG,'sobel'); 
figure(1);
imshow(ISbl);

%估算形心位置
[m,n]=size(ISbl);
X0=0;
Y0=0;
COUNT=0;
for i=3:m-3
    for j=3:n-2
        if(ISbl(i,j))
            X0=X0+j;
            Y0=Y0+i;
            COUNT = COUNT+1;
        end
    end 
end
X0 = round(X0/COUNT);
Y0 = round(Y0/COUNT);

%获取精确边缘
IfineEdge = ISbl;
IfineEdge(:,:)=0;
R1=115;R2=55;
RVector      = zeros(1,366);
Rsearchvalue = zeros(1,60);
tmp = [-1,0,1,-1,1,-1,0,1;-1,-1,-1,0,0,1,1,1];

%从（X0，Y0）向外向内根据变异参数搜索边缘点
for di=1:1:360
    deg = (di-1)*pi/180;
    for R=R1:-1:R2+1
        %disp(R);
        X=R*cos(deg)+X0;
        Y=R*sin(deg)+Y0;
        GMean=0;
        for i=1:8
            GValue=double(IGry(int32(Y+tmp(1,i)),int32(X+tmp(2,i))))-double(IGry(int32(Y),int32(X)));
            GMean=GMean+double(IGry(int32(Y+tmp(1,i)),int32(X+tmp(2,i))))*double(IGry(int32(Y+tmp(1,i)),int32(X+tmp(2,i))));
            Rsearchvalue(1,R-R2)=Rsearchvalue(1,R-R2)+GValue*GValue;
        end
        Rsearchvalue(1,R-R2)=Rsearchvalue(1,R-R2)/(GMean+double(IGry(int32(Y),int32(X)))*double(IGry(int32(Y),int32(X))));
    end  
    [Y,I]=max(Rsearchvalue(:));
    %显示原始搜索边缘
    RVector(di+3)=R2+I;
end
RVector(1)   = RVector(361);
RVector(2)   = RVector(362);
RVector(3)   = RVector(363);
RVector(364) = RVector(4);
RVector(365) = RVector(5);
RVector(366) = RVector(6);
%中值滤波+均值滤波去除噪声点
VectTmp = medfilt1(RVector,7);
RVector = meanfilter(VectTmp, 7);
%显示精确的边缘
for di=1:1:360
    IfineEdge(int32(RVector(di+3)*sin((di-1)*pi/180)+Y0), int32(RVector(di+3)*cos((di-1)*pi/180)+X0)) = 1;
end
%在灰度图像上显示找到的边缘
for i=3:m-3
    for j=3:n-2
    if(IfineEdge(i,j)==1)   
        IClr(i,j,1)=255;
        IClr(i,j,2)=255;
        IClr(i,j,3)=255;
    end 
    end 
end     
figure(2);
imshow(IClr);
figure(3);
imshow(IfineEdge);

%求取精确形心
sigmanmx = 0;
sigmanmy = 0;
sigmadem = 0;
Xk1=RVector(3)*cos((359)*pi/180)+X0;
Yk1=RVector(3)*sin((359)*pi/180)+Y0;
EdgeChain = zeros(2,360);
for di = 1:1:360
    Xk = RVector(di+3)*cos((di-1)*pi/180)+X0;
    Yk = RVector(di+3)*sin((di-1)*pi/180)+Y0;
    sigmanmx = sigmanmx+Yk*(Xk*Xk-Xk1*Xk1)-Xk*Xk*(Yk-Yk1);
    sigmanmy = sigmanmy+Yk*Yk*(Xk-Xk1)-Xk*(Yk*Yk-Yk1*Yk1);
    sigmadem = sigmadem+Yk*(Xk-Xk1)-Xk*(Yk-Yk1);
    Xk1=Xk;
    Yk1=Yk;
    EdgeChain(1,di)=Xk;
    EdgeChain(2,di)=Yk;
end 
X0fine = sigmanmx/sigmadem/2;
Y0fine = sigmanmy/sigmadem/2;
%寻找轴向
RVectorSample  = zeros(1,55);
for di = 6:1:50
    deltax = EdgeChain(1,(di-5)*8)-X0fine;
    deltay = EdgeChain(2,(di-5)*8)-X0fine;
    RVectorSample(di)= sqrt(deltax*deltax+deltay*deltay);
end
RVectorSample(1)=RVectorSample(46);
RVectorSample(2)=RVectorSample(47);
RVectorSample(3)=RVectorSample(48);
RVectorSample(4)=RVectorSample(49);
RVectorSample(5)=RVectorSample(50);
RVectorSample(51)=RVectorSample(6);
RVectorSample(52)=RVectorSample(7);
RVectorSample(53)=RVectorSample(8);
RVectorSample(54)=RVectorSample(9);
RVectorSample(55)=RVectorSample(10);
g=zeros(1,45);
for di = 6:1:50
    for Ri=1:5
        g(di-5)=g(di-5)+abs(RVectorSample(di-Ri)-RVectorSample(di+Ri));
    end 
end
[YR,IR]=min(g(:));
deg1=IR*8-90;
if(deg1<0)
    deg1=deg1+360;
end
deg2=IR*8+90;
if(deg2>360)
    deg2=deg2-360;
end
%在图像上绘出轴向线
figure(4);
%imshow(IfineEdge);
imshow(IClr);
axis on;hold on;
plot([X0fine,EdgeChain(1,IR*8)],[Y0fine,EdgeChain(2,IR*8)],'b','linewidth',2);
plot([EdgeChain(1,deg1),EdgeChain(1,int32(deg2))],[EdgeChain(2,int32(deg1)),EdgeChain(2,int32(deg2))],'b','linewidth',2);
axis off;


       


