检测光纤的同心性

这是matlab自带的一个例程。

问题描述：光纤横截面包括core和cladding两部分，也就是“中心圆”和“漆包线”
希望两者是同心圆，或者说，希望两者圆心差的不要太离谱
原理：先将图像转换为二值图像，得到只有漆包线（圆环）是黑色、其他都是白色的图像
根据漆包线区域计算出图像重心
根据重心得到cladding边界上一个点，进而得到整个边界上的点集
由边界点击得到圆心。使用了最小二乘法来拟合圆形（算法有点复杂。。没看过）
计算圆心和重心之间的偏差，若小于阈值则认为光纤合格

%初始化
NumPts = 250; %边界像素点的最大数目
tolerance = 1; %容忍距离
NumTimes = 8;

%创建一个视频读写系统对象，读取视频文件
hmfr = vision.VideoFileReader('vipconcentricity.avi', 'PlayCount', NumTimes);
%转换RGB到灰度空间。为什么要创建一个转换对象？直接转不好吗？
hcsc = vision.ColorSpaceConverter('Conversion', 'RGB to intensity');

%创建一个光斑分析(BlobAnalysis)系统对象，用来寻找视频中圆形光斑的重心
hblob = vision.BlobAnalysis('AreaOutputPort', false,...
    'BoundingBoxOutputPort', false,...
    'OutputDataType', 'single',...
    'MaximumCount', 1);

%创建一个边界跟踪系统对象，用来寻找包层的中心
%BoundaryTracer函数是在二值图像上找到object的边界，也就是非零的边界点集合
%目测是用dfs来实现的
htracebound = vision.BoundaryTracer(...
    'MaximumPixelCount', NumPts,...
    'NoBoundaryAction', 'Fill with last point found');

%创建两个视频播放对象，用来显示视频的输入和输出
hVideo1 = vision.VideoPlayer('Name', 'Original');
hVideo_gray = vision.VideoPlayer('Name', 'Gray');
hVideo_binary = vision.VideoPlayer('Name', 'Binary');
hVideo2 = vision.VideoPlayer('Name', 'Results');
hVideo2.Position(1) = hVideo1.Position(1) + 450;

%处理视频流
count = 1;
while(count<NumTimes)
    if isDone(hmfr)
        count = count + 1;
    end
    I = step(hmfr);
    image = step(hcsc, I); %转化图片为灰度图
    step(hVideo_gray, image);
    
    BW = image < 0.5; %转化图片为二值图
    step(hVideo_binary, BW);
    centroid = step(hblob, BW); %计算光斑重心
    
    Idx = floor(centroid(1)); %重心x坐标
    max_idx = find(BW(:,Idx), 1); %重心x坐标所在的列中，第一个非零元素的索引.
    %因为图片已经是二值化的，（而且是理想图片），所以第一个非零元素一定是cladding外围圆周上、垂直过圆心的上侧的点
    
    StartPts = [Idx, single(max_idx)]; %计算起始点，也就是cladding外围圆周上垂直过圆心的上侧的点
    
    %寻找外部包层的边界像素
    Pts = step(htracebound, BW, StartPts);
    
    %以下是难点部分：根据边界点，用最小二乘法拟合圆形，得到圆心和半径
    %原理见Paper：http://www.emis.de/journals/BBMS/Bulletin/sup962/gander.pdf
    Row_bound = Pts(:, 1);
    Col_bound = Pts(:, 2);
    t = [Row_bound Col_bound ones(size(Pts, 1), 1)];
    
    X = pinv(t);
    X1 = Row_bound.^2 + Col_bound.^2;
    x2 = X*(-X1);
    
    radius = sqrt((-0.5*x2(1)).^2 + (-0.5*x2(2)).^2 - x2(3));
    center = [(-0.5*x2(1)), (-0.5*x2(2))];
    %拟合计算完毕，得到了圆心和半径。下面，用圆心和重心比较
    %如果两者之间误差在允许范围内，则光纤是可以接受的（合格）

    dist = sqrt(sum(centroid - center).^2);
    disp(dist)

    % 绘制圆形，并标注出圆心和重心
    y1 = insertMarker(I, centroid, '+', 'Color', 'red');
    y2 = insertShape(y1, 'Circle', [center, radius],'Color', 'cyan');
    y3 = insertMarker(y2, center, '*', 'Color', 'green');

    % 插入文本
    qualitySet =  {'Bad', 'Good'};
    textIdx = (dist <= tolerance) + 1;
    textQuality = ['Concentricity: ' qualitySet{textIdx}];
    textDist = sprintf('Distance in pixels: %d', uint8(dist));
    textAll = sprintf([textQuality '
' textDist]);
    image_out = insertText(y3, [1 1], textAll, 'FontSize', 14);
    % 显示结果
    step(hVideo2, image_out);
    step(hVideo1, I);
end

release(hmfr);