rambo

matlab和opencv计算图像梯度

2015年6月10日 Machine Vision 0 Comments

1、Matlab gradient函数

Gradient(F)函数求的是数值上的梯度，假设F为矩阵.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
>> x=[6,9,3,4,0;5,4,1,2,5;6,7,7,8,0;7,8,9,10,0]
x =
   6     9     3     4     0
   5     4     1     2     5
   6     7     7     8     0
   7     8     9    10     0
 
>> [Fx,Fy]=gradient(x)
Fx =
  3.0000   -1.5000   -2.5000   -1.5000   -4.0000
 -1.0000   -2.0000   -1.0000    2.0000    3.0000
  1.0000    0.5000    0.5000   -3.5000   -8.0000
  1.0000    1.0000    1.0000   -4.5000  -10.0000
Fy =
 -1.0000   -5.0000   -2.0000   -2.0000    5.0000
       0   -1.0000    2.0000    2.0000         0
  1.0000    2.0000    4.0000    4.0000   -2.5000
  1.0000    1.0000    2.0000    2.0000         0

计算规则: [Fx,Fy]=gradient(F)，其中Fx为其水平方向上的梯度，Fy为其垂直方向上的梯度，Fx的第一列元素为原矩阵第二列与第一列元素之差，Fx的第二列元素为原矩阵第三列与第一列元素之差除以2，以此类推：Fx(i,j)=(F(i,j+1)-F(i,j-1))/2。最后一列则为最后两列之差。同理，可以得到Fy。
  • 如果F是一维矩阵，则FX=gradient(F，H)返回F的一维数值梯度。H是F中相邻两点间的间距。
  • 如果F是二维矩阵，返回F的二维数值梯度。  [FX，FY]=gradient(F,HX,HY)。HX，HY参数表示各方向相邻两点的距离，默认参数分别为1，1。 [z1,z2]=gradient(x,2,2);
  • 如果F是三维矩阵，返回F的三维数值梯度。[FX，FY，FZ]=gradient(F,HX,HY,HZ)。HX，HY，HZ参数表示各方向相邻两点的距离。

注意：Sobel算子计算的是像素值的二阶微分 (f(i+1,j) + f(i-1,j) - 2f(i,j)) / 2。

2、Matlab conv2函数及opencv实现

matlab conv2函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
I = 'fsrc.jpg';
Img = imread(I);
if size(Img,3) == 3
    G = rgb2gray(Img);
else
    G = Img;
end
 
hf1 = [-1,0,1];
vf1 = [-1,0,1]';
 
G=single(G);
im1 = conv2(G, hf1, 'same');
im11=uint8(im1);//注意，这里需要把single类型转换为uint8类型，这样和opencv得到的函数是一致的
im2= conv2(G, vf1, 'same');

1）前一节我们学习了图像的卷积操作。一个很自然的问题是如何处理卷积边缘。当卷积点在图像边界时会发生什么，如何处理这个问题？

2）大多数用到卷积操作的OpenCV函数都是将给定图像拷贝到另一个轻微变大的图像中，然后自动填充图像边界(通过下面示例代码中的各种方式)。这样卷积操作就可以在边界像素安全执行了(填充边界在操作完成后会自动删除)。

3）图像边界的两种方法:

BORDER_CONSTANT: 使用常数填充边界 (i.e. 黑色或者 0)

BORDER_REPLICATE: 复制原图中最临近的行或者列。

 opencv实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
enum ConvolutionType {  
/* Return the full convolution, including border */
CONVOLUTION_FULL,
 
/* Return only the part that corresponds to the original image */
CONVOLUTION_SAME,
 
/* Return only the submatrix containing elements that were not influenced by the border */
CONVOLUTION_VALID
};
 
void HJconv2(const Mat &img, const Mat& kernel, ConvolutionType type, Mat& dest) {
Mat source = img;
if(CONVOLUTION_FULL == type) {
source = Mat();
const int additionalRows = kernel.rows-1, additionalCols = kernel.cols-1;
copyMakeBorder(img, source, (additionalRows+1)/2, additionalRows/2, (additionalCols+1)/2, additionalCols/2, BORDER_CONSTANT, Scalar(0));
}
 
Point anchor(kernel.cols - kernel.cols/2 - 1, kernel.rows - kernel.rows/2 - 1);
int borderMode = BORDER_CONSTANT;
cv::Mat result;
flip(kernel,result,-1);
//cout<<kernel<<endl;
//cout<<result<<endl;
filter2D(source, dest, img.depth(),result , anchor, 0, borderMode);
 
if(CONVOLUTION_VALID == type) {
dest = dest.colRange((kernel.cols-1)/2, dest.cols - kernel.cols/2)
.rowRange((kernel.rows-1)/2, dest.rows - kernel.rows/2);
}
}

1
2
3
4
5
6
7
8
9
10
11
12
13
//unit test
int _tmain(int argc, _TCHAR* argv[]){
Mat src,fsrc;
src = imread("E:\\fsrc.jpg", CV_LOAD_IMAGE_COLOR); // 这里为矩阵开辟内存
if(src.depth()==0)
cvtColor(src,src,CV_BGR2GRAY);
 
//double a[] = {-1,0,1};
//Mat kernel = Mat(1, 3,CV_8SC1 , a);
Mat kernel = ( Mat_< float>(1,3) << -1, 0 , 1 );
Mat dest;
HJconv2(src,kernel,CONVOLUTION_SAME,dest);
}

输出数据到文件中，来验证是否处理正确

1
2
3
4
5
6
7
8
9
10
11
12
13
14
ofstream angFileimage("D:\\Ang.txt", ios_base::out);  //按新建或覆盖方式写入  
if (! angFileimage.is_open()) {  
cout << "打开文件失败" << endl;  
return -1;
}  
 
for (int r = 0; r < Ang.rows; r++){  
for (int c = 0; c < Ang.cols; c++){  
float data = Ang.at<float>(r,c);  
                        //读取数据，at<type> - type 是矩阵元素的具体数据格式  
angFileimage<< data << "\t" ;   //每列数据用 tab 隔开  
}  
angFileimage << endl;  //换行  
}

参考资料：
1、http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/imgproc/imgtrans/copyMakeBorder/copyMakeBorder.html
2、http://blog.timmlinder.com/2011/07/opencv-equivalent-to-matlabs-conv2-function/
3、http://blog.163.com/yuyang_tech/blog/static/216050083201352594630749/