vcs.maemo.org Git - opencv/blob - src/cv/cvsumpixels.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 "_cv.h"
  44
  45 namespace cv
  46 {
  47
  48 template<typename QT> inline QT sqr(uchar a) { return a*a; }
  49 template<typename QT> inline QT sqr(float a) { return a*a; }
  50 template<typename QT> inline QT sqr(double a) { return a*a; }
  51 template<> inline double sqr(uchar a) { return CV_8TO32F_SQR(a); }
  52
  53
  54 template<typename T, typename ST, typename QT>
  55 void integral_( const Mat& _src, Mat& _sum, Mat& _sqsum, Mat& _tilted )
  56 {
  57     int cn = _src.channels();
  58     Size size = _src.size();
  59     int x, y, k;
  60
  61     const T* src = (const T*)_src.data;
  62     ST* sum = (ST*)_sum.data;
  63     ST* tilted = (ST*)_tilted.data;
  64     QT* sqsum = (QT*)_sqsum.data;
  65
  66     int srcstep = (int)(_src.step/sizeof(T));
  67     int sumstep = (int)(_sum.step/sizeof(ST));
  68     int tiltedstep = (int)(_tilted.step/sizeof(ST));
  69     int sqsumstep = (int)(_sqsum.step/sizeof(QT));
  70
  71     size.width *= cn;
  72
  73     memset( sum, 0, (size.width+cn)*sizeof(sum[0]));
  74     sum += sumstep + cn;
  75
  76     if( sqsum )
  77     {
  78         memset( sqsum, 0, (size.width+cn)*sizeof(sqsum[0]));
  79         sqsum += sqsumstep + cn;
  80     }
  81
  82     if( tilted )
  83     {
  84         memset( tilted, 0, (size.width+cn)*sizeof(tilted[0]));
  85         tilted += tiltedstep + cn;
  86     }
  87
  88     if( sqsum == 0 && tilted == 0 )
  89     {
  90         for( y = 0; y < size.height; y++, src += srcstep - cn, sum += sumstep - cn )
  91         {
  92             for( k = 0; k < cn; k++, src++, sum++ )
  93             {
  94                 ST s = sum[-cn] = 0;
  95                 for( x = 0; x < size.width; x += cn )
  96                 {
  97                     s += src[x];
  98                     sum[x] = sum[x - sumstep] + s;
  99                 }
 100             }
 101         }
 102     }
 103     else if( tilted == 0 )
 104     {
 105         for( y = 0; y < size.height; y++, src += srcstep - cn,
 106                         sum += sumstep - cn, sqsum += sqsumstep - cn )
 107         {
 108             for( k = 0; k < cn; k++, src++, sum++, sqsum++ )
 109             {
 110                 ST s = sum[-cn] = 0;
 111                 QT sq = sqsum[-cn] = 0;
 112                 for( x = 0; x < size.width; x += cn )
 113                 {
 114                     T it = src[x];
 115                     s += it;
 116                     sq += sqr<QT>(it);
 117                     ST t = sum[x - sumstep] + s;
 118                     QT tq = sqsum[x - sqsumstep] + sq;
 119                     sum[x] = t;
 120                     sqsum[x] = tq;
 121                 }
 122             }
 123         }
 124     }
 125     else
 126     {
 127         AutoBuffer<ST> _buf(size.width+cn);
 128         ST* buf = _buf;
 129         ST s;
 130         QT sq;
 131         for( k = 0; k < cn; k++, src++, sum++, tilted++, sqsum++, buf++ )
 132         {
 133             sum[-cn] = tilted[-cn] = 0;
 134             sqsum[-cn] = 0;
 135
 136             for( x = 0, s = 0, sq = 0; x < size.width; x += cn )
 137             {
 138                 T it = src[x];
 139                 buf[x] = tilted[x] = it;
 140                 s += it;
 141                 sq += sqr<QT>(it);
 142                 sum[x] = s;
 143                 sqsum[x] = sq;
 144             }
 145
 146             if( size.width == cn )
 147                 buf[cn] = 0;
 148         }
 149
 150         for( y = 1; y < size.height; y++ )
 151         {
 152             src += srcstep - cn;
 153             sum += sumstep - cn;
 154             sqsum += sqsumstep - cn;
 155             tilted += tiltedstep - cn;
 156             buf += -cn;
 157
 158             for( k = 0; k < cn; k++, src++, sum++, sqsum++, tilted++, buf++ )
 159             {
 160                 T it = src[0];
 161                 ST t0 = s = it;
 162                 QT tq0 = sq = sqr<QT>(it);
 163
 164                 sum[-cn] = 0;
 165                 sqsum[-cn] = 0;
 166                 tilted[-cn] = tilted[-tiltedstep];
 167
 168                 sum[0] = sum[-sumstep] + t0;
 169                 sqsum[0] = sqsum[-sqsumstep] + tq0;
 170                 tilted[0] = tilted[-tiltedstep] + t0 + buf[cn];
 171
 172                 for( x = cn; x < size.width - cn; x += cn )
 173                 {
 174                     ST t1 = buf[x];
 175                     buf[x - cn] = t1 + t0;
 176                     t0 = it = src[x];
 177                     tq0 = sqr<QT>(it);
 178                     s += t0;
 179                     sq += tq0;
 180                     sum[x] = sum[x - sumstep] + s;
 181                     sqsum[x] = sqsum[x - sqsumstep] + sq;
 182                     t1 += buf[x + cn] + t0 + tilted[x - tiltedstep - cn];
 183                     tilted[x] = t1;
 184                 }
 185
 186                 if( size.width > cn )
 187                 {
 188                     ST t1 = buf[x];
 189                     buf[x - cn] = t1 + t0;
 190                     t0 = it = src[x];
 191                     tq0 = sqr<QT>(it);
 192                     s += t0;
 193                     sq += tq0;
 194                     sum[x] = sum[x - sumstep] + s;
 195                     sqsum[x] = sqsum[x - sqsumstep] + sq;
 196                     tilted[x] = t0 + t1 + tilted[x - tiltedstep - cn];
 197                     buf[x] = t0;
 198                 }
 199             }
 200         }
 201     }
 202 }
 203
 204 typedef void (*IntegralFunc)(const Mat& _src, Mat& _sum, Mat& _sqsum, Mat& _tilted );
 205
 206 static void
 207 integral( const Mat& src, Mat& sum, Mat* _sqsum, Mat* _tilted, int sdepth )
 208 {
 209     int depth = src.depth(), cn = src.channels();
 210     Size isize(src.cols + 1, src.rows+1);
 211     Mat sqsum, tilted;
 212
 213     if( sdepth <= 0 )
 214         sdepth = depth == CV_8U ? CV_32S : CV_64F;
 215     sdepth = CV_MAT_DEPTH(sdepth);
 216     sum.create( isize, CV_MAKETYPE(sdepth, cn) );
 217
 218     if( _tilted )
 219         _tilted->create( isize, CV_MAKETYPE(sdepth, cn) );
 220     else
 221         _tilted = &tilted;
 222
 223     if( !_sqsum )
 224         _sqsum = &sqsum;
 225
 226     if( _sqsum != &sqsum || _tilted->data )
 227         _sqsum->create( isize, CV_MAKETYPE(CV_64F, cn) );
 228
 229     IntegralFunc func = 0;
 230
 231     if( depth == CV_8U && sdepth == CV_32S )
 232         func = integral_<uchar, int, double>;
 233     else if( depth == CV_8U && sdepth == CV_32F )
 234         func = integral_<uchar, float, double>;
 235     else if( depth == CV_8U && sdepth == CV_64F )
 236         func = integral_<uchar, double, double>;
 237     else if( depth == CV_32F && sdepth == CV_64F )
 238         func = integral_<float, double, double>;
 239     else if( depth == CV_64F && sdepth == CV_64F )
 240         func = integral_<double, double, double>;
 241     else
 242         CV_Error( CV_StsUnsupportedFormat, "" );
 243
 244     func( src, sum, *_sqsum, *_tilted );
 245 }
 246
 247 void integral( const Mat& src, Mat& sum, int sdepth )
 248 {
 249     integral( src, sum, 0, 0, sdepth );
 250 }
 251
 252 void integral( const Mat& src, Mat& sum, Mat& sqsum, int sdepth )
 253 {
 254     integral( src, sum, &sqsum, 0, sdepth );
 255 }
 256
 257 void integral( const Mat& src, Mat& sum, Mat& sqsum, Mat& tilted, int sdepth )
 258 {
 259     integral( src, sum, &sqsum, &tilted, sdepth );
 260 }
 261
 262 }
 263
 264
 265 CV_IMPL void
 266 cvIntegral( const CvArr* image, CvArr* sumImage,
 267             CvArr* sumSqImage, CvArr* tiltedSumImage )
 268 {
 269     cv::Mat src = cv::cvarrToMat(image), sum = cv::cvarrToMat(sumImage), sum0 = sum;
 270     cv::Mat sqsum0, sqsum, tilted0, tilted;
 271     cv::Mat *psqsum = 0, *ptilted = 0;
 272
 273     if( sumSqImage )
 274     {
 275         sqsum0 = sqsum = cv::cvarrToMat(sumSqImage);
 276         psqsum = &sqsum;
 277     }
 278
 279     if( tiltedSumImage )
 280     {
 281         tilted0 = tilted = cv::cvarrToMat(tiltedSumImage);
 282         ptilted = &tilted;
 283     }
 284     cv::integral( src, sum, psqsum, ptilted, sum.depth() );
 285
 286     CV_Assert( sum.data == sum0.data && sqsum.data == sqsum0.data && tilted.data == tilted0.data );
 287 }
 288
 289 /* End of file. */