modules/imgproc/src/sumpixels.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 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "precomp.hpp"
  44 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
  45 static IppStatus sts = ippInit();
  46 #endif
  47
  48 namespace cv
  49 {
  50
  51 template<typename T, typename ST, typename QT>
  52 void integral_( const T* src, size_t _srcstep, ST* sum, size_t _sumstep,
  53                 QT* sqsum, size_t _sqsumstep, ST* tilted, size_t _tiltedstep,
  54                 Size size, int cn )
  55 {
  56     int x, y, k;
  57
  58     int srcstep = (int)(_srcstep/sizeof(T));
  59     int sumstep = (int)(_sumstep/sizeof(ST));
  60     int tiltedstep = (int)(_tiltedstep/sizeof(ST));
  61     int sqsumstep = (int)(_sqsumstep/sizeof(QT));
  62
  63     size.width *= cn;
  64
  65     memset( sum, 0, (size.width+cn)*sizeof(sum[0]));
  66     sum += sumstep + cn;
  67
  68     if( sqsum )
  69     {
  70         memset( sqsum, 0, (size.width+cn)*sizeof(sqsum[0]));
  71         sqsum += sqsumstep + cn;
  72     }
  73
  74     if( tilted )
  75     {
  76         memset( tilted, 0, (size.width+cn)*sizeof(tilted[0]));
  77         tilted += tiltedstep + cn;
  78     }
  79
  80     if( sqsum == 0 && tilted == 0 )
  81     {
  82         for( y = 0; y < size.height; y++, src += srcstep - cn, sum += sumstep - cn )
  83         {
  84             for( k = 0; k < cn; k++, src++, sum++ )
  85             {
  86                 ST s = sum[-cn] = 0;
  87                 for( x = 0; x < size.width; x += cn )
  88                 {
  89                     s += src[x];
  90                     sum[x] = sum[x - sumstep] + s;
  91                 }
  92             }
  93         }
  94     }
  95     else if( tilted == 0 )
  96     {
  97         for( y = 0; y < size.height; y++, src += srcstep - cn,
  98                         sum += sumstep - cn, sqsum += sqsumstep - cn )
  99         {
 100             for( k = 0; k < cn; k++, src++, sum++, sqsum++ )
 101             {
 102                 ST s = sum[-cn] = 0;
 103                 QT sq = sqsum[-cn] = 0;
 104                 for( x = 0; x < size.width; x += cn )
 105                 {
 106                     T it = src[x];
 107                     s += it;
 108                     sq += (QT)it*it;
 109                     ST t = sum[x - sumstep] + s;
 110                     QT tq = sqsum[x - sqsumstep] + sq;
 111                     sum[x] = t;
 112                     sqsum[x] = tq;
 113                 }
 114             }
 115         }
 116     }
 117     else
 118     {
 119         AutoBuffer<ST> _buf(size.width+cn);
 120         ST* buf = _buf;
 121         ST s;
 122         QT sq;
 123         for( k = 0; k < cn; k++, src++, sum++, tilted++, buf++ )
 124         {
 125             sum[-cn] = tilted[-cn] = 0;
 126
 127             for( x = 0, s = 0, sq = 0; x < size.width; x += cn )
 128             {
 129                 T it = src[x];
 130                 buf[x] = tilted[x] = it;
 131                 s += it;
 132                 sq += (QT)it*it;
 133                 sum[x] = s;
 134                 if( sqsum )
 135                     sqsum[x] = sq;
 136             }
 137
 138             if( size.width == cn )
 139                 buf[cn] = 0;
 140
 141             if( sqsum )
 142             {
 143                 sqsum[-cn] = 0;
 144                 sqsum++;
 145             }
 146         }
 147
 148         for( y = 1; y < size.height; y++ )
 149         {
 150             src += srcstep - cn;
 151             sum += sumstep - cn;
 152             tilted += tiltedstep - cn;
 153             buf += -cn;
 154
 155             if( sqsum )
 156                 sqsum += sqsumstep - cn;
 157
 158             for( k = 0; k < cn; k++, src++, sum++, tilted++, buf++ )
 159             {
 160                 T it = src[0];
 161                 ST t0 = s = it;
 162                 QT tq0 = sq = (QT)it*it;
 163
 164                 sum[-cn] = 0;
 165                 if( sqsum )
 166                     sqsum[-cn] = 0;
 167                 tilted[-cn] = tilted[-tiltedstep];
 168
 169                 sum[0] = sum[-sumstep] + t0;
 170                 if( sqsum )
 171                     sqsum[0] = sqsum[-sqsumstep] + tq0;
 172                 tilted[0] = tilted[-tiltedstep] + t0 + buf[cn];
 173
 174                 for( x = cn; x < size.width - cn; x += cn )
 175                 {
 176                     ST t1 = buf[x];
 177                     buf[x - cn] = t1 + t0;
 178                     t0 = it = src[x];
 179                     tq0 = (QT)it*it;
 180                     s += t0;
 181                     sq += tq0;
 182                     sum[x] = sum[x - sumstep] + s;
 183                     if( sqsum )
 184                         sqsum[x] = sqsum[x - sqsumstep] + sq;
 185                     t1 += buf[x + cn] + t0 + tilted[x - tiltedstep - cn];
 186                     tilted[x] = t1;
 187                 }
 188
 189                 if( size.width > cn )
 190                 {
 191                     ST t1 = buf[x];
 192                     buf[x - cn] = t1 + t0;
 193                     t0 = it = src[x];
 194                     tq0 = (QT)it*it;
 195                     s += t0;
 196                     sq += tq0;
 197                     sum[x] = sum[x - sumstep] + s;
 198                     if( sqsum )
 199                         sqsum[x] = sqsum[x - sqsumstep] + sq;
 200                     tilted[x] = t0 + t1 + tilted[x - tiltedstep - cn];
 201                     buf[x] = t0;
 202                 }
 203
 204                 if( sqsum )
 205                     sqsum++;
 206             }
 207         }
 208     }
 209 }
 210
 211
 212 #define DEF_INTEGRAL_FUNC(suffix, T, ST, QT) \
 213 static void integral_##suffix( T* src, size_t srcstep, ST* sum, size_t sumstep, QT* sqsum, size_t sqsumstep, \
 214                               ST* tilted, size_t tiltedstep, Size size, int cn ) \
 215 { integral_(src, srcstep, sum, sumstep, sqsum, sqsumstep, tilted, tiltedstep, size, cn); }
 216
 217 DEF_INTEGRAL_FUNC(8u32s, uchar, int, double)
 218 DEF_INTEGRAL_FUNC(8u32f, uchar, float, double)
 219 DEF_INTEGRAL_FUNC(8u64f, uchar, double, double)
 220 DEF_INTEGRAL_FUNC(32f, float, float, double)
 221 DEF_INTEGRAL_FUNC(32f64f, float, double, double)
 222 DEF_INTEGRAL_FUNC(64f, double, double, double)
 223
 224 typedef void (*IntegralFunc)(const uchar* src, size_t srcstep, uchar* sum, size_t sumstep,
 225                              uchar* sqsum, size_t sqsumstep, uchar* tilted, size_t tstep,
 226                              Size size, int cn );
 227
 228 }
 229
 230
 231 void cv::integral( InputArray _src, OutputArray _sum, OutputArray _sqsum, OutputArray _tilted, int sdepth )
 232 {
 233     Mat src = _src.getMat(), sum, sqsum, tilted;
 234     int depth = src.depth(), cn = src.channels();
 235     Size isize(src.cols + 1, src.rows+1);
 236
 237     if( sdepth <= 0 )
 238         sdepth = depth == CV_8U ? CV_32S : CV_64F;
 239     sdepth = CV_MAT_DEPTH(sdepth);
 240
 241 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
 242     if( ( depth == CV_8U ) && ( !_tilted.needed() ) )
 243     {
 244         if( sdepth == CV_32F )
 245         {
 246             if( cn == 1 )
 247             {
 248                 IppiSize srcRoiSize = ippiSize( src.cols, src.rows );
 249                 _sum.create( isize, CV_MAKETYPE( sdepth, cn ) );
 250                 sum = _sum.getMat();
 251                 if( _sqsum.needed() )
 252                 {
 253                     _sqsum.create( isize, CV_MAKETYPE( CV_64F, cn ) );
 254                     sqsum = _sqsum.getMat();
 255                     ippiSqrIntegral_8u32f64f_C1R( (const Ipp8u*)src.data, (int)src.step, (Ipp32f*)sum.data, (int)sum.step, (Ipp64f*)sqsum.data, (int)sqsum.step, srcRoiSize, 0, 0 );
 256                 }
 257                 else
 258                 {
 259                     ippiIntegral_8u32f_C1R( (const Ipp8u*)src.data, (int)src.step, (Ipp32f*)sum.data, (int)sum.step, srcRoiSize, 0 );
 260                 }
 261                 return;
 262             }
 263         }
 264         if( sdepth == CV_32S )
 265         {
 266             if( cn == 1 )
 267             {
 268                 IppiSize srcRoiSize = ippiSize( src.cols, src.rows );
 269                 _sum.create( isize, CV_MAKETYPE( sdepth, cn ) );
 270                 sum = _sum.getMat();
 271                 if( _sqsum.needed() )
 272                 {
 273                     _sqsum.create( isize, CV_MAKETYPE( CV_64F, cn ) );
 274                     sqsum = _sqsum.getMat();
 275                     ippiSqrIntegral_8u32s64f_C1R( (const Ipp8u*)src.data, (int)src.step, (Ipp32s*)sum.data, (int)sum.step, (Ipp64f*)sqsum.data, (int)sqsum.step, srcRoiSize, 0, 0 );
 276                 }
 277                 else
 278                 {
 279                     ippiIntegral_8u32s_C1R( (const Ipp8u*)src.data, (int)src.step, (Ipp32s*)sum.data, (int)sum.step, srcRoiSize, 0 );
 280                 }
 281                 return;
 282             }
 283         }
 284     }
 285 #endif
 286
 287     _sum.create( isize, CV_MAKETYPE(sdepth, cn) );
 288     sum = _sum.getMat();
 289
 290     if( _tilted.needed() )
 291     {
 292         _tilted.create( isize, CV_MAKETYPE(sdepth, cn) );
 293         tilted = _tilted.getMat();
 294     }
 295
 296     if( _sqsum.needed() )
 297     {
 298         _sqsum.create( isize, CV_MAKETYPE(CV_64F, cn) );
 299         sqsum = _sqsum.getMat();
 300     }
 301
 302     IntegralFunc func = 0;
 303
 304     if( depth == CV_8U && sdepth == CV_32S )
 305         func = (IntegralFunc)GET_OPTIMIZED(integral_8u32s);
 306     else if( depth == CV_8U && sdepth == CV_32F )
 307         func = (IntegralFunc)integral_8u32f;
 308     else if( depth == CV_8U && sdepth == CV_64F )
 309         func = (IntegralFunc)integral_8u64f;
 310     else if( depth == CV_32F && sdepth == CV_32F )
 311         func = (IntegralFunc)integral_32f;
 312     else if( depth == CV_32F && sdepth == CV_64F )
 313         func = (IntegralFunc)integral_32f64f;
 314     else if( depth == CV_64F && sdepth == CV_64F )
 315         func = (IntegralFunc)integral_64f;
 316     else
 317         CV_Error( CV_StsUnsupportedFormat, "" );
 318
 319     func( src.data, src.step, sum.data, sum.step, sqsum.data, sqsum.step,
 320           tilted.data, tilted.step, src.size(), cn );
 321 }
 322
 323 void cv::integral( InputArray src, OutputArray sum, int sdepth )
 324 {
 325     integral( src, sum, noArray(), noArray(), sdepth );
 326 }
 327
 328 void cv::integral( InputArray src, OutputArray sum, OutputArray sqsum, int sdepth )
 329 {
 330     integral( src, sum, sqsum, noArray(), sdepth );
 331 }
 332
 333
 334 CV_IMPL void
 335 cvIntegral( const CvArr* image, CvArr* sumImage,
 336             CvArr* sumSqImage, CvArr* tiltedSumImage )
 337 {
 338     cv::Mat src = cv::cvarrToMat(image), sum = cv::cvarrToMat(sumImage), sum0 = sum;
 339     cv::Mat sqsum0, sqsum, tilted0, tilted;
 340     cv::Mat *psqsum = 0, *ptilted = 0;
 341
 342     if( sumSqImage )
 343     {
 344         sqsum0 = sqsum = cv::cvarrToMat(sumSqImage);
 345         psqsum = &sqsum;
 346     }
 347
 348     if( tiltedSumImage )
 349     {
 350         tilted0 = tilted = cv::cvarrToMat(tiltedSumImage);
 351         ptilted = &tilted;
 352     }
 353     cv::integral( src, sum, psqsum ? cv::_OutputArray(*psqsum) : cv::_OutputArray(),
 354                   ptilted ? cv::_OutputArray(*ptilted) : cv::_OutputArray(), sum.depth() );
 355
 356     CV_Assert( sum.data == sum0.data && sqsum.data == sqsum0.data && tilted.data == tilted0.data );
 357 }
 358
 359 /* End of file. */