| 2 |
3 |
| 3
| 6 |
6 |
7 |
如果不考慮像素位置的影響，那麼找最小cost的時候會認為左上角的cost=2最小。然而如果考慮到像素間距的影響，我們來看左上角方向，和中心的差異是6-2，做個線性插值的話，則左上角距中心單位距離上的值應該為(6-left( 6- 2 right)times 1/sqrt{2} =3.17>3)，所以正上方的才是最小cost的正確方向。
最短路徑查找
在磁力套索中，一般的用法是先單擊一個點，然後移動滑鼠，在滑鼠和一開始單擊的點之間就會出現自動貼近邊緣的線，這裡我們定義一開始點擊的像素點為種子點(seed)，而磁力套索其實在考慮上部分提到的邊緣相關cost的情況下查找種子點到當前滑鼠的最短路徑。如下圖，紅色的就是種子點，而移動滑鼠時，最貼近邊緣的種子點和滑鼠座標的連線就會即時顯示，這也是為什麼磁力套索也叫Livewire。
實現最短路徑的辦法很多，一般而言就是動態規劃了，這裡介紹的是基於Dijkstra演算法的一種實現，基本思想是，給定種子點後，執行Dijkstra演算法將圖像的所有像素遍歷，得到每個像素到種子點的最短路徑。以下面這張圖為例，在一個cost矩陣中，利用Dijkstra演算法遍歷每一個元素後，每個元素都會指向一個相鄰的元素，這樣任一像素都能找到一條到seed的路徑，例如右上角的42和39對應的像素，沿著箭頭到了0。
演算法如下：
输入:
s // 种子点
l(q,r) // 计算局部cost
数据结构:
L // 当前待处理的像素
N(q) // 当前像素相邻的像素
e(q) // 标记一个像素是否已经做过相邻像素展开的Bool函数
g(q) // 从s到q的总cost
输出:
p // 记录所有路径的map
算法:
g(s)←0; L←s; // 将种子点作为第一点初始化
while L≠Ø: // 遍历尚未结束
q←min(L); // 取出最小cost的像素并从待处理像素中移除
e(q)←TRUE; // 将当前像素记录为已经做过相邻像素展开
for each r∈N(q) and not e(r):
gtemp←g(q)+l(q,r); // 计算相邻像素的总cost
if r∈L and gtemp<g(r): // 找到了更好的路径
r←L; { from list.} // 舍弃较大cost的路径
else if r∉L:
g(r)←gtemp; // 记录当前找到的最小路径
p(r)←q;
L←r; // 加入待处理以试图寻找更短的路径
遍歷的過程會優先經過cost最低的區域，如下圖：
所有像素對應的到種子像素的最短路徑都找到後，移動鼠標時就直接畫出到seed的最短路徑就可以了。
Python實作
演算法部分直接呼叫了OpenCV的Canny函數和Sobel函數(求梯度)，對於RGB的處理也很簡陋，直接用梯度最大的值來近似。另外因為懶，cost map和path map都直接用了字典(dict)，而記錄展開過的像素則直接採用了集合(set)。 GUI部分因為不會用QThread所以用了Python的threading，只有圖像顯示互動區域和狀態欄提示，左鍵點擊設定種子點，右鍵結束，已經提取的邊緣為綠色線，正在提取的為藍色線。
程式碼
演算法部分
from __future__ import division
import cv2
import numpy as np
SQRT_0_5 = 0.70710678118654757
class Livewire():
"""
A simple livewire implementation for verification using
1. Canny edge detector + gradient magnitude + gradient direction
2. Dijkstra algorithm
"""
def __init__(self, image):
self.image = image
self.x_lim = image.shape[0]
self.y_lim = image.shape[1]
# The values in cost matrix ranges from 0~1
self.cost_edges = 1 - cv2.Canny(image, 85, 170)/255.0
self.grad_x, self.grad_y, self.grad_mag = self._get_grad(image)
self.cost_grad_mag = 1 - self.grad_mag/np.max(self.grad_mag)
# Weight for (Canny edges, gradient magnitude, gradient direction)
self.weight = (0.425, 0.425, 0.15)
self.n_pixs = self.x_lim * self.y_lim
self.n_processed = 0
@classmethod
def _get_grad(cls, image):
"""
Return the gradient magnitude of the image using Sobel operator
"""
rgb = True if len(image.shape) > 2 else False
grad_x = cv2.Sobel(image, cv2.CV_64F, 1, 0, ksize=3)
grad_y = cv2.Sobel(image, cv2.CV_64F, 0, 1, ksize=3)
if rgb:
# A very rough approximation for quick verification...
grad_x = np.max(grad_x, axis=2)
grad_y = np.max(grad_y, axis=2)
grad_mag = np.sqrt(grad_x**2+grad_y**2)
grad_x /= grad_mag
grad_y /= grad_mag
return grad_x, grad_y, grad_mag
def _get_neighbors(self, p):
"""
Return 8 neighbors around the pixel p
"""
x, y = p
x0 = 0 if x == 0 else x - 1
x1 = self.x_lim if x == self.x_lim - 1 else x + 2
y0 = 0 if y == 0 else y - 1
y1 = self.y_lim if y == self.y_lim - 1 else y + 2
return [(x, y) for x in xrange(x0, x1) for y in xrange(y0, y1) if (x, y) != p]
def _get_grad_direction_cost(self, p, q):
"""
Calculate the gradient changes refer to the link direction
"""
dp = (self.grad_y[p[0]][p[1]], -self.grad_x[p[0]][p[1]])
dq = (self.grad_y[q[0]][q[1]], -self.grad_x[q[0]][q[1]])
l = np.array([q[0]-p[0], q[1]-p[1]], np.float)
if 0 not in l:
l *= SQRT_0_5
dp_l = np.dot(dp, l)
l_dq = np.dot(l, dq)
if dp_l < 0:
dp_l = -dp_l
l_dq = -l_dq
# 2/3pi * ...
return 0.212206590789 * (np.arccos(dp_l)+np.arccos(l_dq))
def _local_cost(self, p, q):
"""
1. Calculate the Canny edges & gradient magnitude cost taking into account Euclidean distance
2. Combine with gradient direction
Assumption: p & q are neighbors
"""
diagnol = q[0] == p[0] or q[1] == p[1]
# c0, c1 and c2 are costs from Canny operator, gradient magnitude and gradient direction respectively
if diagnol:
c0 = self.cost_edges[p[0]][p[1]]-SQRT_0_5*(self.cost_edges[p[0]][p[1]]-self.cost_edges[q[0]][q[1]])
c1 = self.cost_grad_mag[p[0]][p[1]]-SQRT_0_5*(self.cost_grad_mag[p[0]][p[1]]-self.cost_grad_mag[q[0]][q[1]])
c2 = SQRT_0_5 * self._get_grad_direction_cost(p, q)
else:
c0 = self.cost_edges[q[0]][q[1]]
c1 = self.cost_grad_mag[q[0]][q[1]]
c2 = self._get_grad_direction_cost(p, q)
if np.isnan(c2):
c2 = 0.0
w0, w1, w2 = self.weight
cost_pq = w0*c0 + w1*c1 + w2*c2
return cost_pq * cost_pq
def get_path_matrix(self, seed):
"""
Get the back tracking matrix of the whole image from the cost matrix
"""
neighbors = [] # 8 neighbors of the pixel being processed
processed = set() # Processed point
cost = {seed: 0.0} # Accumulated cost, initialized with seed to itself
paths = {}
self.n_processed = 0
while cost:
# Expand the minimum cost point
p = min(cost, key=cost.get)
neighbors = self._get_neighbors(p)
processed.add(p)
# Record accumulated costs and back tracking point for newly expanded points
for q in [x for x in neighbors if x not in processed]:
temp_cost = cost[p] + self._local_cost(p, q)
if q in cost:
if temp_cost < cost[q]:
cost.pop(q)
else:
cost[q] = temp_cost
processed.add(q)
paths[q] = p
# Pop traversed points
cost.pop(p)
self.n_processed += 1
return paths
livewire.pylivewire
from __future__ import division
import time
import cv2
from PyQt4 import QtGui, QtCore
from threading import Thread
from livewire import Livewire
class ImageWin(QtGui.QWidget):
def __init__(self):
super(ImageWin, self).__init__()
self.setupUi()
self.active = False
self.seed_enabled = True
self.seed = None
self.path_map = {}
self.path = []
def setupUi(self):
self.hbox = QtGui.QVBoxLayout(self)
# Load and initialize image
self.image_path = ''
while self.image_path == '':
self.image_path = QtGui.QFileDialog.getOpenFileName(self, '', '', '(*.bmp *.jpg *.png)')
self.image = QtGui.QPixmap(self.image_path)
self.cv2_image = cv2.imread(str(self.image_path))
self.lw = Livewire(self.cv2_image)
self.w, self.h = self.image.width(), self.image.height()
self.canvas = QtGui.QLabel(self)
self.canvas.setMouseTracking(True)
self.canvas.setPixmap(self.image)
self.status_bar = QtGui.QStatusBar(self)
self.status_bar.showMessage('Left click to set a seed')
self.hbox.addWidget(self.canvas)
self.hbox.addWidget(self.status_bar)
self.setLayout(self.hbox)
def mousePressEvent(self, event):
if self.seed_enabled:
pos = event.pos()
x, y = pos.x()-self.canvas.x(), pos.y()-self.canvas.y()
if x < 0:
x = 0
if x >= self.w:
x = self.w - 1
if y < 0:
y = 0
if y >= self.h:
y = self.h - 1
# Get the mouse cursor position
p = y, x
seed = self.seed
# Export bitmap
if event.buttons() == QtCore.Qt.MidButton:
filepath = QtGui.QFileDialog.getSaveFileName(self, 'Save image audio to', '', '*.bmp\n*.jpg\n*.png')
image = self.image.copy()
draw = QtGui.QPainter()
draw.begin(image)
draw.setPen(QtCore.Qt.blue)
if self.path_map:
while p != seed:
draw.drawPoint(p[1], p[0])
for q in self.lw._get_neighbors(p):
draw.drawPoint(q[1], q[0])
p = self.path_map[p]
if self.path:
draw.setPen(QtCore.Qt.green)
for p in self.path:
draw.drawPoint(p[1], p[0])
for q in self.lw._get_neighbors(p):
draw.drawPoint(q[1], q[0])
draw.end()
image.save(filepath, quality=100)
else:
self.seed = p
if self.path_map:
while p != seed:
p = self.path_map[p]
self.path.append(p)
# Calculate path map
if event.buttons() == QtCore.Qt.LeftButton:
Thread(target=self._cal_path_matrix).start()
Thread(target=self._update_path_map_progress).start()
# Finish current task and reset
elif event.buttons() == QtCore.Qt.RightButton:
self.path_map = {}
self.status_bar.showMessage('Left click to set a seed')
self.active = False
def mouseMoveEvent(self, event):
if self.active and event.buttons() == QtCore.Qt.NoButton:
pos = event.pos()
x, y = pos.x()-self.canvas.x(), pos.y()-self.canvas.y()
if x < 0 or x >= self.w or y < 0 or y >= self.h:
pass
else:
# Draw livewire
p = y, x
path = []
while p != self.seed:
p = self.path_map[p]
path.append(p)
image = self.image.copy()
draw = QtGui.QPainter()
draw.begin(image)
draw.setPen(QtCore.Qt.blue)
for p in path:
draw.drawPoint(p[1], p[0])
if self.path:
draw.setPen(QtCore.Qt.green)
for p in self.path:
draw.drawPoint(p[1], p[0])
draw.end()
self.canvas.setPixmap(image)
def _cal_path_matrix(self):
self.seed_enabled = False
self.active = False
self.status_bar.showMessage('Calculating path map...')
path_matrix = self.lw.get_path_matrix(self.seed)
self.status_bar.showMessage(r'Left: new seed / Right: finish')
self.seed_enabled = True
self.active = True
self.path_map = path_matrix
def _update_path_map_progress(self):
while not self.seed_enabled:
time.sleep(0.1)
message = 'Calculating path map... {:.1f}%'.format(self.lw.n_processed/self.lw.n_pixs*100.0)
self.status_bar.showMessage(message)
self.status_bar.showMessage(r'Left: new seed / Right: finish')
gui.pygui.py
主函數
import sys
from PyQt4 import QtGui
from gui import ImageWin
def main():
app = QtGui.QApplication(sys.argv)
window = ImageWin()
window.setMouseTracking(True)
window.setWindowTitle('Livewire Demo')
window.show()
window.setFixedSize(window.size())
sys.exit(app.exec_())
if __name__ == '__main__':
main()
main.py
main.py蛋疼地上傳到了Github(傳送門)，歡迎fork。
效率的改進因為這個程式碼的原型只是為了用C++開發之前的Python評估和驗證，所以完全沒考慮效率，執行速度是完全不行的，基本上400x400的圖片就不能忍了…至於基於Python版本的效率提升我沒有仔細想過，只是大概看來有這麼多比較明顯的地方：
1) 取出當前最小cost像素操作
p = min(cost, key=cost.get )
這個雖然寫起來很爽但顯然是不行的，至少得用個min heap什麼的。因為我是用dict同時表示待處理像素和cost了，也懶得想一下怎麼和Python的heapq結合起來，所以直接用了粗暴省事的min()。
2) 梯度方向的計算
三角函數的計算應該是盡量避免的，另外原文可能是為了將值域擴展到>π所以把q-p也用上了，其實這一項本來權重就小，那怕直接用兩個像素各自的梯度方向向量做點積然後歸一化結果也是還行的。即使要用arccos，也可以考慮寫look-up table近似。當然我最後想說的是個人覺得其實這項真沒那麼必要，直接自適應spilne或者那怕三點均值平滑去噪效果就不錯了。
3) 計算相鄰像素的位置
如果兩個像素相鄰，則他們各自周圍的8個相鄰像素也會重疊。的我的辦法比較原始，可以考率不用模組化直接計算。
4) 替換部分資料結構
例如path map其實本質是給出每個像素在最短路徑上的上一個像素，是個矩陣。其實可以考慮用線性的資料結構來代替，但如果真這樣做一般來說都是在C/C++程式碼裡了。
5) numpy
我印像中對numpy的調用順序也會影響到效率，連續調用numpy的內置方法似乎會帶來效率的整體提升，不過話還是說回來，實際應用中如果到了這一步，應該也屬於C/C++程式碼範疇了。
6) 演算法層面的改進
這塊沒有深入研究，第一感覺是實際應用中沒必要一上來就計算整幅圖像，可以根據seed位置做一些區塊劃分，滑鼠本身也會留下軌跡，也或許可以考慮只在滑鼠軌跡方向進行啟發式搜尋。另外計算路徑的時候或許可以考慮借鏡有點類似Image Pyramid的思想，沒必要一上來就對全解析度下的路徑進行查找。由於後來做的專案沒有採用這個演算法，所以我也沒有繼續研究，雖然挺好奇的，其實有好多現成的程式碼，比如GIMP，不過沒有精力去看了。
更多的改進
雖然都沒做，大概介紹一下，都是考慮了實用性的改進。
路徑冷卻(Path Cooling)
用過Photoshop和GIMP磁力套索的人都知道，即使滑鼠不點擊圖片，在移動過程中也會自動產生一些將摳圖軌跡固定住的點，這些點其實就是新的種子點，而這個使用過程中自動產生新的種子點的方法叫Path cooling。這個方法的基本想法如下：隨著滑鼠移動過程中如果一定時間內一段路徑都保持固定不變，那麼就把這段路徑中離種子最遠的點設定為新的種子，其實背後隱藏的還是動態規劃的思想，貝爾曼最優。這個名字也是比較形象的，路徑冷卻。
動態訓練(Interactive Dynamic Training)
單純的最短路徑查找在使用的時候常常出現找到的邊緣不是想要的邊緣的問題，比如上圖，綠色的線是上一段提取的邊緣，藍色的是目前正在提取的邊緣。左圖中，鏡子外面Lena的帽子邊緣是我們想要提取的，然而由於鏡子裡的Lena的帽子邊緣的cost更低，所以實際提取出的藍色線段如右圖中貼到右邊了。所以Interactive Dynamic Training的想法是，認為綠色的線段是正確提取的邊緣，然後利用綠色線段作為訓練資料來為目前提取邊緣的cost函數附加一個修正值。
[1]中所採用的方法是統計前一段邊緣上點的梯度強度的直方圖，然後按照梯度出現頻率給目前圖中的像素加權。舉例來說如果綠色線段中的所有像素對應的梯度強度都是在50到100之間的話，那麼可以將50到100以10為單位分為5個bin，統計每個bin裡的出現頻率，也就是直方圖，然後對目前偵測到的梯度強度，做個線性加權。比方說50~60區間內對應的像素最多有10個，那麼把10作為最大值，並且對目前檢測到的梯度強度處於50~60之間的像素均乘上係數1.0；如果訓練資料中70~ 80之間有5個，那麼cost加權係數為5/10=0.5，則對所有目前檢測到的梯度強度處於70~80之間的像素均乘上係數0.5；如果訓練資料中100以上沒有，所以cost附加為0/10=0，則加權係數為0，這樣即使偵測到更強的邊緣也不會偏離前段邊緣了。這是基本思想，當然實際的實作沒有這麼簡單，除了邊緣上的像素還要考慮垂直邊緣上左邊和右邊的兩個像素點，這樣保證了邊緣的pattern。另外隨著滑鼠越來越遠離訓練邊緣，偵測到的邊緣的pattern可能會出現不一樣，所以Training可能會起反作用，所以這種Training的作用範圍也需要考慮到滑鼠離種子點的距離，最後還要有一些平滑去噪的處理，具體都在[1]裡有講到，挺繁瑣的(那會好像還沒有SIFT)，不詳述了。
種子點位置的修正(Cursor Snap)
雖然這個演算法可以自動找出種子點和滑鼠之間最貼近邊緣的路徑，不過，人的手，常常抖，所以種子點未必能很好地設置到邊緣。所以可以在使用者設定完種子點位置之後，自動在其座標周圍小範圍內，例如7x7的區域內搜尋cost最低的像素，作為真正的種子點位置，這個過程叫做Cursor snap。
更多Photoshop中磁力套索的一種簡陋實作(基於Python) 相關文章請關注PHP中文網！
