modules/imgproc/src/templmatch.cpp

   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 "precomp.hpp"
  43
  44 namespace cv
  45 {
  46
  47 void crossCorr( const Mat& img, const Mat& _templ, Mat& corr,
  48                 Size corrsize, int ctype,
  49                 Point anchor, double delta, int borderType )
  50 {
  51     const double blockScale = 4.5;
  52     const int minBlockSize = 256;
  53     std::vector<uchar> buf;
  54
  55     Mat templ = _templ;
  56     int depth = img.depth(), cn = img.channels();
  57     int tdepth = templ.depth(), tcn = templ.channels();
  58     int cdepth = CV_MAT_DEPTH(ctype), ccn = CV_MAT_CN(ctype);
  59
  60     CV_Assert( img.dims <= 2 && templ.dims <= 2 && corr.dims <= 2 );
  61
  62     if( depth != tdepth && tdepth != std::max(CV_32F, depth) )
  63     {
  64         _templ.convertTo(templ, std::max(CV_32F, depth));
  65         tdepth = templ.depth();
  66     }
  67
  68     CV_Assert( depth == tdepth || tdepth == CV_32F);
  69     CV_Assert( corrsize.height <= img.rows + templ.rows - 1 &&
  70                corrsize.width <= img.cols + templ.cols - 1 );
  71
  72     CV_Assert( ccn == 1 || delta == 0 );
  73
  74     corr.create(corrsize, ctype);
  75
  76     int maxDepth = depth > CV_8S ? CV_64F : std::max(std::max(CV_32F, tdepth), cdepth);
  77     Size blocksize, dftsize;
  78
  79     blocksize.width = cvRound(templ.cols*blockScale);
  80     blocksize.width = std::max( blocksize.width, minBlockSize - templ.cols + 1 );
  81     blocksize.width = std::min( blocksize.width, corr.cols );
  82     blocksize.height = cvRound(templ.rows*blockScale);
  83     blocksize.height = std::max( blocksize.height, minBlockSize - templ.rows + 1 );
  84     blocksize.height = std::min( blocksize.height, corr.rows );
  85
  86     dftsize.width = std::max(getOptimalDFTSize(blocksize.width + templ.cols - 1), 2);
  87     dftsize.height = getOptimalDFTSize(blocksize.height + templ.rows - 1);
  88     if( dftsize.width <= 0 || dftsize.height <= 0 )
  89         CV_Error( CV_StsOutOfRange, "the input arrays are too big" );
  90
  91     // recompute block size
  92     blocksize.width = dftsize.width - templ.cols + 1;
  93     blocksize.width = MIN( blocksize.width, corr.cols );
  94     blocksize.height = dftsize.height - templ.rows + 1;
  95     blocksize.height = MIN( blocksize.height, corr.rows );
  96
  97     Mat dftTempl( dftsize.height*tcn, dftsize.width, maxDepth );
  98     Mat dftImg( dftsize, maxDepth );
  99
 100     int i, k, bufSize = 0;
 101     if( tcn > 1 && tdepth != maxDepth )
 102         bufSize = templ.cols*templ.rows*CV_ELEM_SIZE(tdepth);
 103
 104     if( cn > 1 && depth != maxDepth )
 105         bufSize = std::max( bufSize, (blocksize.width + templ.cols - 1)*
 106             (blocksize.height + templ.rows - 1)*CV_ELEM_SIZE(depth));
 107
 108     if( (ccn > 1 || cn > 1) && cdepth != maxDepth )
 109         bufSize = std::max( bufSize, blocksize.width*blocksize.height*CV_ELEM_SIZE(cdepth));
 110
 111     buf.resize(bufSize);
 112
 113     // compute DFT of each template plane
 114     for( k = 0; k < tcn; k++ )
 115     {
 116         int yofs = k*dftsize.height;
 117         Mat src = templ;
 118         Mat dst(dftTempl, Rect(0, yofs, dftsize.width, dftsize.height));
 119         Mat dst1(dftTempl, Rect(0, yofs, templ.cols, templ.rows));
 120
 121         if( tcn > 1 )
 122         {
 123             src = tdepth == maxDepth ? dst1 : Mat(templ.size(), tdepth, &buf[0]);
 124             int pairs[] = {k, 0};
 125             mixChannels(&templ, 1, &src, 1, pairs, 1);
 126         }
 127
 128         if( dst1.data != src.data )
 129             src.convertTo(dst1, dst1.depth());
 130
 131         if( dst.cols > templ.cols )
 132         {
 133             Mat part(dst, Range(0, templ.rows), Range(templ.cols, dst.cols));
 134             part = Scalar::all(0);
 135         }
 136         dft(dst, dst, 0, templ.rows);
 137     }
 138
 139     int tileCountX = (corr.cols + blocksize.width - 1)/blocksize.width;
 140     int tileCountY = (corr.rows + blocksize.height - 1)/blocksize.height;
 141     int tileCount = tileCountX * tileCountY;
 142
 143     Size wholeSize = img.size();
 144     Point roiofs(0,0);
 145     Mat img0 = img;
 146
 147     if( !(borderType & BORDER_ISOLATED) )
 148     {
 149         img.locateROI(wholeSize, roiofs);
 150         img0.adjustROI(roiofs.y, wholeSize.height-img.rows-roiofs.y,
 151                        roiofs.x, wholeSize.width-img.cols-roiofs.x);
 152     }
 153     borderType |= BORDER_ISOLATED;
 154
 155     // calculate correlation by blocks
 156     for( i = 0; i < tileCount; i++ )
 157     {
 158         int x = (i%tileCountX)*blocksize.width;
 159         int y = (i/tileCountX)*blocksize.height;
 160
 161         Size bsz(std::min(blocksize.width, corr.cols - x),
 162                  std::min(blocksize.height, corr.rows - y));
 163         Size dsz(bsz.width + templ.cols - 1, bsz.height + templ.rows - 1);
 164         int x0 = x - anchor.x + roiofs.x, y0 = y - anchor.y + roiofs.y;
 165         int x1 = std::max(0, x0), y1 = std::max(0, y0);
 166         int x2 = std::min(img0.cols, x0 + dsz.width);
 167         int y2 = std::min(img0.rows, y0 + dsz.height);
 168         Mat src0(img0, Range(y1, y2), Range(x1, x2));
 169         Mat dst(dftImg, Rect(0, 0, dsz.width, dsz.height));
 170         Mat dst1(dftImg, Rect(x1-x0, y1-y0, x2-x1, y2-y1));
 171         Mat cdst(corr, Rect(x, y, bsz.width, bsz.height));
 172
 173         for( k = 0; k < cn; k++ )
 174         {
 175             Mat src = src0;
 176             dftImg = Scalar::all(0);
 177
 178             if( cn > 1 )
 179             {
 180                 src = depth == maxDepth ? dst1 : Mat(y2-y1, x2-x1, depth, &buf[0]);
 181                 int pairs[] = {k, 0};
 182                 mixChannels(&src0, 1, &src, 1, pairs, 1);
 183             }
 184
 185             if( dst1.data != src.data )
 186                 src.convertTo(dst1, dst1.depth());
 187
 188             if( x2 - x1 < dsz.width || y2 - y1 < dsz.height )
 189                 copyMakeBorder(dst1, dst, y1-y0, dst.rows-dst1.rows-(y1-y0),
 190                                x1-x0, dst.cols-dst1.cols-(x1-x0), borderType);
 191
 192             dft( dftImg, dftImg, 0, dsz.height );
 193             Mat dftTempl1(dftTempl, Rect(0, tcn > 1 ? k*dftsize.height : 0,
 194                                          dftsize.width, dftsize.height));
 195             mulSpectrums(dftImg, dftTempl1, dftImg, 0, true);
 196             dft( dftImg, dftImg, DFT_INVERSE + DFT_SCALE, bsz.height );
 197
 198             src = dftImg(Rect(0, 0, bsz.width, bsz.height));
 199
 200             if( ccn > 1 )
 201             {
 202                 if( cdepth != maxDepth )
 203                 {
 204                     Mat plane(bsz, cdepth, &buf[0]);
 205                     src.convertTo(plane, cdepth, 1, delta);
 206                     src = plane;
 207                 }
 208                 int pairs[] = {0, k};
 209                 mixChannels(&src, 1, &cdst, 1, pairs, 1);
 210             }
 211             else
 212             {
 213                 if( k == 0 )
 214                     src.convertTo(cdst, cdepth, 1, delta);
 215                 else
 216                 {
 217                     if( maxDepth != cdepth )
 218                     {
 219                         Mat plane(bsz, cdepth, &buf[0]);
 220                         src.convertTo(plane, cdepth);
 221                         src = plane;
 222                     }
 223                     add(src, cdst, cdst);
 224                 }
 225             }
 226         }
 227     }
 228 }
 229
 230 }
 231
 232 /*****************************************************************************************/
 233
 234 void cv::matchTemplate( InputArray _img, InputArray _templ, OutputArray _result, int method )
 235 {
 236     CV_Assert( CV_TM_SQDIFF <= method && method <= CV_TM_CCOEFF_NORMED );
 237
 238     int numType = method == CV_TM_CCORR || method == CV_TM_CCORR_NORMED ? 0 :
 239                   method == CV_TM_CCOEFF || method == CV_TM_CCOEFF_NORMED ? 1 : 2;
 240     bool isNormed = method == CV_TM_CCORR_NORMED ||
 241                     method == CV_TM_SQDIFF_NORMED ||
 242                     method == CV_TM_CCOEFF_NORMED;
 243
 244     Mat img = _img.getMat(), templ = _templ.getMat();
 245     if( img.rows < templ.rows || img.cols < templ.cols )
 246         std::swap(img, templ);
 247
 248     CV_Assert( (img.depth() == CV_8U || img.depth() == CV_32F) &&
 249                img.type() == templ.type() );
 250
 251     Size corrSize(img.cols - templ.cols + 1, img.rows - templ.rows + 1);
 252     _result.create(corrSize, CV_32F);
 253     Mat result = _result.getMat();
 254
 255 #ifdef HAVE_TEGRA_OPTIMIZATION
 256     if (tegra::matchTemplate(img, templ, result, method))
 257         return;
 258 #endif
 259
 260     int cn = img.channels();
 261     crossCorr( img, templ, result, result.size(), result.type(), Point(0,0), 0, 0);
 262
 263     if( method == CV_TM_CCORR )
 264         return;
 265
 266     double invArea = 1./((double)templ.rows * templ.cols);
 267
 268     Mat sum, sqsum;
 269     Scalar templMean, templSdv;
 270     double *q0 = 0, *q1 = 0, *q2 = 0, *q3 = 0;
 271     double templNorm = 0, templSum2 = 0;
 272
 273     if( method == CV_TM_CCOEFF )
 274     {
 275         integral(img, sum, CV_64F);
 276         templMean = mean(templ);
 277     }
 278     else
 279     {
 280         integral(img, sum, sqsum, CV_64F);
 281         meanStdDev( templ, templMean, templSdv );
 282
 283         templNorm = CV_SQR(templSdv[0]) + CV_SQR(templSdv[1]) +
 284                     CV_SQR(templSdv[2]) + CV_SQR(templSdv[3]);
 285
 286         if( templNorm < DBL_EPSILON && method == CV_TM_CCOEFF_NORMED )
 287         {
 288             result = Scalar::all(1);
 289             return;
 290         }
 291
 292         templSum2 = templNorm +
 293                      CV_SQR(templMean[0]) + CV_SQR(templMean[1]) +
 294                      CV_SQR(templMean[2]) + CV_SQR(templMean[3]);
 295
 296         if( numType != 1 )
 297         {
 298             templMean = Scalar::all(0);
 299             templNorm = templSum2;
 300         }
 301
 302         templSum2 /= invArea;
 303         templNorm = sqrt(templNorm);
 304         templNorm /= sqrt(invArea); // care of accuracy here
 305
 306         q0 = (double*)sqsum.data;
 307         q1 = q0 + templ.cols*cn;
 308         q2 = (double*)(sqsum.data + templ.rows*sqsum.step);
 309         q3 = q2 + templ.cols*cn;
 310     }
 311
 312     double* p0 = (double*)sum.data;
 313     double* p1 = p0 + templ.cols*cn;
 314     double* p2 = (double*)(sum.data + templ.rows*sum.step);
 315     double* p3 = p2 + templ.cols*cn;
 316
 317     int sumstep = sum.data ? (int)(sum.step / sizeof(double)) : 0;
 318     int sqstep = sqsum.data ? (int)(sqsum.step / sizeof(double)) : 0;
 319
 320     int i, j, k;
 321
 322     for( i = 0; i < result.rows; i++ )
 323     {
 324         float* rrow = (float*)(result.data + i*result.step);
 325         int idx = i * sumstep;
 326         int idx2 = i * sqstep;
 327
 328         for( j = 0; j < result.cols; j++, idx += cn, idx2 += cn )
 329         {
 330             double num = rrow[j], t;
 331             double wndMean2 = 0, wndSum2 = 0;
 332
 333             if( numType == 1 )
 334             {
 335                 for( k = 0; k < cn; k++ )
 336                 {
 337                     t = p0[idx+k] - p1[idx+k] - p2[idx+k] + p3[idx+k];
 338                     wndMean2 += CV_SQR(t);
 339                     num -= t*templMean[k];
 340                 }
 341
 342                 wndMean2 *= invArea;
 343             }
 344
 345             if( isNormed || numType == 2 )
 346             {
 347                 for( k = 0; k < cn; k++ )
 348                 {
 349                     t = q0[idx2+k] - q1[idx2+k] - q2[idx2+k] + q3[idx2+k];
 350                     wndSum2 += t;
 351                 }
 352
 353                 if( numType == 2 )
 354                 {
 355                     num = wndSum2 - 2*num + templSum2;
 356                     num = MAX(num, 0.);
 357                 }
 358             }
 359
 360             if( isNormed )
 361             {
 362                 t = sqrt(MAX(wndSum2 - wndMean2,0))*templNorm;
 363                 if( fabs(num) < t )
 364                     num /= t;
 365                 else if( fabs(num) < t*1.125 )
 366                     num = num > 0 ? 1 : -1;
 367                 else
 368                     num = method != CV_TM_SQDIFF_NORMED ? 0 : 1;
 369             }
 370
 371             rrow[j] = (float)num;
 372         }
 373     }
 374 }
 375
 376
 377 CV_IMPL void
 378 cvMatchTemplate( const CvArr* _img, const CvArr* _templ, CvArr* _result, int method )
 379 {
 380     cv::Mat img = cv::cvarrToMat(_img), templ = cv::cvarrToMat(_templ),
 381         result = cv::cvarrToMat(_result);
 382     CV_Assert( result.size() == cv::Size(std::abs(img.cols - templ.cols) + 1,
 383                                          std::abs(img.rows - templ.rows) + 1) &&
 384               result.type() == CV_32F );
 385     matchTemplate(img, templ, result, method);
 386 }
 387
 388 /* End of file. */