哈夫变换

基本哈夫变换

点线对偶性

$y=px+q$

$XY$ 空间的一条直线对应着 $PQ$ 空间的一个点。

$q=-px+y$

检测步骤

根据点线对偶性，可以基于投票算法的思想，设计出如下算法以实现检测一系列 $XY$ 空间下的点 $(x_i,y_i)$ 最有可能的共线直线 $y=px+q$ .

其中，参数 $(p,q)$ 正是我们需要确定的参数。

我们可以在参数空间 $PQ$ 中对所以可能的 $p,q$ 取值进行投票，最终得票最高的点对即为所求。具体步骤如下：

对参数空间中的参数 $p,q$ 的取值范围进行量化，根据量化结果，构造一个累加数组 $A[p_{\min}..p_{\max},\;q_{\min}..q_{\max}]$ ，并初始化为0；
对 $XY$ 空间中给定点 $(x_i,y_i)$ 取遍所有可能的 $p$ 值，由 $q=-px_i+y_i$ 计算所有可能的 $q$ 值，再对每一步得出的 $p,q$ 进行累加，即 $A[p,q]\leftarrow A[p,q]+1$ ；
找出 $A$ 中最大值所对应的 $(p,q)$ ，由此即可确立直线。这个最大值表示落在这条直线上的给定点数目最多。

待更

DM码的定位与读取

目标识别

边缘检测

为了能够识别DM码中的定位，我们首先需要对其进行边缘检测，从而将检测结果用于识别。

测试发现，如果先将图片进行滤波降噪再检测，效果会好很多。
于是我们先对原始图像选择 $5\times 5$ 的盒式滤波器进行模糊。
然后我们再利用 $Sobel$ 算子进行边缘检测。

而 $Hough$ 变换通常是对二值化图像进行检测的，因此我们还需要将检测结果二值化。
由于三幅图光线、拍摄环境、二维码位置等因素的影响，三幅图的灰度图像灰度值分布也并不相同。所以应该分别设置不同的阈值进行二值化。这里我们通过不断测试得到了较为理想的阈值选取：

Image	`dm1`	`dm2`	`dm3`
Threshold	0.1	0.5	0.4

圆检测

MATLAB 中并没有提供利用 $Hough$ 变换实现的圆形检测函数，但是理论上是可以直接编程实现的。
幸运地，我们也利用网络资源获取到了由来自哈工大的 CSDN 网友 Weisong Zhao 编写的代码.

他提供的 Hough_circle 函数需要给定二值图、半径步长、最小半径、最大半径、选取阈值的参数。
我们仍然通过不断测试得到了合理的半径选取：

Image	`dm1`	`dm2`	`dm3`
$(R_{min},R_{max})$	(150, 170)	(110, 130)	(150, 170)

此处给出的代码仅给出了对 dm1.bmp 的处理，其他图片的处理是类似的，因此省略。下同

close all;clear;clc;

% 读入图像
dm1 = im2double(imread('DM01.bmp'));

% 盒式滤波去噪
box_H = 1/25*ones(5,5);
dm1 = imfilter(dm1,box_H);

% sobel算子进行边缘检测
dm1_sobel = abs(imfilter(dm1,fspecial('sobel')'))+...
            abs(imfilter(dm1,fspecial('sobel')));
            
% 对结果进行二值化处理
% 三幅图应该分别设置了不同的阈值
dm1_sobel(dm1_sobel>0.1) = 1; dm1_sobel(dm1_sobel<0.1) = 0;

% 调用Hough_circle函数定位圆的位置
[Hough_space,Hough_circle_result,Para] = Hough_circle(dm1_sobel,5,0.1,150,170,0.6);
figure; imshow(dm1_sobel);hold on;
t=0:0.01*pi:2*pi; 
plot(Para(:,2), Para(:,1), 'g+');
plot(cos(t).*Para(1,3)+Para(1,2),sin(t).*Para(1,3)+Para(1,1),'g.');

DM circle检测

直线检测

防止DM码外部的直线对检测进行干扰，我们先根据检测得到的圆的位置信息进行图像过滤，即将圆外的二值图全部置0.

之后，我们再利用自带的 hough、houghpeaks、houghlines 函数进行直线检测。

由于已经提取到了圆内的图像，并且DM码的“最大可投票”的直线就是我们需要的那两条，这里我们选取2作为极值个数正好满足需求。

% 截取圆内图像
for i = 1:size(dm1_sobel,1)
    for j = 1:size(dm1_sobel,2)
        if ((i-Para(1,1))^2+(j-Para(1,2))^2 > Para(1,3)^2)
            dm1_sobel(i,j) = 0;dm1(i,j) = 1;
        end
    end
end
dm1 = imcrop(dm1,[Para(1,2)-Para(1,3) Para(1,1)-Para(1,3) 2*Para(1,3) 2*Para(1,3)]);
dm1_sobel = imcrop(dm1_sobel,[Para(1,2)-Para(1,3) Para(1,1)-Para(1,3) 2*Para(1,3) 2*Para(1,3)]);
figure; imshow(dm1_sobel);hold on;

% 调用自带hough函数定位直线(前2条)
[H,Theta,Rho] = hough(dm1_sobel);
P = houghpeaks(H,2,'threshold',ceil(0.1*max(H(:))));
lines = houghlines(dm1_sobel,Theta,Rho,P,'FillGap',5,'MinLength',7);
for k = 1:length(lines)
 xy = [lines(k).point1; lines(k).point2];
 plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','g');
end

DM 直线检测

几何校正

根据 Data Matrix 码的表示方法，我们需要将图像根据定位到的直线进行旋转，使得直线呈“L”型。

此外，为了获得纯净的二维码图像，我们再根据二维码边长进一步裁剪图像，然后再对原图像进行二值化。

同样地，测试得到良好的阈值：

Image	`dm1`	`dm2`	`dm3`
Threshold	0.23	0.34	0.60

矫正缺陷：由于实验图像基本都是从正面垂直拍摄的，因此获取到的图像除了旋转要求外无需再做处理。但是实际中往往需要进行真正的 几何校正，相当于将倾斜拍摄的二维码也进行转正。这里我们并没有做进一步的优化。

数据读取算法设计

得益于 DM码简单的编码方式，我们只需对获取到的正面二值化 DM码图像中黑色和白色方块进行扫描存储即可。

算法步骤

通过检测到的直线长度计算方格的边长 $a=len/(nums+2)$ ，其中 $len$ 是直线长度， $nums$ 是已知方格的行数；
创建一个DM码的数据存储矩阵 $Data$ ，大小为 $(nums+2)\times (nums+2)$ ，初始化为全0；
创建一个和方格大小（ $a\times a$ ）一样的全1矩阵 $Cube=\begin{bmatrix}1&\cdots&1\\\vdots&&\vdots\\1&\cdots&1\end{bmatrix}$ ；
从图像的左上角 $G_{11}=\begin{bmatrix}f(x_1,y_1)&\cdots&f(x_1,y_1+a)\\\vdots&&\vdots\\f(x_1+a,y_1)&\cdots&f(x_1+a,y_1+a)\end{bmatrix}$ 开始，将其与 $Cube$ 进行卷积（点乘）；
计算出结果矩阵中的“含1量” ： $tmp=\sum (G_{11}*Cube) /a^2$ ；
如果含1量大于给定阈值 $T$ ，则将数据矩阵对于位置填1，即 $Data[1,1] = 1$ ；
滚动到下一个方格 $G_{12},G_{13},\cdots$ 并重复步骤4~5直到所有方格计算完成；
截取 $Data$ 矩阵中的 $nums\times nums$ 数据部分，该内容即为 DM 码内容。

同样地，测试得到良好的阈值：

Image	`dm1`	`dm2`	`dm3`
Threshold	0.400	0.350	0.332

% 对图像进行角度旋转与裁边
len = sqrt((lines(1).point1(1)-lines(1).point2(1))^2+(lines(1).point1(2)-lines(1).point2(2))^2);
dm1 = imcomplement(dm1);
dm1 = imrotate(dm1,lines(1).theta);
dm1 = imcomplement(dm1);
dm1 = im2bw(dm1,0.23);
dm1 = imcrop(dm1,[size(dm1,1)/2-len/2 size(dm1,2)/2-len/2 len len]);
figure; imshow(dm1);

% 根据方格边长遍历读取二维码信息
a = round(len/12);
cube = ones(a,a);
result = zeros(12,12);
for i = 1:a:size(dm1,1)-a
    for j = 1:a:size(dm1,2)-a
        tmp = dm1(i:i+a-1,j:j+a-1).*cube;
        if sum(tmp(:))/(a*a) > 0.4
            result((i-1)/a+1,(j-1)/a+1) = 1;
        end
    end
end
result = result(2:11,2:11);
figure;imagesc(result);axis equal,axis tight,axis off