-
对图像数据进行预处理后,得到平滑的光斑边缘,进而采用曲线拟合对光斑的中心精确定位。圆可以看成是一个特殊的椭圆,本文中研究一般椭圆的光斑中心定位[8-9]。由图像预处理得到的所有的边缘数据对圆进行最小二乘拟合,利用求解出的方程系数可直接计算出光斑的中心位置。
-
在2维坐标系中,平面椭圆方程的隐式方程表示为F(a, x)=0,其中$a = (A, B, C, D, E, F) $为椭圆方程的参量,即椭圆隐式方程为:
$ f\left( {x, y} \right) = A{x^2} + Bxy + C{y^2} + Dx + Ey + F $
(1) 式中,$ f\left( {x, y} \right)$为椭圆方程,其中x和y分别为椭圆上点的横纵坐标值; $ A, B, C, D, E, F$为椭圆方程参量。
比较常见的椭圆拟合方法有基于代数距离的最小二乘法。由最小二乘法可知,椭圆的拟合可以近似看成是求光斑边缘上的i个数据点$\left( {{x_i}, {y_i}} \right) $到曲线$ f\left( {x, y} \right) = 0$的残差平方和的最小值,即求目标函数$f\left( {B, C, D, E, F} \right) = \sum\limits_{i = 1}^m {} {\left( {{x_i}^2 + B{x_i}{y_i} + C{y_i}^2 + D{x_i} + E{y_i} + F} \right)^2} $的最小值。再由极值原理,欲使$f\left( {B, C, D, E, F} \right) $值为最小,必有$ \frac{{\partial f}}{{\partial B}} = \frac{{\partial f}}{{\partial C}} = \frac{{\partial f}}{{\partial D}} = \frac{{\partial f}}{{\partial E}} = \frac{{\partial f}}{{\partial F}} = 0$。由此可得一个线性方程组,将其转换成求线性矩阵方程Ax=b的解。这里A为5×5维实矩阵,x为待求5维列向量,b为n维实向量。
$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^m {} {x_i}^2{y_i}^2}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}^3}&{\sum\limits_{i = 1}^m {} {x_i}^2{y_i}}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}^2}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}}\\ {\sum\limits_{i = 1}^m {} {x_i}{y_i}^3}&{\sum\limits_{i = 1}^m {} {y_i}^4}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}^2}&{\sum\limits_{i = 1}^m {} {y_i}^3}&{\sum\limits_{i = 1}^m {} {y_i}^2}\\ {\sum\limits_{i = 1}^m {} {x_i}^2{y_i}}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}^2}&{\sum\limits_{i = 1}^m {} {x_i}^2}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}}&{\sum\limits_{i = 1}^m {} {x_i}}\\ {\sum\limits_{i = 1}^m {} {x_i}{y_i}^2}&{\sum\limits_{i = 1}^m {} {y_i}^3}&{\sum\limits_{i = 1}^m {} {x_i}{y_i}}&{\sum\limits_{i = 1}^m {} {y_i}^2}&{\sum\limits_{i = 1}^m {} {y_i}}\\ {\sum\limits_{i = 1}^m {} {x_i}{y_i}}&{\sum\limits_{i = 1}^m {} {y_i}^2}&{\sum\limits_{i = 1}^m {} {x_i}}&{\sum\limits_{i = 1}^m {} {y_i}}&{\sum\limits_{i = 1}^m {} 1} \end{array}} \right]\left[ \begin{array}{l} B\\ C\\ D\\ E\\ F \end{array} \right] = \\ \left[ \begin{array}{l} - \sum\limits_{i = 1}^m {} {x_i}^3{y_i}\\ - \sum\limits_{i = 1}^m {} {x_i}^2{y_i}^2\\ - \sum\limits_{i = 1}^m {} {x_i}^3\\ - \sum\limits_{i = 1}^m {} {x_i}^2{y_i}\\ - \sum\limits_{i = 1}^m {} {x_i}^2 \end{array} \right] \end{array} $
(2) 式中,$ \left( {{x_i}, {y_i}} \right)$为光斑边界上点的坐标,$\left( {{x_i}, {y_i}} \right) \subset E $表示光斑图像中边缘上所有的点的集合。通过求解该线性方程组,便可得到椭圆方程参量B, C, D, E和F,则中心坐标$\left( {{x_{\rm{c}}}, {y_{\rm{c}}}} \right) $为:
$ \left\{ \begin{array}{l} {x_{\rm{c}}} = \frac{{BE - 2CD}}{{4C - {B^2}}}\\ {y_{\rm{c}}} = \frac{{BD - 2E}}{{4C - {B^2}}} \end{array} \right. $
(3) -
在数字信号处理(digital signal processing,DSP)系统设计开发工具中以高斯消元法为基础对线性方程组进行求解,鉴于DSP builder的除法模块只能计算出被除数和除数相除以后的商和余数,故而采用移位寄存器将被除数移位倍增,再进行相除,经验证,其精度远远小于单精度浮点数,计算的结果与理论比较误差大于10%,因此选用Quartus Ⅱ的单精度浮点运算知识产权(intellectual property,IP)核,运用IP核设计运算处理单元。
一般地,线性最小二乘问题的求解主要有迭代法和矩阵分解法两类。其中矩阵分解法可以最大限度实现并行计算,从而更适用于FPGA实现[10]。矩阵分解的方法主要有QR分解[11-12]、LU分解[13]、奇异值分解[14]以及Cholesky分解[15]。QR分解的计算过程复杂;而奇异值分解与QR分解的计算难度相当;大规模矩阵一般采用LU分解来进行分块分解,但是LU分解过程中需要大量的乘法和除法运算,会造成一定的延时,并且消耗过多的硬件资源;LDLT分解适用于对称的正定矩阵,分解过程中避免了开方运算,其计算的复杂度仅为LU分解的一半,故而本文中采取LDLT分解算法来实现矩阵的分解。
-
LDLT分解适用于对称的正定矩阵,对任意一个对称正定矩阵A,将其分解成A=LDLT的形式,其中L为对角线为1的下三角矩阵,D为对角阵,LT为L的转置,即:
$ \begin{array}{l} \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{LD}}{\mathit{\boldsymbol{L}}^{\rm{T}}} = \left[ {\begin{array}{*{20}{c}} 1&{}&{}&{}&{}\\ {{l_{21}}}&1&{}&{}&{}\\ {{l_{31}}}&{{l_{32}}}&1&{}&{}\\ \ldots&\ldots&\ldots&\ldots &{}\\ {{l_{n1}}}&{{l_{n2}}}&{{l_{n3}}}& \ldots &1 \end{array}} \right] \cdot \\ \left[ {\begin{array}{*{20}{c}} {{d_1}}&{}&{}&{}&{}\\ {}&{{d_2}}&{}&{}&{}\\ {}&{}&{{d_3}}&{}&{}\\ {}&{}&{}& \ldots &{}\\ {}&{}&{}&{}&{{d_n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&{{l_{21}}}&{{l_{31}}}& \ldots &{{l_{n1}}}\\ {}&1&{{l_{32}}}& \ldots &{{l_{n2}}}\\ {}&{}&1& \ldots&\ldots \\ {}&{}&{}&1& \ldots \\ {}&{}&{}&{}&1 \end{array}} \right] \end{array} $
(4) 矩阵D及L中的元素dj和lij根据下式计算:
$ \left\{ \begin{array}{l} {d_j} = {a_{jj}} - \sum\limits_{k = 1}^{j - 1} {} {l_{jk}}^2{d_k}\\ {l_{ij}} = \frac{{{a_{ij}} - \sum\limits_{k = 1}^{j - 1} {} {l_{ik}}{d_k}{l_{jk}}}}{{{d_j}}} \end{array} \right. $
(5) 式中,$ j = 1, 2, \ldots , n;i = j + 1, j + 2, \ldots , n$。aij和ajj分别为矩阵A中的第i行j列和第j行j列所对应的元素,dj和dk为对角阵D中对角线上所对应的的元素(k=j-1),lij,lik和ljk分别为矩阵L中第i行j列、第i行k列和第j行k列所对应的元素。
分解后即将求解的线性方程组Ax=b转换为求LDLTx=b的解,等同于求解3个三角线性方程组。
$ \left\{ \begin{array}{l} \mathit{\boldsymbol{Lz}} = \mathit{\boldsymbol{b}}\\ \mathit{\boldsymbol{Dr}} = \mathit{\boldsymbol{z}}\\ {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{r}} \end{array} \right. $
(6) 式中, z和r为引入的n维待求列向量。根据以上步骤即可求得方程的解向量x。
-
为了实现在FPGA硬件平台上高效的矩阵分解,必须合理采用流水线技术,设计并行计算[16]。算法并行化的根本是数据之间的依赖关系,下面将分析LDLT分解算法的数据之间的关系。
根据(5)式,每一个内层循环用于计算矩阵L的一列元素和矩阵D对应的一个元素,以第j个内层循环中计算的具体数据依赖关系举例说明。将(5)式展开如下所示:
$ \left\{ \begin{array}{l} {d_j} = {a_{jj}} - \left[ {{l_{i1}}^2{d_1} + {l_{i2}}^2{d_2} + \ldots + {l_{i(j - 1)}}^2{d_{(j - 1)}}} \right]\\ {l_{ij}} = [{a_{ij}} - ({l_{i1}}{d_1}{l_{j1}} + {l_{i2}}{d_2}{l_{j2}} + \ldots + \\ {l_{i(j - 1)}}{d_{(j - 1)}}{l_{j(j - 1)}})]/{d_j} \end{array} \right. $
(7) 由(7)式可知,矩阵D中元素dj的计算依赖于矩阵A中的元素ajj、矩阵L的第j行元素以及矩阵D的前(j-1)个元素。L矩阵的元素lij的计算依赖于矩阵A中的元素aij、矩阵D中的元素dj以及矩阵L的第i行元素。而lij与l(i+1)j之间的计算没有依赖关系,且dj和lij的计算都要计算中间结果ljkdk。
通过分析数据依赖关系可知,LDLT分解算法的外层数据循环计算需要依存已经获得的数据,内层数据循环计算则没有依赖关系,直接并行计算。
-
通过对LDLT分解算法及相关数据分析,可设计如图 1所示的分解结构框图。分解模块主要分为计算对角阵D元素的单元PED、下三角阵L元素的单元PEL和数据选择器,PED与PEL共享PED的计算结果ljkdk,所以可以将此结果缓存于先入先出缓存器(first in, first out, FIFO)中的直接用于lij的计算,从而避免重复运算,节省资源。PED的输入ljk通过四选一数据选择器,按照一定的顺序将4组存储矩阵L中元素的随机存取存储器(random access memory, RAM)中的数据读出,用于计算dj。矩阵A的数据通过数据分配器从外部传输到PEL与PED。其中,矩阵中aij与dj, lij各自关联,但是PED与PEL计算存在延迟,故将数据存于FIFO_a0~FIFO_a4中。当前lij的计算结果与当前dj值关联很大,且lij的计算还会用到计算dj时的中间值ljkdk,所以PED要在PEL之前执行。在理想的情况下,为了实现运算速度最快,PEL的数量越多,效率越高。输入矩阵A为5×5的矩阵,在每一个内层循环中,需要计算的矩阵D的对角线元素dj仅为1个,采用一个PED单元即可;而计算的矩阵L的元素lij的个数最多为4个,且互不依赖,故采用4个PEL并行计算。分解模块具体结构如图 3所示。其中RAMD为存储矩阵D中元素的随机存取存储器。
图 3中,PED和PEL模块均采取流水化设计方案,运算单元主要包括乘法器(multiplication)、除法器(division)、加法器(addition)和减法器(subtraction),均使用Quartus Ⅱ提供的IEEE 754规范的浮点数运算IP核实现。
在设计PED和PEL模块时,鉴于乘法器在计算时存在固有延时,所以将输入的数据ljk写入FIFO1中实现缓存,此时FIFO1就作为第2个乘法器的输入端,与第1个乘法器的计算结果进行第二次乘法运算。由于ljkdk为PED和PEL计算的共有部分,为减少乘法器的使用,故将计算中间结果ljkdk存于FIFO_ld中。累加器是由加法器和锁存器组成,每隔加法器计算的延时时间就将FIFO2中数据读出一个,累加结束后,将结果传递给减法器输入端,并使能减法器。最终的计算结果存入RAM中。内层循环结束后,循环计算5次,即可获得对角阵D和下三角阵L矩阵的元素dj和lij。
-
分解结束后,获得下三角矩阵L和对角阵D,根据(LDLT)x=b,可以等价成3个简单的方程组Lz=b,Dr=z和LTx=r,由此可得求解公式:
$ \left\{ \begin{array}{l} {z_m} = {b_m} - \sum\limits_{p = 1}^{m - 1} {} {l_{mp}}{z_p}, \left( {m = 1, 2, \ldots , q} \right)\\ {r_m} = \frac{{{z_m}}}{{{d_m}}}, \left( {m = 1, 2, \ldots , q} \right)\\ {x_m} = {r_m} - \sum\limits_{u = m + 1}^q {} {l_{um}}{x_u}, \left( {m = q, q - 1, \ldots , 1} \right) \end{array} \right. $
(8) 式中,lmp和lum分别表示矩阵L中第m行p列和第u行m列所对应的元素;dm表示矩阵D对角线上对应的元素;zm, rm, xm, bm分别代表列向量z, r, x, b的第m个元素;zp,xu分别代表列向量z的第p个元素和列向量x的第u个元素;q为待求向量的维数(q=5)。
求解模块要用到PEL模块中计算的数据lij,故每算出一个lij的值,求解模块中就要有一个计算单元与之相对应。求解模块的结构框图与PED和PEL的框图基本相同。
-
本文中基于Quartus Ⅱ 14.0的开发环境,将上述设计在Altera公司的DE1-SOC开发板上进行了验证。鉴于作者无法准确地得到光斑图像真实的中心,因此通过MATLAB程序生成一幅人工光斑图像,并附着一些噪声及干扰,该光斑图像的理论中心坐标为(150,166),分别在FPGA平台和计算机Visual Studio(VS)环境下利用本文中的算法对光斑进行中心定位,根据(3)式计算的最终结果对比如表 1所示。其中,FPGA的工作时钟设定为50MHz,并采取IEEE 754规范下的单精度浮点制IP核来实现设计中的各浮点运算。
Table 1. Deviation between FPGA and PC
coordinate of center absolute deviation visual studio (150.06907, 166.097) (0.06907, 0.097) FPGA (150.06881, 166.048) (0.06881, 0.048) 由表 1可知,本方法对圆的中心定位的精度小于0.1pixel,但FPGA平台与计算机VS环境下实现相同的算法计算的结果有差异,出现此结果的原因是,在FPGA平台使用单精度浮点IP核进行计算时,首先要将数据转换为单精度浮点数(32bit),但对于较大的数据,尾数部分会有部分省略,进而对结果产生一定的误差。不过两个平台计算的中心位置的误差都小于0.1pixel,所以该方法保证了中心定位算法的精度。
-
保证中心定位算法精度的前提下,为测试方法的效率,将实现LDLT解最小二乘问题的方法在FPGA平台和计算机VS下进行了计算效率比较。具体结果如表 2所示。
Table 2. Processing time
FPGA visual studio proportion decomposition/μs 36.38 620 1/17.04 solving/μs 37.06 995 1/26.85 total/μs 73.44 1615 1/21.99 由表 2可知,在FPGA平台下实现解5维最小二乘问题时,本文中方法在保证精度的前提下,可达到相对于计算机21倍以上的速度。
基于FPGA的光纤光斑中心定位算法研究
Locating algorithm of optical fiber spot center based on FPGA
-
摘要: 为了解决传统数字图像处理算法中数据运算量大、复杂度高、耗时长的问题,提出一种基于可编程门阵列(FPGA)光纤光斑中心定位的方法。采用数字信号处理系统,利用开发工具(DSP builder),设计了光斑图像预处理算法和边缘检测算法,用最小二乘法拟合光斑边界,采用流水线设计,增强了数据处理的并行能力,提高了处理速度。在Cyclone V平台上进行理论分析和实验验证,取得了光斑图像边界、中心坐标数据。结果表明,在保证对光斑中心定位的绝对误差小于0.1pixel的条件下,使用FPGA比计算机运算速度能提高21倍以上。该研究能够在FPGA平台上快速准确定位光斑中心。Abstract: In order to solve the problem of large amount of data, high complexity and long time consuming in traditional digital image processing algorithm, a method for center location of optical fiber spot was presented based on field-programmable gate array (FPGA). By using digital signal processing system and development tools (DSP builder), image preprocessing algorithm and edge detection algorithm were designed. The least square method was used to fit the boundary of light spot. The parallel capability for data processing was enhanced and the processing speed was improved by means of pipeline design. Theoretical analysis and experimental verification were carried out on Cyclone V platform. The image boundary and the central coordinate data were obtained. The results show that under the condition that the absolute error of spot center location is less than 0.1pixel, the computer speed can be improved by more than 21 times based on FPGA. The spot center can be located quickly and accurately on FPGA platform.
-
Key words:
- image processing /
- spot center detection /
- least square method /
- optical fiber spot
-
Table 1. Deviation between FPGA and PC
coordinate of center absolute deviation visual studio (150.06907, 166.097) (0.06907, 0.097) FPGA (150.06881, 166.048) (0.06881, 0.048) Table 2. Processing time
FPGA visual studio proportion decomposition/μs 36.38 620 1/17.04 solving/μs 37.06 995 1/26.85 total/μs 73.44 1615 1/21.99 -
[1] FAN Ch L, LIU Ch, DING G, et al.A method of circle curve fitting based on the cumulative error of the radius error[C]//International Conference on Computational Intelligence and Security. New York, USA: IEEE, 2015: 211-214. [2] WANG Q Q, LIU J, PENG Zh, et al.Measurement system for laser divergence angle based on LabView[J]. Chinese Journal of Lasers, 2012, 39(11):1108005(in Chinese). doi: 10.3788/CJL [3] ZHANG L. System of measuring laser spot[D]. Xi'an: Xidian University, 2010: 19-45(in Chinese). [4] WANG L L, HU Zh W, JI H X. Laser spot center location algorithm based on Gaussian fitting[J]. Journal of Applied Optics, 2012, 33(5):985-990(in Chinese). [5] LIU H L, HOU W, FAN Y L, et al. An improved algorithm of laser spot center location[J]. Computer Measurement & Control, 2014, 22(1):139-141(in Chinese). [6] LAN Zh L, YANG X F. Practical improvement of laser spot center location algorithm[J]. Computer Engineering, 2008, 34(6):7-9(in Chinese). [7] NING Sh N, ZHU M, SUN H H, et al. Realization of improved sobel adaptive edge detection algorithm based on FPGA[J]. Chinese Journal of Liquid Crystals and Displays, 2014, 29(3):395-402(in Chinese). doi: 10.3788/YJYXS [8] SHI Y, CHENG X. Laser spot center detection based on the geometric feature[C]//International Symposium on Information Science & Engineering. New York, USA: IEEE, 2010: 322-325. [9] LIU J S, YUAN S C. Subpixel location algorithm for circular targets center based on Zernike moments and curvature[J]. Computer Engineering & Applications, 2010, 46(29):153-156. [10] ZHONG X J. Research on FPGA calculation method of typical matrix decomposition[D]. Harbin: Harbin Institute of Technology, 2012: 14-66(in Chinese). [11] TAI Y G, LO C T D, PSARRIS K. Applying out-of-core QR decomposition algorithms on FPGA-based systems[C]//International Conference on Field Programmable Logic and Applications. New York, USA: IEEE, 2007: 86-91. [12] HAN B, YANG Z, ZHENG Y R. FPGA implementation of QR decomposition for MIMO-OFDM using four CORDIC cores[C]//IEEE International Conference on Communications. New York, USA: IEEE, 2013: 4556-4560. [13] WU G M, DOU Y, PETERSON G D. Blocking LU decomposition for FPGAs[C]//IEEE International Symposium on Field-Programmable Custom Computing Machines. New York, USA: IEEE, 2010: 109-112. [14] TANG B P, JIANG Y H, ZHANG X Ch. Feature extraction method of rolling bearing fault based on singular value decomposition-morphology filter and empirical mode decomposition[J]. Chinese Journal of Mechanical Engineering, 2010, 46(5):37-42(in Chinese). doi: 10.3901/JME.2010.05.037 [15] WU G M, DOU Y, WANG M. A fine-grained parallel algorithm for the Cholesky decomposition[J]. Computer Engineering & Science, 2010, 32(9):102-106(in Chinese). [16] GUO L, TANG Y H, ZHOU J, et al. Research on parallel architecture for LDLT decomposition -processor[J]. Computer Engineering, 2011, 37(21):241-243(in Chinese).