Getting started with Computer Vision on FPGAs Get fami
作者：禅与计算机程序设计艺术
1.简介
近年来，随着越来越多的人将目光转向移动互联网、物联网和智能设备领域，人工智能（AI）成为当前科技热点。而人工智能的一个重要分支——计算机视觉（CV）则是其中的重要组成部分之一。本文将带领读者熟悉并掌握计算机视觉在FPGA上的基础知识、技术要点和基本操作方法，并使用HLS(High-Level Synthesis)工具对其进行实现。最终实现图像边缘检测、图像特征提取等一些基本的图像处理算法的加速计算。本文假设读者了解FPGA硬件、HLS工具和图像处理相关的基本概念。
2.相关概念与术语
首先，我们需要对计算机视觉及其相关术语有所了解。下面是最重要的几个词汇：
① 图像：图像是由像素点组成的矩阵，每一个像素点都可以看作一个二维或三维空间上的点。一般情况下，图像通常是有灰度值或者彩色值的。
② 相机：相机就是一个能够拍摄图像的装置，包括照相机、激光摄影机和红外摄影机等。
③ 显示器：显示器就是用来显示图像的输出设备。
④ 相机矩阵：相机矩阵描述了从世界坐标系到相机坐标系的变换关系。
⑤ 摄像头参数：摄像头参数是指相机内部的某些参数，如焦距、畸变、白平衡等。
⑥ 分辨率：分辨率就是图像在一个像素点的位置所占用的空间大小。分辨率越高，图像细节就越少。
⑦ 分场方式：分场方式也称为图像传感器模式，它是在多光谱范围内同时捕获不同波长的光线，从而获得多通道图像。
⑧ 图像增强：图像增强即对图像进行图像处理，如锐化、去噪、均衡曝光、锯齿纹理消除等。
⑨ 阈值化：阈值化是一种简单有效的方法，通过设定一个阈值将图像划分为不同的区域。
⑩ 轮廓检测：轮廓检测就是识别图像中各个明显的区域和边界。
3.核心算法原理和具体操作步骤
3.1 Sobel算子
Sobel算子是一种对图像进行卷积运算的算法。在Sobel算子中，有一个对角方向上的核函数求导得到另一个方向上的梯度信息。通过求取两个方向上的梯度值，就可以提取出图像边缘信息。下图是Sobel算子示意图：
假设要在图像I上应用Sobel算子，首先需要确定卷积核的大小，一般选取3x3或5x5的矩形核。然后，将核函数横纵拉伸，使得每个元素表示了相邻两行或列像素的差异，再进行离散余弦变换，即可得到边缘强度图。
3.2 Canny边缘检测
Canny边缘检测是一种基于阈值的图像边缘检测算法。首先，利用Sobel算子计算图像梯度，然后根据梯度值、阈值和NMS(Non-Maxima Suppression)操作得到图像边缘。这里的阈值选择是根据经验值，也可采用Otsu方法自动优化。Canny边缘检测的主要步骤如下：
- 低通滤波：通过高斯滤波平滑图像边缘。
- 算子阈值：用低通滤波后的图像生成梯度幅值图。
- NMS：通过极值抑制处理边缘检测结果。
最后，返回原始图像的边缘作为输出。下图是一个Canny边缘检测的示意图：
3.3 Harris角点检测
Harris角点检测是一种检测图像局部特征的方法。该算法利用图像梯度的弯曲程度及其方向信息来检测角点。对于每一个像素，该算法计算周围的16个像素点的梯度幅值，然后计算这些像素点的梯度平方和。如果某个像素的梯度平方和大于某个常数，就认为它是角点。然而，这种方法有一个缺陷就是对噪声敏感。为了缓解这一缺陷，Harris角点检测还引入了一个抗噪声的过程。
Harris角点检测的主要步骤如下：
- 提取特征：对于每一个像素，提取水平和垂直方向上的梯度幅值及其方向。
- 评价特征：对于每个特征点，评估其质量并决定是否保留它。
- 非最大抑制：通过比较其邻域内的其他特征点，消除重复的候选点。
3.4 KLT光流跟踪
KLT光流跟踪是一种快速、准确的用于计算图像移动目标的光流的方法。该算法将图像中的每一个像素看作一个“带状”点，每次迭代时计算这些点的运动向量。通常情况下，KLT算法分成三个阶段：
- 初始化：初始化图像中的所有带状点的运动向量。
- 跟踪：根据前一帧的运动信息更新运动向量。
- 回溯：根据运动向量计算每个像素的位置。
KLT算法的主要优点是速度快、迭代次数少。但是，它要求初始猜测精度较高，而且无法处理目标出现在图像边缘的问题。
4.实践项目
4.1 Sobel算子在FPGA上的实现
下面展示如何在Xilinx Virtex Ultrascale+ FPGA上实现Sobel算子。首先，我们要创建HLS工程并配置好其编译选项。在C++源文件中包含图像处理库和OpenCV API。然后，按照以下步骤进行实现：
- 配置数据类型：由于Xilinx的DSP48E1是双精度浮点型，所以我们设置图像的数据类型为float。
- 将图像输入DDR缓存中：将图像数据存入DDR缓存，以便FPGA进行处理。
- 对图像进行复制：由于Xilinx FPGA的资源限制，我们不能一次性读取整个图像，因此需要对图像进行复制，以便能存放多个图像块。
- 执行卷积核操作：对图像块执行卷积核操作，得到结果存储在DDR缓存中。
- 存储结果：将卷积核操作的结果存入DDR缓存。
- 从DDR缓存中读取结果：从DDR缓存中读取卷积核操作的结果。
- 输出结果：输出卷积核操作的结果。
下面的代码片段展示了在FPGA上实现Sobel算子的代码。
void sobel(cv::Mat &src, cv::Mat &dst) {
const int WIDTH = src.cols;
const int HEIGHT = src.rows;
// Copy the input image to the DDR memory for acceleration
short img[WIDTH*HEIGHT];
for (int i=0; i<WIDTH*HEIGHT; ++i)
img[i] = (short)(src.at<uchar>(i/WIDTH, i%WIDTH));
// Read from global memory into local buffer
short inbuf[KERNEL_SIZE][KERNEL_SIZE], outbuf[KERNEL_SIZE][KERNEL_SIZE];
#pragma HLS ARRAY_PARTITION variable=inbuf dim=1 complete
#pragma HLS ARRAY_PARTITION variable=outbuf dim=1 complete
ap_uint<1> px, py;
int x, y, dx, dy, mag;
float fmag;
// Apply filter to each pixel of the source image
loop_height:for (y=0; y<HEIGHT-KERNEL_SIZE+1; y+=STRIDE) {
loop_width:for (x=0; x<WIDTH-KERNEL_SIZE+1; x+=STRIDE) {
loop_kernel:for (dy=-KERNEL_RADIUS; dy<=KERNEL_RADIUS; ++dy) {
loop_kernel_x:for (dx=-KERNEL_RADIUS; dx<=KERNEL_RADIUS; ++dx) {
if ((x+dx>=0 && x+dx<WIDTH) && (y+dy>=0 && y+dy<HEIGHT)) {
int index = ((y+dy)*WIDTH + (x+dx))*NUM_CHANNELS;
// Read data from global memory to local buffer
loadData(img, index, NUM_CHANNELS, inbuf);
// Compute gradient magnitude
computeGradientMag(&px, &py, inbuf, &fmag);
// Store result in output buffer
storeResult((y/STRIDE)*(WIDTH/STRIDE)+(x/STRIDE),
(px*(py==0))+0*((py==1)<<1)+1*((py==2)<<1),
magToAngle(fmag), &outbuf);
}
}
}
// Write results back to global memory
writeData(img, outbuf, x, y);
}
}
// Convert floating point values to unsigned chars and copy them to destination matrix
dst = cv::Mat(src.size(), CV_8UC1);
for (int i=0; i<WIDTH*HEIGHT; ++i)
dst.at<uchar>(i/WIDTH, i%WIDTH) = (unsigned char)(outbuf[(i/WIDTH)/STRIDE][(i%WIDTH)/STRIDE]);
}
代码解读4.2 Canny边缘检测在FPGA上的实现
下面展示如何在Xilinx Virtex Ultrascale+ FPGA上实现Canny边缘检测。首先，我们要创建HLS工程并配置好其编译选项。在C++源文件中包含图像处理库和OpenCV API。然后，按照以下步骤进行实现：
- 配置数据类型：由于Xilinx的DSP48E1是双精度浮点型，所以我们设置图像的数据类型为float。
- 设置边缘检测的参数：设置经验值、阈值、平滑窗口大小等。
- 创建窗口结构：根据算法参数创建一个窗口结构，其中包含Sobel算子的输出、梯度幅值图、边缘强度图等。
- 图像边缘提取：循环遍历所有的窗口，利用Sobel算子计算每个像素的梯度幅值和方向。如果像素满足阈值条件，就将它添加到边缘强度图中。
- 边缘链接：遍历边缘强度图中的每条边缘，计算该边缘连接的所有连续边缘，并标记它们为边缘点。
- 绘制边缘：根据标记的边缘绘制图像的边缘。
- 返回边缘：将边缘输出给CPU。
下面的代码片段展示了在FPGA上实现Canny边缘检测的代码。
void cannyEdgeDetector(cv::Mat& src, cv::Mat& dst, double threshold1,
double threshold2, int apertureSize) {
const int WIDTH = src.cols;
const int HEIGHT = src.rows;
// Allocate space for window structure
struct WindowStruct {
float gradx[WINDOW_SIZE][WINDOW_SIZE]; // Horizontal gradients
float grady[WINDOW_SIZE][WINDOW_SIZE]; // Vertical gradients
float mgrad[WINDOW_SIZE][WINDOW_SIZE]; // Magnitudes
float orient[WINDOW_SIZE][WINDOW_SIZE]; // Orientations
uchar edge[WINDOW_SIZE][WINDOW_SIZE]; // Edge map
};
WindowStruct winstruc;
// Initialize parameters
const float PI = 3.14159265f;
const float eps = 0.0001f;
int wndRadius = getWindowRadius(apertureSize);
float sigma = 0.3f * ((wndRadius <= 2)? 1 : kGaussianFilterFactor);
const float smoothThreshold = std::max(threshold1 / 255.f,
CANNY_DEFAULT_LOWER_THRESHOLD);
float lowThreshRatio = threshold1 / (eps + maxVal(src));
float highThreshRatio = threshold2 / (eps + maxVal(src));
// Allocate temporary buffers
uchar tmpBuff[MAX_IMG_SIZE][MAX_IMG_SIZE];
// Loop through all pixels and detect edges using Canny algorithm
loop_height:for (int y=wndRadius; y<HEIGHT-wndRadius; ++y) {
loop_width:for (int x=wndRadius; x<WIDTH-wndRadius; ++x) {
bool isBoundaryPoint = false;
// Extract subwindow around current pixel
extractSubwindow(src, x, y, wndRadius, &winstruc);
// Calculate gradients and orientation angle
calculateGradientsAndOrientation(wndRadius, &winstruc);
// Calculate magnitude and thresholding
calculateMagnitudeAndThreshold(smoothThreshold, lowThreshRatio,
highThreshRatio, &winstruc);
// Link connected components
linkEdges(&isBoundaryPoint, &winstruc);
// Mark boundary points
markBoundaries(isBoundaryPoint, &winstruc);
// Draw edges at original resolution
drawEdgesAtOriginalResolution(tmpBuff, x, y, wndRadius, &winstruc);
}
}
// Create final output image by copying from temporal buffer
convertOutputImage(tmpBuff, dst, WIDTH, HEIGHT);
}
代码解读4.3 Harris角点检测在FPGA上的实现
下面展示如何在Xilinx Virtex Ultrascale+ FPGA上实现Harris角点检测。首先，我们要创建HLS工程并配置好其编译选项。在C++源文件中包含图像处理库和OpenCV API。然后，按照以下步骤进行实现：
- 配置数据类型：由于Xilinx的DSP48E1是双精度浮点型，所以我们设置图像的数据类型为float。
- 设置角点检测的参数：设置邻域半径、标准差等。
- 创建窗口结构：根据算法参数创建一个窗口结构，其中包含图像、梯度幅值图、角点强度图等。
- 图像角点检测：循环遍历所有的窗口，利用Harris角点检测算法计算每个窗口的角点强度。
- 过滤角点：根据角点强度筛选角点，并确定它们的方向。
- 返回角点：将角点输出给CPU。
下面的代码片段展示了在FPGA上实现Harris角点检测的代码。
void harrisCornerDetection(cv::Mat& src, cv::Mat& dst,
int blockSize, int apertureSize, double k) {
const int WIDTH = src.cols;
const int HEIGHT = src.rows;
// Define kernel size used in differentiation
const int X_ORDER = 1;
const int Y_ORDER = 0;
// Define number of bins used for computing Harris response function
const int BIN_COUNT = 256;
// Create empty Mat object for storing corner responses
cv::Mat response(cv::Size(WIDTH, HEIGHT), CV_32FC1);
// Set up Gaussian smoothing kernel for derivative calculation
int wndRadius = getDerivativeKernelRadius(blockSize, apertureSize);
cv::Mat xkernel = createDerivativeKernel(wndRadius, BLOCK_BORDER_SIZE, X_ORDER, true);
cv::Mat ykernel = createDerivativeKernel(wndRadius, BLOCK_BORDER_SIZE, Y_ORDER, true);
// Allocate memory for temporary buffers
float Ixx[BLOCK_SIZE+BLOCK_BORDER_SIZE][BLOCK_SIZE+BLOCK_BORDER_SIZE];
float Iyy[BLOCK_SIZE+BLOCK_BORDER_SIZE][BLOCK_SIZE+BLOCK_BORDER_SIZE];
float Ixy[BLOCK_SIZE+BLOCK_BORDER_SIZE][BLOCK_SIZE+BLOCK_BORDER_SIZE];
cv::Mat tempIxx(cv::Size(BLOCK_SIZE+BLOCK_BORDER_SIZE, BLOCK_SIZE+BLOCK_BORDER_SIZE),
CV_32FC1);
cv::Mat tempIyy(cv::Size(BLOCK_SIZE+BLOCK_BORDER_SIZE, BLOCK_SIZE+BLOCK_BORDER_SIZE),
CV_32FC1);
cv::Mat tempIxy(cv::Size(BLOCK_SIZE+BLOCK_BORDER_SIZE, BLOCK_SIZE+BLOCK_BORDER_SIZE),
CV_32FC1);
cv::Mat eigenValues(cv::Size(BLOCK_SIZE+BLOCK_BORDER_SIZE, BLOCK_SIZE+BLOCK_BORDER_SIZE),
CV_32FC1);
cv::Mat eigenVectors(cv::Size(BLOCK_SIZE+BLOCK_BORDER_SIZE, BLOCK_SIZE+BLOCK_BORDER_SIZE),
CV_32FC1);
// Main loop that processes all blocks of the image
loop_height:for (int y=BLOCK_BORDER_SIZE; y<HEIGHT-BLOCK_BORDER_SIZE-BLOCK_SIZE;
y+=BLOCK_STEP_SIZE) {
loop_width:for (int x=BLOCK_BORDER_SIZE; x<WIDTH-BLOCK_BORDER_SIZE-BLOCK_SIZE;
x+=BLOCK_STEP_SIZE) {
// Compute derivatives using gaussian filters
cv::filter2D(src, tempIxx, CV_32F, xkernel, cv::Point(-1,-1));
cv::filter2D(src, tempIyy, CV_32F, ykernel, cv::Point(-1,-1));
cv::multiply(tempIxx, tempIxx, tempIxx);
cv::multiply(tempIyy, tempIyy, tempIyy);
cv::addWeighted(tempIxx,.5, tempIyy,.5, 0, tempIxy);
// Compute trace and determinant of the Hessian matrix
cv::reduce(tempIxx, Ixx, 0, cv::REDUCE_SUM);
cv::reduce(tempIyy, Iyy, 0, cv::REDUCE_SUM);
cv::reduce(tempIxy, Ixy, 0, cv::REDUCE_SUM);
cv::subtract(Ixx, Iyy, Ixy);
// Compute the corner response
double detM = ((double*)Ixx)[0]*((double*)Iyy)[0]-((double*)Ixy)[0]*((double*)Ixy)[0];
double traceM = ((double*)Ixx)[0]+((double*)Iyy)[0];
double R = detM - k*(traceM*traceM);
// Normalize response between [0,BIN_COUNT-1]
R = cv::min(R/(detM*traceM*sqrt(traceM)), 1.);
int cornerResponse = round(R*BIN_COUNT);
// Update response map for this block
response.at<float>(y, x) = (float)cornerResponse;
}
}
// Refine detected corners based on their eigenvectors
cv::cornerEigenValsAndVecs(response, BLOCK_SIZE, 3, eigenValues, eigenVectors);
// Filter out weak corners according to their strength and position
findCorners(eigenValues, eigenVectors, src.channels(), k, blockSize, wndRadius,
dst);
}
代码解读5.总结与展望
本文从相关概念和术语开始介绍，介绍了图像处理的基本概念和相关术语。然后详细介绍了计算机视觉中的几个核心算法——Sobel算子、Canny边缘检测、Harris角点检测和KLT光流跟踪。之后，讲述了FPGA上实现这些算法的步骤，并提供了相应的C++代码。通过这些例子，读者应该能对计算机视觉在FPGA上的实现有初步的认识。除此之外，还有许多计算机视觉的论文、研究报告和项目不断涌现，希望能有更多的同类文章出现。
另外，本文针对一些重难点的地方如NMS（Non Maximum Suppression），也提供了一些解决方案。虽然这些解决方案并不完美，但它们仍然值得参考。在FPGA上实现Sobel算子是非常具有挑战性的任务，但仍然值得尝试。至于后面三个算法，目前还没有看到有大的进展。