2013-02-20 16:36:25 ace12304568 阅读数 3269

头疼了两天总算整明白了:

前言:

       在ocr(数字图像文本识别)过程中,由于图像的不可控性,总会存在一定角度的倾斜。倾斜角度要满足一定范围:   

       θ≤d/L         (θ即倾斜角度,d为文本行距,L为文本行长)

      若超出这个范围,则可能将下(上)一行的文字拼接到本行,替代原有文字,产生误断。

      通常来讲:θ应在2度以内。

      为了避免或者是纠正这种问题,就需要对图像进行预处理,进行倾斜检测:

检测方法

     有如下假设:搜索图像都已进行了二值化处理(非黑即白)(常用方法:中值阈值,最佳阈值,otu算法)。图像大小为M*N,图像只是倾斜,并未上下翻转。

     直线在想素质上的特点:扫描图像中一条直线方程上的坐标点,会发现黑点最多(N<Max(M,N))。

     方法1:投影法。

      思想是:如图所示,扫描有限条直线(m条),统计每条直线上法线方向黑点最多(最少)的一斜行。在所有最值中,最大的那行所对应的斜率mx/n的绝对值即倾斜角度的正切值。

  

       实现方法:输入int[][] pixels;

                        输出int[][]  A;

                        步骤: 1初始化参数int a=n, b=0, count=0, max=0,c=0;

                                     直线方程:y=ax/b;法线方程y=bx/a+c

                                 2     for(;b<m;b++){

                                             max=0;

                                              for(int i=0;i<n;i++){

                                                  //统计对应斜率(a,b)时,过点(i*b/a,i)的法线上的黑点数;

                                                     count=0;  c=(a*a-b*b)*i/a*a 

                                                     for(int j=0;j<m;i++)

                                                        if(pixels[j][b/a*j+c]==1)count++;

                                                    if(count>max){

                                                        max=count;

                                                        A[b][0]=i;//记录每种斜率上黑点最多的斜行,和最大值;

                                                        A[b][1]=max;

                                                     }

                                                }

                                           }

                                    3接下来是对数组A的分析了。这个数组保存了图像各个方向上最长直线的斜率和位置。                                

                               选择合理的其中一个作为倾斜角度,即可。也可简单的只去最大值;因为文本图像通常没有很

                               多的噪声,不会出现大的色块。其实依照数组你甚至可以重绘所有被检测到的直线。



       方法2:hough变换:

   原理:

            图像中任意一条直线的方程可写为Ax+By=C;做归一化处理可得方程Ax+By=1;那么一条直线将唯一的被数对(A,B)确定。

            过任意一点(a,b)的直线方程为Ma+Nb=1;若p1(a,b),p2(c,d)共线(不是几何意义的共线,而是有条线,两点均在此线上)Ax+By=1可得:M1a+N1b=1与M2c+N2d=1交于(A,B)点.

       如图所示:

那么在扫描整幅图像时:可作如下处理:

        0.初始条件:int[][] AB;

        1.i=0->m;j=0->n

                 if(pixels[i][j]==1)

                         a=0->max;b=0->max

                                 if(ai+bj==1)

                                      AB[i][j]++;

         2.寻找AB中的最大值,其对应坐标(i,j);就是所求倾斜角的-i/j就是倾斜角的正切值。

       

     实现:

           上述原理中有一处不好处理就是AB的取值范围,不好确定。

           做如下改进即可,原图像中直线的方程设为:

                  p=xcosa+ysina;(a为图像中直线倾斜角度的余角,p为原点到直线的距离)

           则A->p/sina  b->p/cosa  0<a<180  0<p<sqrt(m*m+n*n);

           即做了一个十字坐标系xoy到极坐标系poa的转换。



PS:由此可以看到直线检测在图像处理中是如何运作的。是跟曲线方程特征而决定的。同理其他曲线,也可根据曲线方程的特点有hough变换得到。投影法就略显困难。因此hough变换被广泛用于图像中图形识别过程中。

但是hough变换也存在缺点:计算量大,优点是图像内容无关。





2018-01-23 17:47:11 LQMIKU 阅读数 1736

文档图像倾斜角检测及校正(二)

996.icu LICENSE

  • 最小距离法直线拟合原理
  • Matlab程序

阅读之前注意:

Hi,你好,我是Cooper Liu,欢迎来到我写的“文档图像校正”系列博客。基于三种原理,我写了四个实验性的Matlab验证程序,以及两个文档校正Matlab程序。在这里你将能够获取所有的源代码以及测试图片,完全可以在你自己的Matlab上跑这些程序。

如果你是学生,请注意不要抄袭,课程设计作业的话,这种程序仅仅只能让你得到80%左右的成绩。
如果你是工作人士或者只是感兴趣的极客,Okay,我想这些程序对于理解原理是如何应用为程序的已经足够 。

最后,请勿将这些资源用于商业用途(如你所见,这些程序都非常的初级)或者是谋取个人利益,知识在传播的过程中能展现更大的价值-

本文阅读建议用时:32min
本文阅读结构如下表:

项目 下属项目 测试用例数量
最小距离法直线拟合原理 0
Matlab程序 1

最小距离法直线拟合原理1

最小距离法拟合直线是指,我们设线条上n个特征点:
(xi,yi),i=1,2,,n(x_i,y_i ), i=1,2,…,n,则第i个特征点到直线 y=kx+by=kx+b 的垂直距离为:
yi(kxi+b)(1+k2)\frac{|y_i-(kx_i+b)|}{\sqrt{(1+k^2 )}}因此n个特征点到直线 y=kx+by=kx+b 的距离平方和为:
D(k,b)=i=1n[yi(kxi+b)]21+k2=11+k2i=1n[yi(kxi+b)]2D(k,b)=\sum_{i=1}^n\frac{[y_i-(kx_i+b)]^2}{1+k^2}= \frac{1}{1+k^2}\sum_{i=1}^n[y_i-(kx_i+b)]^2如果存在kd,bdk_d,b_d使得:
D(kd,bd)=mink,bD(k,b)D(k_d,b_d )=\min_{k,b}D(k,b)
则称kd,bdk_d,b_dk,bk,b的最小距离估计,此方法称为最小距离法。
为了计算和叙述的方便,我们引入以下的记号:
xˉ=1ni=1nxiyˉ=1ni=1nyiLxx=i=1n(xixˉ)2Lyy=i=1n(yiyˉ)2\bar x=\frac{1}{n}\sum_{i=1}^n x_i \quad \bar y=\frac{1}{n}\sum_{i=1}^n y_i\quad L_{xx}=\sum_{i=1}^n (x_i-\bar x)^2\quad L_{yy}=\sum_{i=1}^n (y_i-\bar y)^2 Lxy=i=1n(xixˉ)(yiyˉ) L_{xy}=\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)\quad 利用极值法求得:bd=ykxˉkd=LyyLxx+(LyyLxx)2+4Lxy22Lxyb_d=y-k\bar x,k_d=\frac{L_{yy}-L_{xx}+\sqrt{(L_{yy}-L_{xx} )^2+4L_{xy}^2 }}{2L_{xy}}
所以文档的倾斜角为:αd=arctan(kd)α_d=arctan⁡(k_d)

Matlab程序

基于以上原理,我编写了Matlab程序来进行实验,源代码可以参考angleDetection2.m文件。测试图片为line5.bmp。2
测试图片可以从这里获取,链接:https://pan.baidu.com/s/1dGmmGjn 密码:okt3

以下是实验结果:
实验结果
时间

如果您不想打开新的页面查看matlab源代码,也可以直接参考以下代码:

%%本版基于最小距离法直线拟合原理 
%%2018.01.19 by Cooper Liu
%%Questions? Contact me: angelpoint@foxmail.com
clear;clc; %清空之前的变量
I=imread('line5.bmp'); %读取图像
level=graythresh(I); %使用最大类间方差法找到图片的一个合适的阈值
bw=im2bw(I,level); %根据阈值,使用im2bw函数将灰度图像转换为二值图像时
figure(1);imshow(bw);

[m,n]=size(bw); %获取尺寸
xSum=0;xCount=0;
ySum=0;yCount=0;
tic; %计时开始
for i=1:m
    for j=1:n
        if bw(i,j)==0
           xSum=xSum+i;xCount=xCount+1;
           ySum=ySum+j;yCount=yCount+1;
        end
    end
end
xMean=xSum/xCount;
yMean=ySum/yCount;

Lxx=0;
Lyy=0;
Lxy=0;
for i=1:m
    for j=1:n
        if bw(i,j)==0
          Lxx=Lxx+(i-xMean)^2;
          Lyy=Lyy+(j-yMean)^2;
          Lxy=Lxy+(i-xMean)*(j-yMean);
        end
    end
end
toc %获取计时时间
tmp=((Lyy-Lxx) + nthroot((Lyy-Lxx)^2 + 4*Lxy^2, 2))/(2*Lxy);
if isnan(tmp) 
    tmp=inf; %如果求得tmp是NaN,那么不需要旋转
end
%tmp=atan(2*Lxy/(Lxx-Lyy))/2; %如果倾斜角相对于水平超过正负45
angle=atan(tmp); %用atan求出来的角度在-pi/2到+pi/2之间
angle=angle*180/pi

rot=90-angle;
pic=imrotate(I,rot,'crop'); %旋转图像
figure(2);imshow(pic);

如果本文有帮助到你,不如请我一罐可乐吧 🐬
在这里插入图片描述


  1. 参考文献:[1] 荆雷,张欣,郭金鑫.基于版面的拍照文档图像倾斜校正.激光与红外[J].2010,第10期 ↩︎

  2. 同时感谢愿意在网络上分享自己想法的各位博主。 ↩︎

2018-01-23 14:40:57 LQMIKU 阅读数 3853

文档图像倾斜角检测及校正(一)

996.icu LICENSE

  • 霍夫变换原理
  • Matlab程序

阅读之前注意:

Hi,你好,我是Cooper Liu,欢迎来到我写的“文档图像校正”系列博客。基于三种原理,我写了四个实验性的Matlab验证程序,以及两个文档校正Matlab程序。在这里你将能够获取所有的源代码以及测试图片,完全可以在你自己的Matlab上跑这些程序。

如果你是学生,请注意不要抄袭,课程设计作业的话,这种程序仅仅只能让你得到80%左右的成绩。
如果你是工作人士或者只是感兴趣的极客,Okay,我想这些程序对于理解原理是如何应用为程序的已经足够 。

最后,请勿将这些资源用于商业用途(如你所见,这些程序都非常的初级)或者是谋取个人利益,知识在传播的过程中能展现更大的价值-

本文阅读建议用时:32min
本文阅读结构如下表:

项目 下属项目 测试用例数量
霍夫变换原理 0
Matlab程序 1

霍夫变换原理

所谓霍夫变换,即对于图像平面上的一个点(x , y ),我们采用参数方程p=xcos(θ)+ysin(θ)把这个点映射到参数p-theta平面,那么图像平面上的一个点就对应p-theta平面的一条曲线,其中的p表示图像平面中的这个点所在直线到原点的距离,theta表示这个点所在直线与X轴的夹角。

又因为图像平面上的一个点对应一系列穿过这个点的直线,即有一系列对应的p和theta,所以一个点在参数p-theta平面对应着一条正弦曲线。由此我们可以推导,如果是图像平面上的一条直线,那么直线上的每个点在参数p-theta平面对应的曲线都会相交于同一点,即当前直线的(theta, p)。

基于霍夫变换原理,我们可以在p-theta平面找到最多曲线相交的那点,这一点对应着图像平面最长的直线(可以是连续的也可以是不连续的)。这个点的theta坐标即是我们要寻找的倾斜角。1

Matlab程序

基于以上原理,我编写了Matlab程序来进行实验,源代码可以参考angleDetection1.m文件。测试图片为line5.bmp。2
测试图片可以从这里获取,链接:https://pan.baidu.com/s/1dGmmGjn 密码:okt3

以下是实验结果:
实验结果
时间

如果您不想打开新的页面查看matlab源代码,也可以直接参考以下代码:

%%本版基于霍夫变换原理
%%2018.01.16 by Cooper Liu
%%Questions? Contact me: angelpoint@foxmail.com
clear;clc; %清空之前的变量
I=imread('line5.bmp'); %读取图像
level=graythresh(I); %使用最大类间方差法找到图片的一个合适的阈值
bw=im2bw(I,level); %根据阈值,使用im2bw函数将灰度图像转换为二值图像时
figure(1);imshow(bw);
[m,n]=size(bw); %获取尺寸
pMax=round(sqrt(m^2+n^2)); %计算最大p
thetaMax=180; %设定最大角度
countMatrix=zeros(pMax,thetaMax); %关于p和角度的计数矩阵
tic;
for i=1:m
    for j=1:n
        if bw(i,j)==0
            for theta=1:thetaMax %对theta作循环
                p=floor( abs( i*cos(3.14*theta/180) + j*sin(3.14*theta/180) ) ); %在theta循环过程中计算图像矩阵中一个像素点对应的p值
                countMatrix(p+1,theta)=countMatrix(p+1,theta)+1; %像素点对应的计数矩阵中(p+1,theta)处计数
            end
        end
    end
end
[m,n]=size(countMatrix);
for i=1:m
    for j=1:n
        if countMatrix(i,j)>countMatrix(1,1)
            countMatrix(1,1)=countMatrix(i,j); %获取最多曲线的相交点
            angle=j; %获取相交点处对应的角度
        end
    end
end
toc;
angle %这里得到的角度是 原点到直线的垂线 与X轴的夹角
if angle<=90
    rot=-angle;
else
    rot=180-angle;
end
pic=imrotate(I,rot,'crop'); %旋转图像
figure(2);imshow(pic);

如果本文有帮助到你,不如请我一罐可乐吧
在这里插入图片描述


  1. 参考文献:[1] 荆雷,张欣,郭金鑫.基于版面的拍照文档图像倾斜校正.激光与红外[J].2010,第10期 ↩︎

  2. 同时感谢愿意在网络上分享自己想法的各位博主。 ↩︎

2018-01-25 18:34:26 LQMIKU 阅读数 1552

文档图像倾斜角检测及校正(三)

996.icu LICENSE

  • 斜率空间投票原理
  • Matlab程序

阅读之前注意:

Hi,你好,我是Cooper Liu,欢迎来到我写的“文档图像校正”系列博客。基于三种原理,我写了四个实验性的Matlab验证程序,以及两个文档校正Matlab程序。在这里你将能够获取所有的源代码以及测试图片,完全可以在你自己的Matlab上跑这些程序。

如果你是学生,请注意不要抄袭,课程设计作业的话,这种程序仅仅只能让你得到80%左右的成绩。
如果你是工作人士或者只是感兴趣的极客,Okay,我想这些程序对于理解原理是如何应用为程序的已经足够 。

最后,请勿将这些资源用于商业用途(如你所见,这些程序都非常的初级)或者是谋取个人利益,知识在传播的过程中能展现更大的价值-

本文阅读建议用时:42min
本文阅读结构如下表:

项目 下属项目 测试用例数量
斜率空间投票原理 0
Matlab程序 2

斜率空间投票原理1

斜率空间投票指的是,在文档图像的某一个区域中,如果我们知道有一个点(x0,y0)(x_0,y_0 )是直线上的点,那我们将该点所在的直线斜率定为m,则坐标和斜率的关系可表示为y0=m(xx0)y_0=m(x-x_0)。由此图像区域中的每个目标像素(x,y)(x',y')(x0,y0)(x_0,y_0 )点之间连线的斜率mm'可表示为m=(yy0)/(xx0)m'=(y'-y_0)/(x'-x_0)

我们将斜率值映射到一组累加器上,每求出一个斜率值,对应的累加器加一。因为同一直线上的点求得的斜率相同,所以累加器会出现局部最大值,该值所对应的斜率即是所求直线的斜率,倾斜角也能相应的求出。

值得注意的是,当x=x0x'=x_0时,mm'为无穷大。为了避免这一情况,当x=x0x'=x_0时,我们令m=2m'=2,而当m>1m'>1m<1m'<1时,我们采用公式m=1/m2m''=1/m'-2的计算值代替mm'。这样一来,mm'的范围即被限定在(-1,3)的区间内。

Matlab程序

基于以上原理,我编写了Matlab程序来进行实验,源代码可以参考angleDetection3.m文件和angleDetection4.m文件。测试图片为line5.bmp。2
测试图片可以从这里获取,链接:https://pan.baidu.com/s/1dGmmGjn 密码:okt3

我编写了两个版本的Matlab程序,其中版本一(angleDetection3.m文件)的思想是把斜率区间尽可能的细分,比如分成10000份,再计算对应10000个累加器的最大值,从而反推出斜率。而版本二(angleDetection4.m文件)的思想则是先把斜率空间分为10个子区间,求出累加器的最大值,把最大值对应的子区间再分为10个子区间,如此迭代,直至区间宽度小于我们所控制的迭代精度。

以下是版本一的实验结果:
实验结果
时间

以下是版本二的实验结果:
实验结果
时间

如果您不想打开新的页面查看matlab源代码,也可以直接参考以下代码:
版本一(angleDetection3.m文件):

%%本版基于斜率空间投票原理 
%%2018.01.20_14:09 by Cooper Liu
%%Questions? Contact me: angelpoint@foxmail.com
clear;clc; %清空之前的变量
I=imread('line5.bmp'); %读取图像
level=graythresh(I); %使用最大类间方差法找到图片的一个合适的阈值
bw=im2bw(I,level); %根据阈值,使用im2bw函数将灰度图像转换为二值图像
figure(1);imshow(bw);
[m,n]=size(bw); %获取尺寸
flag=0;
tic; %计时开始
 for i=1:m
    for j=1:n
        if bw(i,j)==0
            if flag==0
                x0=i;y0=j; %获取第一个点
                flag=1;
            end
        end
    end
 end
range=10000; %range越大对角度分得更细,计算量更大
low=-1;
high=3;

a=zeros(range,1);
 for i=1:m
   for j=1:n
        if bw(i,j)==0
            if i==x0
                k=2;
                area=floor((k-(-1))/4*range); %向下取整
                a(area+1,1)=a(area+1,1)+1; %区间计数
            else
                k=(j-y0)/(i-x0);
                k2=k;
                if  k>1 || k<-1
                    k2=2-1/k; %转换斜率,保证斜率范围(-1, 3)
                end
                area=floor((k2-(-1))/4*range); %向下取整
                a(area+1,1)=a(area+1,1)+1;
            end 
        end
    end
 end
 toc %获取时间
[s_i,index] = sort(a,'descend'); %降序排序并获取原来序列的序号至index
max=index(1);
tmp=(-1)+max/range*4;
if tmp>1
    k=1/(2-tmp); %还原原来的斜率
else
    k=tmp;
end
angle=atan(k); %用atan求出来的角度在-pi/2到+pi/2之间
angle=angle*180/pi

rot=90-angle;
pic=imrotate(I,rot,'crop'); %旋转图像
figure(2);imshow(pic);

版本二(angleDetection4.m文件):

%%本版基于斜率空间投票原理,改进自angleDetection3.m 
%%2018.01.20_15:41 by Cooper Liu
%%Questions? Contact me: angelpoint@foxmail.com
clear;clc; %清空之前的变量
I=imread('line5.bmp'); %读取图像
level=graythresh(I); %使用最大类间方差法找到图片的一个合适的阈值
bw=im2bw(I,level); %根据阈值,使用im2bw函数将灰度图像转换为二值图像时
figure(1);imshow(bw);
[m,n]=size(bw); %获取尺寸
flag=0;

 for i=1:m
    for j=1:n
        if bw(i,j)==0
            if flag==0
                x0=i;y0=j;
                flag=1;
            end
        end
    end
 end
range=10; %整个斜率区间分为10个子区间
low=-1;
high=3;
tmpLow=low;
tmpHigh=high;
max=0;
areaAddtion=0; 
tic; %开始计时
while high-low>0.0001
a=zeros(range,1);
    for i=1:m
        for j=1:n
            if bw(i,j)==0
                if i==x0 %如果x坐标与x0相同
                    k=2; %则固定k为2
                    if k>=low && k<=high
                        area=floor((k-low)/(high-low)*range); %向下取整
                        if area==range %如果等于最后一个区间
                            area=area-1; %请考虑到后面区间计数的area+1
                        end
                        a(area+1,1)=a(area+1,1)+1; %区间计数
                    end
                else     %如果x坐标与x0不同
                    k=(j-y0)/(i-x0); %则计算斜率
                    k2=k;
                    if  k>1 || k<-1
                        k2=2-1/k;
                    end
                    if k2>=low && k2<=high
                        area=floor((k2-low)/(high-low)*range); %向下取整
                        if area==range
                            area=area-1;
                        end
                        a(area+1,1)=a(area+1,1)+1;
                    end
                end 
            end
        end
    end
[s_i,index] = sort(a,'descend');
max = index(1);
tmpLow = low;
tmpHigh = high;
if max>=2 && max<=9
    areaAddition=0; %根据需要可以更改子区间的范围
else
    areaAddtion=0;
end
low = tmpLow + (max - 1 - areaAddition)/range*(tmpHigh - tmpLow);
high = tmpLow + (max + areaAddtion)/range*(tmpHigh - tmpLow);
end

toc %获取计时时间

tmp=tmpLow+max/range*(tmpHigh-tmpLow);
if tmp>1
    k=1/(2-tmp);
else
    k=tmp;
end
angle=atan(k); %用atan求出来的角度在-pi/2到+pi/2之间
angle=angle*180/pi

rot=90-angle;
pic=imrotate(I,rot,'crop'); %旋转图像
figure(2);imshow(pic);

如果本文有帮助到你,不如请我一罐可乐吧 🎏
在这里插入图片描述


  1. 参考文献:[1] 荆雷,张欣,郭金鑫.基于版面的拍照文档图像倾斜校正.激光与红外[J].2010,第10期 ↩︎

  2. 同时感谢愿意在网络上分享自己想法的各位博主。 ↩︎

2013-12-08 21:25:26 thnh169 阅读数 5414

       1.图像的几何学变换

       之前的博文里,我简单的介绍了图像的放大与缩小。放大与缩小也算是图像的几何学变换,本文介绍了其他的几何学变换,包括旋转、水平倾斜和垂直倾斜(当然,还有水平移动与垂直移动。这些变换很简单,不需要插值,所以这里就不着重介绍了)。

       假设输入图像为G(u,v),其变换后的图像为F(x,y)。其变化的方法,如下所示。


       图像的几何学变换,主要有两种向前映射与向后映射。

       1.1向前映射

       所谓向前映射,就是从输入图像为g(u,v)的(0,0)点开始,将g(u,v)遍历一遍,依次计算g(u,v)变换后的坐标。当然,计算出来的坐标不会是整数(很大程度上不会是整数),就像下图所示那样。

      借由上图,我们重新理解一些向前映射。假设图像的一部分这5个点,经过变换,我们得到了右边的5个点。其实,这个5个点经过变换之后,计算出来的坐标并不是整数。假设现在遍历到这个点(图中以蓝色表示),通过计算得到的点是。我们将的坐标进行四舍五入,将的值赋值给(这种处理相当于选择最近的点,属于最初级的插值方法)。
       使用上述的向前映射的思想,将图像旋转,我们能得到如下效果。

       可以看出来,所得到的结果非常“斑驳”。其原因是,我们将g(u,v)进行变换,计算之后,所得到的坐标四舍五入之后,有的点没有被赋值,而有点被赋值了多次。这就产生了“空穴”,所以让变换后的图像非常斑驳。
       我看了一些文章,基本到这里都会说,“由于向前映射会产生空穴,所以向前映射一般不使用”。诸如此类的话,其实这样理解不太对,产生空穴的根本原因是由于插值法的不恰当所产生的。使用向前映射,也可以不产生“空穴”,所使用的方法如下所示。

       我们可以通过最接近的四个点,去确定的值。这样的话,就可以不产生空穴。但是,这个方法实现起来很困难,所以,向前映射才一般不被使用。

       1.2 向后映射

       向后映射,就是将输出图像f(x,y)遍历一遍,然后计算输出(x,y)的时候,所需要g(x,y)的坐标。当然,这个数不一定是整数。为了方便理解,还是看下图。

       假设,我们遍历到,通过计算,我们得到这样一个点。也就意味着,如果需要确定的灰度值,那么,我们需要点处的灰度值。这样,几何变换的问题又归结到了插值问题。最方便的办法就是选择最近的点,入上图的例子,的灰度值等于的灰度值。当然,还有更加“高级”点的方法,精度也还算可以。的灰度值的确定方法如下所示。

       如上图,我们使用四个点(),去确定的值。首先之间使用线性插值,确定出的值。同样的方法,确定出的值。之间,进行线性插值,就可以得到的值了。(图画的非常清楚了)
       使用向后映射,将图像旋转,我们可以得到如下结果。

       1.3 上述两部分的Matlab代码

close all;
clear all;
clc;

%% -----------Geometric_spatial_transform------------------
f = imread('./letter_T.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);

%-----by forward mapping-----%
g_fm = zeros(M,N);
seta = -pi/8;

for v = (-M/2):1:(M/2)-1
    for w = (-N/2):1:(N/2)-1
        x = round(v*cos(seta) - w*sin(seta));
        y = round(v*sin(seta) + w*cos(seta));
        if (((y>=(-N/2))&&(y<=(N/2)))&&((x>=(-M/2))&&(x<=(M/2))))
            g_fm(x+(M/2)+1,y+(N/2)+1) = f(v+(M/2)+1,w+(N/2)+1);
        end
    end
end

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(g_fm,[0 1]);
xlabel('b).Ruselt of Geometric spatial transform');

%-----by inverse mapping----
g_im = zeros(M,N);
seta = -pi/8;

for x = (-M/2):1:(M/2)-1
    for y = (-N/2):1:(N/2)-1
        v = x*cos(-seta) - y*sin(-seta);
        w = x*sin(-seta) + y*cos(-seta);
        if (((w>=(-N/2)+1)&&(w<=(N/2)-1))&&((v>=(-M/2)+1)&&(v<=(M/2)-1)))
            %g_im(x+(M/2)+1,y+(N/2)+1) = 1;%f(round(v+(M/2)+1),round(w+(N/2)+1));
            Q_11 = f(floor(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_21 =  f(floor(v)+(M/2)+1,ceil(w)+(N/2)+1);
            Q_12 =  f(ceil(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_22 =   f(ceil(v)+(M/2)+1,ceil(w)+(N/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;

            g_im(x+(M/2)+1,y+(N/2)+1) = (R2-R1)*(v-floor(v)) + R1;
        end
    end
end

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(g_im,[0 1]);
xlabel('b).Ruselt of Geometric spatial transform');

      2.图像的配准

      图像的配准,常常用于超解析领域。当然了,作为基础学习,博主没学那么深。这次的图像配准,试图去还原被倾斜变换的图像。首先,先将图像进行两个方向的倾斜。

       我们这次的目的是,将图像b).还原为a).。这次变换相对简单,我们使用如下模型去拟合变换关系。

我们在原图与变换后的图像上,选择出四个标准点,然后带入方程,并求出系数。使用这个方程,将图像b).还原为a)。这个过程比较简单,下面是结果。另外,所选择的标准点,我已经在上图中用红圈标定出来了。

       可以看到,还原的效果非常的好(额,其实这也是应为畸变比较简单的缘故)。为了看出于原图的区别,我们做出了差分图像。可以看出来,还原不是完美的。
       至此,上个代码作为本文的结尾。额,由于特征点的选择是手工的,这个代码可能没有多少意义。(不要吐槽最后这一句话。一般的,在使用特征点做图象配准的时候,一般选择在某个实际的物体上放入某个实际的标志,然后通过检测这个实际的标志,去进行变换。【这里可以参考《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods的第二章的相关内容,这里有叙述的!!】)
close all;
clear all;
clc;
%% --------------image registration---------------------------
f_Original = imread('./characters_test_pattern.tif');
f_Original = mat2gray(f_Original,[0 255]);
[M,N] = size(f_Original);

f = zeros(M+4,N+4);
for x = 1:M
    f(x+2,:) = [0 0 f_Original(x,:) 0 0];
end

[M,N] = size(f);

s_v = 0.3;
s_w = 0.02;

P = round(M+s_v*N);
Q = round(N+s_w*M);
g_1 = zeros(P,N);

for x = (-P/2):1:(P/2)-1
    for y = (-N/2):1:(N/2)-1
        v = x - s_v * y;
        w = y;
        if (((w>=(-N/2))&&(w<=(N/2)-1))&&((v>=(-M/2))&&(v<=(M/2)-1)))
            Q_11 = f(floor(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_21 =  f(floor(v)+(M/2)+1,ceil(w)+(N/2)+1);
            Q_12 =  f(ceil(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_22 =   f(ceil(v)+(M/2)+1,ceil(w)+(N/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;

            g_1(x+(P/2)+1,y+(N/2)+1) = (R2-R1)*(v-floor(v)) + R1;
        end
    end
end

g = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
    for y = (-Q/2):1:(Q/2)-1
        v = x;
        w = y - s_w * x;
        if (((w>=(-N/2))&&(w<=(N/2)-1))&&((v>=(-P/2))&&(v<=(P/2)-1)))
            Q_11 = g_1(floor(v)+(P/2)+1,floor(w)+(N/2)+1);
            Q_21 =  g_1(floor(v)+(P/2)+1,ceil(w)+(N/2)+1);
            Q_12 =  g_1(ceil(v)+(P/2)+1,floor(w)+(N/2)+1);
            Q_22 =   g_1(ceil(v)+(P/2)+1,ceil(w)+(N/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;

            g(x+(P/2)+1,y+(Q/2)+1) = (R2-R1)*(v-floor(v)) + R1;
        end
    end
end

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');
hold on;
plot(621, 78,'ro','MarkerSize',7);
plot(116,113,'ro','MarkerSize',7);
plot( 85,649,'ro','MarkerSize',7);
plot(624,641,'ro','MarkerSize',7);

subplot(1,2,2);
imshow(g,[0 1]);
xlabel('b).Ruselt of Geometric spatial transform');
hold on;
plot(625,264,'ro','MarkerSize',7);
plot(116,147,'ro','MarkerSize',7);
plot( 96,674,'ro','MarkerSize',7);
plot(638,828,'ro','MarkerSize',7);
%% 
f_reg = zeros(M,N);

Original = [621-(N/2)  78-(M/2)  (621-(N/2))*(78-(M/2)) 1;
            116-(N/2) 113-(M/2) (116-(N/2))*(113-(M/2)) 1;
             85-(N/2) 649-(M/2)  (85-(N/2))*(649-(M/2)) 1;
            624-(N/2) 641-(M/2) (624-(N/2))*(641-(M/2)) 1];
        
Ruselt = [625-(Q/2);116-(Q/2);96-(Q/2);638-(Q/2)];

c = (Original)\Ruselt;  %C1~C4

Ruselt = [264-(P/2);147-(P/2);674-(P/2);828-(P/2)];

c = [c (Original)\Ruselt];  %C5~C8 wwwwww

for y = (-M/2):1:(M/2)-1
    for x = (-N/2):1:(N/2)-1
        w = c(1,1)*y + c(2,1)*x + c(3,1)*x*y + c(4,1);
        v = c(1,2)*y + c(2,2)*x + c(3,2)*x*y + c(4,2);
        if (((w>=(-Q/2))&&(w<=(Q/2)-1))&&((v>=(-P/2))&&(v<=(P/2)-1)))
            
            Q_11 = g(floor(v)+(P/2)+1,floor(w)+(Q/2)+1);
            Q_21 =  g(floor(v)+(P/2)+1,ceil(w)+(Q/2)+1);
            Q_12 =  g(ceil(v)+(P/2)+1,floor(w)+(Q/2)+1);
            Q_22 =   g(ceil(v)+(P/2)+1,ceil(w)+(Q/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;
            
            f_reg(x+(N/2)+1,y+(M/2)+1) = (R2-R1)*(v-floor(v)) + R1;
            %g(round(v+(P/2)+1),round(w+(Q/2)+1)) = 0.5;
        end
    end
end

g_diff = abs(f - f_reg);

figure();
subplot(1,2,1);
imshow(f_reg,[0 1]);
xlabel('c).Ruselt of image registration');

subplot(1,2,2);
imshow(g_diff,[0 1]);
xlabel('d).Difference image');
      原文发于博客:http://blog.csdn.net/thnh169/ 

仪表识别方法汇总

阅读数 1881

没有更多推荐了,返回首页