图像模板匹配算法的研究和实现
这次写一下算法方面的，图像处理中模板匹配算法的研究和实现。
一： 首先我们先上一下模板匹配的理论及其公式描述：
模板匹配是通过在输入图像上滑动模板图像块对实际的图像块和输入图像进行匹配，并且可以利用函数cvMinMaxLoc()找到最佳匹配的位置。例如在工业应用中，可以锁定图像中零部件的位置，并根据具体的位置，进行具体的处理。匹配的过程中可以使用不同的method，通过最合适的method，进行最合适的匹配。
MatchTemplate 比较模板和重叠的图像区域
void cvMatchTemplate( const CvArr* image,const CvArr* templ, CvArr* result, int method );
- image 欲搜索的图像。它应该是单通道、8-比特或32-比特 浮点数图像
- templ 搜索模板，不能大于输入图像，且与输入图像具有一样的数据类型
- result 比较结果的映射图像。单通道、32-比特浮点数. 如果图像是 W×H 而 templ 是 w×h ，则 result 一定是 (W-w+1）×(H-h+1）.
- method 指定匹配方法：
函数 cvMatchTemplate 与函数 cvCalcBackProjectPatch 类似。它滑动过整个图像 image, 用指定方法比较 templ 与图像尺寸为 w×h 的重叠区域，并且将比较结果存到 result 中。 下面是不同的比较方法，可以使用其中的一种 (I 表示图像，T - 模板, R - 结果. 模板与图像重叠区域 x'=0..w-1, y'=0..h-1 之间求和):
OpenCV模板匹配算法.可用的方法有6个:
平方差匹配method=CV_TM_SQDIFF，这类方法利用平方差来进行匹配,最好匹配为0.匹配越差,匹配值越大.
标准平方差匹配method=CV_TM_SQDIFF_NORMED
相关匹配method=CV_TM_CCORR，这类方法采用模板和图像间的乘法操作,所以较大的数表示匹配程度较高,0标识最坏的匹配效果.
标准相关匹配method=CV_TM_CCORR_NORMED
相关系数匹配method=CV_TM_CCOEFF。这类方法将模版对其均值的相对值与图像对其均值的相关值进行匹配,1表示完美匹配,-1表示糟糕的匹配,0表示没有任何相关性(随机序列).
在这里
标准相关系数匹配method=CV_TM_CCOEFF_NORMED
通常,随着从简单的测量(平方差)到更复杂的测量(相关系数),我们可获得越来越准确的匹配(同时也意味着越来越大的计算代价).最好的办法是对所有这些设置多做一些测试实验,以便为自己的应用选择同时兼顾速度和精度的最佳方案.
通过使用 cvMinMaxLoc函数，我们确定结果矩阵 R 的最大值和最小值的位置。函数中的参数有:
二： opencv 实现模板匹配
这里调用 OPENCV 自带的模板匹配函数实现模板匹配 .具体代码如下：
//outDataMat 中相邻扫描线进行匹配，模板矩阵大小为100，目标矩阵大小为200，结果矩阵大小为101
// 参数说明： winSize ：100 ； stepSize ：10 ；outDataMat ：300*8196
void CDisplacement::displacementAlgorithm( CvMat * outDataMat，int winSize ， int stepSize )
{
int multiWin = 2; //大窗口对小窗口的倍数
int winNum = static_cast<int>(outDataMat->cols - multiWin * winSize) / stepSize; //计算需要匹配的段数
CvMat *disMat = cvCreateMat(outDataMat->rows - 1, winNum, CV_32FC1); //用于保存位移的矩阵，匹配得出的整数位移经过插值后应该为浮点数
CvMat *templateMat = cvCreateMatHeader(1, winSize, CV_32FC1); //用于匹配的模板
CvMat *objectMat = cvCreateMatHeader(1, multiWin * winSize, CV_32FC1); //用于匹配的目标
CvMat *resultMat = cvCreateMat(1, (multiWin - 1) * winSize + 1, CV_32FC1); //用于存放匹配结果
double min, max; //匹配结果中的最大值
CvPoint max_loc; //匹配结果最大值所在位置
double pmax, nmax; //最大值的前后两个值
double displacement; //存储位移值
float temp_result[101];
for (int i = 0; i < disMat->rows; ++i) //求位移
{
for (int j = 0; j < disMat->cols; ++j)
{
/* ->data.ptr表示数据的起始位置 */
templateMat->data.ptr = outDataMat->data.ptr + i * outDataMat->step
+ static_cast<int>(sizeof(float) * ((multiWin - 1) * winSize / 2 + j * stepSize));
objectMat->data.ptr = outDataMat->data.ptr + (i + 1) * outDataMat->step + static_cast<int>(sizeof(float) * (j * stepSize)/* + 0.5*/); //相邻扫描线
cvMatchTemplate(objectMat, templateMat, resultMat, CV_TM_CCORR_NORMED); //匹配
cvMinMaxLoc(resultMat, &min, &max, NULL, &max_loc); //查找最大值
pmax = *static_cast<float*>(static_cast<void*>(resultMat->data.ptr + sizeof(float)*max_loc.x - sizeof(float))); //取最大值的前一个值
nmax = *static_cast<float*>(static_cast<void*>(resultMat->data.ptr + sizeof(float)*max_loc.x + sizeof(float))); //取最大值的后一个值
displacement = (multiWin - 1) * winSize / 2 - max_loc.x - (pmax - nmax) / (2 * (pmax - 2 * max + nmax)); //插值得到位移
CV_MAT_ELEM(*disMat,float, i, j) = static_cast<float>(displacement); //结果存放到位移矩阵中
}
//test for wong 2016/06/24
SaveDispDataFileFilter_1("weiyi_cpu.dat", disMat); //保存位移值
}
三： CUDA实现模板匹配
这里采用CUDA实现模板匹配功能
void CudaMain::computeDisplacement_cuda(CvMat* filtOutMat, int multiWin, int winSize, int stepSize, CvMat*outputMat){
int WinNum = (filtOutMat->cols - multiWin*winSize) / stepSize; // 一维位移矩阵
Complex* hInput = (Complex*)filtOutMat->data.fl; // 数据位置；
Complex*hOutput = (Complex*)outputMat->data.fl; // 数据位置；
cudaMemcpy(inputMat, hInput, filtOutMat->cols*filtOutMat->rows*sizeof(Complex), cudaMemcpyHostToDevice); // CPU-GPU
dim3 dBlock;
dim3 dThread;
dBlock.x = filtOutMat->rows - 1; // 输出矩阵行数 ,块数 299
dThread.x = WinNum; // 输出矩阵列数 , 线程数 799
templateMat*templateMatShare; // 模板内存在GPU分配
objectMat* objectMatShare; // 目标内存在GPU分配
resultMat*resultMatShare; // 匹配结果在GPU分配
Complex* min;
Complex* max;
int* max_location;
Complex* displacement;
cudaMalloc(&templateMatShare, dBlock.x*dThread.x* sizeof(templateMat)); //模板矩阵在GPU全局内存分配
cudaMalloc(&objectMatShare, dBlock.x*dThread.x* sizeof(objectMat)); //目标矩阵在GPU全局内存分配
cudaMalloc(&resultMatShare, dBlock.x*dThread.x* sizeof(resultMat)); //结果矩阵在GPU全局内存分配
cudaMalloc(&min, dBlock.x*dThread.x* sizeof(Complex)); // min在GPU全局内存分配
cudaMalloc(&max, dBlock.x*dThread.x* sizeof(Complex)); // max在GPU全局内存分配
cudaMalloc(&max_location, dBlock.x*dThread.x* sizeof(int)); // max_location在GPU全局内存分配
cudaMalloc(&displacement, dBlock.x*dThread.x* sizeof(Complex)); // max_location在GPU全局内存分配
//求位移矩阵
displacement_api_cuda << < dBlock, dThread >> > (inputMat, filtOutMat->rows, filtOutMat->cols, multiWin, winSize, stepSize, templateMatShare, objectMatShare, resultMatShare, min, max, max_location, displacement);
cudaThreadSynchronize();
cudaMemcpy(hOutput, displacement, sizeof(Complex)*outputMat->cols*outputMat->rows, cudaMemcpyDeviceToHost); //拷贝GPU处理后数据到 CPU
SaveDataFile("weiyi_gpu.dat", outputMat); //保存数据
cudaFree(templateMatShare);
cudaFree(objectMatShare);
cudaFree(resultMatShare);
cudaFree(min);
cudaFree(max);
cudaFree(max_location);
}
其中核函数：
__global__ void displacement_api_cuda(Complex*disInputCuda, int rows, int cols, int multiWin, int winSize, int stepSize, templateMat*templateMatShare, objectMat* objectMatShare, resultMat*resultMatShare, Complex*min, Complex*max, int*max_location, Complex* displacement ) {
int out_offset = blockIdx.x *blockDim.x + threadIdx.x;
int bid = blockIdx.x ;
int tid = threadIdx.x;
Complex*templateMatID; //ID
Complex*objectMatID; //ID
templateMatID = (Complex*)(disInputCuda + blockIdx.x*cols + (multiWin - 1) * winSize / 2 + threadIdx.x * stepSize);
for (int i = 0; i < 100;i++) {
if (i < 64) {
templateMatShare[out_offset].tempData.elem[i]= *(templateMatID + i);
}
else
templateMatShare[out_offset].tempData.atom[i-64] = *(templateMatID + i);
}
objectMatID = (Complex*)(disInputCuda + (blockIdx.x + 1)*cols + threadIdx.x * stepSize);
for (int i = 0; i < 200; i++) {
if (i<64)
objectMatShare[out_offset].objData.elem_0[i] = *(objectMatID + i);
else if (i<128)
objectMatShare[out_offset].objData.elem_1[i - 64] = *(objectMatID + i);
else if (i<192)
objectMatShare[out_offset].objData.elem_2[i - 128] = *(objectMatID + i);
else
objectMatShare[out_offset].objData.atom[i - 192] = *(objectMatID + i);
}
for (int i = 0; i < 101; i++) {
if (i<64)
resultMatShare[out_offset].resData.elem[i] = 0;
else
resultMatShare[out_offset].resData.atom[i-64] = 0;
}
for (int i = 0; i < 192; i++) {
if (i<64)
objectMatShare[out_offset].elem_0[i] = *(objectMatID+i);
else if (i<128)
objectMatShare[out_offset].elem_1[i-64] = *(objectMatID + i);
else
objectMatShare[out_offset].elem_2[i - 128] = *(objectMatID + i);
}
for (int j = 0; j < 8; j++) {
objectMatShare[out_offset].atom[j] = *(templateMatID + j + 192);
xcorr_cuda(templateMatShare[out_offset].tempData.elem, objectMatShare[out_offset].objData.elem_0, resultMatShare[out_offset].resData.elem);
minMax_cuda(resultMatShare[out_offset].resData.elem, &min[out_offset], &max[out_offset], &max_location[out_offset]);
interp_cuda(resultMatShare[out_offset].resData.elem, &max_location[out_offset], &max[out_offset], &multiWin, &winSize, &displacement[out_offset]);
__syncthreads();
}
__device__ void xcorr_cuda(const Complex* templateMat_startID, const Complex* objectMat_startID, Complex*resultMat_startID) {
for (int i = 0; i < 101; i++) {
Complex sum_object = 0;
Complex frac_object = 0;
Complex pow_template = 0;
Complex pow_object = 0;
Complex result = 0;
//sum_object
for (int j = 0; j < 100; j++) {
sum_object += *(objectMat_startID + i + j);
}
// ave
Complex ave_object = sum_object / 100;
//fraction
for (int j = 0; j < 100; j++) {
Complex tmp = *(templateMat_startID + j) * (*(objectMat_startID + i + j) - ave_object);
frac_object += tmp;
}
//pow temp
for (int j = 0; j < 100; j++) {
pow_template += *(templateMat_startID + j) * *(templateMat_startID + j);
}
//pow objectMat
for (int j = 0; j < 100; j++) {
pow_object += *(objectMat_startID + i + j)* * (objectMat_startID + i + j);
}
//result
result = sqrt(pow_template*pow_object);
//output
*(resultMat_startID + i) = frac_object / result;
}
}
__device__ void minMax_cuda(Complex*resultMat_startID, Complex* min_value, Complex* max_value, int * max_location) {
//求最大值及位置
for (int i = 0; i < 101; i++) {
if (*(resultMat_startID + i) >= *max_value) {
*max_location = i;
*max_value = *(resultMat_startID + i);
}
}
}
__device__ void interp_cuda(Complex*resultMat_startID, int * max_loc, Complex*max_value, int * multiWin, int * winSize, Complex* displace) {
Complex*pre = (Complex*)resultMat_startID + *max_loc - 1;
Complex*next = (Complex*)resultMat_startID + *max_loc + 1;
*displace = (*multiWin - 1) * *winSize / 2 - *max_loc - (*pre - *next) / (2 * (*pre - 2 * *max_value + *next));
}
四： 结果比较
这里贴一下采用CPU计算的和GPU计算的结果吧：
采用CPU计算的结果：
0.009192 0.008164 0.007095 0.011042 0.011536 0.011247 0.016015 0.010360 0.001595 -0.000845 -0.002765 -0.000039 0.007204 0.005229 0.002261 0.003771 0.002982 0.010227 0.017885 0.016466 0.017094 0.017613 -0.000961 -0.003701 0.015722 0.016839 0.013439 0.007723 0.002980 0.010488 0.012283 0.006072 0.001040 -0.003270 -0.003462 0.004596 0.008017 0.010673 0.020252 0.002328 -0.001016 0.009833 0.011129 0.009137 0.012297 0.014499 0.017363 0.013004 0.003839 0.008114 0.014437 0.007309 0.004949 0.004221 0.002466 -0.001996 -0.015309 -0.020181 -0.013732 -0.012274 0.003202
采用CUDA计算的结果：
0.013173 0.015335 0.011267 0.008420 0.004566 0.006342 0.014983 0.010390 0.001600 -0.000754 -0.002277 0.001073 0.007668 0.005285 0.002635 0.004263 0.003369 0.010357 0.018863 0.017140 0.017338 0.019080 0.000751 -0.004420 0.017763 0.020394 0.015874 0.007988 0.000732 0.010925 0.016618 0.013930 0.008934 -0.004655 -0.013801 -0.003829 0.004890 0.014446 0.028133 -0.001455 -0.005051 0.008987 0.009782 0.003444 0.006536 0.009789 0.012227 0.010090 0.003550 -0.000100 0.000534 -0.005920 -0.008974 -0.008588 -0.002132 -0.002369 -0.017581 -0.032293 -0.021036 -0.009534 0.013279 0.025665 0.012573 0.006627
浙公网安备 33010602011771号