sanweiguandao.m

先把图形改为.bmp形式

%计算100幅图像的管道半径和管道的中轴坐标。
x_Initial_Loc = 0;  %存放前一幅图像管道的中轴横坐标，作为计算本幅图像的管道的中轴横坐标的初始点
y_Initial_Loc = 0;  %存放前一幅图像管道的中轴纵坐标，作为计算本幅图像的管道的中轴纵坐标的初始点

Distance_R_Mat = [];%存放100幅图像管道的半径
x_Center_Loc = [];  %存放100幅图像管道的中轴横坐标
y_Center_Loc = [];  %存放100幅图像管道的中轴纵坐标

for image_no = 0:99
    Edge_x_axis = []; %坐标系调整为原点在图像中心后的边缘点的X坐标
    Edge_y_axis = []; %坐标系调整为原点在图像中心后的边缘点的Y坐标
    Edge_x = [];      %坐标系调整前的边缘点的X坐标
    Edge_y = [];      %坐标系调整前的边缘点的Y坐标
    
    FileName = [];
    FileName = [ FileName, num2str( image_no ), '.bmp' ];

    Data_Bmp = imread( FileName );      %取出图像数据。
    [EdgeDate,Threshold] = edge( uint8(Data_Bmp) ); %取出图像的边缘数据和提取边缘的阈值。
     
    %取出每一幅图像所有边界点的坐标。
    for i = 1:512
        for j = 1:512
            if ( EdgeDate(i, j) ) >= Threshold
                Edge_x = [ Edge_x, i ];
                Edge_y = [ Edge_y, j ];
            end
        end    
    end 
    %图像数据矩阵坐标系是以左上角为原点，右为Y轴正向，下为X轴正向。
    %把图像数据矩阵坐标系调整为以图像的中心为坐标原点，X轴方向修改为竖直向上，Y轴方向不变。
    Edge_x_Axis = -Edge_x + 256;
    Edge_y_Axis =  Edge_y - 256;
    
    if strcmp( FileName, '0.bmp' )
        %求出第一幅图像的边界点的横纵坐标的平均值，以此为管道圆心的初始值。
        x_Center = mean( Edge_x_Axis );%此处也可以考虑使用x或y的最小值和最大值的平均
        y_Center = mean( Edge_y_Axis );         
    else %从第二次循环开始进入此处
        %其他的图像的管道中心都以前一幅图像的管道中心为初始值。
        x_Center = x_Initial_Loc;
        y_Center = y_Initial_Loc;
    end
    
    %取管道中心附近的点为初始点，求出管道半径和管道的中轴线的坐标。具体为:
    %取管道中心附近点周围的n个点，求出这些点到边界距离的最小值，
    %并以其最大值为管道半径，其对应点为管道中心坐标。
    if image_no < 10
        x_Center_Temp = [ x_Center-3, x_Center-2, x_Center-1, x_Center, x_Center+1, x_Center+2, x_Center+3 ];
        y_Center_Temp = [ y_Center-3, y_Center-2, y_Center-1, y_Center, y_Center+1, y_Center+2, y_Center+3 ];
    elseif image_no < 30
        x_Center_Temp = [ x_Center-5, x_Center-4, x_Center-3, x_Center-2, x_Center-1, x_Center, x_Center+1, x_Center+2, x_Center+3, x_Center+4, x_Center+5 ];
        y_Center_Temp = [ y_Center-5, y_Center-4, y_Center-3, y_Center-2, y_Center-1, y_Center, y_Center+1, y_Center+2, y_Center+3, y_Center+4, y_Center+5 ];
    elseif image_no < 50
        x_Center_Temp = [ x_Center-7, x_Center-6, x_Center-5, x_Center-4, x_Center-3, x_Center-2, x_Center-1, x_Center, x_Center+1, x_Center+2, x_Center+3, x_Center+4, x_Center+5, x_Center+6, x_Center+7 ];
        y_Center_Temp = [ y_Center-7, y_Center-6, y_Center-5, y_Center-4, y_Center-3, y_Center-2, y_Center-1, y_Center, y_Center+1, y_Center+2, y_Center+3, y_Center+4, y_Center+5, y_Center+6, y_Center+7 ];
    elseif image_no < 70
        x_Center_Temp = [ x_Center-9, x_Center-8, x_Center-7, x_Center-6, x_Center-5, x_Center-4, x_Center-3, x_Center-2, x_Center-1, x_Center, x_Center+1, x_Center+2, x_Center+3, x_Center+4, x_Center+5, x_Center+6, x_Center+7, x_Center+8, x_Center+9 ];
        y_Center_Temp = [ y_Center-9, y_Center-8, y_Center-7, y_Center-6, y_Center-5, y_Center-4, y_Center-3, y_Center-2, y_Center-1, y_Center, y_Center+1, y_Center+2, y_Center+3, y_Center+4, y_Center+5, y_Center+6, y_Center+7, y_Center+8, y_Center+9 ];
    else
        x_Center_Temp = [ x_Center-11, x_Center-10, x_Center-9, x_Center-8, x_Center-7, x_Center-6, x_Center-5, x_Center-4, x_Center-3, x_Center-2, x_Center-1, x_Center, x_Center+1, x_Center+2, x_Center+3, x_Center+4, x_Center+5, x_Center+6, x_Center+7, x_Center+8, x_Center+9, x_Center+10, x_Center+11 ];
        y_Center_Temp = [ y_Center-11, y_Center-10, y_Center-9, y_Center-8, y_Center-7, y_Center-6, y_Center-5, y_Center-4, y_Center-3, y_Center-2, y_Center-1, y_Center, y_Center+1, y_Center+2, y_Center+3, y_Center+4, y_Center+5, y_Center+6, y_Center+7, y_Center+8, y_Center+9, y_Center+10, y_Center+11 ];
    end
    
    Center_Point_Sum = length( x_Center_Temp );%取临时中心点的点个数，以供后面的循环使用
    Edge_Point_Sum   = length( Edge_x_Axis );  %取边界点的点个数，以供后面的循环使用
    Dist_Center_Mat  = [];                     %边界点到所有中心点的距离
    Distance_Min_Center = [];
    for i = 1:Center_Point_Sum
        for j = 1:Center_Point_Sum
            for k = 1:Edge_Point_Sum
                Dist_Center_Mat(k) = sqrt( (x_Center_Temp(i)-Edge_x_Axis(k))^2 + (y_Center_Temp(j)-Edge_y_Axis(k))^2 );
            end
            Distance_Min_Center = [ Distance_Min_Center, min( Dist_Center_Mat ) ];%所有距离中的最小值是半径
        end    
    end
    %此处产生的误差比较大，可以考虑使用前n个最大距离的平均值对应的点为管道半径。
    [Distance_R,Ind] = max( Distance_Min_Center );%获得一个图层的半径和对应的位置
   
    % 图像0.bmp的管道中轴作为所有图像管道的最初的初始值，其他图像都是以前一幅图像的管道中轴坐标为初始值。
    if strcmp( FileName, '0.bmp' )
        x_Initial_Loc = x_Center;
        y_Initial_Loc = y_Center;         
    else 
        %取出管道中心点作为管道中轴的点的坐标
        if mod( Ind, Center_Point_Sum ) == 0
            x_Initial_Loc = x_Center_Temp( floor( Ind/Center_Point_Sum ) );
            y_Initial_Loc = y_Center_Temp( Center_Point_Sum );
        else
            x_Initial_Loc = x_Center_Temp( floor( Ind/Center_Point_Sum ) + 1 );%行
            y_Initial_Loc = y_Center_Temp( mod( Ind, Center_Point_Sum ) );     %列
        end        
    end
    %保存100幅图像中的管道的半径和中轴坐标
    Distance_R_Mat = [ Distance_R_Mat, Distance_R ];
    x_Center_Loc = [ x_Center_Loc, x_Initial_Loc ];
    y_Center_Loc = [ y_Center_Loc, y_Initial_Loc ];
end
GuanDao_Ridus = mean( Distance_R_Mat )