CN102609943A - 一种基于线性Radon变换算法的图像处理方法 - Google Patents

Abstract

本发明涉及一种基于线性Radon变换算法的图像处理方法，该方法通过归一化线性Radon变换算法，在沿某一方向和位置的直线积分时同时计算积分长度，然后，将积分值除以积分长度得到在变换域反映积分直线上图像灰度的均值，以分辨出图像中的微弱线性特征。本发明提出了归一化线性Radon变换算法，可以在变换域反映积分直线上图像灰度均值，有利于检测出图像中的微弱线性特征。

Description

一种基于线性Radon变换算法的图像处理方法

技术领域

本发明涉及图像处理的技术领域，特别涉及一种基于线性Radon变换算法的图像处理方法。

背景技术

奥地利数学家J.Radon于1917年提出Radon变换及其逆变换以后，在医学、物理学、天文学等许多领域都已得到了广泛的应用。首次将Radon变换理论应用于CT研究的A.M.Cormark(美国)及第一台CT制作者C.N.Hounsfield(英国)因此荣获1979年诺贝尔医学和生理学奖。目前，Radon变换在许多学科领域所取得的研究成果已引起人们的极大关注，是图像处理中的一种重要研究方法。Radon变换分为线性Radon变换和非线性Radon变换。

在图像处理中，线性Radon变换可以将图像域中的一条直线映射为变换域中的一个点，根据这一性质，可以将灰度图像进行Radon变换，使图像中偏亮或偏暗的线性特征转化为变换域中的峰值点或谷值点，因而在图像线性检测中得到了广泛应用。其中，

二维欧氏空间x-y中函数f(x，y)的Radon变换定义为：

$R(\mathrm{\θ},\mathrm{\ρ})=\underset{D}{\∫\∫}f(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx};$

其中，D是整个x-y平面；θ-ρ平面为变换域平面，θ为从所检测直线的法线与x轴的夹角，ρ为从原点到所检测直线的垂直距离；

δ为Dirac函数：

$\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})=\left\{\begin{array}{c}1,\mathrm{\ρ}=x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\\ 0,\mathrm{\ρ}\≠x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\end{array};\right.$

对于图像，f(x，y)＝I(x，y)为(x，y)处的灰度值。

在实际应用时，由于图像往往为一个具有一定宽度和高度长方形(或正方形)区域，因此，根据Radon变换定义，在沿不同方向或位置的直线进行积分时，积分线段的长度会相差较大，积分结果不但与积分直线上图像的灰度值有关，还与积分线段的长度有关，因而变换域中各点值的大小并不能准确反映图像中对应方向和位置直线上的灰度的整体情况。

为解决这一问题，现有技术在一些直线检测的图像处理中，往往先对图像进行边缘检测生成二值图像，然后再对二值图像进行Radon变换。这种作法对直线定位有效，但其检测直线的能力却取决于边缘检测算法，并没有利用Radon变换提取线性特征的功能。

发明内容

本发明的目的在于，提供一种基于线性Radon变换算法的图像处理方法。该方法基于归一化线性Radon变换算法，在沿某一方向和位置的直线积分时，还计算积分长度，然后将积分值除以积分长度，从而可以在变换域反映积分直线上图像灰度均值，有利于检测出图像中的微弱线性特征。

为实现上述发明目的，本发明提供了一种基于线性Radon变换算法的图像处理方法，该方法通过归一化线性Radon变换算法，在沿某一方向和位置的直线积分时同时计算积分长度，然后，将积分值除以积分长度得到在变换域反映积分直线上图像灰度均值，以分辨出图像中的微弱线性特征。

作为上述技术方案的一种改进，所述的积分直线上图像灰度均值通过归一化线性Radon变换得到：

$\stackrel{\‾}{R}(\mathrm{\θ},\mathrm{\ρ})=\frac{1}{K(\mathrm{\θ},\mathrm{\ρ})}\underset{D}{\∫\∫}f(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx};$

其中， $K(\mathrm{\θ},\mathrm{\ρ})=\underset{D}{\∫\∫}\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx},$ 是指积分长度；

式中，D是整个x-y平面；θ-ρ平面为变换域平面，θ为从x轴到所检测直线的法线的角度，ρ为从原点到所检测直线的垂直距离；

δ为Dirac函数：

$\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})=\left\{\begin{array}{c}1,\mathrm{\ρ}=x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\\ 0,\mathrm{\ρ}\≠x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\end{array};\right.$

f(x，y)为图像的元素I(x，y)在(x，y)处的灰度值。

作为上述技术方案的另一种改进，在利用Radon变换检测图像中的线性特征时，所述的图像中如果存在一些明显偏亮或明显偏暗的奇异区，则在剔除奇异区后再对图像进行归一化线性Radon变换：

$\stackrel{\‾}{R}(\mathrm{\θ},\mathrm{\ρ})=\frac{1}{L(\mathrm{\θ},\mathrm{\ρ})}\underset{D}{\∫\∫}f(x,y)M(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx};$

其中， $L(\mathrm{\θ},\mathrm{\ρ})=\∫\underset{D}{\∫}M(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx},$ 是指积分长度；

式中，

作为上述技术方案的再一种改进，所述的归一化线性Radon变换是以图像的中点为图像域卡迪尔坐标系的坐标原点。

作为上述技术方案的再一种改进，所述的图像的数字图像I的大小为H×W个像素，其中，H是指图像的高度，W是指图像的宽度，该图像以二维数组表示，第x行第y列元素表示图像坐标系中位置(x，y)处的像素的灰度值；图像的归一化线性Radon变换结果用二维数组R保存，变换时的积分长度用二维数组IL保存，将二维数组R和对应的积分长度数组IL的对应元素相除，并将结果保存在数组R中，即为归一化线性Radon变换得到图像数据。

作为上述技术方案的再一种改进，所述的二维数组R和IL的大小相同，

它们的行数均为： $L=\mathrm{int}(2\×[0.5+\sqrt{{\left(\frac{H}{2}\right)}^2+{\left(\frac{W}{2}\right)}^2}\left]\right),$ 其中int(·)表示取整；

它们的列数均由变换时积分直线的方位角θ的取值决定，如果θ在[0，π)区间上以等间隔取M个点，

k＝0，1，2，…M-1，则二维数组R和IL的列数为M。

作为上述技术方案的再一种改进，所述的各个积分直线方位角θ[k](k＝0，1，2，…M-1)的正弦和余弦值，分别保存在一维数组CosTheta和SinTheta中，用于直接查值以减少计算量。

作为上述技术方案的再一种改进，所述的图像中的各个元素I(x，y)通过点(x，y)的各方位角θ[k](k＝0，1，2，…M-1)积分直线到坐标原点的垂直距离d为：

$d=\mathrm{int}\left((x-\frac{H}{2})\cos \mathrm{\θ}\right[k]+(y-\frac{W}{2})\sin \mathrm{\θ}\left[k\right]).$

本发明的优点在于，在实际应用时，由于数字图像往往为一个具有一定宽度和高度的长方形(或正方形)区域，在进行Radon变换时，沿不同方向或位置的直线进行积分，积分结果不但与积分直线上图像的灰度值有关，还与积分线段的长度有关，因而变换域中各点值的大小并不能准确反映图像中对应方向和位置直线上的灰度的整体情况。为解决这一问题，本发明提出了归一化线性Radon变换算法，它在沿某一方向和位置的直线积分时，同时计算积分长度，然后将积分值除以积分长度，从而可以在变换域反映积分直线上图像灰度均值，有利于检测出图像中的微弱线性特征。

附图说明

图1是本发明基于线性Radon变换算法的图像处理方法的一实施例的流程图；

图2是本发明基于线性Radon变换算法的图像处理方法的另一实施例的流程图；

图3待处理的原始的包含线性特征的灰度图像；

图4是对图3进行传统Radon变换的结果；

图5是图3进行本发明归一化Radon变换的结果；

图6是采用两种变换方法的结果在相同区域的比较图例，其中，(a)为传统Radon变换结果，(b)为归一化Radon变换结果。

具体实施方式

为了更好地理解本发明的技术方案，以下结合附图对本发明的实施步骤做进一步的描述。

Radon变换在医学、物理学、天文学等许多领域都已得到了广泛的应用。在图像处理中，线性Radon变换可以将图像域中的一条直线映射为变换域中的一个点。根据这一性质，可以将灰度图像进行Radon变换，使图像中偏亮或偏暗的线性特征转化为变换域中的峰值点或谷值点，可用于检测灰度图像中的线性特征。

在实际应用时，由于数字图像往往为一个具有一定宽度和高度的长方形(或正方形)区域，在进行Radon变换时，沿不同方向或位置的直线进行积分，积分结果不但与积分直线上图像的灰度值有关，还与积分线段的长度有关，因而变换域中各点值的大小并不能准确反映图像中对应方向和位置直线上的灰度的整体情况。为解决这一问题，本发明提供的一种基于线性Radon变换算法的图像处理方法的步骤为：

设输入灰度数字图像I的大小为H×W(高度×宽度)个像素，以图像中点为坐标原点，θ在[0，π)区间上以等间隔取M个点，

k＝0，1，2，…M-1；L为图像四角到原点的距离，L＝int(sqrt(H*H/4+W*W/4)+0.5)；输出为输入数字图像I的归一化线性Radon变换R，为二维数组，大小为(2×L)×M；

则归一化Radon变换算法如图1所示：

(1)初始化R数组，将所有元素置为0，即R[i，j]＝0(i＝0，1，2，…2×L-1，j＝0，1，2，…L-1)；

(2)计算各θ[k](k＝0，1，2，…M-1)的正弦和余弦值，分别保存在一维数组CosTheta和SinTheta中；；

(3)初始化积分长度数组IL，将所有元素置为0，即IL[i，j]＝0(i＝0，1，2，…2×L-1，j＝0，1，2，…M-1)；

(4)令整数x＝0；

(5)令整数y＝0；

(6)令整数k＝0；

(7)计算；

d＝int(double(x-H/2)×CosTheta(k)+double(y-W/2)×SinTheta(k)+L)；

R(d，k)+＝double(I(y，x))；

IL(d，k)+＝1；

(8)k+＝1；如果k＜M则返回第(7)步；

(9)y+＝1；如果y＜H则返回第(6)步；

(10)x+＝1；如果x＜W则返回第(5)步；

(11)将二维数组R和IL的对应元素相除并将结果保存在数组R中。

如果需要对多幅相同大小图像进行Radon变换，则可以先计算IL并保存，然后在计算图像的Radon变换时，不再需要进行第(8)步的计算，可以减少计算量。

在利用Radon变换检测图像中的线性特征时，图像中一些明显偏亮或明显偏暗的区域(这里称其为奇异区)会对变换结果影响较大，不利于或明或暗的微弱线性特征的检出。为此，往往需要先检测出这样的区域，然后在对图像进行Radon变换时，跳过这些奇异区，则可以避免其对变换结果的影响。

设剔除奇异区后的有效模板为

则此时的归一化线性Radon变换为：

$\stackrel{\‾}{R}(\mathrm{\θ},\mathrm{\ρ})=\frac{1}{L(\mathrm{\θ},\mathrm{\ρ})}\underset{D}{\∫\∫}f(x,y)M(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx}$

其中， $L(\mathrm{\θ},\mathrm{\ρ})=\underset{D}{\∫\∫}M(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx}$

则此时的归一化Radon变换算法如图2所示：

(1)初始化R数组，将所有元素置为0，即R[i，j]＝0(i＝0，1，2，…2*rho_Max-1，j＝0，1，2，…theta_Max-1)；

(2)计算各theta[i](i＝0，1，2，…theta_Max-1)的正弦和余弦值，分别保存在一维数组CosTheta和SinTheta中；cenX＝W/2；ceny＝H/2；

(3)初始化积分长度数组IL，将所有元素置为0，即IL[i，j]＝0(i＝0，1，2，…2*rho_Max-1，j＝0，1，2，…theta_Max-1)；

(4)令整数x＝0；

(5)令整数y＝0；

(6)如果M[x，y]＝0，则直接跳到第(11)步；

(7)令整数theta＝0；

(8)计算

rho＝int(double(x-cenX)*CosTheta(theta)+double(y-cenY)*SinTheta(theta)+rho_Max)；

R(rho，theta)+＝double(I(y，x))；

IL(rho，theta)+＝1；

(9)theta+＝1；如果theta＜theta_Max则返回第(7)步；

(10)y+＝1；如果y＜H则返回第(6)步；

(11)x+＝1；如果x＜W则返回第(5)步；

(12)将二维数组R和IL的对应元素相除并将结果保存在数组R中。

综上，对于图像I中的各个像素I(x，y)(x＝0，1，2，…H-1，y＝0，1，2，…W-1)，计算沿通过该点的各方位角θ[k](k＝0，1，2，…M-1)直线积分时像素灰度值对积分结果和积分长度的贡献：

(a)计算通过该点(x，y)的各方位角θ[k](k＝0，1，2，…M-1)积分直线到坐标原点的距离d：

$d=(x-\frac{H}{2})\cos \mathrm{\θ}\left[k\right]+(y-\frac{W}{2})\sin \mathrm{\θ}\left[k\right];$

式中，θ[k]的余弦和正弦值可直接从一维数组CosTheta和SinTheta中查值，以减少计算量。

将d加上L/2并取整保存在整型变量r中。；

(b)将数组R的元素R(d，k)加上图像在(x，y)处的像素的灰度值I(x，y)作为R(r，k)的新值；

(c)将积分长度IL的元素IL(r，k)加上1作为IL(r，k)的新值。；

(4)将二维数组R和积分长度二维数组IL中对应的元素相除，并将结果保存在图像灰度值数组R中。

图4和图5分别是采用传统和本发明的方法对图3的原始图像进行处理的结果，由于积分长度不同形成的“XX”形波峰脊影响了图像中线性特征形成的波峰和波谷的检测，而在图5中很容易分辨图像中线性特征形成的波峰和波谷。

图6是两种变换结果在相同区域的比较，(a)为传统Radon变换结果，(b)为归一化Radon变换结果，我们可以更加清楚地看到，在(b)中很容易分辨图像中线性特征形成的波峰和波谷。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (8)

1.一种基于线性Radon变换算法的图像处理方法，该方法通过归一化线性Radon变换算法，在沿某一方向和位置的直线积分时同时计算积分长度，然后，将积分值除以积分长度得到在变换域反映积分直线上图像灰度均值，以分辨出图像中的微弱线性特征。

2.根据权利要求1所述的基于线性Radon变换算法的图像处理方法，其特征在于，所述的积分直线上图像灰度均值通过归一化线性Radon变换得到：

$\stackrel{\‾}{R}(\mathrm{\θ},\mathrm{\ρ})=\frac{1}{K(\mathrm{\θ},\mathrm{\ρ})}\underset{D}{\∫\∫}f(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx};$

其中， $K(\mathrm{\θ},\mathrm{\ρ})=\underset{D}{\∫\∫}\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx},$ 是指积分长度；

式中，D是整个x-y平面；θ-ρ平面为变换域平面，θ为从x轴到所检测直线的法线的角度，ρ为从原点到所检测直线的垂直距离；

δ为Dirac函数：

$\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})=\left\{\begin{array}{c}1,\mathrm{\ρ}=x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\\ 0,\mathrm{\ρ}\≠x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\end{array};\right.$

f(x，y)为图像的元素I(x，y)在(x，y)处的灰度值。

3.根据权利要求1所述的基于线性Radon变换算法的图像处理方法，其特征在于，在利用Radon变换检测图像中的线性特征时，所述的图像中如果存在一些明显偏亮或明显偏暗的奇异区，则在剔除奇异区后再对图像进行归一化线性Radon变换：

$\stackrel{\‾}{R}(\mathrm{\θ},\mathrm{\ρ})=\frac{1}{L(\mathrm{\θ},\mathrm{\ρ})}\underset{D}{\∫\∫}f(x,y)M(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx};$

其中， $L(\mathrm{\θ},\mathrm{\ρ})=\underset{D}{\∫\∫}M(x,y)\mathrm{\δ}(\mathrm{\ρ}-x\cos \mathrm{\θ}-y\sin \mathrm{\θ})\mathrm{dydx},$ 是指积分长度；

式中，

4.根据权利要求1～3中任一项所述的基于线性Radon变换算法的图像处理方法，其特征在于，所述的归一化线性Radon变换是以图像的中点为图像域卡迪尔坐标系的坐标原点。

5.根据权利要求1～3中任一项所述的基于线性Radon变换算法的图像处理方法，其特征在于，所述的图像的数字图像I的大小为H×W个像素，其中，H是指图像的高度，W是指图像的宽度，该图像以二维数组表示，第x行第y列元素表示图像坐标系中位置(x，y)处的像素的灰度值；图像的归一化线性Radon变换结果用二维数组R保存，变换时的积分长度用二维数组IL保存，将二维数组R和对应的积分长度数组IL的对应元素相除，并将结果保存在数组R中，即为归一化线性Radon变换得到结果。

6.根据权利要求5所述的基于线性Radon变换算法的图像处理方法，其特征在于，所述的二维数组R和IL的大小相同，

它们的行数均为： $L=\mathrm{int}(2\×[0.5+\sqrt{{\left(\frac{H}{2}\right)}^2+{\left(\frac{W}{2}\right)}^2}\left]\right),$ 其中int(·)表示取整；

它们的列数均由变换时积分直线的方位角θ的取值决定，如果θ在[0，π)区间上以等间隔取M个点，

k＝0，1，2，…M-1，则二维数组R和IL的列数为M。

7.根据权利要求6所述的基于线性Radon变换算法的图像处理方法，其特征在于，所述的各个积分直线方位角θ[k](k＝0，1，2，…M-1)的正弦和余弦值，分别保存在一维数组CosTheta和SinTheta中，用于直接查值以减少计算量。

8.根据权利要求6所述的基于线性Radon变换算法的图像处理方法，其特征在于，所述的图像中的各个元素I(x，y)通过点(x，y)的各方位角θ[k](k＝0，1，2，…M-1)积分直线到坐标原点的垂直距离d为：

$d=\mathrm{int}\left((x-\frac{H}{2})\cos \mathrm{\θ}\right[k]+(y-\frac{W}{2})\sin \mathrm{\θ}[k\left]\right).$

Images