本來是給一個圖像分割項目算法評估時的Python框架,覺得有點意思,就稍稍拓展了一下,用PyQt加了個殼,在非常簡陋的程度上模擬了一下的磁力套索功能。為什么簡陋:1) 只實現(xiàn)了最基本的邊緣查找。路徑冷卻,動態(tài)訓(xùn)練,鼠標(biāo)位置修正都沒有,更別提曲線閉合,摳圖,Alpha Matting等等;2) 沒考慮性能規(guī)范,只為了書寫方便;3) 我對Qt了解很淺,至今不會寫Signal-Slot,不知道GUI寫得是否合理;4) 沒調(diào)試。
基本算法
相關(guān)算法我并沒有做很深入的調(diào)研,不過相信這類應(yīng)用中影響力最大的算法是來源于[1],也是本文的主要參考,基本思想是把圖片看成是一個無向圖,相鄰像素之間就可以計算出一個局部cost,于是就轉(zhuǎn)化成了最短路徑問題了,接下來就是基于Dijkstra算法產(chǎn)生路徑,就是需要提取的邊緣。主要涉及的算法有兩部分:1) 相鄰像素的cost計算;2) 最短路徑算法。
邊緣檢測
計算相鄰像素cost的最終目的還是為了尋找邊緣,所以本質(zhì)還是邊緣檢測?;舅枷胧牵ㄟ^各種不同手段檢測邊緣,并且根據(jù)檢測到的強度來求加權(quán)值,作為cost。從最短路徑的角度來說,就是邊緣越明顯的地方,cost的值越小。[1]中的建議是用三種指標(biāo)求加權(quán):1) 邊緣檢測算子;2) 梯度強度(Gradient Magnitude);3) 梯度方向(Gradient Direction)。本文的方法和[1]有那么一些不一樣,因為懶,用了OpenCV中的Canny算子檢測邊緣而不是Laplacian Zero-Crossing Operator。表達式如下:
[lleft( p,q ight)={{w}_{E}}{{f}_{E}}left( q ight)+{{w}_{G}}{{f}_{G}}left( q ight)+{{w}_{D}}{{f}_{D}}left( p,q ight)]
Canny算子
基本思想是根據(jù)梯度信息,先檢測出許多連通的像素,然后對于每一坨連通的像素只取其中最大值且連通的部分,將周圍置零,得到初始的邊緣(Edges),這個過程叫做Non-Maximum Suppression。然后用二閾值的辦法將這些檢測到的初始邊緣分為Strong, Weak, and None三個等級,顧名思義,Strong就是很確定一定是邊緣了,None就被舍棄,然后從Weak中挑選和Strong連通的作為保留的邊緣,得到最后的結(jié)果,這個過程叫做Hysteresis Thresholding。這個算法太經(jīng)典,更多細節(jié)一Google出來一大堆,我就不贅述了。公式如下:
[{{f}_{E}}left( q
ight)=left{ egin{matrix}
0; ext{ if }q ext{ is on a edge} \
1; ext{ if }q ext{ is not on a edge} \
end{matrix}
ight.]
其實從權(quán)值的計算上和最大梯度有些重復(fù),因為如果沿著最大梯度方向找出來的路徑基本上就是邊緣,這一項的作用我的理解主要應(yīng)該是1) 避免梯度都很大的區(qū)域出現(xiàn)離明顯邊緣的偏離;2) 保證提取邊緣的連續(xù)性,一定程度上來講也是保證平滑。
梯度強度
就是梯度求模而已,x和y兩個方向的梯度值平方相加在開方,公式如下:
[{{I}_{G}}left( q ight)=sqrt{{{I}_{x}}left( q ight)+{{I}_{y}}left( q ight)}]
因為要求cost,所以反向并歸一化:
[{{f}_{G}}left( q ight)=1-frac{{{I}_{G}}left( q ight)}{max left( {{I}_{G}} ight)}]
梯度方向
這一項其實是個平滑項,會給變化劇烈的邊緣賦一個比較高的cost,讓提取的邊緣避免噪聲的影響。具體公式如下:
[{{f}_{D}}left( p,q ight)=frac{2}{3pi }left( arccos left( {bbtvznv_{p}}left( p,q ight) ight)+arccos left( {f3pdpvb_{q}}left( p,q ight) ight) ight)]
其中,
[{ffhjlrx_{p}}left( p,q ight)=leftlangle {b33399d_{ot }}left( p ight),{{l}_{D}}left( p,q ight) ight angle ]
[{3399vd9_{q}}left( p,q ight)=leftlangle {{l}_{D}}left( p,q ight),{99dpr9b_{ot }}left( q ight) ight angle ]
[{{l}_{D}}left( p,q
ight)=left{ egin{matrix}
q-p; ext{ if }leftlangle {zxjlnd9_{ot }}left( p
ight),q-p
ight
angle ge 0 \
p-q; ext{ if }leftlangle {fdx3zhd_{ot }}left( p
ight),q-p
ight
angle <0 \
end{matrix}
ight.]
({lblfhdl_{ot }}left( p ight))是取p的垂直方向,另外注意上式中符號的判斷會將({3dx3jpf_{ot }}left( p ight))和({{l}_{D}}left( p,q ight))的取值限制在π/2以內(nèi)。
[{jfz3939_{ot }}left( p ight)=left( {{p}_{y}},-{{p}_{x}} ight)]
斜對角方向的cost修正
在二維圖像中,相鄰的像素通常按照間隔歐式距離分為兩種:1) 上下左右相鄰,間隔為像素邊長;2) 斜對角相鄰,間隔為像素邊長的(sqrt{2})倍。在計算局部cost的時候通常要把這種距離差異的影響考慮進去,比如下面這幅圖:
2 | 3 | 4 |
5 | 6 | 6 |
7 | 8 | 9 |
如果不考慮像素位置的影響,那么查找最小cost的時候會認(rèn)為左上角的cost=2最小。然而如果考慮到像素間距的影響,我們來看左上角方向,和中心的差異是6-2,做個線性插值的話,則左上角距中心單位距離上的值應(yīng)該為(6-left( 6-2 ight) imes 1/sqrt{2} =3.17>3),所以正上方的才是最小cost的正確方向。
最短路徑查找
在磁力套索中,一般的用法是先單擊一個點,然后移動鼠標(biāo),在鼠標(biāo)和一開始單擊的點之間就會出現(xiàn)自動貼近邊緣的線,這里我們定義一開始單擊的像素點為種子點(seed),而磁力套索其實在考慮上部分提到的邊緣相關(guān)cost的情況下查找種子點到當(dāng)前鼠標(biāo)的最短路徑。如下圖,紅色的就是種子點,而移動鼠標(biāo)時,最貼近邊緣的種子點和鼠標(biāo)坐標(biāo)的連線就會實時顯示,這也是為什么磁力套索也叫Livewire。
實現(xiàn)最短路徑的辦法很多,一般而言就是動態(tài)規(guī)劃了,這里介紹的是基于Dijkstra算法的一種實現(xiàn),基本思想是,給定種子點后,執(zhí)行Dijkstra算法將圖像的所有像素遍歷,得到每個像素到種子點的最短路徑。以下面這幅圖為例,在一個cost矩陣中,利用Dijkstra算法遍歷每一個元素后,每個元素都會指向一個相鄰的元素,這樣任意一個像素都能找到一條到seed的路徑,比如右上角的42和39對應(yīng)的像素,沿著箭頭到了0。
算法如下:
輸入: s // 種子點 l(q,r) // 計算局部cost 數(shù)據(jù)結(jié)構(gòu): L // 當(dāng)前待處理的像素 N(q) // 當(dāng)前像素相鄰的像素 e(q) // 標(biāo)記一個像素是否已經(jīng)做過相鄰像素展開的Bool函數(shù) g(q) // 從s到q的總cost
遍歷的過程會優(yōu)先經(jīng)過cost最低的區(qū)域,如下圖:
所有像素對應(yīng)的到種子像素的最短路徑都找到后,移動鼠標(biāo)時就直接畫出到seed的最短路徑就可以了。
Python實現(xiàn)
算法部分直接調(diào)用了OpenCV的Canny函數(shù)和Sobel函數(shù)(求梯度),對于RGB的處理也很簡陋,直接用梯度最大的值來近似。另外因為懶,cost map和path map都直接用了字典(dict),而記錄展開過的像素則直接采用了集合(set)。GUI部分因為不會用QThread所以用了Python的threading,只有圖像顯示交互區(qū)域和狀態(tài)欄提示,左鍵點擊設(shè)置種子點,右鍵結(jié)束,已經(jīng)提取的邊緣為綠色線,正在提取的為藍色線。
代碼
算法部分
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.py
livewire.py
GUI部分
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 *.jpg *.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.py
gui.py
主函數(shù)
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++開發(fā)之前的Python評估和驗證,所以完全沒考慮效率,執(zhí)行速度是完全不行的,基本上400x400的圖片就不能忍了……至于基于Python版本的效率提升我沒有仔細想過,只是大概看來有這么幾個比較明顯的地方:
1) 取出當(dāng)前最小cost像素操作
p = min(cost, key=cost.get)
這個雖然寫起來很爽但顯然是不行的,至少得用個min heap什么的。因為我是用dict同時表示待處理像素和cost了,也懶得想一下怎么和Python的heapq結(jié)合起來,所以直接用了粗暴省事的min()。
2) 梯度方向的計算
三角函數(shù)的計算應(yīng)該是盡量避免的,另外原文可能是為了將值域擴展到>π所以把q-p也用上了,其實這一項本來權(quán)重就小,那怕直接用兩個像素各自的梯度方向向量做點積然后歸一化一下結(jié)果也是還行的。即使要用arccos,也可以考慮寫個look-up table近似。當(dāng)然我最后想說的是個人覺得其實這項真沒那么必要,直接自適應(yīng)spilne或者那怕三點均值平滑去噪效果就不錯了。
3) 計算相鄰像素的位置
如果兩個像素相鄰,則他們各自周圍的8個相鄰像素也會重合。的我的辦法比較原始,可以考率不用模塊化直接計算。
4) 替換部分?jǐn)?shù)據(jù)結(jié)構(gòu)
比如path map其實本質(zhì)是給出每個像素在最短路徑上的上一個像素,是個矩陣。其實可以考慮用線性的數(shù)據(jù)結(jié)構(gòu)代替,不過如果真這樣做一般來說都是在C/C++代碼里了。
5) numpy
我印象中對numpy的調(diào)用順序也會影響到效率,連續(xù)調(diào)用numpy的內(nèi)置方法似乎會帶來效率的整體提升,不過話還是說回來,實際應(yīng)用中如果到了這一步,應(yīng)該也屬于C/C++代碼范疇了。
6) 算法層面的改進
這塊沒有深入研究,第一感覺是實際應(yīng)用中沒必要一上來就計算整幅圖像,可以根據(jù)seed位置做一些區(qū)塊劃分,鼠標(biāo)本身也會留下軌跡,也或許可以考慮只在鼠標(biāo)軌跡方向進行啟發(fā)式搜索。另外計算路徑的時候也許可以考慮借鑒有點類似于Image Pyramid的思想,沒必要一上來就對全分辨率下的路徑進行查找。由于后來做的項目沒有采用這個算法,所以我也沒有繼續(xù)研究,雖然挺好奇的,其實有好多現(xiàn)成的代碼,比如GIMP,不過沒有精力去看了。
更多的改進
雖然都沒做,大概介紹一下,都是考慮了實用性的改進。
路徑冷卻(Path Cooling)
用過Photoshop和GIMP磁力套索的人都知道,即使鼠標(biāo)不點擊圖片,在移動過程中也會自動生成一些將摳圖軌跡固定住的點,這些點其實就是新的種子點,而這種使用過程中自動生成新的種子點的方法叫Path cooling。這個方法的基本思路如下:隨著鼠標(biāo)移動過程中如果一定時間內(nèi)一段路徑都保持固定不變,那么就把這段路徑中離種子最遠的點設(shè)置為新的種子,其實背后隱藏的還是動態(tài)規(guī)劃的思想,貝爾曼最優(yōu)。這個名字也是比較形象的,路徑冷卻。
動態(tài)訓(xùn)練(Interactive Dynamic Training)
單純的最短路徑查找在使用的時候常常出現(xiàn)找到的邊緣不是想要的邊緣的問題,比如上圖,綠色的線是上一段提取的邊緣,藍色的是當(dāng)前正在提取的邊緣。左圖中,鏡子外面Lena的帽子邊緣是我們想要提取的,然而由于鏡子里的Lena的帽子邊緣的cost更低,所以實際提取出的藍色線段如右圖中貼到右邊了。所以Interactive Dynamic Training的思想是,認(rèn)為綠色的線段是正確提取的邊緣,然后利用綠色線段作為訓(xùn)練數(shù)據(jù)來給當(dāng)前提取邊緣的cost函數(shù)附加一個修正值。
[1]中采用的方法是統(tǒng)計前一段邊緣上點的梯度強度的直方圖,然后按照梯度出現(xiàn)頻率給當(dāng)前圖中的像素加權(quán)。舉例來說如果綠色線段中的所有像素對應(yīng)的梯度強度都是在50到100之間的話,那么可以將50到100以10為單位分為5個bin,統(tǒng)計每個bin里的出現(xiàn)頻率,也就是直方圖,然后對當(dāng)前檢測到的梯度強度,做個線性加權(quán)。比方說50~60區(qū)間內(nèi)對應(yīng)的像素最多有10個,那么把10作為最大值,并且對當(dāng)前檢測到的梯度強度處于50~60之間的像素均乘上系數(shù)1.0;如果訓(xùn)練數(shù)據(jù)中70~80之間有5個,那么cost加權(quán)系數(shù)為5/10=0.5,則對所有當(dāng)前檢測到的梯度強度處于70~80之間的像素均乘上系數(shù)0.5;如果訓(xùn)練數(shù)據(jù)中100以上沒有,所以cost附加為0/10=0,則加權(quán)系數(shù)為0,這樣即使檢測到更強的邊緣也不會偏離前一段邊緣了。這是基本思想,當(dāng)然實際的實現(xiàn)沒有這么簡單,除了邊緣上的像素還要考慮垂直邊緣上左邊和右邊的兩個像素點,這樣保證了邊緣的pattern。另外隨著鼠標(biāo)越來越遠離訓(xùn)練邊緣,檢測到的邊緣的pattern可能會出現(xiàn)不一樣,所以Training可能會起反作用,所以這種Training的作用范圍也需要考慮到鼠標(biāo)離種子點的距離,最后還要有一些平滑去噪的處理,具體都在[1]里有講到,挺繁瑣的(那會好像還沒有SIFT),不詳述了。
種子點位置的修正(Cursor Snap)
雖然這個算法可以自動找出種子點和鼠標(biāo)之間最貼近邊緣的路徑,不過,人的手,常常抖,所以種子點未必能很好地設(shè)置到邊緣上。所以可以在用戶設(shè)置完種子點位置之后,自動在其坐標(biāo)周圍小范圍內(nèi),比如7x7的區(qū)域內(nèi)搜索cost最低的像素,作為真正的種子點位置,這個過程叫做Cursor snap。
更多Photoshop中磁力套索的一種簡陋實現(xiàn)(基于Python) 相關(guān)文章請關(guān)注PHP中文網(wǎng)!
聲明:本網(wǎng)頁內(nèi)容旨在傳播知識,若有侵權(quán)等問題請及時與本網(wǎng)聯(lián)系,我們將在第一時間刪除處理。TEL:177 7030 7066 E-MAIL:11247931@qq.com