Center Line Extraction
Image Center-Line Extraction
A Center-Line Extraction Method
Recently I came across a method for extracting center-line segments from images. Its extraction result is shown in the figure below:
The algorithm is the Steger method based on the Hessian matrix. Its implementation on the CPU is shown below: The general workflow of the method includes Gaussian filtering, first-order partial derivatives, second-order partial derivatives, solving for the eigenvector corresponding to the largest eigenvalue of the Hessian matrix, and finally placing the image pixels that satisfy the conditions into the set of center-line pixels.
// lane center extraction func
// Steger C. An unbiased detector of curvilinear structures. IEEE Trans Pattern Anal Mach Intell[J].
// IEEE Transactions on Pattern Analysis & Machine Intelligence, 1998, 20(2):113-125.
// https://blog.csdn.net/Dangkie/article/details/78996761
cv::Mat StegerLine(const cv::Mat image) {
cv::Mat img;
//Gaussian filtering
cv::GaussianBlur(image, img, cv::Size(31, 11), 5, 2);
//first-order partial derivatives
cv::Mat m1, m2;
m1 = (cv::Mat_<float>(1, 2) << 1 ,-1); //x partial derivative
m2 = (cv::Mat_<float>(2, 1) << 1 ,-1); //y partial derivative
cv::Mat dx, dy;
filter2D(img, dx, CV_32FC1, m1);
filter2D(img, dy, CV_32FC1, m2);
//second-order partial derivatives
cv::Mat m3, m4, m5;
m3 = (cv::Mat_<float>(1, 3) << 1, -2, 1); //second-order x partial derivative
m4 = (cv::Mat_<float>(3, 1) << 1, -2, 1); //second-order y partial derivative
m5 = (cv::Mat_<float>(2, 2) << 1, -1, -1, 1); //second-order xy partial derivative
cv::Mat dxx, dyy, dxy;
filter2D(img, dxx, CV_32FC1, m3);
filter2D(img, dyy, CV_32FC1, m4);
filter2D(img, dxy, CV_32FC1, m5);
//Hessian matrix
std::vector<double> Pt;
for (int i = 0; i < img.cols; i++) {
for (int j = 0; j < img.rows; j++) {
if (image.at<float>(j, i) > 200 / 255.0) {
cv::Mat hessian(2, 2, CV_32FC1);
hessian.at<float>(0, 0) = dxx.at<float>(j, i);
hessian.at<float>(0, 1) = dxy.at<float>(j, i);
hessian.at<float>(1, 0) = dxy.at<float>(j, i);
hessian.at<float>(1, 1) = dyy.at<float>(j, i);
cv::Mat eValue;
cv::Mat eVectors;
eigen(hessian, eValue, eVectors);
double nx, ny;
if (fabs(eValue.at<float>(0, 0)) >= fabs(eValue.at<float>(1, 0))) //eigenvector corresponding to the largest eigenvalue
{
nx = eVectors.at<float>(0, 0);
ny = eVectors.at<float>(0, 1);
} else {
nx = eVectors.at<float>(1, 0);
ny = eVectors.at<float>(1, 1);
}
// std::cout<< "row, col, a, b, c, nx, ny: " << j << " "
// << i << " " << dxx.at<float>(j, i)<< " "
// << dxy.at<float>(j, i) << " " << dyy.at<float>(j, i)
// << " " << nx <<" " << ny <<" " << std::endl;
double t = -(nx * dx.at<float>(j, i) + ny * dy.at<float>(j, i)) /
(nx * nx * dxx.at<float>(j, i) + 2 * nx * ny * dxy.at<float>(j, i) +
ny * ny * dyy.at<float>(j, i));
if (fabs(t * nx) <= 0.5 && fabs(t * ny) <= 0.5) {
Pt.push_back(i);
Pt.push_back(j);
}
}
}
}
cv::Mat center_line = cv::Mat::zeros(image.rows, image.cols, CV_32FC1);
for (int k = 0; k < static_cast<int>(Pt.size() / 2); k++) {
int i = Pt[2 * k + 0];
int j = Pt[2 * k + 1];
center_line.ptr<float>(j)[i] = 1.0;
}
return center_line;
}
The above method runs slowly on the CPU. Its implementation is ported to run on the GPU, mainly in three steps:
- void GaussianBlur();
- void ComputeDeriv();
- void ComputeHessian();
The first step, Gaussian filtering, can be executed on global memory. Alternatively, the filter kernel parameters can be placed in constant memory to speed up memory access, or shared memory can be used to further accelerate the convolution. In addition, a 2D Gaussian filter is a separable filter that can be decomposed into a combination of a 1D row filter and a 1D column filter. Such an implementation already exists in NVIDIA’s official CUDA sample code; if you have already installed CUDA on your machine, it is located at “/usr/local/cuda/samples/3_Imaging/convolutionSeparable”.
Here we first present the code that performs the convolution using constant memory and global memory:
First, we need to generate the Gaussian kernel:
void GaussianMatrix2d(float *array, int radius_w, int radius_h,
float sigma_x_square, float sigma_y_square) {
int kernel_size_w = radius_w * 2 + 1;
float normalizationFactor = 0.0;
for (int i = -radius_w; i <= radius_w; i++) {
for (int j = -radius_h; j <= radius_h; j++) {
array[(j + radius_h) * kernel_size_w + i + radius_w] =
exp((-1.0 * i * i) / (2 * sigma_x_square) +
(-1.0 * j * j) / (2 * sigma_y_square));
normalizationFactor +=
array[(j + radius_h) * kernel_size_w + i + radius_w];
}
}
// Normalize, since the Gaussian is truncated, we would like to integrate to 1
for (int i = -radius_w; i <= radius_w; i++) {
for (int j = -radius_h; j <= radius_h; j++) {
array[(j + radius_h) * kernel_size_w + i + radius_w] /=
normalizationFactor;
}
}
}
const int kernel_w = 31;
const int kernel_h = 11;
__device__ __constant__ float d_cFilterKernel[kernel_w * kernel_h];
void setConvolutionKernel(float *h_Kernel)
{
cudaMemcpyToSymbol(d_cFilterKernel, h_Kernel, kernel_w * kernel_h * sizeof(float),
0, cudaMemcpyHostToDevice);
}
__global__ void GaussianBlurKernel_Constant(float* src, float* dst,
int width, int height) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
const int kernel_radius_w = (kernel_w - 1) / 2;
const int kernel_radius_h = (kernel_h - 1) / 2;
if(col < width && row < height) {
float convolved = 0.f; // temp variable for result
int locX, locY; // Local x and y
// Main loop. Iterates over all pixels in its domain, and then convolves the
// Gaussian matrix with the submatrix around the pixel.
for (int off_y = -kernel_radius_h ; off_y <= kernel_radius_h; off_y++)
{
for (int off_x = -kernel_radius_w; off_x <= kernel_radius_w; off_x++)
{
locX = min ( width-1, max(0, col + off_x));
locY = min ( height-1,max(0, row + off_y));
//convolved += src[locY* width + locX] *
// device_gaussian[(off_y + kernel_radius_h) * kernel_size_w + kernel_radius_w + off_x];
convolved += src[locY * width + locX] *
d_cFilterKernel[(off_y + kernel_radius_h) * kernel_w + off_x + kernel_radius_w];
//printf("%f ", d_cFilterKernel[(off_y + kernel_radius_h) * kernel_w + off_x + kernel_radius_w]);
}
}
dst[row * width + col] = min(convolved, 1.0);
}
}
For the second step, computing the first- and second-order image derivatives, since the filter kernel is fairly simple, global memory is used directly for the computation:
__global__ void Filter2D_dxy_kernel(float* src, float* dst, int width, int height) {
// dxy kernel
// m5 = (cv::Mat_<float>(2, 2) << 1, -1, -1, 1); //second-order xy partial derivative
int col = threadIdx.x + blockIdx.x * blockDim.x;
int row = threadIdx.y + blockIdx.y * blockDim.y;
if(col < (width - 1) && row < (height - 1)) {
float value = src[row * width + col] -
src[row * width + col + 1] -
src[(row + 1) * width + col] +
src[(row + 1) * width + col + 1];
dst[row * width + col] = value;
}
}
The third step requires computing the eigenvalues and eigenvectors of a 2x2 matrix; the implementation is given directly below:
__global__ void ComputeHessianKernel(float* src, float* dx, float* dy,
float* dxx, float* dyy, float* dxy,
float* result,
int width, int height) {
int col = threadIdx.x + blockIdx.x * blockDim.x;
int row = threadIdx.y + blockIdx.y * blockDim.y;
if(col < width && row < height && src[row * width + col] > 0.75) {
// cal eigen vector and eigen value
/*
* a b
* b c
*/
float a = dxx[row * width + col];
float b = dxy[row * width + col];
float c = dyy[row * width + col];
//eigenvector corresponding to the largest eigenvalue
float lambda1 = (a + c + sqrt((a - c) * (a - c) + 4 * b * b)) / 2.0;
float lambda2 = (a + c - sqrt((a - c) * (a - c) + 4 * b * b)) / 2.0;
double nx, ny;
if(fabs(b) < 1e-15) {
if(fabs(lambda1) >= fabs(lambda2)) {
nx = 1.0;
ny = 0.0;
} else {
nx = 0.0;
ny = 1.0;
}
} else {
if(fabs(lambda1) >= fabs(lambda2)) {
float tmp_nx = 1.0;
float tmp_ny = (lambda1 - a) / b;
nx = tmp_nx / sqrt(1.0 + tmp_ny * tmp_ny);
ny = tmp_ny / sqrt(1.0 + tmp_ny * tmp_ny);
} else {
float tmp_nx = 1.0;
float tmp_ny = (lambda2 - a) / b;
nx = tmp_nx / sqrt(1.0 + tmp_ny * tmp_ny);
ny = tmp_ny / sqrt(1.0 + tmp_ny * tmp_ny);
}
}
// if(row == 240 && col == 639) {
// printf("a, b, c, nx, ny: %f, %f, %f, %f, %f\n", a, b, c, nx, ny);
// }
double t = -(nx * dx[row * width + col] + ny * dy[row * width + col]) /
(nx * nx * dxx[row * width + col] + 2 * nx * ny * dxy[row * width + col] +
ny * ny * dyy[row * width + col]);
if(fabs(t * nx) <= 0.5 && fabs(t * ny) <= 0.5) {
result[row * width + col] = 1.0;
}
}
}
After the above optimization, a speedup of about 7-8x can be achieved, while the center-line extraction result remains consistent with the original.
Comments