1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 5 // By downloading, copying, installing or using the software you agree to this license. 6 // If you do not agree to this license, do not download, install, 7 // copy or use the software. 8 // 9 // 10 // Intel License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000, Intel Corporation, all rights reserved. 14 // Third party copyrights are property of their respective owners. 15 // 16 // Redistribution and use in source and binary forms, with or without modification, 17 // are permitted provided that the following conditions are met: 18 // 19 // * Redistribution's of source code must retain the above copyright notice, 20 // this list of conditions and the following disclaimer. 21 // 22 // * Redistribution's in binary form must reproduce the above copyright notice, 23 // this list of conditions and the following disclaimer in the documentation 24 // and/or other materials provided with the distribution. 25 // 26 // * The name of Intel Corporation may not be used to endorse or promote products 27 // derived from this software without specific prior written permission. 28 // 29 // This software is provided by the copyright holders and contributors "as is" and 30 // any express or implied warranties, including, but not limited to, the implied 31 // warranties of merchantability and fitness for a particular purpose are disclaimed. 32 // In no event shall the Intel Corporation or contributors be liable for any direct, 33 // indirect, incidental, special, exemplary, or consequential damages 34 // (including, but not limited to, procurement of substitute goods or services; 35 // loss of use, data, or profits; or business interruption) however caused 36 // and on any theory of liability, whether in contract, strict liability, 37 // or tort (including negligence or otherwise) arising in any way out of 38 // the use of this software, even if advised of the possibility of such damage. 39 // 40 //M*/ 41 42 #include "_cv.h" 43 44 void 45 icvCrossCorr( const CvArr* _img, const CvArr* _templ, CvArr* _corr, CvPoint anchor ) 46 { 47 const double block_scale = 4.5; 48 const int min_block_size = 256; 49 CvMat* dft_img[CV_MAX_THREADS] = {0}; 50 CvMat* dft_templ = 0; 51 void* buf[CV_MAX_THREADS] = {0}; 52 int k, num_threads = 0; 53 54 CV_FUNCNAME( "icvCrossCorr" ); 55 56 __BEGIN__; 57 58 CvMat istub, *img = (CvMat*)_img; 59 CvMat tstub, *templ = (CvMat*)_templ; 60 CvMat cstub, *corr = (CvMat*)_corr; 61 CvSize dftsize, blocksize; 62 int depth, templ_depth, corr_depth, max_depth = CV_32F, 63 cn, templ_cn, corr_cn, buf_size = 0, 64 tile_count_x, tile_count_y, tile_count; 65 66 CV_CALL( img = cvGetMat( img, &istub )); 67 CV_CALL( templ = cvGetMat( templ, &tstub )); 68 CV_CALL( corr = cvGetMat( corr, &cstub )); 69 70 if( CV_MAT_DEPTH( img->type ) != CV_8U && 71 CV_MAT_DEPTH( img->type ) != CV_16U && 72 CV_MAT_DEPTH( img->type ) != CV_32F ) 73 CV_ERROR( CV_StsUnsupportedFormat, 74 "The function supports only 8u, 16u and 32f data types" ); 75 76 if( !CV_ARE_DEPTHS_EQ( img, templ ) && CV_MAT_DEPTH( templ->type ) != CV_32F ) 77 CV_ERROR( CV_StsUnsupportedFormat, 78 "Template (kernel) must be of the same depth as the input image, or be 32f" ); 79 80 if( !CV_ARE_DEPTHS_EQ( img, corr ) && CV_MAT_DEPTH( corr->type ) != CV_32F && 81 CV_MAT_DEPTH( corr->type ) != CV_64F ) 82 CV_ERROR( CV_StsUnsupportedFormat, 83 "The output image must have the same depth as the input image, or be 32f/64f" ); 84 85 if( (!CV_ARE_CNS_EQ( img, corr ) || CV_MAT_CN(templ->type) > 1) && 86 (CV_MAT_CN( corr->type ) > 1 || !CV_ARE_CNS_EQ( img, templ)) ) 87 CV_ERROR( CV_StsUnsupportedFormat, 88 "The output must have the same number of channels as the input (when the template has 1 channel), " 89 "or the output must have 1 channel when the input and the template have the same number of channels" ); 90 91 depth = CV_MAT_DEPTH(img->type); 92 cn = CV_MAT_CN(img->type); 93 templ_depth = CV_MAT_DEPTH(templ->type); 94 templ_cn = CV_MAT_CN(templ->type); 95 corr_depth = CV_MAT_DEPTH(corr->type); 96 corr_cn = CV_MAT_CN(corr->type); 97 max_depth = MAX( max_depth, templ_depth ); 98 max_depth = MAX( max_depth, depth ); 99 max_depth = MAX( max_depth, corr_depth ); 100 if( depth > CV_8U ) 101 max_depth = CV_64F; 102 103 if( img->cols < templ->cols || img->rows < templ->rows ) 104 CV_ERROR( CV_StsUnmatchedSizes, 105 "Such a combination of image and template/filter size is not supported" ); 106 107 if( corr->rows > img->rows + templ->rows - 1 || 108 corr->cols > img->cols + templ->cols - 1 ) 109 CV_ERROR( CV_StsUnmatchedSizes, 110 "output image should not be greater than (W + w - 1)x(H + h - 1)" ); 111 112 blocksize.width = cvRound(templ->cols*block_scale); 113 blocksize.width = MAX( blocksize.width, min_block_size - templ->cols + 1 ); 114 blocksize.width = MIN( blocksize.width, corr->cols ); 115 blocksize.height = cvRound(templ->rows*block_scale); 116 blocksize.height = MAX( blocksize.height, min_block_size - templ->rows + 1 ); 117 blocksize.height = MIN( blocksize.height, corr->rows ); 118 119 dftsize.width = cvGetOptimalDFTSize(blocksize.width + templ->cols - 1); 120 if( dftsize.width == 1 ) 121 dftsize.width = 2; 122 dftsize.height = cvGetOptimalDFTSize(blocksize.height + templ->rows - 1); 123 if( dftsize.width <= 0 || dftsize.height <= 0 ) 124 CV_ERROR( CV_StsOutOfRange, "the input arrays are too big" ); 125 126 // recompute block size 127 blocksize.width = dftsize.width - templ->cols + 1; 128 blocksize.width = MIN( blocksize.width, corr->cols ); 129 blocksize.height = dftsize.height - templ->rows + 1; 130 blocksize.height = MIN( blocksize.height, corr->rows ); 131 132 CV_CALL( dft_templ = cvCreateMat( dftsize.height*templ_cn, dftsize.width, max_depth )); 133 134 num_threads = cvGetNumThreads(); 135 136 for( k = 0; k < num_threads; k++ ) 137 CV_CALL( dft_img[k] = cvCreateMat( dftsize.height, dftsize.width, max_depth )); 138 139 if( templ_cn > 1 && templ_depth != max_depth ) 140 buf_size = templ->cols*templ->rows*CV_ELEM_SIZE(templ_depth); 141 142 if( cn > 1 && depth != max_depth ) 143 buf_size = MAX( buf_size, (blocksize.width + templ->cols - 1)* 144 (blocksize.height + templ->rows - 1)*CV_ELEM_SIZE(depth)); 145 146 if( (corr_cn > 1 || cn > 1) && corr_depth != max_depth ) 147 buf_size = MAX( buf_size, blocksize.width*blocksize.height*CV_ELEM_SIZE(corr_depth)); 148 149 if( buf_size > 0 ) 150 { 151 for( k = 0; k < num_threads; k++ ) 152 CV_CALL( buf[k] = cvAlloc(buf_size) ); 153 } 154 155 // compute DFT of each template plane 156 for( k = 0; k < templ_cn; k++ ) 157 { 158 CvMat dstub, *src, *dst, temp; 159 CvMat* planes[] = { 0, 0, 0, 0 }; 160 int yofs = k*dftsize.height; 161 162 src = templ; 163 dst = cvGetSubRect( dft_templ, &dstub, cvRect(0,yofs,templ->cols,templ->rows)); 164 165 if( templ_cn > 1 ) 166 { 167 planes[k] = templ_depth == max_depth ? dst : 168 cvInitMatHeader( &temp, templ->rows, templ->cols, templ_depth, buf[0] ); 169 cvSplit( templ, planes[0], planes[1], planes[2], planes[3] ); 170 src = planes[k]; 171 planes[k] = 0; 172 } 173 174 if( dst != src ) 175 cvConvert( src, dst ); 176 177 if( dft_templ->cols > templ->cols ) 178 { 179 cvGetSubRect( dft_templ, dst, cvRect(templ->cols, yofs, 180 dft_templ->cols - templ->cols, templ->rows) ); 181 cvZero( dst ); 182 } 183 cvGetSubRect( dft_templ, dst, cvRect(0,yofs,dftsize.width,dftsize.height) ); 184 cvDFT( dst, dst, CV_DXT_FORWARD + CV_DXT_SCALE, templ->rows ); 185 } 186 187 tile_count_x = (corr->cols + blocksize.width - 1)/blocksize.width; 188 tile_count_y = (corr->rows + blocksize.height - 1)/blocksize.height; 189 tile_count = tile_count_x*tile_count_y; 190 191 { 192 #ifdef _OPENMP 193 #pragma omp parallel for num_threads(num_threads) schedule(dynamic) 194 #endif 195 // calculate correlation by blocks 196 for( k = 0; k < tile_count; k++ ) 197 { 198 int thread_idx = cvGetThreadNum(); 199 int x = (k%tile_count_x)*blocksize.width; 200 int y = (k/tile_count_x)*blocksize.height; 201 int i, yofs; 202 CvMat sstub, dstub, *src, *dst, temp; 203 CvMat* planes[] = { 0, 0, 0, 0 }; 204 CvMat* _dft_img = dft_img[thread_idx]; 205 void* _buf = buf[thread_idx]; 206 CvSize csz = { blocksize.width, blocksize.height }, isz; 207 int x0 = x - anchor.x, y0 = y - anchor.y; 208 int x1 = MAX( 0, x0 ), y1 = MAX( 0, y0 ), x2, y2; 209 csz.width = MIN( csz.width, corr->cols - x ); 210 csz.height = MIN( csz.height, corr->rows - y ); 211 isz.width = csz.width + templ->cols - 1; 212 isz.height = csz.height + templ->rows - 1; 213 x2 = MIN( img->cols, x0 + isz.width ); 214 y2 = MIN( img->rows, y0 + isz.height ); 215 216 for( i = 0; i < cn; i++ ) 217 { 218 CvMat dstub1, *dst1; 219 yofs = i*dftsize.height; 220 221 src = cvGetSubRect( img, &sstub, cvRect(x1,y1,x2-x1,y2-y1) ); 222 dst = cvGetSubRect( _dft_img, &dstub, 223 cvRect(0,0,isz.width,isz.height) ); 224 dst1 = dst; 225 226 if( x2 - x1 < isz.width || y2 - y1 < isz.height ) 227 dst1 = cvGetSubRect( _dft_img, &dstub1, 228 cvRect( x1 - x0, y1 - y0, x2 - x1, y2 - y1 )); 229 230 if( cn > 1 ) 231 { 232 planes[i] = dst1; 233 if( depth != max_depth ) 234 planes[i] = cvInitMatHeader( &temp, y2 - y1, x2 - x1, depth, _buf ); 235 cvSplit( src, planes[0], planes[1], planes[2], planes[3] ); 236 src = planes[i]; 237 planes[i] = 0; 238 } 239 240 if( dst1 != src ) 241 cvConvert( src, dst1 ); 242 243 if( dst != dst1 ) 244 cvCopyMakeBorder( dst1, dst, cvPoint(x1 - x0, y1 - y0), IPL_BORDER_REPLICATE ); 245 246 if( dftsize.width > isz.width ) 247 { 248 cvGetSubRect( _dft_img, dst, cvRect(isz.width, 0, 249 dftsize.width - isz.width,dftsize.height) ); 250 cvZero( dst ); 251 } 252 253 cvDFT( _dft_img, _dft_img, CV_DXT_FORWARD, isz.height ); 254 cvGetSubRect( dft_templ, dst, 255 cvRect(0,(templ_cn>1?yofs:0),dftsize.width,dftsize.height) ); 256 257 cvMulSpectrums( _dft_img, dst, _dft_img, CV_DXT_MUL_CONJ ); 258 cvDFT( _dft_img, _dft_img, CV_DXT_INVERSE, csz.height ); 259 260 src = cvGetSubRect( _dft_img, &sstub, cvRect(0,0,csz.width,csz.height) ); 261 dst = cvGetSubRect( corr, &dstub, cvRect(x,y,csz.width,csz.height) ); 262 263 if( corr_cn > 1 ) 264 { 265 planes[i] = src; 266 if( corr_depth != max_depth ) 267 { 268 planes[i] = cvInitMatHeader( &temp, csz.height, csz.width, 269 corr_depth, _buf ); 270 cvConvert( src, planes[i] ); 271 } 272 cvMerge( planes[0], planes[1], planes[2], planes[3], dst ); 273 planes[i] = 0; 274 } 275 else 276 { 277 if( i == 0 ) 278 cvConvert( src, dst ); 279 else 280 { 281 if( max_depth > corr_depth ) 282 { 283 cvInitMatHeader( &temp, csz.height, csz.width, 284 corr_depth, _buf ); 285 cvConvert( src, &temp ); 286 src = &temp; 287 } 288 cvAcc( src, dst ); 289 } 290 } 291 } 292 } 293 } 294 295 __END__; 296 297 cvReleaseMat( &dft_templ ); 298 299 for( k = 0; k < num_threads; k++ ) 300 { 301 cvReleaseMat( &dft_img[k] ); 302 cvFree( &buf[k] ); 303 } 304 } 305 306 307 /***************************** IPP Match Template Functions ******************************/ 308 309 icvCrossCorrValid_Norm_8u32f_C1R_t icvCrossCorrValid_Norm_8u32f_C1R_p = 0; 310 icvCrossCorrValid_NormLevel_8u32f_C1R_t icvCrossCorrValid_NormLevel_8u32f_C1R_p = 0; 311 icvSqrDistanceValid_Norm_8u32f_C1R_t icvSqrDistanceValid_Norm_8u32f_C1R_p = 0; 312 icvCrossCorrValid_Norm_32f_C1R_t icvCrossCorrValid_Norm_32f_C1R_p = 0; 313 icvCrossCorrValid_NormLevel_32f_C1R_t icvCrossCorrValid_NormLevel_32f_C1R_p = 0; 314 icvSqrDistanceValid_Norm_32f_C1R_t icvSqrDistanceValid_Norm_32f_C1R_p = 0; 315 316 typedef CvStatus (CV_STDCALL * CvTemplMatchIPPFunc) 317 ( const void* img, int imgstep, CvSize imgsize, 318 const void* templ, int templstep, CvSize templsize, 319 void* result, int rstep ); 320 321 /*****************************************************************************************/ 322 323 CV_IMPL void 324 cvMatchTemplate( const CvArr* _img, const CvArr* _templ, CvArr* _result, int method ) 325 { 326 CvMat* sum = 0; 327 CvMat* sqsum = 0; 328 329 CV_FUNCNAME( "cvMatchTemplate" ); 330 331 __BEGIN__; 332 333 int coi1 = 0, coi2 = 0; 334 int depth, cn; 335 int i, j, k; 336 CvMat stub, *img = (CvMat*)_img; 337 CvMat tstub, *templ = (CvMat*)_templ; 338 CvMat rstub, *result = (CvMat*)_result; 339 CvScalar templ_mean = cvScalarAll(0); 340 double templ_norm = 0, templ_sum2 = 0; 341 342 int idx = 0, idx2 = 0; 343 double *p0, *p1, *p2, *p3; 344 double *q0, *q1, *q2, *q3; 345 double inv_area; 346 int sum_step, sqsum_step; 347 int num_type = method == CV_TM_CCORR || method == CV_TM_CCORR_NORMED ? 0 : 348 method == CV_TM_CCOEFF || method == CV_TM_CCOEFF_NORMED ? 1 : 2; 349 int is_normed = method == CV_TM_CCORR_NORMED || 350 method == CV_TM_SQDIFF_NORMED || 351 method == CV_TM_CCOEFF_NORMED; 352 353 CV_CALL( img = cvGetMat( img, &stub, &coi1 )); 354 CV_CALL( templ = cvGetMat( templ, &tstub, &coi2 )); 355 CV_CALL( result = cvGetMat( result, &rstub )); 356 357 if( CV_MAT_DEPTH( img->type ) != CV_8U && 358 CV_MAT_DEPTH( img->type ) != CV_32F ) 359 CV_ERROR( CV_StsUnsupportedFormat, 360 "The function supports only 8u and 32f data types" ); 361 362 if( !CV_ARE_TYPES_EQ( img, templ )) 363 CV_ERROR( CV_StsUnmatchedSizes, "image and template should have the same type" ); 364 365 if( CV_MAT_TYPE( result->type ) != CV_32FC1 ) 366 CV_ERROR( CV_StsUnsupportedFormat, "output image should have 32f type" ); 367 368 if( img->rows < templ->rows || img->cols < templ->cols ) 369 { 370 CvMat* t; 371 CV_SWAP( img, templ, t ); 372 } 373 374 if( result->rows != img->rows - templ->rows + 1 || 375 result->cols != img->cols - templ->cols + 1 ) 376 CV_ERROR( CV_StsUnmatchedSizes, "output image should be (W - w + 1)x(H - h + 1)" ); 377 378 if( method < CV_TM_SQDIFF || method > CV_TM_CCOEFF_NORMED ) 379 CV_ERROR( CV_StsBadArg, "unknown comparison method" ); 380 381 depth = CV_MAT_DEPTH(img->type); 382 cn = CV_MAT_CN(img->type); 383 384 if( is_normed && cn == 1 && templ->rows > 8 && templ->cols > 8 && 385 img->rows > templ->cols && img->cols > templ->cols ) 386 { 387 CvTemplMatchIPPFunc ipp_func = 388 depth == CV_8U ? 389 (method == CV_TM_SQDIFF_NORMED ? (CvTemplMatchIPPFunc)icvSqrDistanceValid_Norm_8u32f_C1R_p : 390 method == CV_TM_CCORR_NORMED ? (CvTemplMatchIPPFunc)icvCrossCorrValid_Norm_8u32f_C1R_p : 391 (CvTemplMatchIPPFunc)icvCrossCorrValid_NormLevel_8u32f_C1R_p) : 392 (method == CV_TM_SQDIFF_NORMED ? (CvTemplMatchIPPFunc)icvSqrDistanceValid_Norm_32f_C1R_p : 393 method == CV_TM_CCORR_NORMED ? (CvTemplMatchIPPFunc)icvCrossCorrValid_Norm_32f_C1R_p : 394 (CvTemplMatchIPPFunc)icvCrossCorrValid_NormLevel_32f_C1R_p); 395 396 if( ipp_func ) 397 { 398 CvSize img_size = cvGetMatSize(img), templ_size = cvGetMatSize(templ); 399 400 IPPI_CALL( ipp_func( img->data.ptr, img->step ? img->step : CV_STUB_STEP, 401 img_size, templ->data.ptr, 402 templ->step ? templ->step : CV_STUB_STEP, 403 templ_size, result->data.ptr, 404 result->step ? result->step : CV_STUB_STEP )); 405 for( i = 0; i < result->rows; i++ ) 406 { 407 float* rrow = (float*)(result->data.ptr + i*result->step); 408 for( j = 0; j < result->cols; j++ ) 409 { 410 if( fabs(rrow[j]) > 1. ) 411 rrow[j] = rrow[j] < 0 ? -1.f : 1.f; 412 } 413 } 414 EXIT; 415 } 416 } 417 418 CV_CALL( icvCrossCorr( img, templ, result )); 419 420 if( method == CV_TM_CCORR ) 421 EXIT; 422 423 inv_area = 1./((double)templ->rows * templ->cols); 424 425 CV_CALL( sum = cvCreateMat( img->rows + 1, img->cols + 1, 426 CV_MAKETYPE( CV_64F, cn ))); 427 if( method == CV_TM_CCOEFF ) 428 { 429 CV_CALL( cvIntegral( img, sum, 0, 0 )); 430 CV_CALL( templ_mean = cvAvg( templ )); 431 q0 = q1 = q2 = q3 = 0; 432 } 433 else 434 { 435 CvScalar _templ_sdv = cvScalarAll(0); 436 CV_CALL( sqsum = cvCreateMat( img->rows + 1, img->cols + 1, 437 CV_MAKETYPE( CV_64F, cn ))); 438 CV_CALL( cvIntegral( img, sum, sqsum, 0 )); 439 CV_CALL( cvAvgSdv( templ, &templ_mean, &_templ_sdv )); 440 441 templ_norm = CV_SQR(_templ_sdv.val[0]) + CV_SQR(_templ_sdv.val[1]) + 442 CV_SQR(_templ_sdv.val[2]) + CV_SQR(_templ_sdv.val[3]); 443 444 if( templ_norm < DBL_EPSILON && method == CV_TM_CCOEFF_NORMED ) 445 { 446 cvSet( result, cvScalarAll(1.) ); 447 EXIT; 448 } 449 450 templ_sum2 = templ_norm + 451 CV_SQR(templ_mean.val[0]) + CV_SQR(templ_mean.val[1]) + 452 CV_SQR(templ_mean.val[2]) + CV_SQR(templ_mean.val[3]); 453 454 if( num_type != 1 ) 455 { 456 templ_mean = cvScalarAll(0); 457 templ_norm = templ_sum2; 458 } 459 460 templ_sum2 /= inv_area; 461 templ_norm = sqrt(templ_norm); 462 templ_norm /= sqrt(inv_area); // care of accuracy here 463 464 q0 = (double*)sqsum->data.ptr; 465 q1 = q0 + templ->cols*cn; 466 q2 = (double*)(sqsum->data.ptr + templ->rows*sqsum->step); 467 q3 = q2 + templ->cols*cn; 468 } 469 470 p0 = (double*)sum->data.ptr; 471 p1 = p0 + templ->cols*cn; 472 p2 = (double*)(sum->data.ptr + templ->rows*sum->step); 473 p3 = p2 + templ->cols*cn; 474 475 sum_step = sum ? sum->step / sizeof(double) : 0; 476 sqsum_step = sqsum ? sqsum->step / sizeof(double) : 0; 477 478 for( i = 0; i < result->rows; i++ ) 479 { 480 float* rrow = (float*)(result->data.ptr + i*result->step); 481 idx = i * sum_step; 482 idx2 = i * sqsum_step; 483 484 for( j = 0; j < result->cols; j++, idx += cn, idx2 += cn ) 485 { 486 double num = rrow[j], t; 487 double wnd_mean2 = 0, wnd_sum2 = 0; 488 489 if( num_type == 1 ) 490 { 491 for( k = 0; k < cn; k++ ) 492 { 493 t = p0[idx+k] - p1[idx+k] - p2[idx+k] + p3[idx+k]; 494 wnd_mean2 += CV_SQR(t); 495 num -= t*templ_mean.val[k]; 496 } 497 498 wnd_mean2 *= inv_area; 499 } 500 501 if( is_normed || num_type == 2 ) 502 { 503 for( k = 0; k < cn; k++ ) 504 { 505 t = q0[idx2+k] - q1[idx2+k] - q2[idx2+k] + q3[idx2+k]; 506 wnd_sum2 += t; 507 } 508 509 if( num_type == 2 ) 510 num = wnd_sum2 - 2*num + templ_sum2; 511 } 512 513 if( is_normed ) 514 { 515 t = sqrt(MAX(wnd_sum2 - wnd_mean2,0))*templ_norm; 516 if( t > DBL_EPSILON ) 517 { 518 num /= t; 519 if( fabs(num) > 1. ) 520 num = num > 0 ? 1 : -1; 521 } 522 else 523 num = method != CV_TM_SQDIFF_NORMED || num < DBL_EPSILON ? 0 : 1; 524 } 525 526 rrow[j] = (float)num; 527 } 528 } 529 530 __END__; 531 532 cvReleaseMat( &sum ); 533 cvReleaseMat( &sqsum ); 534 } 535 536 /* End of file. */ 537
