spatialgradient.cpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                           License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
00014 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
00015 // Third party copyrights are property of their respective owners.
00016 //
00017 // Redistribution and use in source and binary forms, with or without modification,
00018 // are permitted provided that the following conditions are met:
00019 //
00020 //   * Redistribution's of source code must retain the above copyright notice,
00021 //     this list of conditions and the following disclaimer.
00022 //
00023 //   * Redistribution's in binary form must reproduce the above copyright notice,
00024 //     this list of conditions and the following disclaimer in the documentation
00025 //     and/or other materials provided with the distribution.
00026 //
00027 //   * The name of the copyright holders may not be used to endorse or promote products
00028 //     derived from this software without specific prior written permission.
00029 //
00030 // This software is provided by the copyright holders and contributors "as is" and
00031 // any express or implied warranties, including, but not limited to, the implied
00032 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00033 // In no event shall the Intel Corporation or contributors be liable for any direct,
00034 // indirect, incidental, special, exemplary, or consequential damages
00035 // (including, but not limited to, procurement of substitute goods or services;
00036 // loss of use, data, or profits; or business interruption) however caused
00037 // and on any theory of liability, whether in contract, strict liability,
00038 // or tort (including negligence or otherwise) arising in any way out of
00039 // the use of this software, even if advised of the possibility of such damage.
00040 //
00041 //M*/
00042 
00043 #include "precomp.hpp"
00044 #include "opencv2/core/hal/intrin.hpp"
00045 
00046 #include <iostream>
00047 namespace cv
00048 {
00049 
00050 /* NOTE:
00051  *
00052  * Sobel-x: -1  0  1
00053  *          -2  0  2
00054  *          -1  0  1
00055  *
00056  * Sobel-y: -1 -2 -1
00057  *           0  0  0
00058  *           1  2  1
00059  */
00060 template <typename T>
00061 static inline void spatialGradientKernel( T& vx, T& vy,
00062                                           const T& v00, const T& v01, const T& v02,
00063                                           const T& v10,               const T& v12,
00064                                           const T& v20, const T& v21, const T& v22 )
00065 {
00066     // vx = (v22 - v00) + (v02 - v20) + 2 * (v12 - v10)
00067     // vy = (v22 - v00) + (v20 - v02) + 2 * (v21 - v01)
00068 
00069     T tmp_add = v22 - v00,
00070       tmp_sub = v02 - v20,
00071       tmp_x   = v12 - v10,
00072       tmp_y   = v21 - v01;
00073 
00074     vx = tmp_add + tmp_sub + tmp_x + tmp_x;
00075     vy = tmp_add - tmp_sub + tmp_y + tmp_y;
00076 }
00077 
00078 void spatialGradient( InputArray _src, OutputArray _dx, OutputArray _dy,
00079                       int ksize, int borderType )
00080 {
00081 
00082     // Prepare InputArray src
00083     Mat src = _src.getMat();
00084     CV_Assert( !src.empty() );
00085     CV_Assert( src.type() == CV_8UC1 );
00086     CV_Assert( borderType == BORDER_DEFAULT || borderType == BORDER_REPLICATE );
00087 
00088     // Prepare OutputArrays dx, dy
00089     _dx.create( src.size(), CV_16SC1 );
00090     _dy.create( src.size(), CV_16SC1 );
00091     Mat dx = _dx.getMat(),
00092         dy = _dy.getMat();
00093 
00094     // TODO: Allow for other kernel sizes
00095     CV_Assert(ksize == 3);
00096 
00097     // Get dimensions
00098     const int H = src.rows,
00099               W = src.cols;
00100 
00101     // Row, column indices
00102     int i = 0,
00103         j = 0;
00104 
00105     // Handle border types
00106     int i_top    = 0,     // Case for H == 1 && W == 1 && BORDER_REPLICATE
00107         i_bottom = H - 1,
00108         j_offl   = 0,     // j offset from 0th   pixel to reach -1st pixel
00109         j_offr   = 0;     // j offset from W-1th pixel to reach Wth  pixel
00110 
00111     if ( borderType == BORDER_DEFAULT ) // Equiv. to BORDER_REFLECT_101
00112     {
00113         if ( H > 1 )
00114         {
00115             i_top    = 1;
00116             i_bottom = H - 2;
00117         }
00118         if ( W > 1 )
00119         {
00120             j_offl = 1;
00121             j_offr = -1;
00122         }
00123     }
00124 
00125     // Pointer to row vectors
00126     uchar *p_src, *c_src, *n_src; // previous, current, next row
00127     short *c_dx,  *c_dy;
00128 
00129     int i_start = 0;
00130     int j_start = 0;
00131 #if CV_SIMD128 && CV_SSE2
00132     uchar *m_src;
00133     short *n_dx, *n_dy;
00134 
00135     // Characters in variable names have the following meanings:
00136     // u: unsigned char
00137     // s: signed int
00138     //
00139     // [row][column]
00140     // m: offset -1
00141     // n: offset  0
00142     // p: offset  1
00143     // Example: umn is offset -1 in row and offset 0 in column
00144     for ( i = 0; i < H - 1; i += 2 )
00145     {
00146         if   ( i == 0 ) p_src = src.ptr<uchar>(i_top);
00147         else            p_src = src.ptr<uchar>(i-1);
00148 
00149         c_src = src.ptr<uchar>(i);
00150         n_src = src.ptr<uchar>(i+1);
00151 
00152         if ( i == H - 2 ) m_src = src.ptr<uchar>(i_bottom);
00153         else              m_src = src.ptr<uchar>(i+2);
00154 
00155         c_dx = dx.ptr<short>(i);
00156         c_dy = dy.ptr<short>(i);
00157         n_dx = dx.ptr<short>(i+1);
00158         n_dy = dy.ptr<short>(i+1);
00159 
00160         v_uint8x16 v_select_m = v_uint8x16(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
00161                                            0, 0, 0, 0xFF);
00162 
00163         // Process rest of columns 16-column chunks at a time
00164         for ( j = 1; j < W - 16; j += 16 )
00165         {
00166             // Load top row for 3x3 Sobel filter
00167             v_uint8x16 v_um = v_load(&p_src[j-1]);
00168             v_uint8x16 v_up = v_load(&p_src[j+1]);
00169             // TODO: Replace _mm_slli_si128 with hal method
00170             v_uint8x16 v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
00171                                                    v_uint8x16(_mm_srli_si128(v_um.val, 1)));
00172             v_uint16x8 v_um1, v_um2, v_un1, v_un2, v_up1, v_up2;
00173             v_expand(v_um, v_um1, v_um2);
00174             v_expand(v_un, v_un1, v_un2);
00175             v_expand(v_up, v_up1, v_up2);
00176             v_int16x8 v_s1m1 = v_reinterpret_as_s16(v_um1);
00177             v_int16x8 v_s1m2 = v_reinterpret_as_s16(v_um2);
00178             v_int16x8 v_s1n1 = v_reinterpret_as_s16(v_un1);
00179             v_int16x8 v_s1n2 = v_reinterpret_as_s16(v_un2);
00180             v_int16x8 v_s1p1 = v_reinterpret_as_s16(v_up1);
00181             v_int16x8 v_s1p2 = v_reinterpret_as_s16(v_up2);
00182 
00183             // Load second row for 3x3 Sobel filter
00184             v_um = v_load(&c_src[j-1]);
00185             v_up = v_load(&c_src[j+1]);
00186             // TODO: Replace _mm_slli_si128 with hal method
00187             v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
00188                                         v_uint8x16(_mm_srli_si128(v_um.val, 1)));
00189             v_expand(v_um, v_um1, v_um2);
00190             v_expand(v_un, v_un1, v_un2);
00191             v_expand(v_up, v_up1, v_up2);
00192             v_int16x8 v_s2m1 = v_reinterpret_as_s16(v_um1);
00193             v_int16x8 v_s2m2 = v_reinterpret_as_s16(v_um2);
00194             v_int16x8 v_s2n1 = v_reinterpret_as_s16(v_un1);
00195             v_int16x8 v_s2n2 = v_reinterpret_as_s16(v_un2);
00196             v_int16x8 v_s2p1 = v_reinterpret_as_s16(v_up1);
00197             v_int16x8 v_s2p2 = v_reinterpret_as_s16(v_up2);
00198 
00199             // Load third row for 3x3 Sobel filter
00200             v_um = v_load(&n_src[j-1]);
00201             v_up = v_load(&n_src[j+1]);
00202             // TODO: Replace _mm_slli_si128 with hal method
00203             v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
00204                                         v_uint8x16(_mm_srli_si128(v_um.val, 1)));
00205             v_expand(v_um, v_um1, v_um2);
00206             v_expand(v_un, v_un1, v_un2);
00207             v_expand(v_up, v_up1, v_up2);
00208             v_int16x8 v_s3m1 = v_reinterpret_as_s16(v_um1);
00209             v_int16x8 v_s3m2 = v_reinterpret_as_s16(v_um2);
00210             v_int16x8 v_s3n1 = v_reinterpret_as_s16(v_un1);
00211             v_int16x8 v_s3n2 = v_reinterpret_as_s16(v_un2);
00212             v_int16x8 v_s3p1 = v_reinterpret_as_s16(v_up1);
00213             v_int16x8 v_s3p2 = v_reinterpret_as_s16(v_up2);
00214 
00215             // dx & dy for rows 1, 2, 3
00216             v_int16x8 v_sdx1, v_sdy1;
00217             spatialGradientKernel<v_int16x8>( v_sdx1, v_sdy1,
00218                                               v_s1m1, v_s1n1, v_s1p1,
00219                                               v_s2m1,         v_s2p1,
00220                                               v_s3m1, v_s3n1, v_s3p1 );
00221 
00222             v_int16x8 v_sdx2, v_sdy2;
00223             spatialGradientKernel<v_int16x8>( v_sdx2, v_sdy2,
00224                                               v_s1m2, v_s1n2, v_s1p2,
00225                                               v_s2m2,         v_s2p2,
00226                                               v_s3m2, v_s3n2, v_s3p2 );
00227 
00228             // Store
00229             v_store(&c_dx[j],   v_sdx1);
00230             v_store(&c_dx[j+8], v_sdx2);
00231             v_store(&c_dy[j],   v_sdy1);
00232             v_store(&c_dy[j+8], v_sdy2);
00233 
00234             // Load fourth row for 3x3 Sobel filter
00235             v_um = v_load(&m_src[j-1]);
00236             v_up = v_load(&m_src[j+1]);
00237             // TODO: Replace _mm_slli_si128 with hal method
00238             v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
00239                                         v_uint8x16(_mm_srli_si128(v_um.val, 1)));
00240             v_expand(v_um, v_um1, v_um2);
00241             v_expand(v_un, v_un1, v_un2);
00242             v_expand(v_up, v_up1, v_up2);
00243             v_int16x8 v_s4m1 = v_reinterpret_as_s16(v_um1);
00244             v_int16x8 v_s4m2 = v_reinterpret_as_s16(v_um2);
00245             v_int16x8 v_s4n1 = v_reinterpret_as_s16(v_un1);
00246             v_int16x8 v_s4n2 = v_reinterpret_as_s16(v_un2);
00247             v_int16x8 v_s4p1 = v_reinterpret_as_s16(v_up1);
00248             v_int16x8 v_s4p2 = v_reinterpret_as_s16(v_up2);
00249 
00250             // dx & dy for rows 2, 3, 4
00251             spatialGradientKernel<v_int16x8>( v_sdx1, v_sdy1,
00252                                               v_s2m1, v_s2n1, v_s2p1,
00253                                               v_s3m1,         v_s3p1,
00254                                               v_s4m1, v_s4n1, v_s4p1 );
00255 
00256             spatialGradientKernel<v_int16x8>( v_sdx2, v_sdy2,
00257                                               v_s2m2, v_s2n2, v_s2p2,
00258                                               v_s3m2,         v_s3p2,
00259                                               v_s4m2, v_s4n2, v_s4p2 );
00260 
00261             // Store
00262             v_store(&n_dx[j],   v_sdx1);
00263             v_store(&n_dx[j+8], v_sdx2);
00264             v_store(&n_dy[j],   v_sdy1);
00265             v_store(&n_dy[j+8], v_sdy2);
00266         }
00267     }
00268     i_start = i;
00269     j_start = j;
00270 #endif
00271     int j_p, j_n;
00272     uchar v00, v01, v02, v10, v11, v12, v20, v21, v22;
00273     for ( i = 0; i < H; i++ )
00274     {
00275         if   ( i == 0 ) p_src = src.ptr<uchar>(i_top);
00276         else            p_src = src.ptr<uchar>(i-1);
00277 
00278         c_src = src.ptr<uchar>(i);
00279 
00280         if ( i == H - 1 ) n_src = src.ptr<uchar>(i_bottom);
00281         else              n_src = src.ptr<uchar>(i+1);
00282 
00283         c_dx = dx.ptr<short>(i);
00284         c_dy = dy.ptr<short>(i);
00285 
00286         // Process left-most column
00287         j = 0;
00288         j_p = j + j_offl;
00289         j_n = 1;
00290         if ( j_n >= W ) j_n = j + j_offr;
00291         v00 = p_src[j_p]; v01 = p_src[j]; v02 = p_src[j_n];
00292         v10 = c_src[j_p]; v11 = c_src[j]; v12 = c_src[j_n];
00293         v20 = n_src[j_p]; v21 = n_src[j]; v22 = n_src[j_n];
00294         spatialGradientKernel<short>( c_dx[0], c_dy[0], v00, v01, v02, v10,
00295                                       v12, v20, v21, v22 );
00296         v00 = v01; v10 = v11; v20 = v21;
00297         v01 = v02; v11 = v12; v21 = v22;
00298 
00299         // Process middle columns
00300         j = i >= i_start ? 1 : j_start;
00301         j_p = j - 1;
00302         v00 = p_src[j_p]; v01 = p_src[j];
00303         v10 = c_src[j_p]; v11 = c_src[j];
00304         v20 = n_src[j_p]; v21 = n_src[j];
00305 
00306         for ( ; j < W - 1; j++ )
00307         {
00308             // Get values for next column
00309             j_n = j + 1; v02 = p_src[j_n]; v12 = c_src[j_n]; v22 = n_src[j_n];
00310             spatialGradientKernel<short>( c_dx[j], c_dy[j], v00, v01, v02, v10,
00311                                           v12, v20, v21, v22 );
00312 
00313             // Move values back one column for next iteration
00314             v00 = v01; v10 = v11; v20 = v21;
00315             v01 = v02; v11 = v12; v21 = v22;
00316         }
00317 
00318         // Process right-most column
00319         if ( j < W )
00320         {
00321             j_n = j + j_offr; v02 = p_src[j_n]; v12 = c_src[j_n]; v22 = n_src[j_n];
00322             spatialGradientKernel<short>( c_dx[j], c_dy[j], v00, v01, v02, v10,
00323                                           v12, v20, v21, v22 );
00324         }
00325     }
00326 
00327 }
00328 
00329 }
00330