Trích đặc trưng SIFT trong OpenCV
OpenCV - Tut 23: Trích đặc trưng SIFT
Tut 23: Trích đặc trưng SIFT
Hi mọi người, Minh đã trở lại sau ~6 năm kể từ Tut 21: Local Binary Patterns (LBP) và Tut 22: Gabor Filters, series OpenCV của mình tiếp tục với một phương pháp trích đặc trưng kinh điển, đó là: SIFT.
SIFT là viết tắt của Scale-Invariant Feature Transform. Nếu dịch sát nghĩa thì hơi cứng, mình tạm hiểu là: một phép biến đổi giúp tìm ra các điểm đặc trưng cục bộ của ảnh sao cho chúng tương đối bền vững khi ảnh bị thay đổi kích thước, xoay, nhiễu nhẹ hoặc thay đổi sáng ở mức vừa phải.
Nói nôm na: nếu trong ảnh có một góc cửa sổ, logo, họa tiết trên áo, chữ trên bảng hiệu, vân gạch, … thì SIFT cố gắng tìm ra những điểm "đáng nhớ" đó và mô tả mỗi điểm bằng một vector đặc trưng. Sau này ta có thể đem các vector đó đi matching giữa hai ảnh.
Có một slide rất hay giảng dạy về phương pháp trích đặc trưng SIFT: https://homes.cs.washington.edu/~shapiro/EE596/notes/SIFT.pdf. Các hình ảnh minh họa lý thuyết trong bài viết này M sẽ dùng ảnh từ slide để trực quan và dễ hình dung hơn.
Ứng dụng của phương pháp trích đặc trưng SIFT
Trong computer vision truyền thống, SIFT thường được dùng cho:
- Ghép ảnh panorama / image stitching.
- Object matching: tìm cùng một vật thể trong nhiều ảnh khác nhau.
- Image retrieval: tìm ảnh giống nhau trong một database.
- 3D reconstruction / Structure from Motion.
- Theo dõi đặc trưng cục bộ giữa các frame.
- Làm baseline để hiểu các thuật toán feature extraction trước thời deep learning.
Điểm hay của SIFT là nó không cố mô tả toàn bộ bức ảnh bằng một vector duy nhất. Thay vào đó, nó tìm nhiều local keypoints và mô tả từng keypoint. Vì vậy nếu ảnh bị che khuất một phần, còn lại một số keypoint khớp được thì ta vẫn có hy vọng nhận ra vật thể.
Nguồn gốc và paper liên quan
SIFT được đề xuất bởi David G. Lowe. Ý tưởng ban đầu xuất hiện trong paper ICCV 1999:
- David G. Lowe, "Object recognition from local scale-invariant features", International Conference on Computer Vision (ICCV), 1999. Trang tổng hợp của tác giả
Phiên bản hoàn chỉnh và được trích dẫn rất nhiều nằm trong paper IJCV 2004 (được >82,771 bài báo tham chiếu đến theo Google Scholar):
- David G. Lowe, "Distinctive Image Features from Scale-Invariant Keypoints", International Journal of Computer Vision, 60(2), 91–110, 2004. PDF
Trang chủ của tác giả ghi chú rằng SIFT từng có bằng sáng chế tại Mỹ và đã hết hạn vào 07/03/2020. Vì vậy trong các phiên bản OpenCV hiện đại, bạn có thể dùng cv2.SIFT_create() trực tiếp thay vì phải đi qua xfeatures2d như trước đây.
Một tài liệu hiện thực rất đáng đọc thêm là:
- Ives Rey-Otero and Mauricio Delbracio, "Anatomy of the SIFT Method", IPOL, 2014. Bài viết
Trong bài này, mình sẽ giải thích theo tinh thần paper của tác giả Lowe, nhưng code cuối bài sẽ là một bản hiện thực đơn giản từ đầu bằng Python + OpenCV, để các bạn debug và nhìn rõ từng bước. Code này dùng để học thuật toán, không nhằm thay thế bản SIFT tối ưu trong OpenCV nha.
Đặc trưng của SIFT
Giả sử ta chụp cùng một vật thể ở hai ảnh khác nhau:
- Ảnh thứ nhất chụp gần.
- Ảnh thứ hai chụp xa hơn.
- Ảnh thứ ba bị xoay nhẹ.
- Ảnh thứ tư sáng hơn hoặc tối hơn một chút.
Nếu chỉ so sánh pixel trực tiếp thì rất dễ sai biệt. Pixel ở ảnh này không còn nằm đúng tọa độ ở ảnh kia nữa. Do đó SIFT không đi theo hướng "so từng pixel", mà đi theo hướng:
- Tìm những điểm nổi bật cục bộ trong ảnh.
- Ước lượng scale phù hợp cho từng điểm.
- Gán hướng chính cho từng điểm để chống xoay.
- Mô tả vùng lân cận quanh điểm đó bằng histogram theo gradient có hướng.
Kết quả cuối cùng của SIFT gồm:
- Keypoints: vị trí, scale, orientation, response, …
- Descriptors: mỗi keypoint thường có một vector 128 chiều.
Tổng quan các bước của SIFT
Theo paper Lowe 2004, pipeline chính của SIFT gồm 4 bước lớn:
- Scale-space extrema detection: tìm ứng viên keypoint trong không gian scale bằng Difference of Gaussian (DoG).
- Keypoint localization: tinh chỉnh vị trí keypoint, loại điểm yếu và điểm nằm trên cạnh.
- Orientation assignment: gán một hoặc nhiều hướng chính cho mỗi keypoint.
- Keypoint descriptor: tạo vector descriptor 128 chiều từ histogram hướng gradient quanh keypoint.
Nghe hơi căng ha :D. Giờ ta bóc từng lớp một.
1. Xây dựng Gaussian scale-space
Một vật thể trong ảnh có thể xuất hiện ở nhiều kích thước khác nhau. Cùng một cái logo, nếu chụp gần thì nó to, chụp xa thì nó nhỏ. Vậy nếu ta chỉ tìm đặc trưng ở đúng một kích thước ảnh, khả năng cao là sẽ bỏ sót.
SIFT xử lý chuyện này bằng cách xây dựng scale-space. Ý tưởng là tạo nhiều phiên bản làm mờ Gaussian của cùng một ảnh:
\[L(x, y, \sigma) = G(x, y, \sigma) * I(x, y)\]Trong đó:
- \(I(x, y)\) là ảnh đầu vào.
- \(G(x, y, \sigma)\) là Gaussian kernel với độ lệch chuẩn \(\sigma\).
- Dấu \(*\) là convolution.
- \(L(x, y, \sigma)\) là ảnh sau khi làm mờ ở scale \(\sigma\).
Khi \(\sigma\) nhỏ, ảnh chỉ bị làm mờ nhẹ, còn nhiều chi tiết nhỏ. Khi \(\sigma\) lớn, ảnh bị làm mờ mạnh hơn, các chi tiết nhỏ biến mất, chỉ còn cấu trúc lớn.
SIFT tổ chức các ảnh này thành nhiều octave. Mỗi octave là một nhóm ảnh cùng kích thước nhưng khác mức blur. Sang octave tiếp theo, ảnh thường được downsample xuống còn một nửa chiều rộng và chiều cao.
2. Difference of Gaussian (DoG)
Lowe dùng Difference of Gaussian để tìm điểm đặc trưng. Công thức:
\[D(x, y, \sigma) = L(x, y, k\sigma) - L(x, y, \sigma)\]Nói đơn giản: lấy hai ảnh Gaussian blur gần nhau rồi trừ chúng cho nhau.
Tại sao lại làm vậy? Vì DoG là một xấp xỉ hiệu quả của Laplacian of Gaussian (LoG). LoG có khả năng phát hiện blob / cấu trúc nổi bật trong ảnh, nhưng tính trực tiếp thì tốn hơn. DoG rẻ hơn vì ta đã có sẵn các ảnh Gaussian trong pyramid.
Ở bước này, với mỗi octave ta có một chồng ảnh DoG. Mỗi pixel trong chồng DoG không chỉ có tọa độ \(x, y\) mà còn có tọa độ scale \(\sigma\).
3. Tìm extrema trong không gian 3D
Một điểm ứng viên keypoint được chọn nếu nó là cực đại hoặc cực tiểu cục bộ trong không gian 3D: x, y, scale.
Cụ thể, một pixel ở ảnh DoG hiện tại sẽ được so với 26 láng giềng:
- 8 điểm xung quanh trong cùng layer scale.
- 9 điểm ở layer scale phía trên.
- 9 điểm ở layer scale phía dưới.
Nếu nó lớn hơn tất cả 26 điểm lân cận, nó là local maximum. Nếu nó nhỏ hơn tất cả 26 điểm lân cận, nó là local minimum. Cả hai trường hợp đều là ứng viên keypoint.
4. Keypoint localization và loại điểm không ổn định
Không phải extrema nào cũng tốt. Có hai loại điểm ta nên loại bỏ:
4.1. Điểm có contrast thấp
Nếu giá trị DoG tại một điểm quá nhỏ, điểm đó có thể sinh ra bởi nhiễu hoặc vùng ảnh ít thông tin. Lowe đề xuất loại các điểm có độ tương phản thấp sau khi tinh chỉnh vị trí.
Trong code học tập cuối bài, mình dùng một ngưỡng contrast_threshold đơn giản để loại những điểm có giá trị DoG quá nhỏ.
4.2. Điểm nằm dọc theo cạnh
Một cạnh dài cũng có thể tạo response lớn trong DoG, nhưng vị trí dọc theo cạnh lại không ổn định. Trượt một chút dọc theo cạnh thì nhìn vẫn gần như giống nhau, nên matching dễ bị sai.
SIFT dùng ma trận Hessian 2D tại điểm ứng viên:
\[H = \begin{bmatrix} D_{xx} & D_{xy} \\ D_{xy} & D_{yy} \end{bmatrix}\]Sau đó kiểm tra tỉ lệ giữa độ cong chính. Paper Lowe dùng điều kiện dựa trên trace (Tr) và determinant (Det):
\[\frac{Tr(H)^2}{Det(H)} < \frac{(r + 1)^2}{r}\]Trong đó \(r\) thường được chọn là 10. Nếu tỉ lệ quá lớn, điểm đó giống cạnh hơn là góc/đốm (corner/blob) ổn định, ta loại.
Đến đây M xin phép làm rõ thêm chi tiết cách tính toán: ma trận Hessian 2D, trace và det.
Cách xây dựng ma trận Hessian 2D
Ở đây ta đang làm việc trên ảnh DoG \(D(x, y)\). Tại một điểm ứng viên \((x, y)\), ta quan tâm đến việc giá trị DoG cong theo trục \(x\) và \(y\) ra sao. Trong giải tích, điều này được mô tả bằng các đạo hàm bậc hai:
- \(D_{xx}\): độ cong theo trục \(x\)
- \(D_{yy}\): độ cong theo trục \(y\)
- \(D_{xy}\): mức độ hai trục \(x, y\) "dính" với nhau như thế nào
Trên ảnh số, ta thường không có công thức liên tục nên sẽ xấp xỉ các đạo hàm này bằng sai phân hữu hạn quanh một vùng 3x3.
Ta có thể ước lượng:
\[D_{xx} \approx D(x+1, y) + D(x-1, y) - 2D(x, y)\] \[D_{yy} \approx D(x, y+1) + D(x, y-1) - 2D(x, y)\] \[D_{xy} \approx \frac{D(x+1, y+1) + D(x-1, y-1) - D(x+1, y-1) - D(x-1, y+1)}{4}\]Ví dụ ta có một vùng DoG 3x3 quanh điểm đang xét như sau, với tâm là giá trị 10:
\[\begin{bmatrix} 3 & 9 & 5 \\ 4 & 10 & 6 \\ 1 & 14 & 6 \end{bmatrix}\]Dựa theo các công thức ta tính toán được:
\[D_{xx} \approx 4 + 6 - 2 \times 10 = -10\] \[D_{yy} \approx 14 + 9 - 2 \times 10 = 3\] \[D_{xy} \approx \frac{6 + 3 - 5 - 1}{4} = 0.75\]Vậy Hessian 2D tại điểm này là:
\[H = \begin{bmatrix} -10 & 0.75 \\ 0.75 & 3 \end{bmatrix}\]Công thức trace và determinant của Hessian
Với một ma trận 2x2 bất kỳ:
\[H = \begin{bmatrix} a & b \\ b & c \end{bmatrix} = \begin{bmatrix} D_{xx} & D_{xy} \\ D_{xy} & D_{yy} \end{bmatrix}\]thì \(Tr(H)\) và \(Det(H)\) được tính toán như sau:
\[Tr(H) = a + c = D_{xx} + D_{yy}\] \[Det(H) = ac - b^2 = D_{xx}D_{yy} - D_{xy}^2\]Áp dụng vào tính ma trận H ví dụ trên thì:
- \[Tr(H) = -10 + 3 = -7\]
- \[Det(H) = -10 \times 3 - 0.75^2 = -30 - 0.5625 = -30.5625\]
Ý nghĩa cốt lõi của tỉ lệ Tr(H)^2 / Det(H)
Nói chính xác hơn, trong SIFT ta dùng:
\[\frac{Tr(H)^2}{Det(H)}\]chứ không phải chỉ \(Tr(H) / Det(H)\) đơn thuần.
Muốn hiểu sâu bản chất, ta nhìn Hessian qua 2 trị riêng (eigenvalues) của nó, giả sử là \(\alpha\) và \(\beta\). Khi đó:
\[Tr(H) = \alpha + \beta\] \[Det(H) = \alpha\beta\]Suy ra:
\[\frac{Tr(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha\beta}\]Đây mới là ý nghĩa quan trọng nhất:
- Nếu \(\alpha\) và \(\beta\) có độ lớn gần nhau, tức độ cong theo hai hướng đều mạnh và khá cân bằng, thì tỉ lệ này không quá lớn.
- Nếu một trị riêng rất lớn còn trị riêng kia rất nhỏ, tức lệch mạnh theo một hướng nhưng yếu theo hướng còn lại, thì tỉ lệ này sẽ rất lớn.
Ví dụ:
- Nếu \(\alpha = 10\) và \(\beta = 10\) thì tỉ lệ này là \(4\)
- Nếu \(\alpha = 20\) và \(\beta = 1\) thì tỉ lệ này là \(22.05\)
Tỉ lệ nhỏ, nghĩa là hai hướng khá cân bằng -> giống hình ảnh của một góc / đốm (corner/blob).
Tỉ lệ lớn, nghĩa là một hướng áp đảo hướng còn lại -> rất giống cạnh.
Vì vậy, bản chất của việc kiểm tra \(Tr(H)^2 / Det(H)\) là để phát hiện góc lấy làm keypoint và loại bỏ cạnh.
5. Gán hướng chính cho keypoint
Đến đây ta đã có vị trí và scale. Nhưng nếu ảnh bị xoay thì sao?
SIFT xử lý bằng cách gán hướng chính (orientation) cho mỗi keypoint. Cách làm:
- Lấy vùng lân cận quanh keypoint ở đúng scale tương ứng.
- Tính gradient theo trục x và y.
- Từ gradient, tính magnitude và orientation.
- Vote orientation vào histogram 36 bins, mỗi bin ứng với 10 độ.
- Đỉnh histogram cao nhất là hướng chính của keypoint.
Nếu có các đỉnh phụ đạt ít nhất khoảng 80% đỉnh chính, Lowe tạo thêm keypoint mới cùng vị trí và scale nhưng khác orientation. Đây là lý do một vị trí ảnh có thể sinh ra nhiều keypoint.
Sau khi đã có orientation, mọi bước descriptor phía sau sẽ được tính tương đối theo hướng này. Nhờ vậy descriptor bền vững hơn khi ảnh bị xoay.
6. Tạo descriptor 128 chiều
Đây là phần rất đẹp của SIFT → kết xuất keypoint thành vector đặc trưng.
Quanh mỗi keypoint, SIFT lấy một vùng lân cận đã được xoay theo orientation của keypoint. Vùng này được chia thành lưới 4 x 4 cell. Trong mỗi cell, ta tạo histogram hướng gradient gồm 8 bins.
Vậy tổng số chiều descriptor là: \(4 \times 4 \times 8 = 128\)
Mỗi phần tử trong vector 128 chiều là lượng gradient vote vào một hướng nào đó tại một cell nào đó. Descriptor này sau đó được normalize để giảm ảnh hưởng của thay đổi sáng. Paper Lowe còn clamp các giá trị lớn, thường ở mức 0.2, rồi normalize lại lần nữa để tăng độ bền vững.
SIFT bền vững với những biến đổi nào?
SIFT được thiết kế để tương đối bền vững với:
- Thay đổi kích thước (scale): nhờ scale-space và octave.
- Phép xoay (rotation): nhờ gán orientation.
- Thay đổi chiếu sáng nhẹ (illumination): nhờ dùng gradient và normalize descriptor.
- Nhiễu nhẹ (noise): nhờ Gaussian blur và chọn keypoint ổn định.
- Thay đổi góc quay một phần: theo paper Lowe 2004, SIFT vẫn mạnh mẽ trong một khoảng thay đổi "affine distortion" / "3D viewpoint" nhất định, nhưng không phải bất biến tuyệt đối.
Điểm cần nhớ: SIFT rất mạnh, nhưng vẫn "fail" nếu vật thể biến dạng quá nhiều, bị che khuất nặng, mờ do di chuyển, hoặc ảnh có quá ít kết cấu "texture" ─ ảnh trơn.
Code hiện thực SIFT từ đầu bằng Python + OpenCV
Trong thực tế, nếu cần dùng SIFT cho production / thí nghiệm nghiêm túc, bạn nên dùng bản tối ưu của OpenCV:
sift = cv2.SIFT_create()
keypoints, descriptors = sift.detectAndCompute(gray, None)
Nhưng để hiểu thuật toán, ta sẽ tự hiện thực một bản mini SIFT. Bản này:
- Dùng OpenCV để đọc ảnh, resize, Gaussian blur, vẽ keypoint và lưu ảnh.
- Không dùng
cv2.SIFT_create(). - Tự xây Gaussian pyramid.
- Tự tính DoG.
- Tự tìm extrema 3D.
- Tự lọc cạnh (edge).
- Tự gán hướng (orientation).
- Tự tạo vector 128 chiều (descriptor).
- Lưu ảnh visualization ra disk.
Lưu ý: code bên dưới là bản học thuật, đã giản lược "subpixel localization" và "trilinear interpolation" trong descriptor để dễ đọc. Vì vậy kết quả sẽ không giống 100% OpenCV SIFT, nhưng pipeline chính là đúng tinh thần thuật toán.
Ảnh mẫu đầu vào:
linhdan.0712.jpg (Cre: @linhdan.0712.jpg)
sift_from_scratch.py
import argparse
import math
from pathlib import Path
import cv2
import numpy as np
def to_gray_float(image_bgr):
gray = cv2.cvtColor(image_bgr, cv2.COLOR_BGR2GRAY)
return gray.astype(np.float32) / 255.0
def build_gaussian_pyramid(gray, num_octaves=4, num_layers=3, sigma=1.6):
pyramid = []
k = 2 ** (1.0 / num_layers)
base = cv2.GaussianBlur(gray, (0, 0), sigmaX=sigma, sigmaY=sigma)
for _ in range(num_octaves):
gaussian_images = [base]
current_sigma = sigma
for layer in range(1, num_layers + 3):
next_sigma = sigma * (k ** layer)
sigma_to_apply = math.sqrt(max(next_sigma ** 2 - current_sigma ** 2, 1e-8))
blurred = cv2.GaussianBlur(gaussian_images[-1], (0, 0), sigmaX=sigma_to_apply, sigmaY=sigma_to_apply)
gaussian_images.append(blurred)
current_sigma = next_sigma
pyramid.append(gaussian_images)
next_base = gaussian_images[num_layers]
h, w = next_base.shape
if min(h, w) < 32:
break
base = cv2.resize(next_base, (w // 2, h // 2), interpolation=cv2.INTER_NEAREST)
return pyramid, k
def build_dog_pyramid(gaussian_pyramid):
dog_pyramid = []
for gaussian_images in gaussian_pyramid:
dogs = []
for idx in range(1, len(gaussian_images)):
dogs.append(gaussian_images[idx] - gaussian_images[idx - 1])
dog_pyramid.append(dogs)
return dog_pyramid
def is_3d_extremum(cube, threshold):
center = cube[1, 1, 1]
if abs(center) < threshold:
return False
neighbors = np.delete(cube.flatten(), 13)
return center > neighbors.max() or center < neighbors.min()
def pass_edge_response_test(dog, x, y, edge_threshold=10.0):
value = dog[y, x]
dxx = dog[y, x + 1] + dog[y, x - 1] - 2 * value
dyy = dog[y + 1, x] + dog[y - 1, x] - 2 * value
dxy = (dog[y + 1, x + 1] - dog[y + 1, x - 1] - dog[y - 1, x + 1] + dog[y - 1, x - 1]) / 4.0
trace = dxx + dyy
determinant = dxx * dyy - dxy * dxy
if determinant <= 0:
return False
ratio = (trace * trace) / determinant
limit = ((edge_threshold + 1) ** 2) / edge_threshold
return ratio < limit
def orientation_histogram(image, x, y, scale, num_bins=36):
radius = int(round(3 * 1.5 * scale))
weight_sigma = 1.5 * scale
h, w = image.shape
hist = np.zeros(num_bins, dtype=np.float32)
for yy in range(y - radius, y + radius + 1):
if yy <= 0 or yy >= h - 1:
continue
for xx in range(x - radius, x + radius + 1):
if xx <= 0 or xx >= w - 1:
continue
dx = image[yy, xx + 1] - image[yy, xx - 1]
dy = image[yy + 1, xx] - image[yy - 1, xx]
magnitude = math.sqrt(dx * dx + dy * dy)
angle = (math.degrees(math.atan2(dy, dx)) + 360.0) % 360.0
dist2 = (xx - x) ** 2 + (yy - y) ** 2
weight = math.exp(-dist2 / (2 * weight_sigma * weight_sigma))
bin_idx = int(round(angle / 360.0 * num_bins)) % num_bins
hist[bin_idx] += weight * magnitude
for _ in range(2):
hist = (np.roll(hist, 1) + hist + np.roll(hist, -1)) / 3.0
return hist
def detect_keypoints(gaussian_pyramid, dog_pyramid, k, sigma=1.6, contrast_threshold=0.01, edge_threshold=10.0):
keypoints = []
for octave_idx, dogs in enumerate(dog_pyramid):
for layer_idx in range(1, len(dogs) - 1):
dog_prev = dogs[layer_idx - 1]
dog_curr = dogs[layer_idx]
dog_next = dogs[layer_idx + 1]
h, w = dog_curr.shape
for y in range(8, h - 8):
for x in range(8, w - 8):
cube = np.stack([
dog_prev[y - 1:y + 2, x - 1:x + 2],
dog_curr[y - 1:y + 2, x - 1:x + 2],
dog_next[y - 1:y + 2, x - 1:x + 2],
])
if not is_3d_extremum(cube, contrast_threshold):
continue
if not pass_edge_response_test(dog_curr, x, y, edge_threshold):
continue
scale_octave = sigma * (k ** layer_idx)
hist = orientation_histogram(gaussian_pyramid[octave_idx][layer_idx], x, y, scale_octave)
if hist.max() <= 0:
continue
peak_threshold = 0.8 * hist.max()
for bin_idx, value in enumerate(hist):
left = hist[(bin_idx - 1) % len(hist)]
right = hist[(bin_idx + 1) % len(hist)]
if value >= peak_threshold and value >= left and value >= right:
angle = 360.0 * bin_idx / len(hist)
factor = 2 ** octave_idx
keypoints.append({
"x": float(x * factor),
"y": float(y * factor),
"octave_x": x,
"octave_y": y,
"octave": octave_idx,
"layer": layer_idx,
"scale_octave": scale_octave,
"size": float(2 * scale_octave * factor),
"angle": angle,
"response": float(abs(dog_curr[y, x])),
})
return keypoints
def compute_descriptor(image, keypoint, descriptor_width=4, num_bins=8):
x = keypoint["octave_x"]
y = keypoint["octave_y"]
angle = keypoint["angle"]
scale = keypoint["scale_octave"]
h, w = image.shape
angle_rad = math.radians(angle)
cos_t = math.cos(angle_rad)
sin_t = math.sin(angle_rad)
hist_width = 3.0 * scale
radius = int(round(hist_width * math.sqrt(2) * (descriptor_width + 1) / 2))
hist = np.zeros((descriptor_width, descriptor_width, num_bins), dtype=np.float32)
for dy in range(-radius, radius + 1):
yy = y + dy
if yy <= 0 or yy >= h - 1:
continue
for dx in range(-radius, radius + 1):
xx = x + dx
if xx <= 0 or xx >= w - 1:
continue
rot_x = (cos_t * dx + sin_t * dy) / hist_width
rot_y = (-sin_t * dx + cos_t * dy) / hist_width
xbin = rot_x + descriptor_width / 2 - 0.5
ybin = rot_y + descriptor_width / 2 - 0.5
if not (0 <= xbin < descriptor_width and 0 <= ybin < descriptor_width):
continue
gx = image[yy, xx + 1] - image[yy, xx - 1]
gy = image[yy + 1, xx] - image[yy - 1, xx]
magnitude = math.sqrt(gx * gx + gy * gy)
orientation = (math.degrees(math.atan2(gy, gx)) - angle + 360.0) % 360.0
obin = int(orientation / 360.0 * num_bins) % num_bins
gaussian_weight = math.exp(-(rot_x * rot_x + rot_y * rot_y) / (2 * (0.5 * descriptor_width) ** 2))
hist[int(ybin), int(xbin), obin] += gaussian_weight * magnitude
descriptor = hist.flatten()
descriptor /= np.linalg.norm(descriptor) + 1e-7
descriptor = np.clip(descriptor, 0, 0.2)
descriptor /= np.linalg.norm(descriptor) + 1e-7
return descriptor.astype(np.float32)
def sift_from_scratch(image_path, output_path, max_keypoints=300):
image_bgr = cv2.imread(str(image_path))
if image_bgr is None:
raise FileNotFoundError(f"Không đọc được ảnh: {image_path}")
gray = to_gray_float(image_bgr)
min_side = min(gray.shape)
num_octaves = max(1, min(4, int(math.log2(min_side)) - 4))
gaussian_pyramid, k = build_gaussian_pyramid(gray, num_octaves=num_octaves)
dog_pyramid = build_dog_pyramid(gaussian_pyramid)
keypoints = detect_keypoints(gaussian_pyramid, dog_pyramid, k)
keypoints = sorted(keypoints, key=lambda item: item["response"], reverse=True)[:max_keypoints]
descriptors = []
valid_keypoints = []
for kp in keypoints:
octave_image = gaussian_pyramid[kp["octave"]][kp["layer"]]
descriptor = compute_descriptor(octave_image, kp)
if np.linalg.norm(descriptor) > 0:
descriptors.append(descriptor)
valid_keypoints.append(kp)
descriptors = np.vstack(descriptors) if descriptors else np.zeros((0, 128), dtype=np.float32)
cv_keypoints = [
cv2.KeyPoint(kp["x"], kp["y"], kp["size"], kp["angle"], kp["response"])
for kp in valid_keypoints
]
vis = cv2.drawKeypoints(
image_bgr,
cv_keypoints,
None,
flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS,
)
output_path = Path(output_path)
output_path.parent.mkdir(parents=True, exist_ok=True)
cv2.imwrite(str(output_path), vis)
descriptor_path = output_path.with_suffix(".npy")
np.save(str(descriptor_path), descriptors)
print(f"Số keypoint: {len(valid_keypoints)}")
print(f"Descriptor shape: {descriptors.shape}")
print(f"Đã lưu ảnh visualization: {output_path}")
print(f"Đã lưu descriptor: {descriptor_path}")
def main():
parser = argparse.ArgumentParser(description="Mini SIFT from scratch bằng Python + OpenCV")
parser.add_argument("--image", required=True, help="Đường dẫn ảnh đầu vào")
parser.add_argument("--output", default="outputs/sift_keypoints.jpg", help="Đường dẫn ảnh output")
parser.add_argument("--max-keypoints", type=int, default=300, help="Số keypoint tối đa để vẽ")
args = parser.parse_args()
sift_from_scratch(args.image, args.output, args.max_keypoints)
if __name__ == "__main__":
main()
Cách chạy:
python sift_from_scratch.py --image linhdan.0712.jpg --output linhdan.0712-SIFT.jpg
Kết quả ảnh in ra:
Thời gian thực thi 23s cho ảnh đầu vào độ phân giải 640x800:
Read image linhdan.0712.jpg -> resolution: 640x800
Gaussian pyramid octaves: 4
Built Gaussian pyramid
Built DoG pyramid
Detected 877 keypoints (before filtering)
Valid keypoints after response filter: 300 / 300
Saved visualization to linhdan.0712-SIFT.jpg
real 0m22.655s
user 0m23.104s
sys 0m3.757s
Một số điểm quan trọng trong script:
build_gaussian_pyramid(): tạo nhiều ảnh blur trong từng octave.build_dog_pyramid(): trừ hai ảnh Gaussian liên tiếp để tạo DoG.is_3d_extremum(): so một điểm với 26 láng giềng trong DoG volume.pass_edge_response_test(): loại điểm giống cạnh bằng trace / det của Hessian.orientation_histogram(): tạo histogram 36 bins để tìm hướng chính của keypoint.compute_descriptor(): chia vùng quanh keypoint thành lưới 4x4, mỗi cell có histogram 8 hướng, cuối cùng tạo vector 128 chiều.cv2.drawKeypoints(): trực quan hóa keypoint và lưu ảnh ra disk.
Trong bản SIFT đầy đủ, phần localization và descriptor còn tinh tế hơn nữa: nội suy subpixel, nội suy peak orientation, trilinear interpolation khi vote vào descriptor, xử lý duplicate keypoint, nhiều chi tiết tối ưu tốc độ, … Bài này cố tình giữ code vừa đủ để bạn đọc được từ trên xuống dưới.
SIFT trong OpenCV
Sau khi hiểu thuật toán, dùng OpenCV sẽ rất gọn. Thời gian thực thi chưa tới 1s cho ảnh đầu vào độ phân giải 640x800.
sift_opencv.py
import cv2
img = cv2.imread("linhdan.0712.jpg")
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
sift = cv2.SIFT_create()
keypoints, descriptors = sift.detectAndCompute(gray, None)
out = cv2.drawKeypoints(img, keypoints, None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imwrite("linhdan.0712-SIFT-OPENCV.jpg", out)
print(descriptors.shape)
Nhưng nếu chỉ dùng API ngay từ đầu thì ta sẽ dễ bỏ qua các chi tiết của giải thuật (hiểu sâu). Tự hiện thực một bản mini như trên sẽ giúp bạn hiểu và nhớ lâu hơn.
Nếu bạn đã học HOG ở Tut 17 thì sẽ thấy SIFT descriptor cũng có tinh thần tương tự: gom các hướng gradient thành histogram. Khác biệt lớn là SIFT làm điều đó quanh từng keypoint cục bộ, ở scale và orientation riêng của keypoint đó.
Bài viết tiếp theo ✍️: Tut 24 ─ Trích đặc trưng SURF (coming soon)
Cảm ơn bạn đã theo dõi bài viết. Hãy kết nối với tớ nhé!
- Minh: https://www.facebook.com/minhng.info
- Khám phá xử lý ảnh - GVGroup: https://www.facebook.com/groups/ip.gvgroup
Khám phá xử lý ảnh - GVGroup
Danh sách bài viết series OpenCV:
- Hashtag #OpenCV
- Tut 1: Xử lý ảnh - OpenCV đọc ghi hình ảnh (code Python và C++)
- Tut 1.1: Xử lý ảnh - Cấu trúc dữ liệu ảnh trong OpenCV. Pixel là gì?
- Tut 1.2: Xử lý ảnh - Chuyển đổi ảnh OpenCV sang Pillow và ngược lại
- Tut 1.3: Vì sao gốc tọa độ ảnh ở góc trên bên trái?
- Tut 2: Xử lý ảnh - OpenCV resize, crop và padding hình ảnh (code Python và C++)
- Tut 3: Xử lý ảnh - OpenCV biến đổi mức sáng hình ảnh (code Python)
- Tut 4: Xử lý ảnh - OpenCV vùng quan tâm (ROI) là gì? (code Python)
- Tut 4.1: Xử lý ảnh - OpenCV: vẽ văn bản, đường thẳng, mũi tên, hình chữ nhật, hình tròn, ellipse, đa giác
- Tut 4.2: Xử lý ảnh - Pha trộn ảnh trong OpenCV (blending)
- Tut 5: Xử lý ảnh - OpenCV ảnh nhị phân
- Tut 6: Xử lý ảnh - OpenCV cân bằng sáng (histogram equalization)
- Tut 7: Xử lý ảnh - OpenCV kỹ thuật cửa sổ trượt (sliding window)
- Tut 8: Xử lý ảnh - Convolution là gì?
- Tut 9: Xử lý ảnh - Làm mờ ảnh (blur)
- Tut 10: Xử lý ảnh - Gradient của ảnh là gì?
- Tut 11: Xử lý ảnh - Phát hiện cạnh Canny (Canny Edge Detection)
- Tut 12: Xử lý ảnh - Phát hiện đường thẳng bằng Hough Transform (Hough Line)
- Tut 13: Xử lý ảnh - Hiện thực phát hiện đoạn thẳng dùng Hough Transform (Hough Line)
- Tut 14: Xử lý ảnh - Giải thuật phân vùng Region Growing trên ảnh màu
- Tut 15: Xử lý ảnh - Giải thuật Background Subtraction trên ảnh màu
- Tut 16: Xử lý ảnh - Frame Subtraction để phát hiện chuyển động trong video
- Tut 17: Xử lý ảnh - HOG - Histograms of Oriented Gradients
- Tut 18: Xử lý ảnh - HOG - Huấn luyện mô hình phân loại người
- Tut 19: Xử lý ảnh - HOG - Phát hiện người
- Tut 20: Xử lý ảnh - Tổng hợp kinh nghiệm xử lý ảnh (End)
- Tut 21: Xử lý ảnh - Hiện thực trích đặc trưng Local Binary Patterns (LBP)
- Tut 22: Xử lý ảnh - Trích đặc trưng Gabor filters
- Tut 23: Xử lý ảnh - Trích đặc trưng SIFT
