基於OpenCL的影象積分圖演算法實現

積分圖的概念

影象積分圖演算法在影象特徵檢測中有著比較廣泛的應用,主要用於規則區域特徵值的計算。
積分圖的概念可用下圖表示:
這裡寫圖片描述

座標A(x,y)的積分圖是其左上角的所有畫素之和(圖中的陰影部分)。定義為:
這裡寫圖片描述

在上圖中,A(x,y)表示點(x,y)的積分圖;s(x,y)表示點(x,y)的y方向的所有原始影象之和。

積分圖演算法在CPU上的序列實現

在CPU上序列實現積分圖計算的典型程式碼如下:

    /*
* 標準的積分圖演算法(cpu)
* 返回積分圖矩陣物件
* is_square為true時為積方圖物件
*/
template<typename _E=E
,typename _ENABLE=typename std::enable_if<sizeof(_E)<=sizeof(cl_ulong)>::type>
integral_matrix integral_cpu(bool is_square)const {
throw_if(this->get_row_stride()*this->height!=this->v.size())
auto integ_mat = integral_matrix(this->width,this->height);
// 行寬度
const auto row_stride = this->get_row_stride();
if (!integ_mat.v.size())
integ_mat.v = std::vector<integral_matrix::element_type>(row_stride * this->height);
auto last_line  = integ_mat.v.data(); // 積分圖上一行指標
auto cur_line   = last_line;            // 積分圖當前行指標
auto src_line   = this->v.data();       // 原圖當前行指標
// 計算第一行字首和
prefix_sum(src_line, cur_line, this->width,is_square);
src_line  = row_stride; // 積分圖當前行指標步進一行
cur_line  = row_stride;
integral_matrix::element_type line_sum; // 積分圖當前行所有元素之和
typename std::decay<decltype(this->height)>::type y;
typename std::decay<decltype(this->width )>::type x;
// 從第二行開始計算積分圖
if(is_square){
for (y = 1; y < this->height;     y,
src_line     = row_stride,
cur_line     = row_stride,
last_line = row_stride) {
line_sum = 0;
for ( x = 0; x < this->width;   x) {
line_sum  = src_line[x]*src_line[x];
cur_line[x] = line_sum   last_line[x];
}
}
}else{
for (y = 1; y < this->height;     y,
src_line     = row_stride,
cur_line     = row_stride,
last_line = row_stride) {
line_sum = 0;
for ( x = 0; x < this->width;   x) {
line_sum  = src_line[x];
cur_line[x] = line_sum   last_line[x];
}
}
}
return std::move(integ_mat);//返回積分圖物件
}

字首和(prefix sum)

說到積分圖,就得引入字首和(prefix sum)的概念:
給定一個陣列A[1..n],字首和陣列prefix_sum[1..n]定義為:prefix_sum[i] = A[0] A1 … A[i-1];
例如:A[5,6,7,8] –> prefix_sum[5,11,18,26]
prefix_sum[0] =A[0] ;
prefix_sum1 =A[0] A1 ;
prefix_sum2 =A[0] A1 A2 ;
prefix_sum3 =A[0] A1 A2 A3 ;

下面是字首和陣列計算的典型程式碼,非常簡單

/* 字首和計算模板函式 is_square為true時計算平方和*/
template<typename _T1,typename _T2,
typename _ENABLE=typename std::enable_if<sizeof(_T1)<=sizeof(_T2)>::type>
inline static void prefix_sum(const _T1 *const src,_T2 *const dst,size_t size,bool is_square){
assert(nullptr!=src&&nullptr!=dst);
if(is_square){
dst[0] = src[0]*src[0];
for( size_t i=1; i<size;   i)   dst[i]=(_T2)(src[i])*(_T2)(src[i]) dst[i-1];
}
else{
dst[0] = src[0];
for( size_t i=1; i<size;   i)   dst[i]=(_T2)(src[i]) dst[i-1];
}
}

OpenCL並行實現

積分圖也可以用下面的公式(2)和公式(3)得出:
這裡寫圖片描述

從公式(2)和公式(3)可以看出,積分圖的演算法類似於字首和計算(prefix sum)

對於只有一行的畫素的影象,它的積分圖就是其字首和陣列

所以,如果要用OpenCL平行計算影象矩陣A的積分圖,可以把積分圖演算法分拆成兩個步驟:

  1. 首先計算矩陣A在x方向的字首和矩陣A1
  2. 然後再在計算矩陣A1在y方向字首和矩陣A2,A2就是影象矩陣A的積分圖矩陣

在OpenCL實現中為了提高記憶體訪問效能,計算矩陣A1在y方向字首和矩陣的時候,通常先將矩陣A1轉置,然後再進行計算x方向的字首和。

所以OpenCL具體實現的時候,分為下面4步

  1. 計算矩陣A在x方向的字首和矩陣A1
  2. A1轉置
  3. 計算矩陣A1在x方向的字首和矩陣A2
  4. A2轉置
    也就是說,基於OpenCL的積分圖演算法最終被分解為兩次x方向字首和計算和2次矩陣轉置

下面是主機端的部分實現程式碼:

/*
* 計算影象的積分圖/積方圖,
* 返回積分圖矩陣物件
* is_square為true時為積方圖物件
*/
gray_matrix_cl::integral_matrix gray_matrix_cl::integral(bool is_square) const {
auto integral_mat1=to_matrix_cl().prefix_sum_line<cl_uint>(std::string(KERNEL_NAME(prefix_sum_line)"_uchar_uint"),is_square);//執行kernel,計算x方向字首和
//原圖為灰度影象矩陣,所以原圖元素型別為uchar,生成的字首和矩陣integral_mat1的元素型別為uint
integral_mat1= integral_mat1.transpose(KERNEL_NAME(matrix_transpose)"_uint");//執行kernel,integral_mat1矩陣轉置
auto integral_mat2=integral_mat1.prefix_sum_line<cl_ulong>(std::string(KERNEL_NAME(prefix_sum_line)"_uint_ulong"),false);//執行kernel,計算x方向字首和
// 源矩陣integral_mat1的元素型別為uint,生成的字首和矩陣元素型別為ulong
integral_mat2= integral_mat2.transpose(KERNEL_NAME(matrix_transpose)"_ulong");//執行kernel,integral_mat2矩陣轉置
return std::move(integral_mat2);//返回積分圖物件
}
/*
* 計算矩陣中的每一行的字首和(prefix_sum)
* 結果輸出到int_mat
*/
template< typename DST_E,
typename _CL_TYPE=CL_TYPE,typename _E=E,
typename _ENABLE=typename   std::enable_if<std::is_scalar<_E>::value
&&std::is_scalar<DST_E>::value
&&sizeof(DST_E)>=sizeof(_E)
&&std::is_base_of<cl::Buffer,_CL_TYPE>::value
>::type>
matrix_cl<DST_E,cl::Buffer>
prefix_sum_line(const std::string &kernel_name,bool is_square)const{
this->check_cl_mem_obj(SOURCE_AT);
auto context=dynamic_cast<const cl::Memory&>(this->cl_mem_obj).getInfo<CL_MEM_CONTEXT>();
matrix_cl<DST_E,cl::Buffer> integral_mat(width, height,0,context);
cl::CommandQueue q = cl::CommandQueue(context);
run_kernel(global_facecl_context.getKernel(kernel_name)
, cl::EnqueueArgs(q, { 1, height })//每個work-iteam處理1行
, false
, *this
, integral_mat
, width
, get_row_stride()
, is_square?CL_TRUE:CL_FALSE
);
return std::move(integral_mat);
}
/*
* 矩陣轉置
*/
template<typename _CL_TYPE=CL_TYPE>
typename std::enable_if<std::is_base_of<cl::Buffer, _CL_TYPE>::value,self_type>::type
transpose(const std::string &kernel_name)const{
this->check_cl_mem_obj(SOURCE_AT);
auto context = dynamic_cast<const cl::Memory&>(this->cl_mem_obj).getInfo<CL_MEM_CONTEXT>();
self_type dst_mat(height, width, align_v, context);
dst_mat.align_v = this->align; // 記錄轉置前的水平對齊值 
cl::CommandQueue queue = cl::CommandQueue(context);
run_kernel(global_facecl_context.getKernel(kernel_name)
, cl::EnqueueArgs(queue,{ 1, height })//每個work-iteam處理1行
, false
, *this
, dst_mat
, width
, this->get_row_stride()
, dst_mat.get_row_stride()
);
return std::move(dst_mat);
}

上面程式碼中用到的run_kernel函式參見我的部落格《opencl:cl::make_kernel的進化》

下面是上面程式碼中執行的kernel函式prefix_sum_line的程式碼,每個work-item處理一行資料,實現的功能很簡單,就是計算矩陣中一行資料的字首和(prefix sum),
為減少對global記憶體的訪問,kernel函式中用到了local memory(程式碼中的local_block陣列)來暫存每行的部分資料。local_block陣列的大小在編譯內kernel程式碼時由編譯器提供,參見我的部落格《opencl::kernel中獲取local memory size》

 ///////////////////////////////////////////////////////////////////////////////
//! @file   : prefix_sum_line.cl
//! @date   : 2016/03/04
//! @author: guyadong
//! @brief  : Calculates the integral sum scan of an image
////////////////////////////////////////////////////////////////////////////////
#ifndef CL_DEVICE_LOCAL_MEM_SIZE //local memory的大小,由編譯器提供
#error not defined CL_DEVICE_LOCAL_MEM_SIZE by complier with options -D
#endif
#ifndef SRC_TYPE //源矩陣資料型別 uchar,uinit,ulong.....
#error not defined SRC_TYPE by complier with options -D
#endif
#ifndef DST_TYPE //目標矩陣資料型別 uchar,uinit,ulong.....
#error not defined DST_TYPE by complier with options -D
#endif
#define LOCAL_BUFFER_SIZE (CL_DEVICE_LOCAL_MEM_SIZE/sizeof(DST_TYPE))//編譯時確定local buffer陣列的大小
#define _KERNEL_NAME(s,d) prefix_sum_line_##s##_##d
#define KERNEL_NAME(s,d) _KERNEL_NAME(s,d)
// kernel function的名字在編譯期根據SRC_TYPE 和DST_TYPE新增型別字尾
///////////////////////////////////////////////////////////////////////////////
//! @brief  :   Calculates the prefix sum for each line of an image if is_square is CL_FALSE,
//                  Calculates the prefix sum of is_square if is_square is CL_TRUE, 
////////////////////////////////////////////////////////////////////////////////
__kernel void KERNEL_NAME(SRC_TYPE,DST_TYPE)( __global SRC_TYPE *sourceImage, __global DST_TYPE * dest, int width, int width_step,int is_square ){
__local DST_TYPE local_block[ LOCAL_BUFFER_SIZE ];
const int line_index = get_global_id(1)*width_step;// 計算當前行的起始位置
__global SRC_TYPE * const src_ptr   = line_index   sourceImage;// 源矩陣的起始指標
__global DST_TYPE * const dst_ptr   = line_index   dest;    // 目標矩陣的起始指標
__global SRC_TYPE * block_src_ptr   = src_ptr;  
__global DST_TYPE * block_dst_ptr   = dst_ptr;  
int block_size = 0; // 塊大小
DST_TYPE last_sum=0;// 上一塊陣列的字首和
// 將一行資料按local_block陣列的大小來分塊處理
for( int start_x = 0 ; start_x < width ; 
start_x                  = LOCAL_BUFFER_SIZE,
block_src_ptr        = LOCAL_BUFFER_SIZE,
block_dst_ptr    = LOCAL_BUFFER_SIZE,
last_sum             = local_block[block_size -1] ){
block_size = min( (int)LOCAL_BUFFER_SIZE, width - start_x );
// compute prefix sum of a block with local memory      
if(is_square){
local_block[0] = last_sum   ((DST_TYPE)block_src_ptr[0])*((DST_TYPE)block_src_ptr[0]);
for( int i=1; i<block_size;   i)    local_block[i]=((DST_TYPE)block_src_ptr[i])*((DST_TYPE)block_src_ptr[i]) local_block[i-1];
}else{
local_block[0] = last_sum   (DST_TYPE)block_src_ptr[0];
for( int i=1; i<block_size;   i)    local_block[i]=block_src_ptr[i] local_block[i-1];
}       
// copy local_block to dest
for(int i = 0 ; i < block_size;   i){
block_dst_ptr[i]=local_block[i];
}       
}
}

矩陣轉置的kernel程式碼,每個work-item處理一行資料

 ///////////////////////////////////////////////////////////////////////////////
//! @file   : matrix_transpose.cl
//! @date   : 2016/03/04
//! @author: guyadong
//! @brief  : matrix transpose
////////////////////////////////////////////////////////////////////////////////
#ifndef SRC_TYPE 
#error not defined SRC_TYPE by complier with options -D
#endif
#define _KERNEL_NAME(s) matrix_transpose_##s
#define KERNEL_NAME(s) _KERNEL_NAME(s) 
// kernel function的名字在編譯期根據SRC_TYPE 加一個型別字尾
__kernel void KERNEL_NAME(SRC_TYPE)( __global SRC_TYPE *matrix_src,__global SRC_TYPE *matrix_dst, int width, int src_width_step,int dst_width_step ){
const int y                 = get_global_id(1);
__global SRC_TYPE * src_ptr = matrix_src   y*src_width_step;
for( int x = 0; x < width;   x,  src_ptr ){
matrix_dst[  x*dst_width_step   y ] = *src_ptr;
} 
}

補充:
後來我對這個演算法進行了改進,參見我的後續博文《基於OpenCL的影象積分圖演算法改進》

參考文章

《AdaBoost人臉檢測演算法1(轉)》
《基於OpenCL的影象積分圖演算法優化研究》