Update to 2.0.0 tree from current Fremantle build
[opencv] / src / cv / cvsumpixels.cpp
diff --git a/src/cv/cvsumpixels.cpp b/src/cv/cvsumpixels.cpp
new file mode 100644 (file)
index 0000000..bb69e7e
--- /dev/null
@@ -0,0 +1,289 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "_cv.h"
+
+namespace cv
+{
+
+template<typename QT> inline QT sqr(uchar a) { return a*a; }
+template<typename QT> inline QT sqr(float a) { return a*a; }
+template<typename QT> inline QT sqr(double a) { return a*a; }
+template<> inline double sqr(uchar a) { return CV_8TO32F_SQR(a); }
+
+
+template<typename T, typename ST, typename QT>
+void integral_( const Mat& _src, Mat& _sum, Mat& _sqsum, Mat& _tilted )
+{
+    int cn = _src.channels();
+    Size size = _src.size();
+    int x, y, k;
+
+    const T* src = (const T*)_src.data;
+    ST* sum = (ST*)_sum.data;
+    ST* tilted = (ST*)_tilted.data;
+    QT* sqsum = (QT*)_sqsum.data;
+
+    int srcstep = (int)(_src.step/sizeof(T));
+    int sumstep = (int)(_sum.step/sizeof(ST));
+    int tiltedstep = (int)(_tilted.step/sizeof(ST));
+    int sqsumstep = (int)(_sqsum.step/sizeof(QT));
+
+    size.width *= cn;
+
+    memset( sum, 0, (size.width+cn)*sizeof(sum[0]));
+    sum += sumstep + cn;
+
+    if( sqsum )
+    {
+        memset( sqsum, 0, (size.width+cn)*sizeof(sqsum[0]));
+        sqsum += sqsumstep + cn;
+    }
+
+    if( tilted )
+    {
+        memset( tilted, 0, (size.width+cn)*sizeof(tilted[0]));
+        tilted += tiltedstep + cn;
+    }
+
+    if( sqsum == 0 && tilted == 0 )
+    {
+        for( y = 0; y < size.height; y++, src += srcstep - cn, sum += sumstep - cn )
+        {
+            for( k = 0; k < cn; k++, src++, sum++ )
+            {
+                ST s = sum[-cn] = 0;
+                for( x = 0; x < size.width; x += cn )
+                {
+                    s += src[x];
+                    sum[x] = sum[x - sumstep] + s;
+                }
+            }
+        }
+    }
+    else if( tilted == 0 )
+    {
+        for( y = 0; y < size.height; y++, src += srcstep - cn,
+                        sum += sumstep - cn, sqsum += sqsumstep - cn )
+        {
+            for( k = 0; k < cn; k++, src++, sum++, sqsum++ )
+            {
+                ST s = sum[-cn] = 0;
+                QT sq = sqsum[-cn] = 0;
+                for( x = 0; x < size.width; x += cn )
+                {
+                    T it = src[x];
+                    s += it;
+                    sq += sqr<QT>(it);
+                    ST t = sum[x - sumstep] + s;
+                    QT tq = sqsum[x - sqsumstep] + sq;
+                    sum[x] = t;
+                    sqsum[x] = tq;
+                }
+            }
+        }
+    }
+    else
+    {
+        AutoBuffer<ST> _buf(size.width+cn);
+        ST* buf = _buf;
+        ST s;
+        QT sq;
+        for( k = 0; k < cn; k++, src++, sum++, tilted++, sqsum++, buf++ )
+        {
+            sum[-cn] = tilted[-cn] = 0;
+            sqsum[-cn] = 0;
+
+            for( x = 0, s = 0, sq = 0; x < size.width; x += cn )
+            {
+                T it = src[x];
+                buf[x] = tilted[x] = it;
+                s += it;
+                sq += sqr<QT>(it);
+                sum[x] = s;
+                sqsum[x] = sq;
+            }
+
+            if( size.width == cn )
+                buf[cn] = 0;
+        }
+
+        for( y = 1; y < size.height; y++ )
+        {
+            src += srcstep - cn;
+            sum += sumstep - cn;
+            sqsum += sqsumstep - cn;
+            tilted += tiltedstep - cn;
+            buf += -cn;
+
+            for( k = 0; k < cn; k++, src++, sum++, sqsum++, tilted++, buf++ )
+            {
+                T it = src[0];
+                ST t0 = s = it;
+                QT tq0 = sq = sqr<QT>(it);
+
+                sum[-cn] = 0;
+                sqsum[-cn] = 0;
+                tilted[-cn] = tilted[-tiltedstep];
+
+                sum[0] = sum[-sumstep] + t0;
+                sqsum[0] = sqsum[-sqsumstep] + tq0;
+                tilted[0] = tilted[-tiltedstep] + t0 + buf[cn];
+
+                for( x = cn; x < size.width - cn; x += cn )
+                {
+                    ST t1 = buf[x];
+                    buf[x - cn] = t1 + t0;
+                    t0 = it = src[x];
+                    tq0 = sqr<QT>(it);
+                    s += t0;
+                    sq += tq0;
+                    sum[x] = sum[x - sumstep] + s;
+                    sqsum[x] = sqsum[x - sqsumstep] + sq;
+                    t1 += buf[x + cn] + t0 + tilted[x - tiltedstep - cn];
+                    tilted[x] = t1;
+                }
+
+                if( size.width > cn )
+                {
+                    ST t1 = buf[x];
+                    buf[x - cn] = t1 + t0;
+                    t0 = it = src[x];
+                    tq0 = sqr<QT>(it);
+                    s += t0;
+                    sq += tq0;
+                    sum[x] = sum[x - sumstep] + s;
+                    sqsum[x] = sqsum[x - sqsumstep] + sq;
+                    tilted[x] = t0 + t1 + tilted[x - tiltedstep - cn];
+                    buf[x] = t0;
+                }
+            }
+        }
+    }
+}
+
+typedef void (*IntegralFunc)(const Mat& _src, Mat& _sum, Mat& _sqsum, Mat& _tilted );
+
+static void
+integral( const Mat& src, Mat& sum, Mat* _sqsum, Mat* _tilted, int sdepth )
+{
+    int depth = src.depth(), cn = src.channels();
+    Size isize(src.cols + 1, src.rows+1);
+    Mat sqsum, tilted;
+
+    if( sdepth <= 0 )
+        sdepth = depth == CV_8U ? CV_32S : CV_64F;
+    sdepth = CV_MAT_DEPTH(sdepth);
+    sum.create( isize, CV_MAKETYPE(sdepth, cn) );
+    
+    if( _tilted )
+        _tilted->create( isize, CV_MAKETYPE(sdepth, cn) );
+    else
+        _tilted = &tilted;
+    
+    if( !_sqsum )
+        _sqsum = &sqsum;
+    
+    if( _sqsum != &sqsum || _tilted->data )
+        _sqsum->create( isize, CV_MAKETYPE(CV_64F, cn) );
+
+    IntegralFunc func = 0;
+
+    if( depth == CV_8U && sdepth == CV_32S )
+        func = integral_<uchar, int, double>;
+    else if( depth == CV_8U && sdepth == CV_32F )
+        func = integral_<uchar, float, double>;
+    else if( depth == CV_8U && sdepth == CV_64F )
+        func = integral_<uchar, double, double>;
+    else if( depth == CV_32F && sdepth == CV_64F )
+        func = integral_<float, double, double>;
+    else if( depth == CV_64F && sdepth == CV_64F )
+        func = integral_<double, double, double>;
+    else
+        CV_Error( CV_StsUnsupportedFormat, "" );
+
+    func( src, sum, *_sqsum, *_tilted );
+}
+
+void integral( const Mat& src, Mat& sum, int sdepth )
+{
+    integral( src, sum, 0, 0, sdepth );
+}
+
+void integral( const Mat& src, Mat& sum, Mat& sqsum, int sdepth )
+{
+    integral( src, sum, &sqsum, 0, sdepth );
+}
+
+void integral( const Mat& src, Mat& sum, Mat& sqsum, Mat& tilted, int sdepth )
+{
+    integral( src, sum, &sqsum, &tilted, sdepth );
+}
+
+}
+
+
+CV_IMPL void
+cvIntegral( const CvArr* image, CvArr* sumImage,
+            CvArr* sumSqImage, CvArr* tiltedSumImage )
+{
+    cv::Mat src = cv::cvarrToMat(image), sum = cv::cvarrToMat(sumImage), sum0 = sum;
+    cv::Mat sqsum0, sqsum, tilted0, tilted;
+    cv::Mat *psqsum = 0, *ptilted = 0;
+
+    if( sumSqImage )
+    {
+        sqsum0 = sqsum = cv::cvarrToMat(sumSqImage);
+        psqsum = &sqsum;
+    }
+
+    if( tiltedSumImage )
+    {
+        tilted0 = tilted = cv::cvarrToMat(tiltedSumImage);
+        ptilted = &tilted;
+    }
+    cv::integral( src, sum, psqsum, ptilted, sum.depth() );
+
+    CV_Assert( sum.data == sum0.data && sqsum.data == sqsum0.data && tilted.data == tilted0.data );
+}
+
+/* End of file. */