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 "_cv.h"
43
44#define ICV_DEF_INTEGRAL_OP_C1( flavor, arrtype, sumtype, sqsumtype, worktype,  \
45                                cast_macro, cast_sqr_macro )    \
46static CvStatus CV_STDCALL                                      \
47icvIntegralImage_##flavor##_C1R( const arrtype* src, int srcstep,\
48                                 sumtype* sum, int sumstep,     \
49                                 sqsumtype* sqsum, int sqsumstep,\
50                                 sumtype* tilted, int tiltedstep,\
51                                 CvSize size )                  \
52{                                                               \
53    int x, y;                                                   \
54    sumtype s;                                                  \
55    sqsumtype sq;                                               \
56    sumtype* buf = 0;                                           \
57                                                                \
58    srcstep /= sizeof(src[0]);                                  \
59                                                                \
60    memset( sum, 0, (size.width+1)*sizeof(sum[0]));             \
61    sumstep /= sizeof(sum[0]);                                  \
62    sum += sumstep + 1;                                         \
63                                                                \
64    if( sqsum )                                                 \
65    {                                                           \
66        memset( sqsum, 0, (size.width+1)*sizeof(sqsum[0]));     \
67        sqsumstep /= sizeof(sqsum[0]);                          \
68        sqsum += sqsumstep + 1;                                 \
69    }                                                           \
70                                                                \
71    if( tilted )                                                \
72    {                                                           \
73        memset( tilted, 0, (size.width+1)*sizeof(tilted[0]));   \
74        tiltedstep /= sizeof(tilted[0]);                        \
75        tilted += tiltedstep + 1;                               \
76    }                                                           \
77                                                                \
78    if( sqsum == 0 && tilted == 0 )                             \
79    {                                                           \
80        for( y = 0; y < size.height; y++, src += srcstep,       \
81                                          sum += sumstep )      \
82        {                                                       \
83            sum[-1] = 0;                                        \
84            for( x = 0, s = 0; x < size.width; x++ )            \
85            {                                                   \
86                sumtype t = cast_macro(src[x]);                 \
87                s += t;                                         \
88                sum[x] = sum[x - sumstep] + s;                  \
89            }                                                   \
90        }                                                       \
91    }                                                           \
92    else if( tilted == 0 )                                      \
93    {                                                           \
94        for( y = 0; y < size.height; y++, src += srcstep,       \
95                        sum += sumstep, sqsum += sqsumstep )    \
96        {                                                       \
97            sum[-1] = 0;                                        \
98            sqsum[-1] = 0;                                      \
99                                                                \
100            for( x = 0, s = 0, sq = 0; x < size.width; x++ )    \
101            {                                                   \
102                worktype it = src[x];                           \
103                sumtype t = cast_macro(it);                     \
104                sqsumtype tq = cast_sqr_macro(it);              \
105                s += t;                                         \
106                sq += tq;                                       \
107                t = sum[x - sumstep] + s;                       \
108                tq = sqsum[x - sqsumstep] + sq;                 \
109                sum[x] = t;                                     \
110                sqsum[x] = tq;                                  \
111            }                                                   \
112        }                                                       \
113    }                                                           \
114    else                                                        \
115    {                                                           \
116        if( sqsum == 0 )                                        \
117        {                                                       \
118            assert(0);                                          \
119            return CV_NULLPTR_ERR;                              \
120        }                                                       \
121                                                                \
122        buf = (sumtype*)cvStackAlloc((size.width + 1 )* sizeof(buf[0]));\
123        sum[-1] = tilted[-1] = 0;                               \
124        sqsum[-1] = 0;                                          \
125                                                                \
126        for( x = 0, s = 0, sq = 0; x < size.width; x++ )        \
127        {                                                       \
128            worktype it = src[x];                               \
129            sumtype t = cast_macro(it);                         \
130            sqsumtype tq = cast_sqr_macro(it);                  \
131            buf[x] = tilted[x] = t;                             \
132            s += t;                                             \
133            sq += tq;                                           \
134            sum[x] = s;                                         \
135            sqsum[x] = sq;                                      \
136        }                                                       \
137                                                                \
138        if( size.width == 1 )                                   \
139            buf[1] = 0;                                         \
140                                                                \
141        for( y = 1; y < size.height; y++ )                      \
142        {                                                       \
143            worktype it;                                        \
144            sumtype t0;                                         \
145            sqsumtype tq0;                                      \
146                                                                \
147            src += srcstep;                                     \
148            sum += sumstep;                                     \
149            sqsum += sqsumstep;                                 \
150            tilted += tiltedstep;                               \
151                                                                \
152            it = src[0/*x*/];                                   \
153            s = t0 = cast_macro(it);                            \
154            sq = tq0 = cast_sqr_macro(it);                      \
155                                                                \
156            sum[-1] = 0;                                        \
157            sqsum[-1] = 0;                                      \
158            /*tilted[-1] = buf[0];*/                            \
159            tilted[-1] = tilted[-tiltedstep];                   \
160                                                                \
161            sum[0] = sum[-sumstep] + t0;                        \
162            sqsum[0] = sqsum[-sqsumstep] + tq0;                 \
163            tilted[0] = tilted[-tiltedstep] + t0 + buf[1];      \
164                                                                \
165            for( x = 1; x < size.width - 1; x++ )               \
166            {                                                   \
167                sumtype t1 = buf[x];                            \
168                buf[x-1] = t1 + t0;                             \
169                it = src[x];                                    \
170                t0 = cast_macro(it);                            \
171                tq0 = cast_sqr_macro(it);                       \
172                s += t0;                                        \
173                sq += tq0;                                      \
174                sum[x] = sum[x - sumstep] + s;                  \
175                sqsum[x] = sqsum[x - sqsumstep] + sq;           \
176                t1 += buf[x+1] + t0 + tilted[x - tiltedstep - 1];\
177                tilted[x] = t1;                                 \
178            }                                                   \
179                                                                \
180            if( size.width > 1 )                                \
181            {                                                   \
182                sumtype t1 = buf[x];                            \
183                buf[x-1] = t1 + t0;                             \
184                it = src[x];    /*+*/                           \
185                t0 = cast_macro(it);                            \
186                tq0 = cast_sqr_macro(it);                       \
187                s += t0;                                        \
188                sq += tq0;                                      \
189                sum[x] = sum[x - sumstep] + s;                  \
190                sqsum[x] = sqsum[x - sqsumstep] + sq;           \
191                tilted[x] = t0 + t1 + tilted[x - tiltedstep - 1];\
192                buf[x] = t0;                                    \
193            }                                                   \
194        }                                                       \
195    }                                                           \
196                                                                \
197    return CV_OK;                                               \
198}
199
200
201ICV_DEF_INTEGRAL_OP_C1( 8u32s, uchar, int, double, int, CV_NOP, CV_8TO32F_SQR )
202ICV_DEF_INTEGRAL_OP_C1( 8u64f, uchar, double, double, int, CV_8TO32F, CV_8TO32F_SQR )
203ICV_DEF_INTEGRAL_OP_C1( 32f64f, float, double, double, double, CV_NOP, CV_SQR )
204ICV_DEF_INTEGRAL_OP_C1( 64f, double, double, double, double, CV_NOP, CV_SQR )
205
206
207#define ICV_DEF_INTEGRAL_OP_CN( flavor, arrtype, sumtype, sqsumtype,    \
208                                worktype, cast_macro, cast_sqr_macro )  \
209static CvStatus CV_STDCALL                                      \
210icvIntegralImage_##flavor##_CnR( const arrtype* src, int srcstep,\
211                                 sumtype* sum, int sumstep,     \
212                                 sqsumtype* sqsum, int sqsumstep,\
213                                 CvSize size, int cn )          \
214{                                                               \
215    int x, y;                                                   \
216    srcstep /= sizeof(src[0]);                                  \
217                                                                \
218    memset( sum, 0, (size.width+1)*cn*sizeof(sum[0]));          \
219    sumstep /= sizeof(sum[0]);                                  \
220    sum += sumstep + cn;                                        \
221                                                                \
222    if( sqsum )                                                 \
223    {                                                           \
224        memset( sqsum, 0, (size.width+1)*cn*sizeof(sqsum[0]));  \
225        sqsumstep /= sizeof(sqsum[0]);                          \
226        sqsum += sqsumstep + cn;                                \
227    }                                                           \
228                                                                \
229    size.width *= cn;                                           \
230                                                                \
231    if( sqsum == 0 )                                            \
232    {                                                           \
233        for( y = 0; y < size.height; y++, src += srcstep,       \
234                                          sum += sumstep )      \
235        {                                                       \
236            for( x = -cn; x < 0; x++ )                          \
237                sum[x] = 0;                                     \
238                                                                \
239            for( x = 0; x < size.width; x++ )                   \
240                sum[x] = cast_macro(src[x]) + sum[x - cn];      \
241                                                                \
242            for( x = 0; x < size.width; x++ )                   \
243                sum[x] = sum[x] + sum[x - sumstep];             \
244        }                                                       \
245    }                                                           \
246    else                                                        \
247    {                                                           \
248        for( y = 0; y < size.height; y++, src += srcstep,       \
249                        sum += sumstep, sqsum += sqsumstep )    \
250        {                                                       \
251            for( x = -cn; x < 0; x++ )                          \
252            {                                                   \
253                sum[x] = 0;                                     \
254                sqsum[x] = 0;                                   \
255            }                                                   \
256                                                                \
257            for( x = 0; x < size.width; x++ )                   \
258            {                                                   \
259                worktype it = src[x];                           \
260                sumtype t = cast_macro(it) + sum[x-cn];         \
261                sqsumtype tq = cast_sqr_macro(it) + sqsum[x-cn];\
262                sum[x] = t;                                     \
263                sqsum[x] = tq;                                  \
264            }                                                   \
265                                                                \
266            for( x = 0; x < size.width; x++ )                   \
267            {                                                   \
268                sumtype t = sum[x] + sum[x - sumstep];          \
269                sqsumtype tq = sqsum[x] + sqsum[x - sqsumstep]; \
270                sum[x] = t;                                     \
271                sqsum[x] = tq;                                  \
272            }                                                   \
273        }                                                       \
274    }                                                           \
275                                                                \
276    return CV_OK;                                               \
277}
278
279
280ICV_DEF_INTEGRAL_OP_CN( 8u32s, uchar, int, double, int, CV_NOP, CV_8TO32F_SQR )
281ICV_DEF_INTEGRAL_OP_CN( 8u64f, uchar, double, double, int, CV_8TO32F, CV_8TO32F_SQR )
282ICV_DEF_INTEGRAL_OP_CN( 32f64f, float, double, double, double, CV_NOP, CV_SQR )
283ICV_DEF_INTEGRAL_OP_CN( 64f, double, double, double, double, CV_NOP, CV_SQR )
284
285
286static void icvInitIntegralImageTable( CvFuncTable* table_c1, CvFuncTable* table_cn )
287{
288    table_c1->fn_2d[CV_8U] = (void*)icvIntegralImage_8u64f_C1R;
289    table_c1->fn_2d[CV_32F] = (void*)icvIntegralImage_32f64f_C1R;
290    table_c1->fn_2d[CV_64F] = (void*)icvIntegralImage_64f_C1R;
291
292    table_cn->fn_2d[CV_8U] = (void*)icvIntegralImage_8u64f_CnR;
293    table_cn->fn_2d[CV_32F] = (void*)icvIntegralImage_32f64f_CnR;
294    table_cn->fn_2d[CV_64F] = (void*)icvIntegralImage_64f_CnR;
295}
296
297
298typedef CvStatus (CV_STDCALL * CvIntegralImageFuncC1)(
299    const void* src, int srcstep, void* sum, int sumstep,
300    void* sqsum, int sqsumstep, void* tilted, int tiltedstep,
301    CvSize size );
302
303typedef CvStatus (CV_STDCALL * CvIntegralImageFuncCn)(
304    const void* src, int srcstep, void* sum, int sumstep,
305    void* sqsum, int sqsumstep, CvSize size, int cn );
306
307icvIntegral_8u32s_C1R_t icvIntegral_8u32s_C1R_p = 0;
308icvSqrIntegral_8u32s64f_C1R_t icvSqrIntegral_8u32s64f_C1R_p = 0;
309
310CV_IMPL void
311cvIntegral( const CvArr* image, CvArr* sumImage,
312            CvArr* sumSqImage, CvArr* tiltedSumImage )
313{
314    static CvFuncTable tab_c1, tab_cn;
315    static int inittab = 0;
316
317    CV_FUNCNAME( "cvIntegralImage" );
318
319    __BEGIN__;
320
321    CvMat src_stub, *src = (CvMat*)image;
322    CvMat sum_stub, *sum = (CvMat*)sumImage;
323    CvMat sqsum_stub, *sqsum = (CvMat*)sumSqImage;
324    CvMat tilted_stub, *tilted = (CvMat*)tiltedSumImage;
325    int coi0 = 0, coi1 = 0, coi2 = 0, coi3 = 0;
326    int depth, cn;
327    int src_step, sum_step, sqsum_step, tilted_step;
328    CvIntegralImageFuncC1 func_c1 = 0;
329    CvIntegralImageFuncCn func_cn = 0;
330    CvSize size;
331
332    if( !inittab )
333    {
334        icvInitIntegralImageTable( &tab_c1, &tab_cn );
335        inittab = 1;
336    }
337
338    CV_CALL( src = cvGetMat( src, &src_stub, &coi0 ));
339    CV_CALL( sum = cvGetMat( sum, &sum_stub, &coi1 ));
340
341    if( sum->width != src->width + 1 ||
342        sum->height != src->height + 1 )
343        CV_ERROR( CV_StsUnmatchedSizes, "" );
344
345    if( (CV_MAT_DEPTH( sum->type ) != CV_64F &&
346        (CV_MAT_DEPTH( src->type ) != CV_8U ||
347         CV_MAT_DEPTH( sum->type ) != CV_32S )) ||
348        !CV_ARE_CNS_EQ( src, sum ))
349        CV_ERROR( CV_StsUnsupportedFormat,
350        "Sum array must have 64f type (or 32s type in case of 8u source array) "
351        "and the same number of channels as the source array" );
352
353    if( sqsum )
354    {
355        CV_CALL( sqsum = cvGetMat( sqsum, &sqsum_stub, &coi2 ));
356        if( !CV_ARE_SIZES_EQ( sum, sqsum ) )
357            CV_ERROR( CV_StsUnmatchedSizes, "" );
358        if( CV_MAT_DEPTH( sqsum->type ) != CV_64F || !CV_ARE_CNS_EQ( src, sqsum ))
359            CV_ERROR( CV_StsUnsupportedFormat,
360                      "Squares sum array must be 64f "
361                      "and the same number of channels as the source array" );
362    }
363
364    if( tilted )
365    {
366        if( !sqsum )
367            CV_ERROR( CV_StsNullPtr,
368            "Squared sum array must be passed if tilted sum array is passed" );
369
370        CV_CALL( tilted = cvGetMat( tilted, &tilted_stub, &coi3 ));
371        if( !CV_ARE_SIZES_EQ( sum, tilted ) )
372            CV_ERROR( CV_StsUnmatchedSizes, "" );
373        if( !CV_ARE_TYPES_EQ( sum, tilted ) )
374            CV_ERROR( CV_StsUnmatchedFormats,
375                      "Sum and tilted sum must have the same types" );
376        if( CV_MAT_CN(tilted->type) != 1 )
377            CV_ERROR( CV_StsNotImplemented,
378                      "Tilted sum can not be computed for multi-channel arrays" );
379    }
380
381    if( coi0 || coi1 || coi2 || coi3 )
382        CV_ERROR( CV_BadCOI, "COI is not supported by the function" );
383
384    depth = CV_MAT_DEPTH(src->type);
385    cn = CV_MAT_CN(src->type);
386
387    if( CV_MAT_DEPTH( sum->type ) == CV_32S )
388    {
389        func_c1 = (CvIntegralImageFuncC1)icvIntegralImage_8u32s_C1R;
390        func_cn = (CvIntegralImageFuncCn)icvIntegralImage_8u32s_CnR;
391    }
392    else
393    {
394        func_c1 = (CvIntegralImageFuncC1)tab_c1.fn_2d[depth];
395        func_cn = (CvIntegralImageFuncCn)tab_cn.fn_2d[depth];
396        if( !func_c1 && !func_cn )
397            CV_ERROR( CV_StsUnsupportedFormat, "This source image format is unsupported" );
398    }
399
400    size = cvGetMatSize(src);
401    src_step = src->step ? src->step : CV_STUB_STEP;
402    sum_step = sum->step ? sum->step : CV_STUB_STEP;
403    sqsum_step = !sqsum ? 0 : sqsum->step ? sqsum->step : CV_STUB_STEP;
404    tilted_step = !tilted ? 0 : tilted->step ? tilted->step : CV_STUB_STEP;
405
406    if( cn == 1 )
407    {
408        if( depth == CV_8U && !tilted && CV_MAT_DEPTH(sum->type) == CV_32S )
409        {
410            if( !sqsum && icvIntegral_8u32s_C1R_p &&
411                icvIntegral_8u32s_C1R_p( src->data.ptr, src_step,
412                            sum->data.i, sum_step, size, 0 ) >= 0 )
413                EXIT;
414
415            if( sqsum && icvSqrIntegral_8u32s64f_C1R_p &&
416                icvSqrIntegral_8u32s64f_C1R_p( src->data.ptr, src_step, sum->data.i,
417                            sum_step, sqsum->data.db, sqsum_step, size, 0, 0 ) >= 0 )
418                EXIT;
419        }
420
421        IPPI_CALL( func_c1( src->data.ptr, src_step, sum->data.ptr, sum_step,
422                        sqsum ? sqsum->data.ptr : 0, sqsum_step,
423                        tilted ? tilted->data.ptr : 0, tilted_step, size ));
424    }
425    else
426    {
427        IPPI_CALL( func_cn( src->data.ptr, src_step, sum->data.ptr, sum_step,
428                        sqsum ? sqsum->data.ptr : 0, sqsum_step, size, cn ));
429    }
430
431    __END__;
432}
433
434
435/* End of file. */
436