zoukankan      html  css  js  c++  java
  • 使用索贝尔(Sobel)进行梯度运算时的数学意义和代码实现研究

    对于做图像处理的工程师来说,Sobel非常熟悉且常用。但是当我们需要使用Sobel进行梯度运算,且希望得到“数学结果”(作为下一步运算的基础)而不是“图片效果”的时候,就必须深入了解Sobel的知识原理和OpenCV实现的细节(当然我们是OpenCV支持则)。这里对具体内容进行研究。

    一、基本原理
    一般来说,用来表示微分的最常用的算子是索贝尔(Sobel)算子,它可以实现任意阶导数和混合偏导数(例如: ∂2/∂x∂y)。

    在x方向上用Sobel算子进行近似一阶求导的结果

    Sobel算子效果,y方向近似一阶导数

    OpenCV中给出了函数使用的定义

    void cv::Sobel(
      cv::InputArray  src,                 // 源图像
      cv::OutputArray dst,                 // 目标图像
      int             ddepth,              // 像素深度 (如CV_8U)
      int             xorder,              // x方向对应的倒数顺序
      int             yorder,              // y方向对应的倒数顺序
      cv::Size        ksize      = 3,      // 核大小
      double          scale      = 1,      // 阈值
      double          delta      = 0,      // 偏移
      int             borderType = cv::BORDER_DEFAULT  // 边框外推方法
    );

    1、其中src和dst是源图像和目标图像,可以通过指明参数ddepth来确定目标图像的深度或类型(如CV_32F)。举个例子,如果src是一幅8位图像,那么dst需要至少CV_16S的深度保证不出现溢出;

    2、xorder和yorder是求导顺序,其取值范围为0、1和2。0表示在这个方向上不进行求导,那2代表什么?

    3、ksize是一个奇数,表示了调用的滤波器的宽和高,目前最大支持到31;

    4、阈值和偏移将在把结果存入dst前调用,这有助于你将求导结果可视化.borderType参数的功能与其他卷积操作完全一样。

    上面有一个遗留问题,就是xorder(yorder)取2的时候代表什么?为此翻阅OpenCV源码

    step1

     

     

    step 2

    step3

     

    static void getSobelKernels( OutputArray _kx, OutputArray _ky,
                                 int dx, int dy, int _ksize, bool normalize, int ktype )
    {
        int i, j, ksizeX = _ksize, ksizeY = _ksize;
        if( ksizeX == 1 && dx > 0 )
            ksizeX = 3;
        if( ksizeY == 1 && dy > 0 )
            ksizeY = 3;
        CV_Assert( ktype == CV_32F || ktype == CV_64F );
        _kx.create(ksizeX, 1, ktype, -1true);
        _ky.create(ksizeY, 1, ktype, -1true);
        Mat kx = _kx.getMat();
        Mat ky = _ky.getMat();
        if( _ksize % 2 == 0 || _ksize > 31 )
            CV_Error( CV_StsOutOfRange, "The kernel size must be odd and not larger than 31" );
        std::vector<int> kerI(std::max(ksizeX, ksizeY) + 1);
        CV_Assert( dx >= 0 && dy >= 0 && dx+dy > 0 );
        forint k = 0; k < 2; k++ )
        {
            Mat* kernel = k == 0 ? &kx : &ky;
            int order = k == 0 ? dx : dy;
            int ksize = k == 0 ? ksizeX : ksizeY;
            CV_Assert( ksize > order );
            if( ksize == 1 )
                kerI[0= 1;
            else if( ksize == 3 )
            {
                if( order == 0 )
                    kerI[0= 1, kerI[1= 2, kerI[2= 1;
                else if( order == 1 )
                    kerI[0= -1, kerI[1= 0, kerI[2= 1;
                else
                    kerI[0= 1, kerI[1= -2, kerI[2= 1;
            }
            else
            {
                int oldval, newval;
                kerI[0= 1;
                for( i = 0; i < ksize; i++ )
                    kerI[i+1= 0;
                for( i = 0; i < ksize - order - 1; i++ )
                {
                    oldval = kerI[0];
                    for( j = 1; j <= ksize; j++ )
                    {
                        newval = kerI[j]+kerI[j-1];
                        kerI[j-1= oldval;
                        oldval = newval;
                    }
                }
                for( i = 0; i < order; i++ )
                {
                    oldval = -kerI[0];
                    for( j = 1; j <= ksize; j++ )
                    {
                        newval = kerI[j-1- kerI[j];
                        kerI[j-1= oldval;
                        oldval = newval;
                    }
                }
            }
            Mat temp(kernel->rows, kernel->cols, CV_32S, &kerI[0]);
            double scale = !normalize ? 1: 1./(1 << (ksize-order-1));
            temp.convertTo(*kernel, ktype, scale);
        }
    }
    其中
     

    那么,可以看见,当order ==2 时候,生成了[1 -2 1]作为类似的模板,不管是什么,这个不是我想要的。

         除了上面看到的,还可以发现同时设置xorder 和 yorder的时候,最后并没有看到相加的动作。而如果我们计算的结果是梯度场的时候,就不仅要算xorder,而且要算yorder,并且最后要把这两个结果求和。如果自己编码,那么可能如下:

    //求x方向偏导
       vx = *(lpSrc + IMGW + 1- *(lpSrc + IMGW - 1+
       *(lpSrc + 1)*2 - *(lpSrc - 1)*2 +
       *(lpSrc - IMGW + 1- *(lpSrc - IMGW - 1);
       //求y方向偏导
       vy = *(lpSrc + IMGW - 1- *(lpSrc - IMGW - 1+
       *(lpSrc + IMGW)*2 - *(lpSrc - IMGW)*2 +
       *(lpSrc + IMGW + 1- *(lpSrc - IMGW + 1);
       gradSum += (abs(vx)+abs(vy));        
     
    二、定量研究
        如果直接用自然图片(比如lena),sobel计算之后可能啥都看不出来。为了定量研究,必须自己做图片。
        搞成这样,看的清楚
    //自己生成实验图片
        Mat matTst = Mat(Size(11,11),CV_8UC1,Scalar(0));
        line(matTst,Point(5,0),Point(5,11),Scalar(255));
        line(matTst,Point(0,5),Point(11,5),Scalar(255));
    求一遍X方向的导数
    //x方向求导
        Sobel(matTst,matX,CV_16SC1,1,0);
    为什么要用16s?首先,所有的模式为CV_{8U,16S,16U,32S,32F,64F}C{1,2,3},其中U代表char,S代表short,F代表float,这个和c++里面一样;8只能是正数,16/32都可以是负数。所有的C1和C都是一样的,比如cv_8uc1 == cv_8u。那么原始图像是CV_8UC1,也就是CV_8U了,那么进行卷积计算,也就是原图像和
    -1 0 1
    -2 0 2
    -1 0 1
    或者类似的东西进行卷积计算,那么结果得到最大为 1*255+2*255+1*255 - 0  > 10000,最小的结果为 0 - 1*255+2*255+1*255  < -10000,所以要用CV_16SC1来保存,才合适。那么得到的结果为
    ​首先是有正有负,集中在中间区域值变化非常的大。取其中一个像素来看
    中心位置 使用
    -1 0 1
    -2 0 2
    -1 0 1
    进行卷积,那么结果为
    0*-1 + 255*0 + 0*1 +255*-2 +255*0 + 255*2 + 0 *-1 +255*0+0*1 = 0+0+0-255+0+255*2+0+0+0 = 0
    其他结果也是,那么结果是符合卷积运算的。
    计算Y方向卷积
    //y方向求导
        Sobel(matTst,matY,CV_16SC1,0,1);
      
    计算XY的卷积
      //xy方向
        Sobel(matTst,matXY,CV_16SC1,1,1);
    这个结果不是我想要的,我想要的是 (abs(vx)+abs(vy));    
    那么这样实现
    //求和
        matXY = abs(matX) + abs(matY);
    甚至可以进一步帮助显示
    //方便显示
     normalize(matXY,matXY,0,255,NORM_MINMAX);
        matXY.convertTo(matXY,CV_8UC1);
    无论如何,这个结果和直接Sobel(matTst,matXY,CV_16SC1,1,1);都是不一样的
    三、小结反思
    应该说,Sobel大概能够做什么?这个很早之前就已经知道了。但是为什么能够达到这样的效果?这个问题,一直到我需要使用Sobel进行相关的数学计算的时候才能够搞明白。
    掌握知识最后要归结到数学抽象层次,才能够算是彻底掌握。这是Sobel之外的获得。
    感谢阅读至此,希望共同进步!





  • 相关阅读:
    /etc/sysctl.conf 控制内核相关配置文件
    python 并发编程 非阻塞IO模型
    python 并发编程 多路复用IO模型
    python 并发编程 异步IO模型
    python 并发编程 阻塞IO模型
    python 并发编程 基于gevent模块 协程池 实现并发的套接字通信
    python 并发编程 基于gevent模块实现并发的套接字通信
    python 并发编程 io模型 目录
    python 并发编程 socket 服务端 客户端 阻塞io行为
    python 并发编程 IO模型介绍
  • 原文地址:https://www.cnblogs.com/jsxyhelu/p/7435680.html
Copyright © 2011-2022 走看看