设计任务

这是一幅表盘图像，试采用图像处理和分析技术，设计适当的算法和程序，计算出表盘上的指示时间（要求精确到秒）。请按统一要求写出算法原理、设计流程，并完成测试分析等报告内容。

算法设计

解题思路：使用霍夫变换检测出表针对应的直线，通过表针的长度区分出时针、分针、和秒针。通过各直线间的角度确定指示时间。

设计方案说明：

1. 为了使霍夫变换只检测出三条直线，把原图像转化为灰度图，并使用灰度范围提取和开闭运算得到只有表针和表盘中心的图像；
2. 然后使用霍夫变换。得到了直线两端的端点。为了确定哪个点是指针末端端点，计算两个点的纵坐标到表盘中心纵坐标（311）的距离，距离大的是指针末端的端点。
3. 使用霍夫变换后得到的直线不是在表针中心的，故使用三个大小和原图像相等的逻辑值矩阵掩膜（除了感兴趣区其他值为0），其中分别包含以三条直线末端点（表针末端）为中心，大小为13*13（使用imtool函数，发现指针的宽度小于6个像素，因此使用13*13的感兴趣区能包含完整的指针末端）的感兴趣区（值为1），与原图像进行乘法运算，提取目标区域；
4. 对于每个目标区域，行遍历的方向是从没有白色行到有白色行，逐行向左遍历，找到碰到值从0变为1的点，以及边缘对应的斜率值突变（由正变负或由负变正）的点，记录下来，而后又逐行向右遍历，记录下点。若在向左或向右遍历中得到两个点，那么这两个点就是指针末端边的相聚最远的边缘点。
5. 对于指针的端边，由于其纵坐标随横坐标的变化是单调的，可通过其上所有点中横坐标的最大值加上横坐标的最小值除以2，再向下取整，得到指末端针中心对应最接近的横坐标点；纵坐标的算法同理。进而得到指针末端中心坐标。
6. 还需要确定表盘中心的坐标。使用imtool函数和scatter函数，得到表盘中心的坐标为(329,311)。
7. 用指针末端中心坐标为尾端，表盘中心坐标为首端做一个向量，把它变成单位向量，求出它与和指向12点指针同方向的单位向量（1,0）的夹角。通过指针末端中心坐标和表盘中心位置的相对关系与夹角，即可算出表盘上显示的时间。

关键算法的设计原理：

形态学开闭运算：消除白色小边缘。

乘法运算：0乘以任何数都是0。

从左到右遍历和从右到左遍历可以筛选出指针末端相聚最远的边缘点原因：设置条件为只有当边缘斜率值突变（由正变负或由负变正）才可以统计此点。从两边都遍历，其中有一边遍历斜率值为0，即被掩模挡住不能统计，那就会有一个点不能被统计，应该选择另外一边统计到的点（可看图5理解）。

说明7：有两种时钟，第一种时钟是时针、分针、秒针不停转的；第二种时钟是等一种针转过一圈，另一种（或两种）针才会转一个格。由图像可知，分针在两个时针刻度之间，因此此时钟是第一种时钟。因此有时针：分针：秒针的角速度关系为1:60:3600。由此出现有两种解法，第一种解法是通过时针和经过数字12的刻度和表盘中心的直线间的角度，直接算出现在的时间，但由于识别直线存在系统误差，这种方法读数不准确；因此选用第二种解法。

其他原理：观察+经验+高等数学。

算法步骤：

灰度图→灰度范围提取→开闭运算→霍夫变换→确定指针末端区域→掩模提取→求指针末端距离最远两点→指针末端中心坐标→算指针方向的单位向量与12点种方向的夹角→求时间

程序设计

算法名称：灰度特征提取、灰度范围提取、二值化图像、形态学开闭运算、掩模、乘法运算等。

工具函数：霍夫变换：hough、houghpeaks、houghlines等。

设计分析：完成任务。

测试分析

 图1

图1为原图→灰度图→提取灰度范围为[20,120)为亮部的二值图→开闭运算后图。可以看到，经过提取灰度范围一步，表面上数字和刻度剩下边缘；经过开闭运算后，边缘去除。

 图2

图2为“开闭运算后图”经过MATLAB内置霍夫变换后得到的线段和端点图。可以发现，线段上对应指针的末端的端点并不与指针末端中心重合。

 图3

图3为使用掩模提取后图。

图4

图4为感兴趣区放大图。

图5

（已手动放大）在图5中，红*为从左到右遍历找到的边缘点，绿*为从右到左遍历找到的边缘点。可以印证若一种遍历方法能找到两个点，那这两个点将是指针末端相聚最远的边缘点。

 图6

（已手动放大）图6，红*为优化找到的指针末端中点，绿*为仅通过霍夫变换得到直线末端点。对比发现，优化有效。

 图8中的蓝线为求出的指针中点和表盘中心的连线，基本与指针重合。绿*为仅通过霍夫变换得到直线末端点，有一定偏移。

结果：时钟上显示10时7分39秒 

特点与优势：算法有效，综合运用高数、经验、观察、图像处理知识

存在的问题和不足：1.代码过长，应使用函数式编程 2.秒针的角度偏大，需优化

MATLAB代码

clc,clear
dial_center=[329,311]; % 表盘中心

%% step1 -- 读入图像
I=imread('ipa01.jpg');
%% step2 -- 转化为灰度图Gray
Gray=rgb2gray(I);
sz=size(Gray);

%% step3 -- 提取Gray中灰度范围为[20,120)作为二值图的亮部，得到图像A
A=(Gray<120)-(Gray<20);

%% step4 -- 开闭运算，去掉图像A中的数字轮廓、钟表刻度轮廓
se=strel('disk',2);
B=imopen(A,se);
C=imclose(B,se);

%% step5 -- 霍夫变换，得到三条与指针中线不重合的直线
[H,theta,rho] = hough(C);
peaks = houghpeaks(H,3);
lines = houghlines(C,theta,rho,peaks);

%% step6 -- 确定指针末端端点
hms_point=zeros(3,3);% 指针末端端点
for i=1:3
    if abs(dial_center(2)-lines(i).point1(2))<abs(dial_center(2)-lines(i).point2(2))
        hms_point(i,1)=lines(i).point2(1);
        hms_point(i,2)=lines(i).point2(2);
    else
        hms_point(i,1)=lines(i).point1(1);
        hms_point(i,2)=lines(i).point1(2);
    end
    hms_point(i,3)=norm(lines(i).point1 - lines(i).point2);% 直线长度
end

%% step7 -- 通过长度排序，确定三条直线对应哪根指针
% 时针最短，分针最长，秒针长度居中
hms_point=sortrows(hms_point,3);
hms_point=[hms_point(1,:);hms_point(3,:);...
    hms_point(2,:)];

%% step7 -- 以各直线末端点为中心，造三种掩模
hms_mask=zeros(sz(1),sz(2),3); % ROI为每个指针末端的掩模
len=6; % （掩模边长-中心点（1像素））/2的值
for i=1:3
    x=hms_point(i,1);
    y=hms_point(i,2);
    hms_mask(y-len:y+len,x-len:x+len,i)=C(y-len:y+len,x-len:x+len);
end

%% step8 -- 求边缘点，当出现斜率突变时
endpoints=zeros(2,2,3);
lpoints=zeros(2,2,3);
rpoints=zeros(2,2,3);
for i=1:3
    st=hms_point(i,2)-len+1;
    ed=hms_point(i,2)+len;
    dir=1;
    temp1=find(hms_mask(st-1,:,i),1);
    if ~isempty(temp1)
        dir=-1;
        st=hms_point(i,2)+len;
        ed=hms_point(i,2)-len+1;
    end
    lflag=1000;
    rflag=1000;
    
    ltemp1=double.empty(1,0);
    rtemp1=double.empty(1,0);
    
    lis_vertical=0;
    ris_vertical=0;
    for j=st:dir:ed
        ltemp2=find(hms_mask(j,:,i),1);
        rtemp2=find(hms_mask(j,:,i),1,'last');
        
        % 记录从黑→白点
        if isempty(ltemp1) && ~isempty(ltemp2)
            lpoints(1,1,i)=ltemp2;
            lpoints(1,2,i)=j;
            rpoints(1,1,i)=rtemp2;
            rpoints(1,2,i)=j;
        end
        
        % 边缘斜率突变
        if ~isempty(ltemp1) && ~isempty(ltemp2)
            if lflag~=1000 && sign(lflag)~=sign(ltemp2-ltemp1)&&...
                    (ltemp2-ltemp1)~=0
                lpoints(2,1,i)=ltemp1;
                lpoints(2,2,i)=j-dir;
                if lis_vertical==1
                    lpoints(2,2,i)=j-2*dir;
                end
            end
            lflag=ltemp2-ltemp1;
        end
        
        if ~isempty(rtemp1) && ~isempty(rtemp2)
            if rflag~=1000 && sign(rflag)~=sign(rtemp2-rtemp1)&&...
                    (rtemp2-rtemp1)~=0
                rpoints(2,1,i)=rtemp1;
                rpoints(2,2,i)=j-dir;
                if ris_vertical==1
                    rpoints(2,2,i)=j-2*dir;
                end
            end
            rflag=rtemp2-rtemp1;
            
        end
        lis_vertical=(ltemp1-ltemp2)==0;
        ris_vertical=(rtemp1-rtemp2)==0;
        ltemp1=ltemp2;
        rtemp1=rtemp2;
    end
end
%% step9 -- 计算指针末端点坐标
old_hms_point=hms_point(:,1:2);
for i=1:3
    if sum(rpoints(:,:,i)==0,'all')==0
        hms_point(i,1)=floor((rpoints(1,1,i)+rpoints(2,1,i))/2);
        hms_point(i,2)=floor((rpoints(1,2,i)+rpoints(2,2,i))/2);
    else
        hms_point(i,1)=floor((lpoints(1,1,i)+lpoints(2,1,i))/2);
        hms_point(i,2)=floor((lpoints(1,2,i)+lpoints(2,2,i))/2);
    end
end
%% step10 --计算夹角
e=[1,0];% 取12点方向的单位向量
hms_vector=zeros(3,2);%算以时针、分针、秒针末端中点为尾端，以表盘中心为首端的向量的单位向量
theta_v=[];
for i =1:3
    hms_vector(i,1)=hms_point(i,1)-dial_center(1);
    hms_vector(i,2)=hms_point(i,2)-dial_center(2);
    len_v=norm(hms_vector(i,:));
    hms_vector(i,1)=hms_vector(i,1)/len_v;
    hms_vector(i,2)=hms_vector(i,2)/len_v;
    theta_v=[theta_v acos(abs(dot(hms_vector(i,:),e))...
        /(norm(hms_vector(i,:))*norm(e)))*180/pi];
end

%% step11 --结合指针点和夹角，计算时间并显示
refer=[30 6 6]; % 变化多少角度等于变化一个单位
add=[3 15 15]; % 每个象限
hms_times=zeros(3,1);
for i=1:3
    if hms_point(i,1)>=dial_center(1) && hms_point(i,2)<=dial_center(2)
        hms_times(i)=floor(theta_v(i)/refer(i));
    end
    if hms_point(i,1)>dial_center(1) && hms_point(i,2)>dial_center(2)
        hms_times(i)=floor(theta_v(i)/refer(i))+add(i);
    end
    if hms_point(i,1)<dial_center(1) && hms_point(i,2)>dial_center(2)
        hms_times(i)=floor(theta_v(i)/refer(i))+2*add(i);
    end
    if hms_point(i,1)<dial_center(1) && hms_point(i,2)<dial_center(2)
        hms_times(i)=floor(theta_v(i)/refer(i))+3*add(i);
    end
end

disp("时钟上显示"+num2str(hms_times(1))+"时"+num2str(hms_times(2))+"分"...
    +num2str(hms_times(3))+"秒");

%% 测试图展示
figure(1)
subplot(141),imshow(I);title("原图",'FontWeight','bold');
subplot(142),imshow(Gray);title("灰度图",'FontWeight','bold');
subplot(143),imshow(A);title("提取灰度范围为[20,120)为亮部的二值图",...
    'FontWeight','bold');
subplot(144),imshow(C);title("开闭运算后图",'FontWeight','bold');

% 注：figure(2)内容来自MATLAB官方文档
figure(2), imshow(C), hold on
max_len = 0;
for k = 1:length(lines)
    xy = [lines(k).point1; lines(k).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

    % Plot beginnings and ends of lines
    plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
    plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');

    % Determine the endpoints of the longest line segment
    len_temp = norm(lines(k).point1 - lines(k).point2);
    if ( len_temp > max_len)
        max_len = len_temp;
        xy_long = xy;
    end
    title("霍夫变换",'FontWeight','bold');
end

figure(3)
for i=1:3
    eval(['subplot(13',num2str(i),')'])
    imshow(hms_mask(:,:,i));
    title("掩模"+num2str(i)+"图",'FontWeight','bold');
end

figure(4)
for i=1:3
    x=hms_point(i,1);
    y=hms_point(i,2);
    eval(['subplot(13',num2str(i),')'])
    imshow(hms_mask(y-len:y+len,x-len:x+len,i));
    title("掩模"+num2str(i)+"局部放大图",'FontWeight','bold');
end

figure(5)
for i=1:3
    eval(['subplot(13',num2str(i),')'])
    imshow(hms_mask(:,:,i));
    hold on;
    plot(lpoints(:,1,i),lpoints(:,2,i),'*r');
    plot(rpoints(:,1,i),rpoints(:,2,i),'*g');
end
sgtitle("指针末端边的相聚最远的边缘点（红*为从左到右遍历，绿*相反）",...
        'FontWeight','bold');

figure(6)
for i=1:3
    eval(['subplot(13',num2str(i),')'])
    imshow(I);
    hold on;
    plot(hms_point(i,1),hms_point(i,2),'*r');
    plot(old_hms_point(i,1),old_hms_point(i,2),'*g');
    plot(dial_center(1),dial_center(2),'*b');
end
sgtitle("定位指针末端中点(红*是优化后的结果，绿*是霍夫变换得到直线上点的结果）",...
        'FontWeight','bold');

figure(7)
for i=1:3
    eval(['subplot(13',num2str(i),')'])
    imshow(I);
    hold on;
    plot(hms_point(i,1),hms_point(i,2),'*r');
    plot(old_hms_point(i,1),old_hms_point(i,2),'*g');
    plot([hms_point(i,1),dial_center(1)],[hms_point(i,2),dial_center(2)],'b','LineWidth',2);
end
sgtitle("定位指针末端点(蓝线为求出的指针中点与表盘中心的连线，绿*为仅使用霍夫变换得到的端点）",...
        'FontWeight','bold');