Update to 2.0.0 tree from current Fremantle build
[opencv] / src / cv / cvpyramids.cpp
diff --git a/src/cv/cvpyramids.cpp b/src/cv/cvpyramids.cpp
new file mode 100644 (file)
index 0000000..ab7c4d0
--- /dev/null
@@ -0,0 +1,586 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////\r
+//\r
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
+//\r
+//  By downloading, copying, installing or using the software you agree to this license.\r
+//  If you do not agree to this license, do not download, install,\r
+//  copy or use the software.\r
+//\r
+//\r
+//                           License Agreement\r
+//                For Open Source Computer Vision Library\r
+//\r
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.\r
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.\r
+// Third party copyrights are property of their respective owners.\r
+//\r
+// Redistribution and use in source and binary forms, with or without modification,\r
+// are permitted provided that the following conditions are met:\r
+//\r
+//   * Redistribution's of source code must retain the above copyright notice,\r
+//     this list of conditions and the following disclaimer.\r
+//\r
+//   * Redistribution's in binary form must reproduce the above copyright notice,\r
+//     this list of conditions and the following disclaimer in the documentation\r
+//     and/or other materials provided with the distribution.\r
+//\r
+//   * The name of the copyright holders may not be used to endorse or promote products\r
+//     derived from this software without specific prior written permission.\r
+//\r
+// This software is provided by the copyright holders and contributors "as is" and\r
+// any express or implied warranties, including, but not limited to, the implied\r
+// warranties of merchantability and fitness for a particular purpose are disclaimed.\r
+// In no event shall the Intel Corporation or contributors be liable for any direct,\r
+// indirect, incidental, special, exemplary, or consequential damages\r
+// (including, but not limited to, procurement of substitute goods or services;\r
+// loss of use, data, or profits; or business interruption) however caused\r
+// and on any theory of liability, whether in contract, strict liability,\r
+// or tort (including negligence or otherwise) arising in any way out of\r
+// the use of this software, even if advised of the possibility of such damage.\r
+//\r
+//M*/\r
+\r
+#include "_cv.h"\r
+\r
+namespace cv\r
+{\r
+\r
+//#undef CV_SSE2\r
+\r
+template<typename T, int shift> struct FixPtCast\r
+{\r
+    typedef int type1;\r
+    typedef T rtype;\r
+    rtype operator ()(type1 arg) const { return (T)((arg + (1 << (shift-1))) >> shift); }\r
+};\r
+\r
+template<typename T, int shift> struct FltCast\r
+{\r
+    typedef T type1;\r
+    typedef T rtype;\r
+    rtype operator ()(type1 arg) const { return arg*(T)(1./(1 << shift)); }\r
+};\r
+\r
+struct NoVec\r
+{\r
+    int operator()(const uchar**, uchar*, int, int) const { return 0; }\r
+};\r
+\r
+#if CV_SSE2\r
+\r
+struct PyrDownVec_32s8u\r
+{\r
+    int operator()(const uchar** _src, uchar* dst, int, int width) const\r
+    {\r
+        int x = 0;\r
+        const int** src = (const int**)_src;\r
+        const int *row0 = src[0], *row1 = src[1], *row2 = src[2], *row3 = src[3], *row4 = src[4];\r
+        __m128i delta = _mm_set1_epi16(128);\r
+\r
+        for( ; x <= width - 16; x += 16 )\r
+        {\r
+            __m128i r0, r1, r2, r3, r4, t0, t1;\r
+            r0 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row0 + x)),\r
+                                 _mm_load_si128((const __m128i*)(row0 + x + 4)));\r
+            r1 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row1 + x)),\r
+                                 _mm_load_si128((const __m128i*)(row1 + x + 4)));\r
+            r2 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row2 + x)),\r
+                                 _mm_load_si128((const __m128i*)(row2 + x + 4)));\r
+            r3 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row3 + x)),\r
+                                 _mm_load_si128((const __m128i*)(row3 + x + 4)));\r
+            r4 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row4 + x)),\r
+                                 _mm_load_si128((const __m128i*)(row4 + x + 4)));\r
+            r0 = _mm_add_epi16(r0, r4);\r
+            r1 = _mm_add_epi16(_mm_add_epi16(r1, r3), r2);\r
+            r0 = _mm_add_epi16(r0, _mm_add_epi16(r2, r2));\r
+            t0 = _mm_add_epi16(r0, _mm_slli_epi16(r1, 2));\r
+            r0 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row0 + x + 8)),\r
+                                 _mm_load_si128((const __m128i*)(row0 + x + 12)));\r
+            r1 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row1 + x + 8)),\r
+                                 _mm_load_si128((const __m128i*)(row1 + x + 12)));\r
+            r2 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row2 + x + 8)),\r
+                                 _mm_load_si128((const __m128i*)(row2 + x + 12)));\r
+            r3 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row3 + x + 8)),\r
+                                 _mm_load_si128((const __m128i*)(row3 + x + 12)));\r
+            r4 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row4 + x + 8)),\r
+                                 _mm_load_si128((const __m128i*)(row4 + x + 12)));\r
+            r0 = _mm_add_epi16(r0, r4);\r
+            r1 = _mm_add_epi16(_mm_add_epi16(r1, r3), r2);\r
+            r0 = _mm_add_epi16(r0, _mm_add_epi16(r2, r2));\r
+            t1 = _mm_add_epi16(r0, _mm_slli_epi16(r1, 2));\r
+            t0 = _mm_srli_epi16(_mm_add_epi16(t0, delta), 8);\r
+            t1 = _mm_srli_epi16(_mm_add_epi16(t1, delta), 8);\r
+            _mm_storeu_si128((__m128i*)(dst + x), _mm_packus_epi16(t0, t1));\r
+        }\r
+\r
+        for( ; x <= width - 4; x += 4 )\r
+        {\r
+            __m128i r0, r1, r2, r3, r4, z = _mm_setzero_si128();\r
+            r0 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row0 + x)), z);\r
+            r1 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row1 + x)), z);\r
+            r2 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row2 + x)), z);\r
+            r3 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row3 + x)), z);\r
+            r4 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(row4 + x)), z);\r
+            r0 = _mm_add_epi16(r0, r4);\r
+            r1 = _mm_add_epi16(_mm_add_epi16(r1, r3), r2);\r
+            r0 = _mm_add_epi16(r0, _mm_add_epi16(r2, r2));\r
+            r0 = _mm_add_epi16(r0, _mm_slli_epi16(r1, 2));\r
+            r0 = _mm_srli_epi16(_mm_add_epi16(r0, delta), 8);\r
+            *(int*)(dst + x) = _mm_cvtsi128_si32(_mm_packus_epi16(r0, r0));\r
+        }\r
+\r
+        return x;\r
+    }\r
+};\r
+\r
+struct PyrDownVec_32f\r
+{\r
+    int operator()(const uchar** _src, uchar* _dst, int, int width) const\r
+    {\r
+        int x = 0;\r
+        const float** src = (const float**)_src;\r
+        float* dst = (float*)_dst;\r
+        const float *row0 = src[0], *row1 = src[1], *row2 = src[2], *row3 = src[3], *row4 = src[4];\r
+        __m128 _4 = _mm_set1_ps(4.f), _scale = _mm_set1_ps(1.f/256);\r
+        for( ; x <= width - 8; x += 8 )\r
+        {\r
+            __m128 r0, r1, r2, r3, r4, t0, t1;\r
+            r0 = _mm_load_ps(row0 + x);\r
+            r1 = _mm_load_ps(row1 + x);\r
+            r2 = _mm_load_ps(row2 + x);\r
+            r3 = _mm_load_ps(row3 + x);\r
+            r4 = _mm_load_ps(row4 + x);\r
+            r0 = _mm_add_ps(r0, r4);\r
+            r1 = _mm_add_ps(_mm_add_ps(r1, r3), r2);\r
+            r0 = _mm_add_ps(r0, _mm_add_ps(r2, r2));\r
+            t0 = _mm_add_ps(r0, _mm_mul_ps(r1, _4));\r
+\r
+            r0 = _mm_load_ps(row0 + x + 4);\r
+            r1 = _mm_load_ps(row1 + x + 4);\r
+            r2 = _mm_load_ps(row2 + x + 4);\r
+            r3 = _mm_load_ps(row3 + x + 4);\r
+            r4 = _mm_load_ps(row4 + x + 4);\r
+            r0 = _mm_add_ps(r0, r4);\r
+            r1 = _mm_add_ps(_mm_add_ps(r1, r3), r2);\r
+            r0 = _mm_add_ps(r0, _mm_add_ps(r2, r2));\r
+            t1 = _mm_add_ps(r0, _mm_mul_ps(r1, _4));\r
+\r
+            t0 = _mm_mul_ps(t0, _scale);\r
+            t1 = _mm_mul_ps(t1, _scale);\r
+\r
+            _mm_storeu_ps(dst + x, t0);\r
+            _mm_storeu_ps(dst + x + 4, t1);\r
+        }\r
+\r
+        return x;\r
+    }\r
+};\r
+\r
+#else\r
+\r
+typedef NoVec PyrDownVec_32s8u;\r
+typedef NoVec PyrDownVec_32f;\r
+\r
+#endif\r
+\r
+template<class CastOp, class VecOp> void\r
+pyrDown_( const Mat& _src, Mat& _dst )\r
+{\r
+    const int PD_SZ = 5;\r
+    typedef typename CastOp::type1 WT;\r
+    typedef typename CastOp::rtype T;\r
+\r
+    Size ssize = _src.size(), dsize = _dst.size();\r
+    int cn = _src.channels();\r
+    int bufstep = (int)alignSize(dsize.width*cn, 16);\r
+    AutoBuffer<WT> _buf(bufstep*PD_SZ + 16);\r
+    WT* buf = alignPtr((WT*)_buf, 16);\r
+    int tabL[CV_CN_MAX*(PD_SZ+2)], tabR[CV_CN_MAX*(PD_SZ+2)];\r
+    AutoBuffer<int> _tabM(dsize.width*cn);\r
+    int* tabM = _tabM;\r
+    WT* rows[PD_SZ];\r
+    CastOp castOp;\r
+    VecOp vecOp;\r
+\r
+    CV_Assert( std::abs(dsize.width*2 - ssize.width) <= 2 &&\r
+               std::abs(dsize.height*2 - ssize.height) <= 2 );\r
+    int k, x, sy0 = -PD_SZ/2, sy = sy0, width0 = std::min((ssize.width-PD_SZ/2-1)/2 + 1, dsize.width);\r
+\r
+    for( x = 0; x <= PD_SZ+1; x++ )\r
+    {\r
+        int sx0 = borderInterpolate(x - PD_SZ/2, ssize.width, BORDER_REFLECT_101)*cn;\r
+        int sx1 = borderInterpolate(x + width0*2 - PD_SZ/2, ssize.width, BORDER_REFLECT_101)*cn;\r
+        for( k = 0; k < cn; k++ )\r
+        {\r
+            tabL[x*cn + k] = sx0 + k;\r
+            tabR[x*cn + k] = sx1 + k;\r
+        }\r
+    }\r
+    \r
+    ssize.width *= cn;\r
+    dsize.width *= cn;\r
+    width0 *= cn;\r
+\r
+    for( x = 0; x < dsize.width; x++ )\r
+        tabM[x] = (x/cn)*2*cn + x % cn;\r
+\r
+    for( int y = 0; y < dsize.height; y++ )\r
+    {\r
+        T* dst = (T*)(_dst.data + _dst.step*y);\r
+        WT *row0, *row1, *row2, *row3, *row4;\r
+\r
+        // fill the ring buffer (horizontal convolution and decimation)\r
+        for( ; sy <= y*2 + 2; sy++ )\r
+        {\r
+            WT* row = buf + ((sy - sy0) % PD_SZ)*bufstep;\r
+            int _sy = borderInterpolate(sy, ssize.height, BORDER_REFLECT_101);\r
+            const T* src = (const T*)(_src.data + _src.step*_sy);\r
+            int limit = cn;\r
+            const int* tab = tabL;\r
+\r
+            for( x = 0;;)\r
+            {\r
+                for( ; x < limit; x++ )\r
+                {\r
+                    row[x] = src[tab[x+cn*2]]*6 + (src[tab[x+cn]] + src[tab[x+cn*3]])*4 +\r
+                        src[tab[x]] + src[tab[x+cn*4]];\r
+                }\r
+\r
+                if( x == dsize.width )\r
+                    break;\r
+\r
+                if( cn == 1 )\r
+                {\r
+                    for( ; x < width0; x++ )\r
+                        row[x] = src[x*2]*6 + (src[x*2 - 1] + src[x*2 + 1])*4 +\r
+                            src[x*2 - 2] + src[x*2 + 2];\r
+                }\r
+                else if( cn == 3 )\r
+                {\r
+                    for( ; x < width0; x += 3 )\r
+                    {\r
+                        const T* s = src + x*2;\r
+                        WT t0 = s[0]*6 + (s[-3] + s[3])*4 + s[-6] + s[6];\r
+                        WT t1 = s[1]*6 + (s[-2] + s[4])*4 + s[-5] + s[7];\r
+                        WT t2 = s[2]*6 + (s[-1] + s[5])*4 + s[-4] + s[8];\r
+                        row[x] = t0; row[x+1] = t1; row[x+2] = t2;\r
+                    }\r
+                }\r
+                else if( cn == 4 )\r
+                {\r
+                    for( ; x < width0; x += 4 )\r
+                    {\r
+                        const T* s = src + x*2;\r
+                        WT t0 = s[0]*6 + (s[-4] + s[4])*4 + s[-8] + s[8];\r
+                        WT t1 = s[1]*6 + (s[-3] + s[5])*4 + s[-7] + s[9];\r
+                        row[x] = t0; row[x+1] = t1;\r
+                        t0 = s[2]*6 + (s[-2] + s[6])*4 + s[-6] + s[10];\r
+                        t1 = s[3]*6 + (s[-1] + s[7])*4 + s[-5] + s[11];\r
+                        row[x+2] = t0; row[x+3] = t1;\r
+                    }\r
+                }\r
+                else\r
+                {\r
+                    for( ; x < width0; x++ )\r
+                    {\r
+                        int sx = tabM[x];\r
+                        row[x] = src[sx]*6 + (src[sx - cn] + src[sx + cn])*4 +\r
+                            src[sx - cn*2] + src[sx + cn*2];\r
+                    }\r
+                }\r
+\r
+                limit = dsize.width;\r
+                tab = tabR - x;\r
+            }\r
+        }\r
+\r
+        // do vertical convolution and decimation and write the result to the destination image\r
+        for( k = 0; k < PD_SZ; k++ )\r
+            rows[k] = buf + ((y*2 - PD_SZ/2 + k - sy0) % PD_SZ)*bufstep;\r
+        row0 = rows[0]; row1 = rows[1]; row2 = rows[2]; row3 = rows[3]; row4 = rows[4];\r
+\r
+        x = vecOp((const uchar**)rows, (uchar*)dst, (int)_dst.step, dsize.width);\r
+        for( ; x < dsize.width; x++ )\r
+            dst[x] = castOp(row2[x]*6 + (row1[x] + row3[x])*4 + row0[x] + row4[x]);\r
+    }\r
+}\r
+\r
+\r
+template<class CastOp, class VecOp> void\r
+pyrUp_( const Mat& _src, Mat& _dst )\r
+{\r
+    const int PU_SZ = 3;\r
+    typedef typename CastOp::type1 WT;\r
+    typedef typename CastOp::rtype T;\r
+\r
+    Size ssize = _src.size(), dsize = _dst.size();\r
+    int cn = _src.channels();\r
+    int bufstep = (int)alignSize((dsize.width+1)*cn, 16);\r
+    AutoBuffer<WT> _buf(bufstep*PU_SZ + 16);\r
+    WT* buf = alignPtr((WT*)_buf, 16);\r
+    AutoBuffer<int> _dtab(ssize.width*cn);\r
+    int* dtab = _dtab;\r
+    WT* rows[PU_SZ];\r
+    CastOp castOp;\r
+    VecOp vecOp;\r
+\r
+    CV_Assert( std::abs(dsize.width - ssize.width*2) == dsize.width % 2 &&\r
+               std::abs(dsize.height - ssize.height*2) == dsize.height % 2);\r
+    int k, x, sy0 = -PU_SZ/2, sy = sy0, width0 = ssize.width - 1;\r
+\r
+    ssize.width *= cn;\r
+    dsize.width *= cn;\r
+    width0 *= cn;\r
+\r
+    for( x = 0; x < ssize.width; x++ )\r
+        dtab[x] = (x/cn)*2*cn + x % cn;\r
+\r
+    for( int y = 0; y < ssize.height; y++ )\r
+    {\r
+        T* dst0 = (T*)(_dst.data + _dst.step*y*2);\r
+        T* dst1 = (T*)(_dst.data + _dst.step*(y*2+1));\r
+        WT *row0, *row1, *row2;\r
+\r
+        if( y*2+1 >= dsize.height )\r
+            dst1 = dst0;\r
+\r
+        // fill the ring buffer (horizontal convolution and decimation)\r
+        for( ; sy <= y + 1; sy++ )\r
+        {\r
+            WT* row = buf + ((sy - sy0) % PU_SZ)*bufstep;\r
+            int _sy = borderInterpolate(sy, ssize.height, BORDER_REFLECT_101);\r
+            const T* src = (const T*)(_src.data + _src.step*_sy);\r
+\r
+            if( ssize.width == cn )\r
+            {\r
+                for( x = 0; x < cn; x++ )\r
+                    row[x] = row[x + cn] = src[x]*8;\r
+                continue;\r
+            }\r
+\r
+            for( x = 0; x < cn; x++ )\r
+            {\r
+                int dx = dtab[x];\r
+                WT t0 = src[x]*6 + src[x + cn]*2;\r
+                WT t1 = (src[x] + src[x + cn])*4;\r
+                row[dx] = t0; row[dx + cn] = t1;\r
+                dx = dtab[ssize.width - cn + x];\r
+                int sx = ssize.width - cn + x;\r
+                t0 = src[sx - cn] + src[sx]*7;\r
+                t1 = src[sx]*8;\r
+                row[dx] = t0; row[dx + cn] = t1;\r
+            }\r
+\r
+            for( x = cn; x < ssize.width - cn; x++ )\r
+            {\r
+                int dx = dtab[x];\r
+                WT t0 = src[x-cn] + src[x]*6 + src[x+cn];\r
+                WT t1 = (src[x] + src[x+cn])*4;\r
+                row[dx] = t0;\r
+                row[dx+cn] = t1;\r
+            }\r
+        }\r
+\r
+        // do vertical convolution and decimation and write the result to the destination image\r
+        for( k = 0; k < PU_SZ; k++ )\r
+            rows[k] = buf + ((y - PU_SZ/2 + k - sy0) % PU_SZ)*bufstep;\r
+        row0 = rows[0]; row1 = rows[1]; row2 = rows[2];\r
+\r
+        x = vecOp((const uchar**)rows, (uchar*)dst0, (int)_dst.step, dsize.width);\r
+        for( ; x < dsize.width; x++ )\r
+        {\r
+            T t1 = castOp((row1[x] + row2[x])*4);\r
+            T t0 = castOp(row0[x] + row1[x]*6 + row2[x]);\r
+            dst1[x] = t1; dst0[x] = t0;\r
+        }\r
+    }\r
+}\r
+\r
+typedef void (*PyrFunc)(const Mat&, Mat&);\r
+\r
+void pyrDown( const Mat& _src, Mat& _dst, const Size& _dsz )\r
+{\r
+    Size dsz = _dsz == Size() ? Size((_src.cols + 1)/2, (_src.rows + 1)/2) : _dsz;\r
+    _dst.create( dsz, _src.type() );\r
+    int depth = _src.depth();\r
+    PyrFunc func = 0;\r
+    if( depth == CV_8U )\r
+        func = pyrDown_<FixPtCast<uchar, 8>, PyrDownVec_32s8u>;\r
+    else if( depth == CV_16U )\r
+        func = pyrDown_<FixPtCast<ushort, 8>, NoVec>;\r
+    else if( depth == CV_32F )\r
+        func = pyrDown_<FltCast<float, 8>, PyrDownVec_32f>;\r
+    else if( depth == CV_64F )\r
+        func = pyrDown_<FltCast<double, 8>, NoVec>;\r
+    else\r
+        CV_Error( CV_StsUnsupportedFormat, "" );\r
+\r
+    func( _src, _dst );\r
+}\r
+\r
+void pyrUp( const Mat& _src, Mat& _dst, const Size& _dsz )\r
+{\r
+    Size dsz = _dsz == Size() ? Size(_src.cols*2, _src.rows*2) : _dsz;\r
+    _dst.create( dsz, _src.type() );\r
+    int depth = _src.depth();\r
+    PyrFunc func = 0;\r
+    if( depth == CV_8U )\r
+        func = pyrUp_<FixPtCast<uchar, 6>, NoVec>;\r
+    else if( depth == CV_16U )\r
+        func = pyrUp_<FixPtCast<ushort, 6>, NoVec>;\r
+    else if( depth == CV_32F )\r
+        func = pyrUp_<FltCast<float, 6>, NoVec>;\r
+    else if( depth == CV_64F )\r
+        func = pyrUp_<FltCast<double, 6>, NoVec>;\r
+    else\r
+        CV_Error( CV_StsUnsupportedFormat, "" );\r
+\r
+    func( _src, _dst );\r
+}\r
+\r
+void buildPyramid( const Mat& _src, vector<Mat>& _dst, int maxlevel )\r
+{\r
+    _dst.resize( maxlevel + 1 );\r
+    _dst[0] = _src;\r
+    for( int i = 1; i <= maxlevel; i++ )\r
+        pyrDown( _dst[i-1], _dst[i] );\r
+}\r
+\r
+}\r
+\r
+CV_IMPL void cvPyrDown( const void* srcarr, void* dstarr, int _filter )\r
+{\r
+    cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);\r
+\r
+    CV_Assert( _filter == CV_GAUSSIAN_5x5 && src.type() == dst.type());\r
+    cv::pyrDown( src, dst, dst.size() );\r
+}\r
+\r
+CV_IMPL void cvPyrUp( const void* srcarr, void* dstarr, int _filter )\r
+{\r
+    cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);\r
+\r
+    CV_Assert( _filter == CV_GAUSSIAN_5x5 && src.type() == dst.type());\r
+    cv::pyrUp( src, dst, dst.size() );\r
+}\r
+\r
+\r
+CV_IMPL void\r
+cvReleasePyramid( CvMat*** _pyramid, int extra_layers )\r
+{\r
+    CV_FUNCNAME( "cvReleasePyramid" );\r
+\r
+    __BEGIN__;\r
+\r
+    CvMat** pyramid;\r
+    int i;\r
+\r
+    if( !_pyramid )\r
+        CV_ERROR( CV_StsNullPtr, "" );\r
+\r
+    pyramid = *_pyramid;\r
+    \r
+    if( pyramid )\r
+    {\r
+        for( i = 0; i <= extra_layers; i++ )\r
+            cvReleaseMat( &pyramid[i] );\r
+    }\r
+    \r
+    cvFree( _pyramid );\r
+\r
+    __END__;\r
+}\r
+\r
+\r
+CV_IMPL CvMat**\r
+cvCreatePyramid( const CvArr* srcarr, int extra_layers, double rate,\r
+                 const CvSize* layer_sizes, CvArr* bufarr,\r
+                 int calc, int filter )\r
+{\r
+    CvMat** pyramid = 0;\r
+    const float eps = 0.1f;\r
+\r
+    CV_FUNCNAME( "cvCreatePyramid" );\r
+\r
+    __BEGIN__;\r
+    \r
+    int i, elem_size, layer_step;\r
+    CvMat stub, *src;\r
+    CvSize size, layer_size;\r
+    uchar* ptr = 0;\r
+\r
+    CV_CALL( src = cvGetMat( srcarr, &stub ));\r
+\r
+    if( extra_layers < 0 )\r
+        CV_ERROR( CV_StsOutOfRange, "The number of extra layers must be non negative" );\r
+\r
+    elem_size = CV_ELEM_SIZE(src->type);\r
+    size = cvGetMatSize(src);\r
+\r
+    if( bufarr )\r
+    {\r
+        CvMat bstub, *buf;\r
+        int bufsize = 0;\r
+\r
+        CV_CALL( buf = cvGetMat( bufarr, &bstub ));\r
+        bufsize = buf->rows*buf->cols*CV_ELEM_SIZE(buf->type);\r
+        layer_size = size;\r
+        for( i = 1; i <= extra_layers; i++ )\r
+        {\r
+            if( !layer_sizes )\r
+            {\r
+                layer_size.width = cvRound(layer_size.width*rate+eps);\r
+                layer_size.height = cvRound(layer_size.height*rate+eps);\r
+            }\r
+            else\r
+                layer_size = layer_sizes[i-1];\r
+            layer_step = layer_size.width*elem_size;\r
+            bufsize -= layer_step*layer_size.height;\r
+        }\r
+\r
+        if( bufsize < 0 )\r
+            CV_ERROR( CV_StsOutOfRange, "The buffer is too small to fit the pyramid" );\r
+        ptr = buf->data.ptr;\r
+    }\r
+\r
+    CV_CALL( pyramid = (CvMat**)cvAlloc( (extra_layers+1)*sizeof(pyramid[0]) ));\r
+    memset( pyramid, 0, (extra_layers+1)*sizeof(pyramid[0]) );\r
+\r
+    pyramid[0] = cvCreateMatHeader( size.height, size.width, src->type );\r
+    cvSetData( pyramid[0], src->data.ptr, src->step );\r
+    layer_size = size;\r
+\r
+    for( i = 1; i <= extra_layers; i++ )\r
+    {\r
+        if( !layer_sizes )\r
+        {\r
+            layer_size.width = cvRound(layer_size.width*rate + eps);\r
+            layer_size.height = cvRound(layer_size.height*rate + eps);\r
+        }\r
+        else\r
+            layer_size = layer_sizes[i];\r
+\r
+        if( bufarr )\r
+        {\r
+            pyramid[i] = cvCreateMatHeader( layer_size.height, layer_size.width, src->type );\r
+            layer_step = layer_size.width*elem_size;\r
+            cvSetData( pyramid[i], ptr, layer_step );\r
+            ptr += layer_step*layer_size.height;\r
+        }\r
+        else\r
+            pyramid[i] = cvCreateMat( layer_size.height, layer_size.width, src->type );\r
+\r
+        if( calc )\r
+            cvPyrDown( pyramid[i-1], pyramid[i], filter );\r
+            //cvResize( pyramid[i-1], pyramid[i], CV_INTER_LINEAR );\r
+    }\r
+    \r
+    __END__;\r
+\r
+    if( cvGetErrStatus() < 0 )\r
+        cvReleasePyramid( &pyramid, extra_layers );\r
+\r
+    return pyramid;\r
+}\r
+\r
+/* End of file. */\r