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#include <stdio.h>
44
45
46static void
47icvCalcMinEigenVal( const float* cov, int cov_step, float* dst,
48                    int dst_step, CvSize size, CvMat* buffer )
49{
50    int j;
51    float* buf = buffer->data.fl;
52    cov_step /= sizeof(cov[0]);
53    dst_step /= sizeof(dst[0]);
54    buffer->rows = 1;
55
56    for( ; size.height--; cov += cov_step, dst += dst_step )
57    {
58        for( j = 0; j < size.width; j++ )
59        {
60            double a = cov[j*3]*0.5;
61            double b = cov[j*3+1];
62            double c = cov[j*3+2]*0.5;
63
64            buf[j + size.width] = (float)(a + c);
65            buf[j] = (float)((a - c)*(a - c) + b*b);
66        }
67
68        cvPow( buffer, buffer, 0.5 );
69
70        for( j = 0; j < size.width ; j++ )
71            dst[j] = (float)(buf[j + size.width] - buf[j]);
72    }
73}
74
75
76static void
77icvCalcHarris( const float* cov, int cov_step, float* dst,
78               int dst_step, CvSize size, CvMat* /*buffer*/, double k )
79{
80    int j;
81    cov_step /= sizeof(cov[0]);
82    dst_step /= sizeof(dst[0]);
83
84    for( ; size.height--; cov += cov_step, dst += dst_step )
85    {
86        for( j = 0; j < size.width; j++ )
87        {
88            double a = cov[j*3];
89            double b = cov[j*3+1];
90            double c = cov[j*3+2];
91            dst[j] = (float)(a*c - b*b - k*(a + c)*(a + c));
92        }
93    }
94}
95
96
97static void
98icvCalcEigenValsVecs( const float* cov, int cov_step, float* dst,
99                      int dst_step, CvSize size, CvMat* buffer )
100{
101    static int y0 = 0;
102
103    int j;
104    float* buf = buffer->data.fl;
105    cov_step /= sizeof(cov[0]);
106    dst_step /= sizeof(dst[0]);
107
108    for( ; size.height--; cov += cov_step, dst += dst_step )
109    {
110        for( j = 0; j < size.width; j++ )
111        {
112            double a = cov[j*3]*0.5;
113            double b = cov[j*3+1];
114            double c = cov[j*3+2]*0.5;
115
116            buf[j + size.width] = (float)(a + c);
117            buf[j] = (float)((a - c)*(a - c) + b*b);
118        }
119
120        buffer->rows = 1;
121        cvPow( buffer, buffer, 0.5 );
122
123        for( j = 0; j < size.width; j++ )
124        {
125            double a = cov[j*3];
126            double b = cov[j*3+1];
127            double c = cov[j*3+2];
128
129            double l1 = buf[j + size.width] + buf[j];
130            double l2 = buf[j + size.width] - buf[j];
131
132            double x = b;
133            double y = l1 - a;
134            double e = fabs(x);
135
136            if( e + fabs(y) < 1e-4 )
137            {
138                y = b;
139                x = l1 - c;
140                e = fabs(x);
141                if( e + fabs(y) < 1e-4 )
142                {
143                    e = 1./(e + fabs(y) + FLT_EPSILON);
144                    x *= e, y *= e;
145                }
146            }
147
148            buf[j] = (float)(x*x + y*y + DBL_EPSILON);
149            dst[6*j] = (float)l1;
150            dst[6*j + 2] = (float)x;
151            dst[6*j + 3] = (float)y;
152
153            x = b;
154            y = l2 - a;
155            e = fabs(x);
156
157            if( e + fabs(y) < 1e-4 )
158            {
159                y = b;
160                x = l2 - c;
161                e = fabs(x);
162                if( e + fabs(y) < 1e-4 )
163                {
164                    e = 1./(e + fabs(y) + FLT_EPSILON);
165                    x *= e, y *= e;
166                }
167            }
168
169            buf[j + size.width] = (float)(x*x + y*y + DBL_EPSILON);
170            dst[6*j + 1] = (float)l2;
171            dst[6*j + 4] = (float)x;
172            dst[6*j + 5] = (float)y;
173        }
174
175        buffer->rows = 2;
176        cvPow( buffer, buffer, -0.5 );
177
178        for( j = 0; j < size.width; j++ )
179        {
180            double t0 = buf[j]*dst[6*j + 2];
181            double t1 = buf[j]*dst[6*j + 3];
182
183            dst[6*j + 2] = (float)t0;
184            dst[6*j + 3] = (float)t1;
185
186            t0 = buf[j + size.width]*dst[6*j + 4];
187            t1 = buf[j + size.width]*dst[6*j + 5];
188
189            dst[6*j + 4] = (float)t0;
190            dst[6*j + 5] = (float)t1;
191        }
192
193        y0++;
194    }
195}
196
197
198#define ICV_MINEIGENVAL     0
199#define ICV_HARRIS          1
200#define ICV_EIGENVALSVECS   2
201
202static void
203icvCornerEigenValsVecs( const CvMat* src, CvMat* eigenv, int block_size,
204                        int aperture_size, int op_type, double k=0. )
205{
206    CvSepFilter dx_filter, dy_filter;
207    CvBoxFilter blur_filter;
208    CvMat *tempsrc = 0;
209    CvMat *Dx = 0, *Dy = 0, *cov = 0;
210    CvMat *sqrt_buf = 0;
211
212    int buf_size = 1 << 12;
213
214    CV_FUNCNAME( "icvCornerEigenValsVecs" );
215
216    __BEGIN__;
217
218    int i, j, y, dst_y = 0, max_dy, delta = 0;
219    int aperture_size0 = aperture_size;
220    int temp_step = 0, d_step;
221    uchar* shifted_ptr = 0;
222    int depth, d_depth;
223    int stage = CV_START;
224    CvSobelFixedIPPFunc ipp_sobel_vert = 0, ipp_sobel_horiz = 0;
225    CvFilterFixedIPPFunc ipp_scharr_vert = 0, ipp_scharr_horiz = 0;
226    CvSize el_size, size, stripe_size;
227    int aligned_width;
228    CvPoint el_anchor;
229    double factorx, factory;
230    bool use_ipp = false;
231
232    if( block_size < 3 || !(block_size & 1) )
233        CV_ERROR( CV_StsOutOfRange, "averaging window size must be an odd number >= 3" );
234
235    if( (aperture_size < 3 && aperture_size != CV_SCHARR) || !(aperture_size & 1) )
236        CV_ERROR( CV_StsOutOfRange,
237        "Derivative filter aperture size must be a positive odd number >=3 or CV_SCHARR" );
238
239    depth = CV_MAT_DEPTH(src->type);
240    d_depth = depth == CV_8U ? CV_16S : CV_32F;
241
242    size = cvGetMatSize(src);
243    aligned_width = cvAlign(size.width, 4);
244
245    aperture_size = aperture_size == CV_SCHARR ? 3 : aperture_size;
246    el_size = cvSize( aperture_size, aperture_size );
247    el_anchor = cvPoint( aperture_size/2, aperture_size/2 );
248
249    if( aperture_size <= 5 && icvFilterSobelVert_8u16s_C1R_p )
250    {
251        if( depth == CV_8U && aperture_size0 == CV_SCHARR )
252        {
253            ipp_scharr_vert = icvFilterScharrVert_8u16s_C1R_p;
254            ipp_scharr_horiz = icvFilterScharrHoriz_8u16s_C1R_p;
255        }
256        else if( depth == CV_32F && aperture_size0 == CV_SCHARR )
257        {
258            ipp_scharr_vert = icvFilterScharrVert_32f_C1R_p;
259            ipp_scharr_horiz = icvFilterScharrHoriz_32f_C1R_p;
260        }
261        else if( depth == CV_8U )
262        {
263            ipp_sobel_vert = icvFilterSobelVert_8u16s_C1R_p;
264            ipp_sobel_horiz = icvFilterSobelHoriz_8u16s_C1R_p;
265        }
266        else if( depth == CV_32F )
267        {
268            ipp_sobel_vert = icvFilterSobelVert_32f_C1R_p;
269            ipp_sobel_horiz = icvFilterSobelHoriz_32f_C1R_p;
270        }
271    }
272
273    if( (ipp_sobel_vert && ipp_sobel_horiz) ||
274        (ipp_scharr_vert && ipp_scharr_horiz) )
275    {
276        CV_CALL( tempsrc = icvIPPFilterInit( src, buf_size,
277            cvSize(el_size.width,el_size.height + block_size)));
278        shifted_ptr = tempsrc->data.ptr + el_anchor.y*tempsrc->step +
279                      el_anchor.x*CV_ELEM_SIZE(depth);
280        temp_step = tempsrc->step ? tempsrc->step : CV_STUB_STEP;
281        max_dy = tempsrc->rows - aperture_size + 1;
282        use_ipp = true;
283    }
284    else
285    {
286        ipp_sobel_vert = ipp_sobel_horiz = 0;
287        ipp_scharr_vert = ipp_scharr_horiz = 0;
288
289        CV_CALL( dx_filter.init_deriv( size.width, depth, d_depth, 1, 0, aperture_size0 ));
290        CV_CALL( dy_filter.init_deriv( size.width, depth, d_depth, 0, 1, aperture_size0 ));
291        max_dy = buf_size / src->cols;
292        max_dy = MAX( max_dy, aperture_size + block_size );
293    }
294
295    CV_CALL( Dx = cvCreateMat( max_dy, aligned_width, d_depth ));
296    CV_CALL( Dy = cvCreateMat( max_dy, aligned_width, d_depth ));
297    CV_CALL( cov = cvCreateMat( max_dy + block_size + 1, size.width, CV_32FC3 ));
298    CV_CALL( sqrt_buf = cvCreateMat( 2, size.width, CV_32F ));
299    Dx->cols = Dy->cols = size.width;
300
301    if( !use_ipp )
302        max_dy -= aperture_size - 1;
303    d_step = Dx->step ? Dx->step : CV_STUB_STEP;
304
305    CV_CALL(blur_filter.init(size.width, CV_32FC3, CV_32FC3, 0, cvSize(block_size,block_size)));
306    stripe_size = size;
307
308    factorx = (double)(1 << (aperture_size - 1)) * block_size;
309    if( aperture_size0 == CV_SCHARR )
310        factorx *= 2;
311    if( depth == CV_8U )
312        factorx *= 255.;
313    factory = factorx = 1./factorx;
314    if( ipp_sobel_vert )
315        factory = -factory;
316
317    for( y = 0; y < size.height; y += delta )
318    {
319        if( !use_ipp )
320        {
321            delta = MIN( size.height - y, max_dy );
322            if( y + delta == size.height )
323                stage = stage & CV_START ? CV_START + CV_END : CV_END;
324            dx_filter.process( src, Dx, cvRect(0,y,-1,delta), cvPoint(0,0), stage );
325            stripe_size.height = dy_filter.process( src, Dy, cvRect(0,y,-1,delta),
326                                                    cvPoint(0,0), stage );
327        }
328        else
329        {
330            delta = icvIPPFilterNextStripe( src, tempsrc, y, el_size, el_anchor );
331            stripe_size.height = delta;
332
333            if( ipp_sobel_vert )
334            {
335                IPPI_CALL( ipp_sobel_vert( shifted_ptr, temp_step,
336                        Dx->data.ptr, d_step, stripe_size,
337                        aperture_size*10 + aperture_size ));
338                IPPI_CALL( ipp_sobel_horiz( shifted_ptr, temp_step,
339                        Dy->data.ptr, d_step, stripe_size,
340                        aperture_size*10 + aperture_size ));
341            }
342            else /*if( ipp_scharr_vert )*/
343            {
344                IPPI_CALL( ipp_scharr_vert( shifted_ptr, temp_step,
345                        Dx->data.ptr, d_step, stripe_size ));
346                IPPI_CALL( ipp_scharr_horiz( shifted_ptr, temp_step,
347                        Dy->data.ptr, d_step, stripe_size ));
348            }
349        }
350
351        for( i = 0; i < stripe_size.height; i++ )
352        {
353            float* cov_data = (float*)(cov->data.ptr + i*cov->step);
354            if( d_depth == CV_16S )
355            {
356                const short* dxdata = (const short*)(Dx->data.ptr + i*Dx->step);
357                const short* dydata = (const short*)(Dy->data.ptr + i*Dy->step);
358
359                for( j = 0; j < size.width; j++ )
360                {
361                    double dx = dxdata[j]*factorx;
362                    double dy = dydata[j]*factory;
363
364                    cov_data[j*3] = (float)(dx*dx);
365                    cov_data[j*3+1] = (float)(dx*dy);
366                    cov_data[j*3+2] = (float)(dy*dy);
367                }
368            }
369            else
370            {
371                const float* dxdata = (const float*)(Dx->data.ptr + i*Dx->step);
372                const float* dydata = (const float*)(Dy->data.ptr + i*Dy->step);
373
374                for( j = 0; j < size.width; j++ )
375                {
376                    double dx = dxdata[j]*factorx;
377                    double dy = dydata[j]*factory;
378
379                    cov_data[j*3] = (float)(dx*dx);
380                    cov_data[j*3+1] = (float)(dx*dy);
381                    cov_data[j*3+2] = (float)(dy*dy);
382                }
383            }
384        }
385
386        if( y + stripe_size.height >= size.height )
387            stage = stage & CV_START ? CV_START + CV_END : CV_END;
388
389        stripe_size.height = blur_filter.process(cov,cov,
390            cvRect(0,0,-1,stripe_size.height),cvPoint(0,0),stage+CV_ISOLATED_ROI);
391
392        if( op_type == ICV_MINEIGENVAL )
393            icvCalcMinEigenVal( cov->data.fl, cov->step,
394                (float*)(eigenv->data.ptr + dst_y*eigenv->step), eigenv->step,
395                stripe_size, sqrt_buf );
396        else if( op_type == ICV_HARRIS )
397            icvCalcHarris( cov->data.fl, cov->step,
398                (float*)(eigenv->data.ptr + dst_y*eigenv->step), eigenv->step,
399                stripe_size, sqrt_buf, k );
400        else if( op_type == ICV_EIGENVALSVECS )
401            icvCalcEigenValsVecs( cov->data.fl, cov->step,
402                (float*)(eigenv->data.ptr + dst_y*eigenv->step), eigenv->step,
403                stripe_size, sqrt_buf );
404
405        dst_y += stripe_size.height;
406        stage = CV_MIDDLE;
407    }
408
409    __END__;
410
411    cvReleaseMat( &Dx );
412    cvReleaseMat( &Dy );
413    cvReleaseMat( &cov );
414    cvReleaseMat( &sqrt_buf );
415    cvReleaseMat( &tempsrc );
416}
417
418
419CV_IMPL void
420cvCornerMinEigenVal( const void* srcarr, void* eigenvarr,
421                     int block_size, int aperture_size )
422{
423    CV_FUNCNAME( "cvCornerMinEigenVal" );
424
425    __BEGIN__;
426
427    CvMat stub, *src = (CvMat*)srcarr;
428    CvMat eigstub, *eigenv = (CvMat*)eigenvarr;
429
430    CV_CALL( src = cvGetMat( srcarr, &stub ));
431    CV_CALL( eigenv = cvGetMat( eigenv, &eigstub ));
432
433    if( (CV_MAT_TYPE(src->type) != CV_8UC1 && CV_MAT_TYPE(src->type) != CV_32FC1) ||
434        CV_MAT_TYPE(eigenv->type) != CV_32FC1 )
435        CV_ERROR( CV_StsUnsupportedFormat, "Input must be 8uC1 or 32fC1, output must be 32fC1" );
436
437    if( !CV_ARE_SIZES_EQ( src, eigenv ))
438        CV_ERROR( CV_StsUnmatchedSizes, "" );
439
440    CV_CALL( icvCornerEigenValsVecs( src, eigenv, block_size, aperture_size, ICV_MINEIGENVAL ));
441
442    __END__;
443}
444
445
446CV_IMPL void
447cvCornerHarris( const CvArr* srcarr, CvArr* harris_responce,
448                int block_size, int aperture_size, double k )
449{
450    CV_FUNCNAME( "cvCornerHarris" );
451
452    __BEGIN__;
453
454    CvMat stub, *src = (CvMat*)srcarr;
455    CvMat eigstub, *eigenv = (CvMat*)harris_responce;
456
457    CV_CALL( src = cvGetMat( srcarr, &stub ));
458    CV_CALL( eigenv = cvGetMat( eigenv, &eigstub ));
459
460    if( (CV_MAT_TYPE(src->type) != CV_8UC1 && CV_MAT_TYPE(src->type) != CV_32FC1) ||
461        CV_MAT_TYPE(eigenv->type) != CV_32FC1 )
462        CV_ERROR( CV_StsUnsupportedFormat, "Input must be 8uC1 or 32fC1, output must be 32fC1" );
463
464    if( !CV_ARE_SIZES_EQ( src, eigenv ))
465        CV_ERROR( CV_StsUnmatchedSizes, "" );
466
467    CV_CALL( icvCornerEigenValsVecs( src, eigenv, block_size, aperture_size, ICV_HARRIS, k ));
468
469    __END__;
470}
471
472
473CV_IMPL void
474cvCornerEigenValsAndVecs( const void* srcarr, void* eigenvarr,
475                          int block_size, int aperture_size )
476{
477    CV_FUNCNAME( "cvCornerEigenValsAndVecs" );
478
479    __BEGIN__;
480
481    CvMat stub, *src = (CvMat*)srcarr;
482    CvMat eigstub, *eigenv = (CvMat*)eigenvarr;
483
484    CV_CALL( src = cvGetMat( srcarr, &stub ));
485    CV_CALL( eigenv = cvGetMat( eigenv, &eigstub ));
486
487    if( CV_MAT_CN(eigenv->type)*eigenv->cols != src->cols*6 ||
488        eigenv->rows != src->rows )
489        CV_ERROR( CV_StsUnmatchedSizes, "Output array should be 6 times "
490            "wider than the input array and they should have the same height");
491
492    if( (CV_MAT_TYPE(src->type) != CV_8UC1 && CV_MAT_TYPE(src->type) != CV_32FC1) ||
493        CV_MAT_TYPE(eigenv->type) != CV_32FC1 )
494        CV_ERROR( CV_StsUnsupportedFormat, "Input must be 8uC1 or 32fC1, output must be 32fC1" );
495
496    CV_CALL( icvCornerEigenValsVecs( src, eigenv, block_size, aperture_size, ICV_EIGENVALSVECS ));
497
498    __END__;
499}
500
501
502CV_IMPL void
503cvPreCornerDetect( const void* srcarr, void* dstarr, int aperture_size )
504{
505    CvSepFilter dx_filter, dy_filter, d2x_filter, d2y_filter, dxy_filter;
506    CvMat *Dx = 0, *Dy = 0, *D2x = 0, *D2y = 0, *Dxy = 0;
507    CvMat *tempsrc = 0;
508
509    int buf_size = 1 << 12;
510
511    CV_FUNCNAME( "cvPreCornerDetect" );
512
513    __BEGIN__;
514
515    int i, j, y, dst_y = 0, max_dy, delta = 0;
516    int temp_step = 0, d_step;
517    uchar* shifted_ptr = 0;
518    int depth, d_depth;
519    int stage = CV_START;
520    CvSobelFixedIPPFunc ipp_sobel_vert = 0, ipp_sobel_horiz = 0,
521                        ipp_sobel_vert_second = 0, ipp_sobel_horiz_second = 0,
522                        ipp_sobel_cross = 0;
523    CvSize el_size, size, stripe_size;
524    int aligned_width;
525    CvPoint el_anchor;
526    double factor;
527    CvMat stub, *src = (CvMat*)srcarr;
528    CvMat dststub, *dst = (CvMat*)dstarr;
529    bool use_ipp = false;
530
531    CV_CALL( src = cvGetMat( srcarr, &stub ));
532    CV_CALL( dst = cvGetMat( dst, &dststub ));
533
534    if( (CV_MAT_TYPE(src->type) != CV_8UC1 && CV_MAT_TYPE(src->type) != CV_32FC1) ||
535        CV_MAT_TYPE(dst->type) != CV_32FC1 )
536        CV_ERROR( CV_StsUnsupportedFormat, "Input must be 8uC1 or 32fC1, output must be 32fC1" );
537
538    if( !CV_ARE_SIZES_EQ( src, dst ))
539        CV_ERROR( CV_StsUnmatchedSizes, "" );
540
541    if( aperture_size == CV_SCHARR )
542        CV_ERROR( CV_StsOutOfRange, "CV_SCHARR is not supported by this function" );
543
544    if( aperture_size < 3 || aperture_size > 7 || !(aperture_size & 1) )
545        CV_ERROR( CV_StsOutOfRange,
546        "Derivative filter aperture size must be 3, 5 or 7" );
547
548    depth = CV_MAT_DEPTH(src->type);
549    d_depth = depth == CV_8U ? CV_16S : CV_32F;
550
551    size = cvGetMatSize(src);
552    aligned_width = cvAlign(size.width, 4);
553
554    el_size = cvSize( aperture_size, aperture_size );
555    el_anchor = cvPoint( aperture_size/2, aperture_size/2 );
556
557    if( aperture_size <= 5 && icvFilterSobelVert_8u16s_C1R_p )
558    {
559        if( depth == CV_8U )
560        {
561            ipp_sobel_vert = icvFilterSobelVert_8u16s_C1R_p;
562            ipp_sobel_horiz = icvFilterSobelHoriz_8u16s_C1R_p;
563            ipp_sobel_vert_second = icvFilterSobelVertSecond_8u16s_C1R_p;
564            ipp_sobel_horiz_second = icvFilterSobelHorizSecond_8u16s_C1R_p;
565            ipp_sobel_cross = icvFilterSobelCross_8u16s_C1R_p;
566        }
567        else if( depth == CV_32F )
568        {
569            ipp_sobel_vert = icvFilterSobelVert_32f_C1R_p;
570            ipp_sobel_horiz = icvFilterSobelHoriz_32f_C1R_p;
571            ipp_sobel_vert_second = icvFilterSobelVertSecond_32f_C1R_p;
572            ipp_sobel_horiz_second = icvFilterSobelHorizSecond_32f_C1R_p;
573            ipp_sobel_cross = icvFilterSobelCross_32f_C1R_p;
574        }
575    }
576
577    if( ipp_sobel_vert && ipp_sobel_horiz && ipp_sobel_vert_second &&
578        ipp_sobel_horiz_second && ipp_sobel_cross )
579    {
580        CV_CALL( tempsrc = icvIPPFilterInit( src, buf_size, el_size ));
581        shifted_ptr = tempsrc->data.ptr + el_anchor.y*tempsrc->step +
582                      el_anchor.x*CV_ELEM_SIZE(depth);
583        temp_step = tempsrc->step ? tempsrc->step : CV_STUB_STEP;
584        max_dy = tempsrc->rows - aperture_size + 1;
585        use_ipp = true;
586    }
587    else
588    {
589        ipp_sobel_vert = ipp_sobel_horiz = 0;
590        ipp_sobel_vert_second = ipp_sobel_horiz_second = ipp_sobel_cross = 0;
591        dx_filter.init_deriv( size.width, depth, d_depth, 1, 0, aperture_size );
592        dy_filter.init_deriv( size.width, depth, d_depth, 0, 1, aperture_size );
593        d2x_filter.init_deriv( size.width, depth, d_depth, 2, 0, aperture_size );
594        d2y_filter.init_deriv( size.width, depth, d_depth, 0, 2, aperture_size );
595        dxy_filter.init_deriv( size.width, depth, d_depth, 1, 1, aperture_size );
596        max_dy = buf_size / src->cols;
597        max_dy = MAX( max_dy, aperture_size );
598    }
599
600    CV_CALL( Dx = cvCreateMat( max_dy, aligned_width, d_depth ));
601    CV_CALL( Dy = cvCreateMat( max_dy, aligned_width, d_depth ));
602    CV_CALL( D2x = cvCreateMat( max_dy, aligned_width, d_depth ));
603    CV_CALL( D2y = cvCreateMat( max_dy, aligned_width, d_depth ));
604    CV_CALL( Dxy = cvCreateMat( max_dy, aligned_width, d_depth ));
605    Dx->cols = Dy->cols = D2x->cols = D2y->cols = Dxy->cols = size.width;
606
607    if( !use_ipp )
608        max_dy -= aperture_size - 1;
609    d_step = Dx->step ? Dx->step : CV_STUB_STEP;
610
611    stripe_size = size;
612
613    factor = 1 << (aperture_size - 1);
614    if( depth == CV_8U )
615        factor *= 255;
616    factor = 1./(factor * factor * factor);
617
618    aperture_size = aperture_size * 10 + aperture_size;
619
620    for( y = 0; y < size.height; y += delta )
621    {
622        if( !use_ipp )
623        {
624            delta = MIN( size.height - y, max_dy );
625            CvRect roi = cvRect(0,y,size.width,delta);
626            CvPoint origin=cvPoint(0,0);
627
628            if( y + delta == size.height )
629                stage = stage & CV_START ? CV_START + CV_END : CV_END;
630
631            dx_filter.process(src,Dx,roi,origin,stage);
632            dy_filter.process(src,Dy,roi,origin,stage);
633            d2x_filter.process(src,D2x,roi,origin,stage);
634            d2y_filter.process(src,D2y,roi,origin,stage);
635            stripe_size.height = dxy_filter.process(src,Dxy,roi,origin,stage);
636        }
637        else
638        {
639            delta = icvIPPFilterNextStripe( src, tempsrc, y, el_size, el_anchor );
640            stripe_size.height = delta;
641
642            IPPI_CALL( ipp_sobel_vert( shifted_ptr, temp_step,
643                Dx->data.ptr, d_step, stripe_size, aperture_size ));
644            IPPI_CALL( ipp_sobel_horiz( shifted_ptr, temp_step,
645                Dy->data.ptr, d_step, stripe_size, aperture_size ));
646            IPPI_CALL( ipp_sobel_vert_second( shifted_ptr, temp_step,
647                D2x->data.ptr, d_step, stripe_size, aperture_size ));
648            IPPI_CALL( ipp_sobel_horiz_second( shifted_ptr, temp_step,
649                D2y->data.ptr, d_step, stripe_size, aperture_size ));
650            IPPI_CALL( ipp_sobel_cross( shifted_ptr, temp_step,
651                Dxy->data.ptr, d_step, stripe_size, aperture_size ));
652        }
653
654        for( i = 0; i < stripe_size.height; i++, dst_y++ )
655        {
656            float* dstdata = (float*)(dst->data.ptr + dst_y*dst->step);
657
658            if( d_depth == CV_16S )
659            {
660                const short* dxdata = (const short*)(Dx->data.ptr + i*Dx->step);
661                const short* dydata = (const short*)(Dy->data.ptr + i*Dy->step);
662                const short* d2xdata = (const short*)(D2x->data.ptr + i*D2x->step);
663                const short* d2ydata = (const short*)(D2y->data.ptr + i*D2y->step);
664                const short* dxydata = (const short*)(Dxy->data.ptr + i*Dxy->step);
665
666                for( j = 0; j < stripe_size.width; j++ )
667                {
668                    double dx = dxdata[j];
669                    double dx2 = dx * dx;
670                    double dy = dydata[j];
671                    double dy2 = dy * dy;
672
673                    dstdata[j] = (float)(factor*(dx2*d2ydata[j] + dy2*d2xdata[j] - 2*dx*dy*dxydata[j]));
674                }
675            }
676            else
677            {
678                const float* dxdata = (const float*)(Dx->data.ptr + i*Dx->step);
679                const float* dydata = (const float*)(Dy->data.ptr + i*Dy->step);
680                const float* d2xdata = (const float*)(D2x->data.ptr + i*D2x->step);
681                const float* d2ydata = (const float*)(D2y->data.ptr + i*D2y->step);
682                const float* dxydata = (const float*)(Dxy->data.ptr + i*Dxy->step);
683
684                for( j = 0; j < stripe_size.width; j++ )
685                {
686                    double dx = dxdata[j];
687                    double dy = dydata[j];
688                    dstdata[j] = (float)(factor*(dx*dx*d2ydata[j] + dy*dy*d2xdata[j] - 2*dx*dy*dxydata[j]));
689                }
690            }
691        }
692
693        stage = CV_MIDDLE;
694    }
695
696    __END__;
697
698    cvReleaseMat( &Dx );
699    cvReleaseMat( &Dy );
700    cvReleaseMat( &D2x );
701    cvReleaseMat( &D2y );
702    cvReleaseMat( &Dxy );
703    cvReleaseMat( &tempsrc );
704}
705
706/* End of file */
707