###### tags: `Learning Opencv`

# 投影轉換

---

## 原理

在前一張所提及的方法皆應用於二維空間的轉換，如果物體是在3D空間發生旋轉，一樣可以透過方程式解法(最小平方法擬合)，OpenCv 提供函式進行解析，由於多一個Z方向的解，因此需要四組方程組才有辦法得到解析解。

cv2.getPerspectiveTransform(src, dst)

$\bigl(\begin{smallmatrix}
\tilde{x}\\\\
\tilde{y}\\\\
\tilde{z}
\end{smallmatrix}\bigr) = \bigl(\begin{smallmatrix}
a_{11} &a_{12} & a_{13}\\\\
a_{21} &a_{22} & a_{23}\\\\
a_{31} &a_{32} & a_{33}
\end{smallmatrix}\bigr)
\bigl(\begin{smallmatrix}
x\\\\
y\\\\
z
\end{smallmatrix}\bigr)$

### Example

假設(0, 0)、(200, 0)、(0, 200)、(200, 200)是原座標，透過某投影轉換依序為(100, 20)、(200, 20)、(50, 70)、(250, 70)

sol:

import cv2
import numpy as np
src = np.array([[0, 0], [200, 0], [0, 200], [200, 200]], np.float32)
dst = np.array([[100, 20], [200, 20], [50, 70], [250, 70]], np.float32)
p = cv2.getPerspectiveTransform(src, dst)

傳回的投影矩陣為：

array([[ 5.00e-01, -3.75e-01,  1.00e+02],
       [ 0.00e+00,  7.50e-02,  2.00e+01],
       [-0.00e+00, -2.50e-03,  1.00e+00]])

## 圖片投影轉換

跟2D的類似方法，不過函式名稱改為

cv2.warpPerspective(src, M, dsize)

import numpy as np
import cv2
import sys

if __name__ == "__main__":
    if len(sys.argv) > 1:
        image = cv2.imread(sys.argv[1], 0)
    else:
        print("Usage: python warpPercpective.py image")
        exit()
        
    h, w = image.shape[:2]
    src = np.array([[0, 0], [w - 1, 0], [0, h - 1], [w - 1, h - 1]])
    dst = np.array([[50, 50], [w / 3, 50], [50, h - 1], [w - 1, h - 1]])
    
    p = cv2.getPerspectiveTransform(src, dst)
    
    r = cv2.warpPerspective(image, p, (w, h), borderValue=125)
    
    cv2.imshow("warpPerspective", r)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

## 極座標轉換

通常利用極座標轉換來校正影像中的圓形物體或被包含在圓環中的物體。

### 將笛卡兒座標轉為極座標

將從(x, y)到($\theta , r$)表示，公式如下：

$r = \sqrt{(x - \bar{x})^2 + (y - \bar{y})^2}$

$\theta = \left\{\begin{matrix}
2 \pi + \arctan2(y - \bar{y}, x - \bar{x}) &, y - \bar{y}\leq 0\\
\arctan2(y - \bar{y}, x - \bar{x}) &, y - \bar{y}\geq 0
\end{matrix}\right.$

OpenCv 提供函式

import cv2
cv2.cartToPolar(x, y)

### Example

import cv2
import numpy as np

x = np.array([[0, 1, 2], [0, 1, 2], [0, 1, 2]], np.float64)
y = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], np.float64)
r, theta = cv2.cartToPolar(x, y, angleInDegrees=True)

>>> r
array([[0.        , 1.        , 2.        ],
       [1.        , 1.41421356, 2.23606798],
       [2.        , 2.23606798, 2.82842712]])
>>> theta
array([[ 0.        ,  0.        ,  0.        ],
       [90.        , 44.99045563, 26.56710434],
       [90.        , 63.43289566, 44.99045563]])

### 將極座標轉為笛卡兒座標

極座標轉換是可逆的，再以知道極座標與笛卡兒座標的條件下可以透過以下式子得到笛卡兒座標

$x = \bar{x} + r\cos \theta$

$y = \bar{y} + r \sin \theta$

cv2.polarToCart(magnitude, angle)

### Example

import cv2
import numpy as np
angle = np.array([[30, 31], [30, 31]], np.float32)
r = np.array([[10, 10], [11, 11]], np.float32)
# 這樣得到的x, y 是針對(0, 0) 進行旋轉
x, y = cv2.polarToCart(r, angle, angleInDegrees=True)

# 只要加上 -12 與 15，就可以針對不同中心進行轉換
x += -12
y += 15

>>> x
array([[8.660255, 8.571674],
       [9.52628 , 9.428843]], dtype=float32)
>>> y
array([[5.0000005, 5.150382 ],
       [5.5000005, 5.66542  ]], dtype=float32)

### 極座標對影像進行轉換

假設輸入影像矩陣為I，$(\bar{x}, \bar{y})$代表極座標空間轉換的中心，輸出矩陣為O，直觀的寫法為下式表示

$O(r, \theta) = f_1 (\bar{x}+r\cos\theta, \bar{y} + r\sin\theta)$

上述的$\theta$與$r$都以1為單位進行離散化，由於轉換步距較大，輸出的影像O可能損失原圖許多資訊，可以透過以下方法進行改進

1. 設定範圍 $r_{min}, r_{max}$、$\theta_{min}, \theta_{max}$ 內進行座標轉換，因此步距值可以用$r_{step}$進行表示

公式最後可以變成

$O(i, j) = f_1 (\bar{x} + r_{min} + r_{step}i * \cos(\theta_{min} + \theta_{step}j), \, \bar{y} + (r_{min} + r_{step}i) * \cos(\theta_{min} + \theta_{step}j))$

#### tips: numpy tile(a, (m, n))

import numpy as np
a = np.array([[1, 2], [3, 4]])
b = np.tile(a, (2, 3)) # 將A矩陣分別在垂直方向和水平方向複製2次與3次

array([[1, 2, 1, 2, 1, 2],
       [3, 4, 3, 4, 3, 4],
       [1, 2, 1, 2, 1, 2],
       [3, 4, 3, 4, 3, 4]])

### polar 函數來實現影像極座標轉換

1. I 代表輸入影像
2. center 代表極座標轉換中心
3. r 代表最小距離與最大距離
4. theta 代表角度範圍 (0, 360)
5. rstep 代表r的轉換步進值
6. thetastep 代表角度的轉換步進值

def polar(I, center, r, theta=(0, 360), rstep=1.0, thetastep=360.0 / (180 * 8)):
    minr, maxr = r
    mintheta, maxtheta = theta
    H = int((maxr - minr) / rstep) + 1
    W = int((maxtheta - mintheta) / thetastep) + 1
    O = 125 * np.ones((H, W), I.dtype)
    
    r = np.linspace(minr, maxr, H)
    r = tile(r, (W, 1))
    r = np.transpose(r)
    theta = np.linspace(mintheta, maxtheta, W)
    theta = np.tile(theta, (H, 1))
    x, y = cv2.polarToCart(r, theta, angleInDegrees=True)
    
    # 最近鄰插值
    for i in range(H):
        for j in range(W):
            px = int(round(x[i][j] + cx))
            py = int(round(y[i][j] + cy))
            if ((px >= 0 and px <= w - 1) and (py >= 0 and py <= h - 1)):
                O[i][j] = I[py][px]
    return O

import cv2
import numpy as np
import sys

if __name__ == "__main__":
if len(sys.argv) > 1:
image = cv2.imread(sys.argv[1], 0)
else:
print("Usage: python warpPercpective.py image")
exit()
h, w = I.sha[:2]
cx, cy = 508, 503
cv2.circle(I, (int(cx), int(cy)), 10, (255.0, 0, 0), 3)
O = polar(I, (cx, cy), (200, 550))
O = cv2.flip(O, 0)
cv2.imshow("O", O)
cv2.waitKey(0)
cv2.destroyAllWindows()