圖像矯正原理說明

        在進行光學掃描時，會因為客觀原因，導致掃描的圖像位置不正，影響后期的圖像處理，因此需對圖像進行圖像矯正工作。

1、圖像傾斜矯正基礎

        圖像傾斜矯正關鍵在於根據圖像特征自動檢測出圖像傾斜方向和傾斜角度。目前常用的傾斜角度方法有：基於投影的方法、基於Hough變換、基於線性擬合，還有進行傅里葉變換到頻域來進行檢測的方法。

1.1 基於投影的傾斜矯正

      基於投影的矯正方法，利用圖像水平方向和垂直方向的投影特征來進行傾斜判斷。主要用到投影的方差和均方差、投影特征矢量、梯度方向場等統計特性。

      在圖像投影中，一條直線沿着它的發現方向投影最長，沿着水平方向投影最短，此稱之為Radon變換。定義：二元函數f（x,y)的投影是在某一方向上的線積分，例如f(x,y)在垂直方向上的線積分是f(x,y)在x方向上的投影，在水平方向上的線積分是在y方向上的投影，沿y'方向的線積分是沿x'方向上的投影。

圖1 矩形函數在水平垂直方向和沿θ角方向的投影

投影可沿任意角度進行，通常f(x,y)的Radon變換是f(x,y)平行於y’軸的線積分，格式如下：

其中

1.2 Hough變換法

        Hough 變換是數字圖象技術中一種有效的發現直線的算法 .它是先把直角坐標系的目標點映射到極坐標系上進行累積 ,即它是先使直角坐標系平面上任一直線上的所有點均累積到極坐標系的同一點集中去 ,然后通過尋找極坐標系中點集的峰值 ,來發現長的直線特征 .由於這種點集是通過累積統計得到的 ,因而能夠容忍直線的間斷 。

       原理很簡單：假設有一條與原點距離為S，方向角為0的一條直線，如下圖所示：

直線上的每一點滿足：

證明過程如下：

在上圖中有效根據三角形相似原理進行證明，找出線段與x，y，θ之間的關系：

所以：

然后我們再在直線上找一一點x1,y1只要證明也滿足這個公式即可。在小三角形中：

所以：

進而：

 

所以直線上任意一點滿足上述公式。

在x-o-y平面內可以看出（x1,y1)即屬於直線L1又屬於直線L2，且滿足：

 

所以在x-o-y平面上的一點（x1,y1)對應於s-o-θ面上的的一條曲線，也就意味着在x-o-y平面內的直線對應s-o-θ存在很多條曲線上。同時由於一條直線上的點都滿足公式：

，即上述公式中我們知道s,θ為一個定值，所以意味着是個定點，所以所以x-o-y平面上處在一條直線上的點經過變換在S-O-θ平面上所得的曲線相交於一點。如下圖所示

 

因此可以把x-o-y平面內直線的問題轉化為S-O-θ平面內點的問題。

1.3、線性回歸

如果我們有一系列相互獨立的點，其近似分布在一條直線附近，我們可以通過線性回歸擬合這條直線。

 

這條直線的一元線性回歸模型為：

設上圖中點的坐標分別是(x₁,y₁) , (x₂,y₂) , …… , (x_n,y_n) 。用最小二乘算法來估計和。則有：

I = 1 , 2 , 3 , … ,n ，用最小二乘算法來計算，找到准則函數，記為

求Q的最小值：

由上式可得：

其中：

 

估計出了直線的斜率和截距，我們就可以估計出這條直線的方程了。

2、圖像矯正算法實現

 

2.1 Rando變換算法實現

利用Radon變換檢測直線傾斜角度的具體步驟為：

（1）用edge函數計算圖像的邊緣二值圖像，檢測出原始圖像中的直線。

（2）計算邊緣圖像的Radon變換，對每一個象素為1的點進行運算（0-179度方向上分別做投影）

（3）檢測出Radon變換矩陣中的峰值，這些峰值對應原始圖像中的直線（上圖中的四個亮點對應圖3.9中的四條直線）。Radon變換矩陣中的這些峰值的列坐標θ就是與原始圖像中的直線垂直的直線的傾斜角度，所以圖像中直線的傾角為90-θ。

bw1=edge(l1,'sobel', 'horizontal');%用Sobel水平算子對圖像邊緣化
bw1=imcrop(bw1,[0 0 500 100]);%對圖像進行剪切，保留圖像中的一條直線，減小運算量
theta=0:179;%定義theta角度范圍
r=radon(bw1,theta);%對圖像進行Radon變換
[m,n]=size(r);
c=1;　　　
for i=1:m
    for j=1:n
        if  r(1,1)<r(i,j)
           r(1,1)=r(i,j);
            c=j;
        end
    end
end      %檢測Radon變換矩陣中的峰值所對應的列坐標
rot=90-c;%確定旋轉角度
pic=imrotate(l1,rot,'crop');%對圖像進行旋轉矯正

2.2 Hough算法實現

（1）對圖像進行邊緣檢測，這里選用了Sobel算子檢測圖像中水平方向的直線。

（2）假設圖像對應於x-o-y空間，定義一個S-o-θ（θ角的范圍為1-180）空間，對圖像中象素為1的每一個點進行計算（應用公式（3.10）），做出每一個象素為1的點的曲線，同時把S-θ平面分成等間隔(1×1)的小網格，這個小網格對應一個記數矩陣。如圖3.5所示，凡是曲線所經過的網格，對應的記數矩陣元素值加1，所以對原圖像中的每一點進行計算以后記數矩陣元素的值等於共線的點數。我們可以認為記數矩陣中元素的最大值對應原始圖像中最長的直線。

（3）檢測出記數矩陣的最大的元素所對應的列坐標θ，θ即為這條直線的法線與X軸的夾角。因此我們可以通過θ角來確定直線的傾斜角度，進而對圖像進行矯正。

 

Hough變換法矯正圖像程序實現如下：

bw=edge(l,'sobel','horizontal');%檢測圖像邊緣直線
[m,n]=size(bw);%計算圖像大小
S=round(sqrt(m^2+n^2));%S可以取到的最大值
ma=180;%θ角最大值
r=zeros(md,ma);%產生初值為零的計數矩陣
for i=1:m
    for j=1:n
        if bw(i,j)==1
            for k=1:ma
                ru=round(abs(i*cos(k*3.14/180)+j*sin(k*3.14/180)));
                r(ru+1,k)=r(ru+1,k)+1;%對矩陣記數
            end
        end
    end
end
[m,n]=size(r);
for i=1:m
    for j=1:n
        if r(i,j)>r(1,1)
r(1,1)=r(i,j);
            c=j;%把矩陣元素最大值所對應的列坐標送給c。
        end
    end
end
if c<=90
rot=-c;    %確定旋轉角度
else
    rot=180-c;
end
pic=imrotate(l,rot,'crop'); %對圖片進行旋轉，矯正圖像

下面是python版本的，里面講解的比較詳細

"""
Created on Tue May 30 13:59:20 2017

@author: Francesco Lacriola
"""
import numpy as np
from PIL import ImageDraw,Image

class Hough(object):
    
    def __init__(self,image, n_theta, n_rho, plot_point):
        self.image= image     #immagine da analizzare
        step=0                #
        self.n_theta= n_theta #numero di angoli da analizzare per ogni punto
        self.n_rho= n_rho     #numero di intervalli di rho
        thetas = np.deg2rad(np.linspace(-90.0, 90.0,self.n_theta)) #array contenente gli angoli da analizzare
        width, height = self.image.shape
        diag_len = np.ceil(np.sqrt(width * width + height * height))   # il rho massimo (la diagonale dell'immagine)
        rhos = np.linspace(-diag_len, diag_len, self.n_rho+1)          # i limiti degli intervalli
        # Cache some resuable values
        cos_t = np.cos(thetas)  #calcola una volta tutti i cos degli angoli da analizzare e li salva in un array
        sin_t = np.sin(thetas)  #calcola una volta tutti i sin degli angoli da analizzare e li salva in un array
        num_thetas = len(thetas)
        
        accumulator = np.zeros((self.n_rho, num_thetas), dtype=np.uint64) #accumulatore
        y_idxs, x_idxs = np.nonzero(self.image)  # una lista di x ed y contenenti le cordinate dei punti bianchi (non zero)
       # Vote in the hough accumulator
        for i in range(len(x_idxs)): #cicla su tutti i punti bianchi (x è la x del punto corrente, y la sua y)
           x = x_idxs[i]
           y = y_idxs[i]
    
           for t_idx in range(num_thetas): #cicla su tutti gli angoli da analizzare
               # Calculate rho. diag_len is added for a positive index 
               rho = x * cos_t[t_idx] + y * sin_t[t_idx] #calcola rho
               ind_rho = Hough.__binary_int__(rhos,rho) #approssima rho ad un indice all'interno dell'array degli intervalli
               accumulator[ind_rho, t_idx] += 1 #assegna il voto
               if plot_point is not None: #un controllo limite
                   plot_point(rho, np.rad2deg(thetas[t_idx%accumulator.shape[1]]), accumulator[ind_rho, t_idx], step) #aggiunge un punto al grafico
                   step+=1
        self.accumulator=accumulator
        self.thetas= thetas
        self.rhos= rhos
    
    def getHoughImage(self):
        accmax= np.argmax(self.accumulator)
        redacc=self.accumulator/self.accumulator[int(accmax/self.accumulator.shape[1])][accmax%self.accumulator.shape[1]]
        return redacc
    
    def getPatternImage(self):
        image=np.zeros(np.shape(self.image))
        image=Image.fromarray((image).astype('uint8'))
        accmax= np.argmax(self.accumulator)
        maxacc=self.accumulator[int(accmax/self.accumulator.shape[1])][accmax%self.accumulator.shape[1]]
        draw = ImageDraw.Draw(image)
        ord_accumulator= self.accumulator
        ind_min=np.argmin(ord_accumulator)
        i=int(ind_min/self.accumulator.shape[1])
        j=ind_min%self.accumulator.shape[1]
        while ord_accumulator[i][j]<maxacc:
             rho = self.rhos[i]
             theta = self.thetas[j]
             y1=(rho/np.sin(theta))-(np.cos(theta)/np.sin(theta)*(np.shape(self.image)[1]-1))
             y2=(rho/np.sin(theta))-(np.cos(theta)/np.sin(theta)*0)
             #per trovare linee orizzontali if (y2-y1)/(-(np.shape(self.image)[1]-1)) >= -0.1 and (y2-y1)/(-(np.shape(self.image)[1]-1)) <=0.1:
             draw.line((np.shape(self.image)[1]-1,y1, 0,y2), fill=int((ord_accumulator[i][j]/maxacc)*255))
             ord_accumulator[i][j]=maxacc
             ind_min=np.argmin(ord_accumulator)
             i=int(ind_min/self.accumulator.shape[1])
             j=ind_min%self.accumulator.shape[1]
        i=int(accmax/self.accumulator.shape[1])
        j=accmax%self.accumulator.shape[1]
        rho = self.rhos[i]
        theta = self.thetas[j]
        y1=(rho/np.sin(theta))-(np.cos(theta)/np.sin(theta)*(np.shape(self.image)[1]-1))
        y2=(rho/np.sin(theta))-(np.cos(theta)/np.sin(theta)*0)
        #per trovare linee orizzontali if (y2-y1)/(-(np.shape(self.image)[1]-1)) >= -0.1 and (y2-y1)/(-(np.shape(self.image)[1]-1)) <=0.1:
        draw.line((np.shape(self.image)[1]-1,y1, 0,y2), fill=int((ord_accumulator[i][j]/maxacc)*255))
        return image
    
    def __binary_int__(array,val):
        assert val<= array[len(array)-1] and val>=array[0]
        found=False
        maxm=len(array)
        start = int((maxm/2))
        minm=0
        if val == array[0]:
            return 0
        while not found:
            if val<=array[start]:
                maxm=start
                start=int(((minm+start)/2))
                
            elif val>array[start+1]:
                minm=start
                start=int(((maxm+start)/2))
                
            else:
                found=True
        return start

2.3 線性回歸算法實現

能夠檢測到圖像上邊緣的一系列的邊界點，我們就可以通過最小二乘法擬合這條邊界直線，從而確定圖像的傾角。

具體方法如下：

（1）找出邊界直線上的點（每列第一次由黑變白的點，且這一列的下兩點還是白的話就可以判為邊界點^[6]），將其行坐標存入數組a即，列坐標存入數組b。

（2）通過最小二乘法擬和這條邊界直線，計算出其斜率L。

（3）通過rot=atan（L），計算直線的傾斜角度，然后對其矯正。

[m,n]=size(l);
bw=im2bw(l,0.3);
%將圖像二值化
t=1;s=1;
for j=144:1: n-144
    for i=1:fix(m/4)
if bw(i,j)==0&bw(i+1,j)==1&bw(i+2,j)==1&bw(i+3,j)==1
%檢測邊緣點
            c(t)=i;%邊緣橫坐標存入數組C
            b1(s)=j;邊緣縱坐標存入數組b1
          break;
       end
    end 
       t=t+1;s=s+1;
end
x=0;y=0;x1=0;
for i=1:length(c)
x=x+c(i);
 x1=x1+c(i)^2;
   end 
   for i=1:length(b1)
      y=y+b1(i);
  end
      y=y/length(b1);
      x=x/length(c);
      c1=x;
      x=length(c)*x*x;
      lxy=x1-x;
      lxx=0;
for i=1:length(b1)
    lxx=lxx+(c(i)-c1)*(b1(i)-y);
end
r=lxy/lxx;
%以上為計算直線參數
rot=atan(r);
%取余切
theta=rot*180/3.142;
%將弧度轉換為角度
pic=imrotate(l,theta,'crop');